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