##// END OF EJS Templates
schain-jc
RSS
Now there are two receiver units one for data and one for plots
Juan C. Espinoza -
r957:d3acc9060c1d
parent child
Show More
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
@@ -1,1228 +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 def __init__(self):
118
119 raise NotImplementedError
120
121 117 def copy(self, inputObj=None):
122 118
123 119 if inputObj == None:
124 120 return copy.deepcopy(self)
125 121
126 122 for key in inputObj.__dict__.keys():
127 123
128 124 attribute = inputObj.__dict__[key]
129 125
130 126 #If this attribute is a tuple or list
131 127 if type(inputObj.__dict__[key]) in (tuple, list):
132 128 self.__dict__[key] = attribute[:]
133 129 continue
134 130
135 131 #If this attribute is another object or instance
136 132 if hasattr(attribute, '__dict__'):
137 133 self.__dict__[key] = attribute.copy()
138 134 continue
139 135
140 136 self.__dict__[key] = inputObj.__dict__[key]
141 137
142 138 def deepcopy(self):
143 139
144 140 return copy.deepcopy(self)
145 141
146 142 def isEmpty(self):
147 143
148 144 return self.flagNoData
149 145
150 146 class JROData(GenericData):
151 147
152 148 # m_BasicHeader = BasicHeader()
153 149 # m_ProcessingHeader = ProcessingHeader()
154 150
155 151 systemHeaderObj = SystemHeader()
156 152
157 153 radarControllerHeaderObj = RadarControllerHeader()
158 154
159 155 # data = None
160 156
161 157 type = None
162 158
163 159 datatype = None #dtype but in string
164 160
165 161 # dtype = None
166 162
167 163 # nChannels = None
168 164
169 165 # nHeights = None
170 166
171 167 nProfiles = None
172 168
173 169 heightList = None
174 170
175 171 channelList = None
176 172
177 173 flagDiscontinuousBlock = False
178 174
179 175 useLocalTime = False
180 176
181 177 utctime = None
182 178
183 179 timeZone = None
184 180
185 181 dstFlag = None
186 182
187 183 errorCount = None
188 184
189 185 blocksize = None
190 186
191 187 # nCode = None
192 188 #
193 189 # nBaud = None
194 190 #
195 191 # code = None
196 192
197 193 flagDecodeData = False #asumo q la data no esta decodificada
198 194
199 195 flagDeflipData = False #asumo q la data no esta sin flip
200 196
201 197 flagShiftFFT = False
202 198
203 199 # ippSeconds = None
204 200
205 201 # timeInterval = None
206 202
207 203 nCohInt = None
208 204
209 205 # noise = None
210 206
211 207 windowOfFilter = 1
212 208
213 209 #Speed of ligth
214 210 C = 3e8
215 211
216 212 frequency = 49.92e6
217 213
218 214 realtime = False
219 215
220 216 beacon_heiIndexList = None
221 217
222 218 last_block = None
223 219
224 220 blocknow = None
225 221
226 222 azimuth = None
227 223
228 224 zenith = None
229 225
230 226 beam = Beam()
231 227
232 228 profileIndex = None
233 229
234 def __init__(self):
235
236 raise NotImplementedError
237
238 230 def getNoise(self):
239 231
240 232 raise NotImplementedError
241 233
242 234 def getNChannels(self):
243 235
244 236 return len(self.channelList)
245 237
246 238 def getChannelIndexList(self):
247 239
248 240 return range(self.nChannels)
249 241
250 242 def getNHeights(self):
251 243
252 244 return len(self.heightList)
253 245
254 246 def getHeiRange(self, extrapoints=0):
255 247
256 248 heis = self.heightList
257 249 # deltah = self.heightList[1] - self.heightList[0]
258 250 #
259 251 # heis.append(self.heightList[-1])
260 252
261 253 return heis
262 254
263 255 def getDeltaH(self):
264 256
265 257 delta = self.heightList[1] - self.heightList[0]
266 258
267 259 return delta
268 260
269 261 def getltctime(self):
270 262
271 263 if self.useLocalTime:
272 264 return self.utctime - self.timeZone*60
273 265
274 266 return self.utctime
275 267
276 268 def getDatatime(self):
277 269
278 270 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
279 271 return datatimeValue
280 272
281 273 def getTimeRange(self):
282 274
283 275 datatime = []
284 276
285 277 datatime.append(self.ltctime)
286 278 datatime.append(self.ltctime + self.timeInterval+1)
287 279
288 280 datatime = numpy.array(datatime)
289 281
290 282 return datatime
291 283
292 284 def getFmaxTimeResponse(self):
293 285
294 286 period = (10**-6)*self.getDeltaH()/(0.15)
295 287
296 288 PRF = 1./(period * self.nCohInt)
297 289
298 290 fmax = PRF
299 291
300 292 return fmax
301 293
302 294 def getFmax(self):
303 295
304 296 PRF = 1./(self.ippSeconds * self.nCohInt)
305 297
306 298 fmax = PRF
307 299
308 300 return fmax
309 301
310 302 def getVmax(self):
311 303
312 304 _lambda = self.C/self.frequency
313 305
314 306 vmax = self.getFmax() * _lambda/2
315 307
316 308 return vmax
317 309
318 310 def get_ippSeconds(self):
319 311 '''
320 312 '''
321 313 return self.radarControllerHeaderObj.ippSeconds
322 314
323 315 def set_ippSeconds(self, ippSeconds):
324 316 '''
325 317 '''
326 318
327 319 self.radarControllerHeaderObj.ippSeconds = ippSeconds
328 320
329 321 return
330 322
331 323 def get_dtype(self):
332 324 '''
333 325 '''
334 326 return getNumpyDtype(self.datatype)
335 327
336 328 def set_dtype(self, numpyDtype):
337 329 '''
338 330 '''
339 331
340 332 self.datatype = getDataTypeCode(numpyDtype)
341 333
342 334 def get_code(self):
343 335 '''
344 336 '''
345 337 return self.radarControllerHeaderObj.code
346 338
347 339 def set_code(self, code):
348 340 '''
349 341 '''
350 342 self.radarControllerHeaderObj.code = code
351 343
352 344 return
353 345
354 346 def get_ncode(self):
355 347 '''
356 348 '''
357 349 return self.radarControllerHeaderObj.nCode
358 350
359 351 def set_ncode(self, nCode):
360 352 '''
361 353 '''
362 354 self.radarControllerHeaderObj.nCode = nCode
363 355
364 356 return
365 357
366 358 def get_nbaud(self):
367 359 '''
368 360 '''
369 361 return self.radarControllerHeaderObj.nBaud
370 362
371 363 def set_nbaud(self, nBaud):
372 364 '''
373 365 '''
374 366 self.radarControllerHeaderObj.nBaud = nBaud
375 367
376 368 return
377 369
378 370 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
379 371 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
380 372 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
381 373 #noise = property(getNoise, "I'm the 'nHeights' property.")
382 374 datatime = property(getDatatime, "I'm the 'datatime' property")
383 375 ltctime = property(getltctime, "I'm the 'ltctime' property")
384 376 ippSeconds = property(get_ippSeconds, set_ippSeconds)
385 377 dtype = property(get_dtype, set_dtype)
386 378 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
387 379 code = property(get_code, set_code)
388 380 nCode = property(get_ncode, set_ncode)
389 381 nBaud = property(get_nbaud, set_nbaud)
390 382
391 383 class Voltage(JROData):
392 384
393 385 #data es un numpy array de 2 dmensiones (canales, alturas)
394 386 data = None
395 387
396 388 def __init__(self):
397 389 '''
398 390 Constructor
399 391 '''
400 392
401 393 self.useLocalTime = True
402 394
403 395 self.radarControllerHeaderObj = RadarControllerHeader()
404 396
405 397 self.systemHeaderObj = SystemHeader()
406 398
407 399 self.type = "Voltage"
408 400
409 401 self.data = None
410 402
411 403 # self.dtype = None
412 404
413 405 # self.nChannels = 0
414 406
415 407 # self.nHeights = 0
416 408
417 409 self.nProfiles = None
418 410
419 411 self.heightList = None
420 412
421 413 self.channelList = None
422 414
423 415 # self.channelIndexList = None
424 416
425 417 self.flagNoData = True
426 418
427 419 self.flagDiscontinuousBlock = False
428 420
429 421 self.utctime = None
430 422
431 423 self.timeZone = None
432 424
433 425 self.dstFlag = None
434 426
435 427 self.errorCount = None
436 428
437 429 self.nCohInt = None
438 430
439 431 self.blocksize = None
440 432
441 433 self.flagDecodeData = False #asumo q la data no esta decodificada
442 434
443 435 self.flagDeflipData = False #asumo q la data no esta sin flip
444 436
445 437 self.flagShiftFFT = False
446 438
447 439 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
448 440
449 441 self.profileIndex = 0
450 442
451 443 def getNoisebyHildebrand(self, channel = None):
452 444 """
453 445 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
454 446
455 447 Return:
456 448 noiselevel
457 449 """
458 450
459 451 if channel != None:
460 452 data = self.data[channel]
461 453 nChannels = 1
462 454 else:
463 455 data = self.data
464 456 nChannels = self.nChannels
465 457
466 458 noise = numpy.zeros(nChannels)
467 459 power = data * numpy.conjugate(data)
468 460
469 461 for thisChannel in range(nChannels):
470 462 if nChannels == 1:
471 463 daux = power[:].real
472 464 else:
473 465 daux = power[thisChannel,:].real
474 466 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
475 467
476 468 return noise
477 469
478 470 def getNoise(self, type = 1, channel = None):
479 471
480 472 if type == 1:
481 473 noise = self.getNoisebyHildebrand(channel)
482 474
483 475 return noise
484 476
485 477 def getPower(self, channel = None):
486 478
487 479 if channel != None:
488 480 data = self.data[channel]
489 481 else:
490 482 data = self.data
491 483
492 484 power = data * numpy.conjugate(data)
493 485 powerdB = 10*numpy.log10(power.real)
494 486 powerdB = numpy.squeeze(powerdB)
495 487
496 488 return powerdB
497 489
498 490 def getTimeInterval(self):
499 491
500 492 timeInterval = self.ippSeconds * self.nCohInt
501 493
502 494 return timeInterval
503 495
504 496 noise = property(getNoise, "I'm the 'nHeights' property.")
505 497 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
506 498
507 499 class Spectra(JROData):
508 500
509 501 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
510 502 data_spc = None
511 503
512 504 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
513 505 data_cspc = None
514 506
515 507 #data dc es un numpy array de 2 dmensiones (canales, alturas)
516 508 data_dc = None
517 509
518 510 #data power
519 511 data_pwr = None
520 512
521 513 nFFTPoints = None
522 514
523 515 # nPairs = None
524 516
525 517 pairsList = None
526 518
527 519 nIncohInt = None
528 520
529 521 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
530 522
531 523 nCohInt = None #se requiere para determinar el valor de timeInterval
532 524
533 525 ippFactor = None
534 526
535 527 profileIndex = 0
536 528
537 529 plotting = "spectra"
538 530
539 531 def __init__(self):
540 532 '''
541 533 Constructor
542 534 '''
543 535
544 536 self.useLocalTime = True
545 537
546 538 self.radarControllerHeaderObj = RadarControllerHeader()
547 539
548 540 self.systemHeaderObj = SystemHeader()
549 541
550 542 self.type = "Spectra"
551 543
552 544 # self.data = None
553 545
554 546 # self.dtype = None
555 547
556 548 # self.nChannels = 0
557 549
558 550 # self.nHeights = 0
559 551
560 552 self.nProfiles = None
561 553
562 554 self.heightList = None
563 555
564 556 self.channelList = None
565 557
566 558 # self.channelIndexList = None
567 559
568 560 self.pairsList = None
569 561
570 562 self.flagNoData = True
571 563
572 564 self.flagDiscontinuousBlock = False
573 565
574 566 self.utctime = None
575 567
576 568 self.nCohInt = None
577 569
578 570 self.nIncohInt = None
579 571
580 572 self.blocksize = None
581 573
582 574 self.nFFTPoints = None
583 575
584 576 self.wavelength = None
585 577
586 578 self.flagDecodeData = False #asumo q la data no esta decodificada
587 579
588 580 self.flagDeflipData = False #asumo q la data no esta sin flip
589 581
590 582 self.flagShiftFFT = False
591 583
592 584 self.ippFactor = 1
593 585
594 586 #self.noise = None
595 587
596 588 self.beacon_heiIndexList = []
597 589
598 590 self.noise_estimation = None
599 591
600 592
601 593 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
602 594 """
603 595 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
604 596
605 597 Return:
606 598 noiselevel
607 599 """
608 600
609 601 noise = numpy.zeros(self.nChannels)
610 602
611 603 for channel in range(self.nChannels):
612 604 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
613 605 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
614 606
615 607 return noise
616 608
617 609 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
618 610
619 611 if self.noise_estimation is not None:
620 612 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
621 613 else:
622 614 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
623 615 return noise
624 616
625 617 def getFreqRangeTimeResponse(self, extrapoints=0):
626 618
627 619 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
628 620 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
629 621
630 622 return freqrange
631 623
632 624 def getAcfRange(self, extrapoints=0):
633 625
634 626 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
635 627 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
636 628
637 629 return freqrange
638 630
639 631 def getFreqRange(self, extrapoints=0):
640 632
641 633 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
642 634 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
643 635
644 636 return freqrange
645 637
646 638 def getVelRange(self, extrapoints=0):
647 639
648 640 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
649 641 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
650 642
651 643 return velrange
652 644
653 645 def getNPairs(self):
654 646
655 647 return len(self.pairsList)
656 648
657 649 def getPairsIndexList(self):
658 650
659 651 return range(self.nPairs)
660 652
661 653 def getNormFactor(self):
662 654
663 655 pwcode = 1
664 656
665 657 if self.flagDecodeData:
666 658 pwcode = numpy.sum(self.code[0]**2)
667 659 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
668 660 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
669 661
670 662 return normFactor
671 663
672 664 def getFlagCspc(self):
673 665
674 666 if self.data_cspc is None:
675 667 return True
676 668
677 669 return False
678 670
679 671 def getFlagDc(self):
680 672
681 673 if self.data_dc is None:
682 674 return True
683 675
684 676 return False
685 677
686 678 def getTimeInterval(self):
687 679
688 680 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
689 681
690 682 return timeInterval
691 683
692 684 def getPower(self):
693 685
694 686 factor = self.normFactor
695 687 z = self.data_spc/factor
696 688 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
697 689 avg = numpy.average(z, axis=1)
698 690
699 691 return 10*numpy.log10(avg)
700 692
701 693 def getCoherence(self, pairsList=None, phase=False):
702 694
703 695 z = []
704 696 if pairsList is None:
705 697 pairsIndexList = self.pairsIndexList
706 698 else:
707 699 pairsIndexList = []
708 700 for pair in pairsList:
709 701 if pair not in self.pairsList:
710 702 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
711 703 pairsIndexList.append(self.pairsList.index(pair))
712 704 for i in range(len(pairsIndexList)):
713 705 pair = self.pairsList[pairsIndexList[i]]
714 706 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
715 707 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
716 708 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
717 709 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
718 710 if phase:
719 711 data = numpy.arctan2(avgcoherenceComplex.imag,
720 712 avgcoherenceComplex.real)*180/numpy.pi
721 713 else:
722 714 data = numpy.abs(avgcoherenceComplex)
723 715
724 716 z.append(data)
725 717
726 718 return numpy.array(z)
727 719
728 720 def setValue(self, value):
729 721
730 722 print "This property should not be initialized"
731 723
732 724 return
733 725
734 726 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
735 727 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
736 728 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
737 729 flag_cspc = property(getFlagCspc, setValue)
738 730 flag_dc = property(getFlagDc, setValue)
739 731 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
740 732 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
741 733
742 734 class SpectraHeis(Spectra):
743 735
744 736 data_spc = None
745 737
746 738 data_cspc = None
747 739
748 740 data_dc = None
749 741
750 742 nFFTPoints = None
751 743
752 744 # nPairs = None
753 745
754 746 pairsList = None
755 747
756 748 nCohInt = None
757 749
758 750 nIncohInt = None
759 751
760 752 def __init__(self):
761 753
762 754 self.radarControllerHeaderObj = RadarControllerHeader()
763 755
764 756 self.systemHeaderObj = SystemHeader()
765 757
766 758 self.type = "SpectraHeis"
767 759
768 760 # self.dtype = None
769 761
770 762 # self.nChannels = 0
771 763
772 764 # self.nHeights = 0
773 765
774 766 self.nProfiles = None
775 767
776 768 self.heightList = None
777 769
778 770 self.channelList = None
779 771
780 772 # self.channelIndexList = None
781 773
782 774 self.flagNoData = True
783 775
784 776 self.flagDiscontinuousBlock = False
785 777
786 778 # self.nPairs = 0
787 779
788 780 self.utctime = None
789 781
790 782 self.blocksize = None
791 783
792 784 self.profileIndex = 0
793 785
794 786 self.nCohInt = 1
795 787
796 788 self.nIncohInt = 1
797 789
798 790 def getNormFactor(self):
799 791 pwcode = 1
800 792 if self.flagDecodeData:
801 793 pwcode = numpy.sum(self.code[0]**2)
802 794
803 795 normFactor = self.nIncohInt*self.nCohInt*pwcode
804 796
805 797 return normFactor
806 798
807 799 def getTimeInterval(self):
808 800
809 801 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
810 802
811 803 return timeInterval
812 804
813 805 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
814 806 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
815 807
816 808 class Fits(JROData):
817 809
818 810 heightList = None
819 811
820 812 channelList = None
821 813
822 814 flagNoData = True
823 815
824 816 flagDiscontinuousBlock = False
825 817
826 818 useLocalTime = False
827 819
828 820 utctime = None
829 821
830 822 timeZone = None
831 823
832 824 # ippSeconds = None
833 825
834 826 # timeInterval = None
835 827
836 828 nCohInt = None
837 829
838 830 nIncohInt = None
839 831
840 832 noise = None
841 833
842 834 windowOfFilter = 1
843 835
844 836 #Speed of ligth
845 837 C = 3e8
846 838
847 839 frequency = 49.92e6
848 840
849 841 realtime = False
850 842
851 843
852 844 def __init__(self):
853 845
854 846 self.type = "Fits"
855 847
856 848 self.nProfiles = None
857 849
858 850 self.heightList = None
859 851
860 852 self.channelList = None
861 853
862 854 # self.channelIndexList = None
863 855
864 856 self.flagNoData = True
865 857
866 858 self.utctime = None
867 859
868 860 self.nCohInt = 1
869 861
870 862 self.nIncohInt = 1
871 863
872 864 self.useLocalTime = True
873 865
874 866 self.profileIndex = 0
875 867
876 868 # self.utctime = None
877 869 # self.timeZone = None
878 870 # self.ltctime = None
879 871 # self.timeInterval = None
880 872 # self.header = None
881 873 # self.data_header = None
882 874 # self.data = None
883 875 # self.datatime = None
884 876 # self.flagNoData = False
885 877 # self.expName = ''
886 878 # self.nChannels = None
887 879 # self.nSamples = None
888 880 # self.dataBlocksPerFile = None
889 881 # self.comments = ''
890 882 #
891 883
892 884
893 885 def getltctime(self):
894 886
895 887 if self.useLocalTime:
896 888 return self.utctime - self.timeZone*60
897 889
898 890 return self.utctime
899 891
900 892 def getDatatime(self):
901 893
902 894 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
903 895 return datatime
904 896
905 897 def getTimeRange(self):
906 898
907 899 datatime = []
908 900
909 901 datatime.append(self.ltctime)
910 902 datatime.append(self.ltctime + self.timeInterval)
911 903
912 904 datatime = numpy.array(datatime)
913 905
914 906 return datatime
915 907
916 908 def getHeiRange(self):
917 909
918 910 heis = self.heightList
919 911
920 912 return heis
921 913
922 914 def getNHeights(self):
923 915
924 916 return len(self.heightList)
925 917
926 918 def getNChannels(self):
927 919
928 920 return len(self.channelList)
929 921
930 922 def getChannelIndexList(self):
931 923
932 924 return range(self.nChannels)
933 925
934 926 def getNoise(self, type = 1):
935 927
936 928 #noise = numpy.zeros(self.nChannels)
937 929
938 930 if type == 1:
939 931 noise = self.getNoisebyHildebrand()
940 932
941 933 if type == 2:
942 934 noise = self.getNoisebySort()
943 935
944 936 if type == 3:
945 937 noise = self.getNoisebyWindow()
946 938
947 939 return noise
948 940
949 941 def getTimeInterval(self):
950 942
951 943 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
952 944
953 945 return timeInterval
954 946
955 947 datatime = property(getDatatime, "I'm the 'datatime' property")
956 948 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
957 949 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
958 950 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
959 951 noise = property(getNoise, "I'm the 'nHeights' property.")
960 952
961 953 ltctime = property(getltctime, "I'm the 'ltctime' property")
962 954 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
963 955
964 956
965 957 class Correlation(JROData):
966 958
967 959 noise = None
968 960
969 961 SNR = None
970 962
971 963 #--------------------------------------------------
972 964
973 965 mode = None
974 966
975 967 split = False
976 968
977 969 data_cf = None
978 970
979 971 lags = None
980 972
981 973 lagRange = None
982 974
983 975 pairsList = None
984 976
985 977 normFactor = None
986 978
987 979 #--------------------------------------------------
988 980
989 981 # calculateVelocity = None
990 982
991 983 nLags = None
992 984
993 985 nPairs = None
994 986
995 987 nAvg = None
996 988
997 989
998 990 def __init__(self):
999 991 '''
1000 992 Constructor
1001 993 '''
1002 994 self.radarControllerHeaderObj = RadarControllerHeader()
1003 995
1004 996 self.systemHeaderObj = SystemHeader()
1005 997
1006 998 self.type = "Correlation"
1007 999
1008 1000 self.data = None
1009 1001
1010 1002 self.dtype = None
1011 1003
1012 1004 self.nProfiles = None
1013 1005
1014 1006 self.heightList = None
1015 1007
1016 1008 self.channelList = None
1017 1009
1018 1010 self.flagNoData = True
1019 1011
1020 1012 self.flagDiscontinuousBlock = False
1021 1013
1022 1014 self.utctime = None
1023 1015
1024 1016 self.timeZone = None
1025 1017
1026 1018 self.dstFlag = None
1027 1019
1028 1020 self.errorCount = None
1029 1021
1030 1022 self.blocksize = None
1031 1023
1032 1024 self.flagDecodeData = False #asumo q la data no esta decodificada
1033 1025
1034 1026 self.flagDeflipData = False #asumo q la data no esta sin flip
1035 1027
1036 1028 self.pairsList = None
1037 1029
1038 1030 self.nPoints = None
1039 1031
1040 1032 def getPairsList(self):
1041 1033
1042 1034 return self.pairsList
1043 1035
1044 1036 def getNoise(self, mode = 2):
1045 1037
1046 1038 indR = numpy.where(self.lagR == 0)[0][0]
1047 1039 indT = numpy.where(self.lagT == 0)[0][0]
1048 1040
1049 1041 jspectra0 = self.data_corr[:,:,indR,:]
1050 1042 jspectra = copy.copy(jspectra0)
1051 1043
1052 1044 num_chan = jspectra.shape[0]
1053 1045 num_hei = jspectra.shape[2]
1054 1046
1055 1047 freq_dc = jspectra.shape[1]/2
1056 1048 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1057 1049
1058 1050 if ind_vel[0]<0:
1059 1051 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1060 1052
1061 1053 if mode == 1:
1062 1054 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1063 1055
1064 1056 if mode == 2:
1065 1057
1066 1058 vel = numpy.array([-2,-1,1,2])
1067 1059 xx = numpy.zeros([4,4])
1068 1060
1069 1061 for fil in range(4):
1070 1062 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1071 1063
1072 1064 xx_inv = numpy.linalg.inv(xx)
1073 1065 xx_aux = xx_inv[0,:]
1074 1066
1075 1067 for ich in range(num_chan):
1076 1068 yy = jspectra[ich,ind_vel,:]
1077 1069 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1078 1070
1079 1071 junkid = jspectra[ich,freq_dc,:]<=0
1080 1072 cjunkid = sum(junkid)
1081 1073
1082 1074 if cjunkid.any():
1083 1075 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1084 1076
1085 1077 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1086 1078
1087 1079 return noise
1088 1080
1089 1081 def getTimeInterval(self):
1090 1082
1091 1083 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1092 1084
1093 1085 return timeInterval
1094 1086
1095 1087 def splitFunctions(self):
1096 1088
1097 1089 pairsList = self.pairsList
1098 1090 ccf_pairs = []
1099 1091 acf_pairs = []
1100 1092 ccf_ind = []
1101 1093 acf_ind = []
1102 1094 for l in range(len(pairsList)):
1103 1095 chan0 = pairsList[l][0]
1104 1096 chan1 = pairsList[l][1]
1105 1097
1106 1098 #Obteniendo pares de Autocorrelacion
1107 1099 if chan0 == chan1:
1108 1100 acf_pairs.append(chan0)
1109 1101 acf_ind.append(l)
1110 1102 else:
1111 1103 ccf_pairs.append(pairsList[l])
1112 1104 ccf_ind.append(l)
1113 1105
1114 1106 data_acf = self.data_cf[acf_ind]
1115 1107 data_ccf = self.data_cf[ccf_ind]
1116 1108
1117 1109 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1118 1110
1119 1111 def getNormFactor(self):
1120 1112 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1121 1113 acf_pairs = numpy.array(acf_pairs)
1122 1114 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1123 1115
1124 1116 for p in range(self.nPairs):
1125 1117 pair = self.pairsList[p]
1126 1118
1127 1119 ch0 = pair[0]
1128 1120 ch1 = pair[1]
1129 1121
1130 1122 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1131 1123 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1132 1124 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1133 1125
1134 1126 return normFactor
1135 1127
1136 1128 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1137 1129 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1138 1130
1139 1131 class Parameters(Spectra):
1140 1132
1141 1133 experimentInfo = None #Information about the experiment
1142 1134
1143 1135 #Information from previous data
1144 1136
1145 1137 inputUnit = None #Type of data to be processed
1146 1138
1147 1139 operation = None #Type of operation to parametrize
1148 1140
1149 1141 #normFactor = None #Normalization Factor
1150 1142
1151 1143 groupList = None #List of Pairs, Groups, etc
1152 1144
1153 1145 #Parameters
1154 1146
1155 1147 data_param = None #Parameters obtained
1156 1148
1157 1149 data_pre = None #Data Pre Parametrization
1158 1150
1159 1151 data_SNR = None #Signal to Noise Ratio
1160 1152
1161 1153 # heightRange = None #Heights
1162 1154
1163 1155 abscissaList = None #Abscissa, can be velocities, lags or time
1164 1156
1165 1157 # noise = None #Noise Potency
1166 1158
1167 1159 utctimeInit = None #Initial UTC time
1168 1160
1169 1161 paramInterval = None #Time interval to calculate Parameters in seconds
1170 1162
1171 1163 useLocalTime = True
1172 1164
1173 1165 #Fitting
1174 1166
1175 1167 data_error = None #Error of the estimation
1176 1168
1177 1169 constants = None
1178 1170
1179 1171 library = None
1180 1172
1181 1173 #Output signal
1182 1174
1183 1175 outputInterval = None #Time interval to calculate output signal in seconds
1184 1176
1185 1177 data_output = None #Out signal
1186 1178
1187 1179 nAvg = None
1188 1180
1189 1181 noise_estimation = None
1190 1182
1191 1183
1192 1184 def __init__(self):
1193 1185 '''
1194 1186 Constructor
1195 1187 '''
1196 1188 self.radarControllerHeaderObj = RadarControllerHeader()
1197 1189
1198 1190 self.systemHeaderObj = SystemHeader()
1199 1191
1200 1192 self.type = "Parameters"
1201 1193
1202 1194 def getTimeRange1(self, interval):
1203 1195
1204 1196 datatime = []
1205 1197
1206 1198 if self.useLocalTime:
1207 1199 time1 = self.utctimeInit - self.timeZone*60
1208 1200 else:
1209 1201 time1 = self.utctimeInit
1210 1202
1211 1203 datatime.append(time1)
1212 1204 datatime.append(time1 + interval)
1213 1205 datatime = numpy.array(datatime)
1214 1206
1215 1207 return datatime
1216 1208
1217 1209 def getTimeInterval(self):
1218 1210
1219 1211 if hasattr(self, 'timeInterval1'):
1220 1212 return self.timeInterval1
1221 1213 else:
1222 1214 return self.paramInterval
1223 1215
1224 1216 def getNoise(self):
1225 1217
1226 1218 return self.spc_noise
1227 1219
1228 1220 timeInterval = property(getTimeInterval)
Show file before | Show file after | Hide comments
@@ -1,464 +1,502
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 from profilehooks import profile
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 from schainpy.model.data.jrodata import JROData
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 34 # dx = int(len(self.x)/self.__MAXNUMX) + 1
34 35
35 36 dy = int(len(z[0])/MAXNUMY) + 1
36 37
37 38 return z[::, ::dy]
38 39
39 40 class throttle(object):
40 41 """Decorator that prevents a function from being called more than once every
41 42 time period.
42 43 To create a function that cannot be called more than once a minute, but
43 44 will sleep until it can be called:
44 45 @throttle(minutes=1)
45 46 def foo():
46 47 pass
47 48
48 49 for i in range(10):
49 50 foo()
50 51 print "This function has run %s times." % i
51 52 """
52 53
53 54 def __init__(self, seconds=0, minutes=0, hours=0):
54 55 self.throttle_period = datetime.timedelta(
55 56 seconds=seconds, minutes=minutes, hours=hours
56 57 )
57 58
58 59 self.time_of_last_call = datetime.datetime.min
59 60
60 61 def __call__(self, fn):
61 62 @wraps(fn)
62 63 def wrapper(*args, **kwargs):
63 64 now = datetime.datetime.now()
64 65 time_since_last_call = now - self.time_of_last_call
65 66 time_left = self.throttle_period - time_since_last_call
66 67
67 68 if time_left > datetime.timedelta(seconds=0):
68 69 return
69 70
70 71 self.time_of_last_call = datetime.datetime.now()
71 72 return fn(*args, **kwargs)
72 73
73 74 return wrapper
74 75
75 76
76 77 class PublishData(Operation):
77 78 """Clase publish."""
78 79
79 80 def __init__(self, **kwargs):
80 81 """Inicio."""
81 82 Operation.__init__(self, **kwargs)
82 83 self.isConfig = False
83 84 self.client = None
84 85 self.zeromq = None
85 86 self.mqtt = None
86 87
87 88 def on_disconnect(self, client, userdata, rc):
88 89 if rc != 0:
89 90 print("Unexpected disconnection.")
90 91 self.connect()
91 92
92 93 def connect(self):
93 94 print 'trying to connect'
94 95 try:
95 96 self.client.connect(
96 97 host=self.host,
97 98 port=self.port,
98 99 keepalive=60*10,
99 100 bind_address='')
100 101 self.client.loop_start()
101 102 # self.client.publish(
102 103 # self.topic + 'SETUP',
103 104 # json.dumps(setup),
104 105 # retain=True
105 106 # )
106 107 except:
107 108 print "MQTT Conection error."
108 109 self.client = False
109 110
110 111 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, verbose=True, **kwargs):
111 112 self.counter = 0
112 113 self.topic = kwargs.get('topic', 'schain')
113 114 self.delay = kwargs.get('delay', 0)
114 115 self.plottype = kwargs.get('plottype', 'spectra')
115 116 self.host = kwargs.get('host', "10.10.10.82")
116 117 self.port = kwargs.get('port', 3000)
117 118 self.clientId = clientId
118 119 self.cnt = 0
119 120 self.zeromq = zeromq
120 121 self.mqtt = kwargs.get('plottype', 0)
121 122 self.client = None
122 123 self.verbose = verbose
123 124 self.dataOut.firstdata = True
124 125 setup = []
125 126 if mqtt is 1:
126 127 self.client = mqtt.Client(
127 128 client_id=self.clientId + self.topic + 'SCHAIN',
128 129 clean_session=True)
129 130 self.client.on_disconnect = self.on_disconnect
130 131 self.connect()
131 132 for plot in self.plottype:
132 133 setup.append({
133 134 'plot': plot,
134 135 'topic': self.topic + plot,
135 136 'title': getattr(self, plot + '_' + 'title', False),
136 137 'xlabel': getattr(self, plot + '_' + 'xlabel', False),
137 138 'ylabel': getattr(self, plot + '_' + 'ylabel', False),
138 139 'xrange': getattr(self, plot + '_' + 'xrange', False),
139 140 'yrange': getattr(self, plot + '_' + 'yrange', False),
140 141 'zrange': getattr(self, plot + '_' + 'zrange', False),
141 142 })
142 143 if zeromq is 1:
143 144 context = zmq.Context()
144 145 self.zmq_socket = context.socket(zmq.PUSH)
145 146 server = kwargs.get('server', 'zmq.pipe')
146 147
147 148 if 'tcp://' in server:
148 149 address = server
149 150 else:
150 151 address = 'ipc:///tmp/%s' % server
151 152
152 153 self.zmq_socket.connect(address)
153 154 time.sleep(1)
154 155
155 156
156 157 def publish_data(self):
157 158 self.dataOut.finished = False
158 159 if self.mqtt is 1:
159 160 yData = self.dataOut.heightList[:2].tolist()
160 161 if self.plottype == 'spectra':
161 162 data = getattr(self.dataOut, 'data_spc')
162 163 z = data/self.dataOut.normFactor
163 164 zdB = 10*numpy.log10(z)
164 165 xlen, ylen = zdB[0].shape
165 166 dx = int(xlen/MAXNUMX) + 1
166 167 dy = int(ylen/MAXNUMY) + 1
167 168 Z = [0 for i in self.dataOut.channelList]
168 169 for i in self.dataOut.channelList:
169 170 Z[i] = zdB[i][::dx, ::dy].tolist()
170 171 payload = {
171 172 'timestamp': self.dataOut.utctime,
172 173 'data': roundFloats(Z),
173 174 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
174 175 'interval': self.dataOut.getTimeInterval(),
175 176 'type': self.plottype,
176 177 'yData': yData
177 178 }
178 179 # print payload
179 180
180 181 elif self.plottype in ('rti', 'power'):
181 182 data = getattr(self.dataOut, 'data_spc')
182 183 z = data/self.dataOut.normFactor
183 184 avg = numpy.average(z, axis=1)
184 185 avgdB = 10*numpy.log10(avg)
185 186 xlen, ylen = z[0].shape
186 187 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
187 188 AVG = [0 for i in self.dataOut.channelList]
188 189 for i in self.dataOut.channelList:
189 190 AVG[i] = avgdB[i][::dy].tolist()
190 191 payload = {
191 192 'timestamp': self.dataOut.utctime,
192 193 'data': roundFloats(AVG),
193 194 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
194 195 'interval': self.dataOut.getTimeInterval(),
195 196 'type': self.plottype,
196 197 'yData': yData
197 198 }
198 199 elif self.plottype == 'noise':
199 200 noise = self.dataOut.getNoise()/self.dataOut.normFactor
200 201 noisedB = 10*numpy.log10(noise)
201 202 payload = {
202 203 'timestamp': self.dataOut.utctime,
203 204 'data': roundFloats(noisedB.reshape(-1, 1).tolist()),
204 205 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
205 206 'interval': self.dataOut.getTimeInterval(),
206 207 'type': self.plottype,
207 208 'yData': yData
208 209 }
209 210 elif self.plottype == 'snr':
210 211 data = getattr(self.dataOut, 'data_SNR')
211 212 avgdB = 10*numpy.log10(data)
212 213
213 214 ylen = data[0].size
214 215 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
215 216 AVG = [0 for i in self.dataOut.channelList]
216 217 for i in self.dataOut.channelList:
217 218 AVG[i] = avgdB[i][::dy].tolist()
218 219 payload = {
219 220 'timestamp': self.dataOut.utctime,
220 221 'data': roundFloats(AVG),
221 222 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
222 223 'type': self.plottype,
223 224 'yData': yData
224 225 }
225 226 else:
226 227 print "Tipo de grafico invalido"
227 228 payload = {
228 229 'data': 'None',
229 230 'timestamp': 'None',
230 231 'type': None
231 232 }
232 233 # print 'Publishing data to {}'.format(self.host)
233 234 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
234 235
235 236 if self.zeromq is 1:
236 237 if self.verbose:
237 238 print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime)
238 239 self.zmq_socket.send_pyobj(self.dataOut)
239 240 self.dataOut.firstdata = False
240 241
241 242
242 243 def run(self, dataOut, **kwargs):
243 244 self.dataOut = dataOut
244 245 if not self.isConfig:
245 246 self.setup(**kwargs)
246 247 self.isConfig = True
247 248
248 249 self.publish_data()
249 250 time.sleep(self.delay)
250 251
251 252 def close(self):
252 253 if self.zeromq is 1:
253 254 self.dataOut.finished = True
254 255 self.zmq_socket.send_pyobj(self.dataOut)
255 256 self.zmq_socket.close()
256 257 if self.client:
257 258 self.client.loop_stop()
258 259 self.client.disconnect()
259 260
260 class ReceiverData(ProcessingUnit, Process):
261
262 class ReceiverData(ProcessingUnit):
263
264 def __init__(self, **kwargs):
265
266 ProcessingUnit.__init__(self, **kwargs)
267
268 self.isConfig = False
269 server = kwargs.get('server', 'zmq.pipe')
270 if 'tcp://' in server:
271 address = server
272 else:
273 address = 'ipc:///tmp/%s' % server
274
275 self.address = address
276 self.dataOut = JROData()
277
278 def setup(self):
279
280 self.context = zmq.Context()
281 self.receiver = self.context.socket(zmq.PULL)
282 self.receiver.bind(self.address)
283 time.sleep(0.5)
284 print '[Starting] ReceiverData from {}'.format(self.address)
285
286
287 def run(self):
288
289 if not self.isConfig:
290 self.setup()
291 self.isConfig = True
292
293 self.dataOut = self.receiver.recv_pyobj()
294 print '[Receiving] {} - {}'.format(self.dataOut.type,
295 self.dataOut.datatime.ctime())
296
297
298 class PlotterReceiver(ProcessingUnit, Process):
261 299
262 300 throttle_value = 5
263 301
264 302 def __init__(self, **kwargs):
265 303
266 304 ProcessingUnit.__init__(self, **kwargs)
267 305 Process.__init__(self)
268 306 self.mp = False
269 307 self.isConfig = False
270 308 self.isWebConfig = False
271 309 self.plottypes =[]
272 310 self.connections = 0
273 311 server = kwargs.get('server', 'zmq.pipe')
274 312 plot_server = kwargs.get('plot_server', 'zmq.web')
275 313 if 'tcp://' in server:
276 314 address = server
277 315 else:
278 316 address = 'ipc:///tmp/%s' % server
279 317
280 318 if 'tcp://' in plot_server:
281 319 plot_address = plot_server
282 320 else:
283 321 plot_address = 'ipc:///tmp/%s' % plot_server
284 322
285 323 self.address = address
286 324 self.plot_address = plot_address
287 325 self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')]
288 326 self.realtime = kwargs.get('realtime', False)
289 327 self.throttle_value = kwargs.get('throttle', 5)
290 328 self.sendData = self.initThrottle(self.throttle_value)
291 329 self.setup()
292 330
293 331 def setup(self):
294 332
295 333 self.data = {}
296 334 self.data['times'] = []
297 335 for plottype in self.plottypes:
298 336 self.data[plottype] = {}
299 337 self.data['noise'] = {}
300 338 self.data['throttle'] = self.throttle_value
301 339 self.data['ENDED'] = False
302 340 self.isConfig = True
303 341 self.data_web = {}
304 342
305 343 def event_monitor(self, monitor):
306 344
307 345 events = {}
308 346
309 347 for name in dir(zmq):
310 348 if name.startswith('EVENT_'):
311 349 value = getattr(zmq, name)
312 350 events[value] = name
313 351
314 352 while monitor.poll():
315 353 evt = recv_monitor_message(monitor)
316 354 if evt['event'] == 32:
317 355 self.connections += 1
318 356 if evt['event'] == 512:
319 357 pass
320 358 if self.connections == 0 and self.started is True:
321 359 self.ended = True
322 360
323 361 evt.update({'description': events[evt['event']]})
324 362
325 363 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
326 364 break
327 365 monitor.close()
328 366 print("event monitor thread done!")
329 367
330 368 def initThrottle(self, throttle_value):
331 369
332 370 @throttle(seconds=throttle_value)
333 371 def sendDataThrottled(fn_sender, data):
334 372 fn_sender(data)
335 373
336 374 return sendDataThrottled
337 375
338 376
339 377 def send(self, data):
340 378 # print '[sending] data=%s size=%s' % (data.keys(), len(data['times']))
341 379 self.sender.send_pyobj(data)
342 380
343 381
344 382 def update(self):
345 383 t = self.dataOut.utctime
346 384
347 385 if t in self.data['times']:
348 386 return
349 387
350 388 self.data['times'].append(t)
351 389 self.data['dataOut'] = self.dataOut
352 390
353 391 for plottype in self.plottypes:
354 392 if plottype == 'spc':
355 393 z = self.dataOut.data_spc/self.dataOut.normFactor
356 394 self.data[plottype] = 10*numpy.log10(z)
357 395 self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor)
358 396 if plottype == 'cspc':
359 397 jcoherence = self.dataOut.data_cspc/numpy.sqrt(self.dataOut.data_spc*self.dataOut.data_spc)
360 398 self.data['cspc_coh'] = numpy.abs(jcoherence)
361 399 self.data['cspc_phase'] = numpy.arctan2(jcoherence.imag, jcoherence.real)*180/numpy.pi
362 400 if plottype == 'rti':
363 401 self.data[plottype][t] = self.dataOut.getPower()
364 402 if plottype == 'snr':
365 403 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR)
366 404 if plottype == 'dop':
367 405 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP)
368 406 if plottype == 'mean':
369 407 self.data[plottype][t] = self.dataOut.data_MEAN
370 408 if plottype == 'std':
371 409 self.data[plottype][t] = self.dataOut.data_STD
372 410 if plottype == 'coh':
373 411 self.data[plottype][t] = self.dataOut.getCoherence()
374 412 if plottype == 'phase':
375 413 self.data[plottype][t] = self.dataOut.getCoherence(phase=True)
376 414 if plottype == 'output':
377 415 self.data[plottype][t] = self.dataOut.data_output
378 416 if plottype == 'param':
379 417 self.data[plottype][t] = self.dataOut.data_param
380 418 if self.realtime:
381 419 self.data_web['timestamp'] = t
382 420 if plottype == 'spc':
383 421 self.data_web[plottype] = roundFloats(decimate(self.data[plottype]).tolist())
384 422 elif plottype == 'cspc':
385 423 self.data_web['cspc_coh'] = roundFloats(decimate(self.data['cspc_coh']).tolist())
386 424 self.data_web['cspc_phase'] = roundFloats(decimate(self.data['cspc_phase']).tolist())
387 425 elif plottype == 'noise':
388 426 self.data_web['noise'] = roundFloats(self.data['noise'][t].tolist())
389 427 else:
390 428 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
391 429 self.data_web['interval'] = self.dataOut.getTimeInterval()
392 430 self.data_web['type'] = plottype
393 431
394 432 def run(self):
395 433
396 434 print '[Starting] {} from {}'.format(self.name, self.address)
397 435
398 436 self.context = zmq.Context()
399 437 self.receiver = self.context.socket(zmq.PULL)
400 438 self.receiver.bind(self.address)
401 439 monitor = self.receiver.get_monitor_socket()
402 440 self.sender = self.context.socket(zmq.PUB)
403 441 if self.realtime:
404 442 self.sender_web = self.context.socket(zmq.PUB)
405 443 self.sender_web.connect(self.plot_address)
406 444 time.sleep(1)
407 445
408 446 if 'server' in self.kwargs:
409 447 self.sender.bind("ipc:///tmp/{}.plots".format(self.kwargs['server']))
410 448 else:
411 449 self.sender.bind("ipc:///tmp/zmq.plots")
412 450
413 451 time.sleep(3)
414 452
415 453 t = Thread(target=self.event_monitor, args=(monitor,))
416 454 t.start()
417 455
418 456 while True:
419 457 self.dataOut = self.receiver.recv_pyobj()
420 458 # print '[Receiving] {} - {}'.format(self.dataOut.type,
421 459 # self.dataOut.datatime.ctime())
422 460
423 461 self.update()
424 462
425 463 if self.dataOut.firstdata is True:
426 464 self.data['STARTED'] = True
427 465
428 466 if self.dataOut.finished is True:
429 467 self.send(self.data)
430 468 self.connections -= 1
431 469 if self.connections == 0 and self.started:
432 470 self.ended = True
433 471 self.data['ENDED'] = True
434 472 self.send(self.data)
435 473 self.setup()
436 474 self.started = False
437 475 else:
438 476 if self.realtime:
439 477 self.send(self.data)
440 478 self.sender_web.send_string(json.dumps(self.data_web))
441 479 else:
442 480 self.sendData(self.send, self.data)
443 481 self.started = True
444 482
445 483 self.data['STARTED'] = False
446 484 return
447 485
448 486 def sendToWeb(self):
449 487
450 488 if not self.isWebConfig:
451 489 context = zmq.Context()
452 490 sender_web_config = context.socket(zmq.PUB)
453 491 if 'tcp://' in self.plot_address:
454 492 dum, address, port = self.plot_address.split(':')
455 493 conf_address = '{}:{}:{}'.format(dum, address, int(port)+1)
456 494 else:
457 495 conf_address = self.plot_address + '.config'
458 496 sender_web_config.bind(conf_address)
459 497 time.sleep(1)
460 498 for kwargs in self.operationKwargs.values():
461 499 if 'plot' in kwargs:
462 500 print '[Sending] Config data to web for {}'.format(kwargs['code'].upper())
463 501 sender_web_config.send_string(json.dumps(kwargs))
464 502 self.isWebConfig = True
General Comments 0
You need to be logged in to leave comments. Login now