##// END OF EJS Templates
schain
RSS
DP+LP Join Spc+Voltage
rflores -
r1549:c70288c93a07
parent child
Show More
Side by Side Unified
Context file:
r1549:c70288c93a07 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,1079 +1,1079
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Definition of diferent Data objects for different types of data
6 6
7 7 Here you will find the diferent data objects for the different types
8 8 of data, this data objects must be used as dataIn or dataOut objects in
9 9 processing units and operations. Currently the supported data objects are:
10 10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 11 """
12 12
13 13 import copy
14 14 import numpy
15 15 import datetime
16 16 import json
17 17
18 18 import schainpy.admin
19 19 from schainpy.utils import log
20 20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 21 from schainpy.model.data import _noise
22 22
23 23
24 24 def getNumpyDtype(dataTypeCode):
25 25
26 26 if dataTypeCode == 0:
27 27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 28 elif dataTypeCode == 1:
29 29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 30 elif dataTypeCode == 2:
31 31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 32 elif dataTypeCode == 3:
33 33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 34 elif dataTypeCode == 4:
35 35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 36 elif dataTypeCode == 5:
37 37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 38 else:
39 39 raise ValueError('dataTypeCode was not defined')
40 40
41 41 return numpyDtype
42 42
43 43
44 44 def getDataTypeCode(numpyDtype):
45 45
46 46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 47 datatype = 0
48 48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 49 datatype = 1
50 50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 51 datatype = 2
52 52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 53 datatype = 3
54 54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 55 datatype = 4
56 56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 57 datatype = 5
58 58 else:
59 59 datatype = None
60 60
61 61 return datatype
62 62
63 63
64 64 def hildebrand_sekhon(data, navg):
65 65 """
66 66 This method is for the objective determination of the noise level in Doppler spectra. This
67 67 implementation technique is based on the fact that the standard deviation of the spectral
68 68 densities is equal to the mean spectral density for white Gaussian noise
69 69
70 70 Inputs:
71 71 Data : heights
72 72 navg : numbers of averages
73 73
74 74 Return:
75 75 mean : noise's level
76 76 """
77 77
78 78 sortdata = numpy.sort(data, axis=None)
79 79 #print(numpy.shape(data))
80 80 #exit()
81 81 '''
82 82 lenOfData = len(sortdata)
83 83 nums_min = lenOfData*0.2
84 84
85 85 if nums_min <= 5:
86 86
87 87 nums_min = 5
88 88
89 89 sump = 0.
90 90 sumq = 0.
91 91
92 92 j = 0
93 93 cont = 1
94 94
95 95 while((cont == 1)and(j < lenOfData)):
96 96
97 97 sump += sortdata[j]
98 98 sumq += sortdata[j]**2
99 99
100 100 if j > nums_min:
101 101 rtest = float(j)/(j-1) + 1.0/navg
102 102 if ((sumq*j) > (rtest*sump**2)):
103 103 j = j - 1
104 104 sump = sump - sortdata[j]
105 105 sumq = sumq - sortdata[j]**2
106 106 cont = 0
107 107
108 108 j += 1
109 109
110 110 lnoise = sump / j
111 111
112 112 return lnoise
113 113 '''
114 114 return _noise.hildebrand_sekhon(sortdata, navg)
115 115
116 116
117 117 class Beam:
118 118
119 119 def __init__(self):
120 120 self.codeList = []
121 121 self.azimuthList = []
122 122 self.zenithList = []
123 123
124 124
125 125 class GenericData(object):
126 126
127 127 flagNoData = True
128 128
129 129 def copy(self, inputObj=None):
130 130
131 131 if inputObj == None:
132 132 return copy.deepcopy(self)
133 133
134 134 for key in list(inputObj.__dict__.keys()):
135 135
136 136 attribute = inputObj.__dict__[key]
137 137
138 138 # If this attribute is a tuple or list
139 139 if type(inputObj.__dict__[key]) in (tuple, list):
140 140 self.__dict__[key] = attribute[:]
141 141 continue
142 142
143 143 # If this attribute is another object or instance
144 144 if hasattr(attribute, '__dict__'):
145 145 self.__dict__[key] = attribute.copy()
146 146 continue
147 147
148 148 self.__dict__[key] = inputObj.__dict__[key]
149 149
150 150 def deepcopy(self):
151 151
152 152 return copy.deepcopy(self)
153 153
154 154 def isEmpty(self):
155 155
156 156 return self.flagNoData
157 157
158 158 def isReady(self):
159 159
160 160 return not self.flagNoData
161 161
162 162
163 163 class JROData(GenericData):
164 164
165 165 systemHeaderObj = SystemHeader()
166 166 radarControllerHeaderObj = RadarControllerHeader()
167 167 type = None
168 168 datatype = None # dtype but in string
169 169 nProfiles = None
170 170 heightList = None
171 171 channelList = None
172 172 flagDiscontinuousBlock = False
173 173 useLocalTime = False
174 174 utctime = None
175 175 timeZone = None
176 176 dstFlag = None
177 177 errorCount = None
178 178 blocksize = None
179 179 flagDecodeData = False # asumo q la data no esta decodificada
180 180 flagDeflipData = False # asumo q la data no esta sin flip
181 181 flagShiftFFT = False
182 182 nCohInt = None
183 183 windowOfFilter = 1
184 184 C = 3e8
185 185 frequency = 49.92e6
186 186 realtime = False
187 187 beacon_heiIndexList = None
188 188 last_block = None
189 189 blocknow = None
190 190 azimuth = None
191 191 zenith = None
192 192 beam = Beam()
193 193 profileIndex = None
194 194 error = None
195 195 data = None
196 196 nmodes = None
197 197 metadata_list = ['heightList', 'timeZone', 'type']
198 198
199 199 def __str__(self):
200 200
201 201 return '{} - {}'.format(self.type, self.datatime())
202 202
203 203 def getNoise(self):
204 204
205 205 raise NotImplementedError
206 206
207 207 @property
208 208 def nChannels(self):
209 209
210 210 return len(self.channelList)
211 211
212 212 @property
213 213 def channelIndexList(self):
214 214
215 215 return list(range(self.nChannels))
216 216
217 217 @property
218 218 def nHeights(self):
219 219
220 220 return len(self.heightList)
221 221
222 222 def getDeltaH(self):
223 223
224 224 return self.heightList[1] - self.heightList[0]
225 225
226 226 @property
227 227 def ltctime(self):
228 228
229 229 if self.useLocalTime:
230 230 return self.utctime - self.timeZone * 60
231 231
232 232 return self.utctime
233 233
234 234 @property
235 235 def datatime(self):
236 236
237 237 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
238 238 return datatimeValue
239 239
240 240 def getTimeRange(self):
241 241
242 242 datatime = []
243 243
244 244 datatime.append(self.ltctime)
245 245 datatime.append(self.ltctime + self.timeInterval + 1)
246 246
247 247 datatime = numpy.array(datatime)
248 248
249 249 return datatime
250 250
251 251 def getFmaxTimeResponse(self):
252 252
253 253 period = (10**-6) * self.getDeltaH() / (0.15)
254 254
255 255 PRF = 1. / (period * self.nCohInt)
256 256
257 257 fmax = PRF
258 258
259 259 return fmax
260 260
261 261 def getFmax(self):
262 262 PRF = 1. / (self.ippSeconds * self.nCohInt)
263
263 #print("ippsec",self.ippSeconds)
264 264 fmax = PRF
265 265 return fmax
266 266
267 267 def getVmax(self):
268 268
269 269 _lambda = self.C / self.frequency
270 270
271 271 vmax = self.getFmax() * _lambda / 2
272 272
273 273 return vmax
274 274
275 275 @property
276 276 def ippSeconds(self):
277 277 '''
278 278 '''
279 279 return self.radarControllerHeaderObj.ippSeconds
280 280
281 281 @ippSeconds.setter
282 282 def ippSeconds(self, ippSeconds):
283 283 '''
284 284 '''
285 285 self.radarControllerHeaderObj.ippSeconds = ippSeconds
286 286
287 287 @property
288 288 def code(self):
289 289 '''
290 290 '''
291 291 return self.radarControllerHeaderObj.code
292 292
293 293 @code.setter
294 294 def code(self, code):
295 295 '''
296 296 '''
297 297 self.radarControllerHeaderObj.code = code
298 298
299 299 @property
300 300 def nCode(self):
301 301 '''
302 302 '''
303 303 return self.radarControllerHeaderObj.nCode
304 304
305 305 @nCode.setter
306 306 def nCode(self, ncode):
307 307 '''
308 308 '''
309 309 self.radarControllerHeaderObj.nCode = ncode
310 310
311 311 @property
312 312 def nBaud(self):
313 313 '''
314 314 '''
315 315 return self.radarControllerHeaderObj.nBaud
316 316
317 317 @nBaud.setter
318 318 def nBaud(self, nbaud):
319 319 '''
320 320 '''
321 321 self.radarControllerHeaderObj.nBaud = nbaud
322 322
323 323 @property
324 324 def ipp(self):
325 325 '''
326 326 '''
327 327 return self.radarControllerHeaderObj.ipp
328 328
329 329 @ipp.setter
330 330 def ipp(self, ipp):
331 331 '''
332 332 '''
333 333 self.radarControllerHeaderObj.ipp = ipp
334 334
335 335 @property
336 336 def metadata(self):
337 337 '''
338 338 '''
339 339
340 340 return {attr: getattr(self, attr) for attr in self.metadata_list}
341 341
342 342
343 343 class Voltage(JROData):
344 344
345 345 dataPP_POW = None
346 346 dataPP_DOP = None
347 347 dataPP_WIDTH = None
348 348 dataPP_SNR = None
349 349
350 350 def __init__(self):
351 351 '''
352 352 Constructor
353 353 '''
354 354
355 355 self.useLocalTime = True
356 356 self.radarControllerHeaderObj = RadarControllerHeader()
357 357 self.systemHeaderObj = SystemHeader()
358 358 self.type = "Voltage"
359 359 self.data = None
360 360 self.nProfiles = None
361 361 self.heightList = None
362 362 self.channelList = None
363 363 self.flagNoData = True
364 364 self.flagDiscontinuousBlock = False
365 365 self.utctime = None
366 366 self.timeZone = 0
367 367 self.dstFlag = None
368 368 self.errorCount = None
369 369 self.nCohInt = None
370 370 self.blocksize = None
371 371 self.flagCohInt = False
372 372 self.flagDecodeData = False # asumo q la data no esta decodificada
373 373 self.flagDeflipData = False # asumo q la data no esta sin flip
374 374 self.flagShiftFFT = False
375 375 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
376 376 self.profileIndex = 0
377 377 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
378 378 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
379 379
380 380 def getNoisebyHildebrand(self, channel=None, Profmin_index=None, Profmax_index=None):
381 381 """
382 382 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
383 383
384 384 Return:
385 385 noiselevel
386 386 """
387 387
388 388 if channel != None:
389 389 data = self.data[channel]
390 390 nChannels = 1
391 391 else:
392 392 data = self.data
393 393 nChannels = self.nChannels
394 394
395 395 noise = numpy.zeros(nChannels)
396 396 power = data * numpy.conjugate(data)
397 397
398 398 for thisChannel in range(nChannels):
399 399 if nChannels == 1:
400 400 daux = power[:].real
401 401 else:
402 402 #print(power.shape)
403 403 daux = power[thisChannel, Profmin_index:Profmax_index, :].real
404 404 #print(daux.shape)
405 405 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
406 406
407 407 return noise
408 408
409 409 def getNoise(self, type=1, channel=None, Profmin_index=None, Profmax_index=None):
410 410
411 411 if type == 1:
412 412 noise = self.getNoisebyHildebrand(channel, Profmin_index, Profmax_index)
413 413
414 414 return noise
415 415
416 416 def getPower(self, channel=None):
417 417
418 418 if channel != None:
419 419 data = self.data[channel]
420 420 else:
421 421 data = self.data
422 422
423 423 power = data * numpy.conjugate(data)
424 424 powerdB = 10 * numpy.log10(power.real)
425 425 powerdB = numpy.squeeze(powerdB)
426 426
427 427 return powerdB
428 428
429 429 @property
430 430 def timeInterval(self):
431 431
432 432 return self.ippSeconds * self.nCohInt
433 433
434 434 noise = property(getNoise, "I'm the 'nHeights' property.")
435 435
436 436
437 437 class Spectra(JROData):
438 438
439 439 def __init__(self):
440 440 '''
441 441 Constructor
442 442 '''
443 443
444 444 self.data_dc = None
445 445 self.data_spc = None
446 446 self.data_cspc = None
447 447 self.useLocalTime = True
448 448 self.radarControllerHeaderObj = RadarControllerHeader()
449 449 self.systemHeaderObj = SystemHeader()
450 450 self.type = "Spectra"
451 451 self.timeZone = 0
452 452 self.nProfiles = None
453 453 self.heightList = None
454 454 self.channelList = None
455 455 self.pairsList = None
456 456 self.flagNoData = True
457 457 self.flagDiscontinuousBlock = False
458 458 self.utctime = None
459 459 self.nCohInt = None
460 460 self.nIncohInt = None
461 461 self.blocksize = None
462 462 self.nFFTPoints = None
463 463 self.wavelength = None
464 464 self.flagDecodeData = False # asumo q la data no esta decodificada
465 465 self.flagDeflipData = False # asumo q la data no esta sin flip
466 466 self.flagShiftFFT = False
467 467 self.ippFactor = 1
468 468 self.beacon_heiIndexList = []
469 469 self.noise_estimation = None
470 470 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
471 471 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
472 472
473 473 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
474 474 """
475 475 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
476 476
477 477 Return:
478 478 noiselevel
479 479 """
480 480
481 481 noise = numpy.zeros(self.nChannels)
482 482
483 483 for channel in range(self.nChannels):
484 484 #print(self.data_spc[0])
485 485 #exit(1)
486 486 daux = self.data_spc[channel,
487 487 xmin_index:xmax_index, ymin_index:ymax_index]
488 488 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
489 489
490 490 return noise
491 491
492 492 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
493 493
494 494 if self.noise_estimation is not None:
495 495 # this was estimated by getNoise Operation defined in jroproc_spectra.py
496 496 return self.noise_estimation
497 497 else:
498 498
499 499 noise = self.getNoisebyHildebrand(
500 500 xmin_index, xmax_index, ymin_index, ymax_index)
501 501 return noise
502 502
503 503 def getFreqRangeTimeResponse(self, extrapoints=0):
504 504
505 505 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
506 506 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
507 507
508 508 return freqrange
509 509
510 510 def getAcfRange(self, extrapoints=0):
511 511
512 512 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
513 513 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
514 514
515 515 return freqrange
516 516
517 517 def getFreqRange(self, extrapoints=0):
518 518
519 519 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
520 520 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
521 521
522 522 return freqrange
523 523
524 524 def getVelRange(self, extrapoints=0):
525 525
526 526 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
527 527 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
528 528
529 529 if self.nmodes:
530 530 return velrange/self.nmodes
531 531 else:
532 532 return velrange
533 533
534 534 @property
535 535 def nPairs(self):
536 536
537 537 return len(self.pairsList)
538 538
539 539 @property
540 540 def pairsIndexList(self):
541 541
542 542 return list(range(self.nPairs))
543 543
544 544 @property
545 545 def normFactor(self):
546 546
547 547 pwcode = 1
548 548
549 549 if self.flagDecodeData:
550 550 pwcode = numpy.sum(self.code[0]**2)
551 551 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
552 552 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
553 553
554 554 return normFactor
555 555
556 556 @property
557 557 def flag_cspc(self):
558 558
559 559 if self.data_cspc is None:
560 560 return True
561 561
562 562 return False
563 563
564 564 @property
565 565 def flag_dc(self):
566 566
567 567 if self.data_dc is None:
568 568 return True
569 569
570 570 return False
571 571
572 572 @property
573 573 def timeInterval(self):
574 574
575 575 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
576 576 if self.nmodes:
577 577 return self.nmodes*timeInterval
578 578 else:
579 579 return timeInterval
580 580
581 581 def getPower(self):
582 582
583 583 factor = self.normFactor
584 584 z = self.data_spc / factor
585 585 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
586 586 avg = numpy.average(z, axis=1)
587 587
588 588 return 10 * numpy.log10(avg)
589 589
590 590 def getCoherence(self, pairsList=None, phase=False):
591 591
592 592 z = []
593 593 if pairsList is None:
594 594 pairsIndexList = self.pairsIndexList
595 595 else:
596 596 pairsIndexList = []
597 597 for pair in pairsList:
598 598 if pair not in self.pairsList:
599 599 raise ValueError("Pair %s is not in dataOut.pairsList" % (
600 600 pair))
601 601 pairsIndexList.append(self.pairsList.index(pair))
602 602 for i in range(len(pairsIndexList)):
603 603 pair = self.pairsList[pairsIndexList[i]]
604 604 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
605 605 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
606 606 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
607 607 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
608 608 if phase:
609 609 data = numpy.arctan2(avgcoherenceComplex.imag,
610 610 avgcoherenceComplex.real) * 180 / numpy.pi
611 611 else:
612 612 data = numpy.abs(avgcoherenceComplex)
613 613
614 614 z.append(data)
615 615
616 616 return numpy.array(z)
617 617
618 618 def setValue(self, value):
619 619
620 620 print("This property should not be initialized")
621 621
622 622 return
623 623
624 624 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
625 625
626 626
627 627 class SpectraHeis(Spectra):
628 628
629 629 def __init__(self):
630 630
631 631 self.radarControllerHeaderObj = RadarControllerHeader()
632 632 self.systemHeaderObj = SystemHeader()
633 633 self.type = "SpectraHeis"
634 634 self.nProfiles = None
635 635 self.heightList = None
636 636 self.channelList = None
637 637 self.flagNoData = True
638 638 self.flagDiscontinuousBlock = False
639 639 self.utctime = None
640 640 self.blocksize = None
641 641 self.profileIndex = 0
642 642 self.nCohInt = 1
643 643 self.nIncohInt = 1
644 644
645 645 @property
646 646 def normFactor(self):
647 647 pwcode = 1
648 648 if self.flagDecodeData:
649 649 pwcode = numpy.sum(self.code[0]**2)
650 650
651 651 normFactor = self.nIncohInt * self.nCohInt * pwcode
652 652
653 653 return normFactor
654 654
655 655 @property
656 656 def timeInterval(self):
657 657
658 658 return self.ippSeconds * self.nCohInt * self.nIncohInt
659 659
660 660
661 661 class Fits(JROData):
662 662
663 663 def __init__(self):
664 664
665 665 self.type = "Fits"
666 666 self.nProfiles = None
667 667 self.heightList = None
668 668 self.channelList = None
669 669 self.flagNoData = True
670 670 self.utctime = None
671 671 self.nCohInt = 1
672 672 self.nIncohInt = 1
673 673 self.useLocalTime = True
674 674 self.profileIndex = 0
675 675 self.timeZone = 0
676 676
677 677 def getTimeRange(self):
678 678
679 679 datatime = []
680 680
681 681 datatime.append(self.ltctime)
682 682 datatime.append(self.ltctime + self.timeInterval)
683 683
684 684 datatime = numpy.array(datatime)
685 685
686 686 return datatime
687 687
688 688 def getChannelIndexList(self):
689 689
690 690 return list(range(self.nChannels))
691 691
692 692 def getNoise(self, type=1):
693 693
694 694
695 695 if type == 1:
696 696 noise = self.getNoisebyHildebrand()
697 697
698 698 if type == 2:
699 699 noise = self.getNoisebySort()
700 700
701 701 if type == 3:
702 702 noise = self.getNoisebyWindow()
703 703
704 704 return noise
705 705
706 706 @property
707 707 def timeInterval(self):
708 708
709 709 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
710 710
711 711 return timeInterval
712 712
713 713 @property
714 714 def ippSeconds(self):
715 715 '''
716 716 '''
717 717 return self.ipp_sec
718 718
719 719 noise = property(getNoise, "I'm the 'nHeights' property.")
720 720
721 721
722 722 class Correlation(JROData):
723 723
724 724 def __init__(self):
725 725 '''
726 726 Constructor
727 727 '''
728 728 self.radarControllerHeaderObj = RadarControllerHeader()
729 729 self.systemHeaderObj = SystemHeader()
730 730 self.type = "Correlation"
731 731 self.data = None
732 732 self.dtype = None
733 733 self.nProfiles = None
734 734 self.heightList = None
735 735 self.channelList = None
736 736 self.flagNoData = True
737 737 self.flagDiscontinuousBlock = False
738 738 self.utctime = None
739 739 self.timeZone = 0
740 740 self.dstFlag = None
741 741 self.errorCount = None
742 742 self.blocksize = None
743 743 self.flagDecodeData = False # asumo q la data no esta decodificada
744 744 self.flagDeflipData = False # asumo q la data no esta sin flip
745 745 self.pairsList = None
746 746 self.nPoints = None
747 747
748 748 def getPairsList(self):
749 749
750 750 return self.pairsList
751 751
752 752 def getNoise(self, mode=2):
753 753
754 754 indR = numpy.where(self.lagR == 0)[0][0]
755 755 indT = numpy.where(self.lagT == 0)[0][0]
756 756
757 757 jspectra0 = self.data_corr[:, :, indR, :]
758 758 jspectra = copy.copy(jspectra0)
759 759
760 760 num_chan = jspectra.shape[0]
761 761 num_hei = jspectra.shape[2]
762 762
763 763 freq_dc = jspectra.shape[1] / 2
764 764 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
765 765
766 766 if ind_vel[0] < 0:
767 767 ind_vel[list(range(0, 1))] = ind_vel[list(
768 768 range(0, 1))] + self.num_prof
769 769
770 770 if mode == 1:
771 771 jspectra[:, freq_dc, :] = (
772 772 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
773 773
774 774 if mode == 2:
775 775
776 776 vel = numpy.array([-2, -1, 1, 2])
777 777 xx = numpy.zeros([4, 4])
778 778
779 779 for fil in range(4):
780 780 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
781 781
782 782 xx_inv = numpy.linalg.inv(xx)
783 783 xx_aux = xx_inv[0, :]
784 784
785 785 for ich in range(num_chan):
786 786 yy = jspectra[ich, ind_vel, :]
787 787 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
788 788
789 789 junkid = jspectra[ich, freq_dc, :] <= 0
790 790 cjunkid = sum(junkid)
791 791
792 792 if cjunkid.any():
793 793 jspectra[ich, freq_dc, junkid.nonzero()] = (
794 794 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
795 795
796 796 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
797 797
798 798 return noise
799 799
800 800 @property
801 801 def timeInterval(self):
802 802
803 803 return self.ippSeconds * self.nCohInt * self.nProfiles
804 804
805 805 def splitFunctions(self):
806 806
807 807 pairsList = self.pairsList
808 808 ccf_pairs = []
809 809 acf_pairs = []
810 810 ccf_ind = []
811 811 acf_ind = []
812 812 for l in range(len(pairsList)):
813 813 chan0 = pairsList[l][0]
814 814 chan1 = pairsList[l][1]
815 815
816 816 # Obteniendo pares de Autocorrelacion
817 817 if chan0 == chan1:
818 818 acf_pairs.append(chan0)
819 819 acf_ind.append(l)
820 820 else:
821 821 ccf_pairs.append(pairsList[l])
822 822 ccf_ind.append(l)
823 823
824 824 data_acf = self.data_cf[acf_ind]
825 825 data_ccf = self.data_cf[ccf_ind]
826 826
827 827 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
828 828
829 829 @property
830 830 def normFactor(self):
831 831 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
832 832 acf_pairs = numpy.array(acf_pairs)
833 833 normFactor = numpy.zeros((self.nPairs, self.nHeights))
834 834
835 835 for p in range(self.nPairs):
836 836 pair = self.pairsList[p]
837 837
838 838 ch0 = pair[0]
839 839 ch1 = pair[1]
840 840
841 841 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
842 842 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
843 843 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
844 844
845 845 return normFactor
846 846
847 847
848 848 class Parameters(Spectra):
849 849
850 850 groupList = None # List of Pairs, Groups, etc
851 851 data_param = None # Parameters obtained
852 852 data_pre = None # Data Pre Parametrization
853 853 data_SNR = None # Signal to Noise Ratio
854 854 abscissaList = None # Abscissa, can be velocities, lags or time
855 855 utctimeInit = None # Initial UTC time
856 856 paramInterval = None # Time interval to calculate Parameters in seconds
857 857 useLocalTime = True
858 858 # Fitting
859 859 data_error = None # Error of the estimation
860 860 constants = None
861 861 library = None
862 862 # Output signal
863 863 outputInterval = None # Time interval to calculate output signal in seconds
864 864 data_output = None # Out signal
865 865 nAvg = None
866 866 noise_estimation = None
867 867 GauSPC = None # Fit gaussian SPC
868 868
869 869 def __init__(self):
870 870 '''
871 871 Constructor
872 872 '''
873 873 self.radarControllerHeaderObj = RadarControllerHeader()
874 874 self.systemHeaderObj = SystemHeader()
875 875 self.type = "Parameters"
876 876 self.timeZone = 0
877 877
878 878 def getTimeRange1(self, interval):
879 879
880 880 datatime = []
881 881
882 882 if self.useLocalTime:
883 883 time1 = self.utctimeInit - self.timeZone * 60
884 884 else:
885 885 time1 = self.utctimeInit
886 886
887 887 datatime.append(time1)
888 888 datatime.append(time1 + interval)
889 889 datatime = numpy.array(datatime)
890 890
891 891 return datatime
892 892
893 893 @property
894 894 def timeInterval(self):
895 895
896 896 if hasattr(self, 'timeInterval1'):
897 897 return self.timeInterval1
898 898 else:
899 899 return self.paramInterval
900 900
901 901
902 902 def setValue(self, value):
903 903
904 904 print("This property should not be initialized")
905 905
906 906 return
907 907
908 908 def getNoise(self):
909 909
910 910 return self.spc_noise
911 911
912 912 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
913 913
914 914
915 915 class PlotterData(object):
916 916 '''
917 917 Object to hold data to be plotted
918 918 '''
919 919
920 920 MAXNUMX = 200
921 921 MAXNUMY = 200
922 922
923 923 def __init__(self, code, exp_code, localtime=True):
924 924
925 925 self.key = code
926 926 self.exp_code = exp_code
927 927 self.ready = False
928 928 self.flagNoData = False
929 929 self.localtime = localtime
930 930 self.data = {}
931 931 self.meta = {}
932 932 self.__heights = []
933 933
934 934 def __str__(self):
935 935 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
936 936 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
937 937
938 938 def __len__(self):
939 939 return len(self.data)
940 940
941 941 def __getitem__(self, key):
942 942 if isinstance(key, int):
943 943 return self.data[self.times[key]]
944 944 elif isinstance(key, str):
945 945 ret = numpy.array([self.data[x][key] for x in self.times])
946 946 if ret.ndim > 1:
947 947 ret = numpy.swapaxes(ret, 0, 1)
948 948 return ret
949 949
950 950 def __contains__(self, key):
951 951 return key in self.data[self.min_time]
952 952
953 953 def setup(self):
954 954 '''
955 955 Configure object
956 956 '''
957 957 self.type = ''
958 958 self.ready = False
959 959 del self.data
960 960 self.data = {}
961 961 self.__heights = []
962 962 self.__all_heights = set()
963 963
964 964 def shape(self, key):
965 965 '''
966 966 Get the shape of the one-element data for the given key
967 967 '''
968 968
969 969 if len(self.data[self.min_time][key]):
970 970 return self.data[self.min_time][key].shape
971 971 return (0,)
972 972
973 973 def update(self, data, tm, meta={}):
974 974 '''
975 975 Update data object with new dataOut
976 976 '''
977 977
978 978 self.data[tm] = data
979 979
980 980 for key, value in meta.items():
981 981 setattr(self, key, value)
982 982
983 983 def normalize_heights(self):
984 984 '''
985 985 Ensure same-dimension of the data for different heighList
986 986 '''
987 987
988 988 H = numpy.array(list(self.__all_heights))
989 989 H.sort()
990 990 for key in self.data:
991 991 shape = self.shape(key)[:-1] + H.shape
992 992 for tm, obj in list(self.data[key].items()):
993 993 h = self.__heights[self.times.tolist().index(tm)]
994 994 if H.size == h.size:
995 995 continue
996 996 index = numpy.where(numpy.in1d(H, h))[0]
997 997 dummy = numpy.zeros(shape) + numpy.nan
998 998 if len(shape) == 2:
999 999 dummy[:, index] = obj
1000 1000 else:
1001 1001 dummy[index] = obj
1002 1002 self.data[key][tm] = dummy
1003 1003
1004 1004 self.__heights = [H for tm in self.times]
1005 1005
1006 1006 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1007 1007 '''
1008 1008 Convert data to json
1009 1009 '''
1010 1010
1011 1011 meta = {}
1012 1012 meta['xrange'] = []
1013 1013 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1014 1014 tmp = self.data[tm][self.key]
1015 1015 shape = tmp.shape
1016 1016 if len(shape) == 2:
1017 1017 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1018 1018 elif len(shape) == 3:
1019 1019 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1020 1020 data = self.roundFloats(
1021 1021 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1022 1022 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1023 1023 else:
1024 1024 data = self.roundFloats(self.data[tm][self.key].tolist())
1025 1025
1026 1026 ret = {
1027 1027 'plot': plot_name,
1028 1028 'code': self.exp_code,
1029 1029 'time': float(tm),
1030 1030 'data': data,
1031 1031 }
1032 1032 meta['type'] = plot_type
1033 1033 meta['interval'] = float(self.interval)
1034 1034 meta['localtime'] = self.localtime
1035 1035 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1036 1036 meta.update(self.meta)
1037 1037 ret['metadata'] = meta
1038 1038 return json.dumps(ret)
1039 1039
1040 1040 @property
1041 1041 def times(self):
1042 1042 '''
1043 1043 Return the list of times of the current data
1044 1044 '''
1045 1045
1046 1046 ret = [t for t in self.data]
1047 1047 ret.sort()
1048 1048 return numpy.array(ret)
1049 1049
1050 1050 @property
1051 1051 def min_time(self):
1052 1052 '''
1053 1053 Return the minimun time value
1054 1054 '''
1055 1055
1056 1056 return self.times[0]
1057 1057
1058 1058 @property
1059 1059 def max_time(self):
1060 1060 '''
1061 1061 Return the maximun time value
1062 1062 '''
1063 1063
1064 1064 return self.times[-1]
1065 1065
1066 1066 # @property
1067 1067 # def heights(self):
1068 1068 # '''
1069 1069 # Return the list of heights of the current data
1070 1070 # '''
1071 1071
1072 1072 # return numpy.array(self.__heights[-1])
1073 1073
1074 1074 @staticmethod
1075 1075 def roundFloats(obj):
1076 1076 if isinstance(obj, list):
1077 1077 return list(map(PlotterData.roundFloats, obj))
1078 1078 elif isinstance(obj, float):
1079 1079 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,697 +1,698
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Base class to create plot operations
6 6
7 7 """
8 8
9 9 import os
10 10 import sys
11 11 import zmq
12 12 import time
13 13 import numpy
14 14 import datetime
15 15 from collections import deque
16 16 from functools import wraps
17 17 from threading import Thread
18 18 import matplotlib
19 19
20 20 if 'BACKEND' in os.environ:
21 21 matplotlib.use(os.environ['BACKEND'])
22 22 elif 'linux' in sys.platform:
23 23 matplotlib.use("TkAgg")
24 24 elif 'darwin' in sys.platform:
25 25 matplotlib.use('MacOSX')
26 26 else:
27 27 from schainpy.utils import log
28 28 log.warning('Using default Backend="Agg"', 'INFO')
29 29 matplotlib.use('Agg')
30 30
31 31 import matplotlib.pyplot as plt
32 32 from matplotlib.patches import Polygon
33 33 from mpl_toolkits.axes_grid1 import make_axes_locatable
34 34 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
35 35
36 36 from schainpy.model.data.jrodata import PlotterData
37 37 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
38 38 from schainpy.utils import log
39 39
40 40 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
41 41 blu_values = matplotlib.pyplot.get_cmap(
42 42 'seismic_r', 20)(numpy.arange(20))[10:15]
43 43 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
44 44 'jro', numpy.vstack((blu_values, jet_values)))
45 45 matplotlib.pyplot.register_cmap(cmap=ncmap)
46 46
47 47 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
48 48 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
49 49
50 50 EARTH_RADIUS = 6.3710e3
51 51
52 52 def ll2xy(lat1, lon1, lat2, lon2):
53 53
54 54 p = 0.017453292519943295
55 55 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
56 56 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
57 57 r = 12742 * numpy.arcsin(numpy.sqrt(a))
58 58 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
59 59 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
60 60 theta = -theta + numpy.pi/2
61 61 return r*numpy.cos(theta), r*numpy.sin(theta)
62 62
63 63
64 64 def km2deg(km):
65 65 '''
66 66 Convert distance in km to degrees
67 67 '''
68 68
69 69 return numpy.rad2deg(km/EARTH_RADIUS)
70 70
71 71
72 72 def figpause(interval):
73 73 backend = plt.rcParams['backend']
74 74 if backend in matplotlib.rcsetup.interactive_bk:
75 75 figManager = matplotlib._pylab_helpers.Gcf.get_active()
76 76 if figManager is not None:
77 77 canvas = figManager.canvas
78 78 if canvas.figure.stale:
79 79 canvas.draw()
80 80 try:
81 81 canvas.start_event_loop(interval)
82 82 except:
83 83 pass
84 84 return
85 85
86 86 def popup(message):
87 87 '''
88 88 '''
89 89
90 90 fig = plt.figure(figsize=(12, 8), facecolor='r')
91 91 text = '\n'.join([s.strip() for s in message.split(':')])
92 92 fig.text(0.01, 0.5, text, ha='left', va='center',
93 93 size='20', weight='heavy', color='w')
94 94 fig.show()
95 95 figpause(1000)
96 96
97 97
98 98 class Throttle(object):
99 99 '''
100 100 Decorator that prevents a function from being called more than once every
101 101 time period.
102 102 To create a function that cannot be called more than once a minute, but
103 103 will sleep until it can be called:
104 104 @Throttle(minutes=1)
105 105 def foo():
106 106 pass
107 107
108 108 for i in range(10):
109 109 foo()
110 110 print "This function has run %s times." % i
111 111 '''
112 112
113 113 def __init__(self, seconds=0, minutes=0, hours=0):
114 114 self.throttle_period = datetime.timedelta(
115 115 seconds=seconds, minutes=minutes, hours=hours
116 116 )
117 117
118 118 self.time_of_last_call = datetime.datetime.min
119 119
120 120 def __call__(self, fn):
121 121 @wraps(fn)
122 122 def wrapper(*args, **kwargs):
123 123 coerce = kwargs.pop('coerce', None)
124 124 if coerce:
125 125 self.time_of_last_call = datetime.datetime.now()
126 126 return fn(*args, **kwargs)
127 127 else:
128 128 now = datetime.datetime.now()
129 129 time_since_last_call = now - self.time_of_last_call
130 130 time_left = self.throttle_period - time_since_last_call
131 131
132 132 if time_left > datetime.timedelta(seconds=0):
133 133 return
134 134
135 135 self.time_of_last_call = datetime.datetime.now()
136 136 return fn(*args, **kwargs)
137 137
138 138 return wrapper
139 139
140 140 def apply_throttle(value):
141 141
142 142 @Throttle(seconds=value)
143 143 def fnThrottled(fn):
144 144 fn()
145 145
146 146 return fnThrottled
147 147
148 148
149 149 @MPDecorator
150 150 class Plot(Operation):
151 151 """Base class for Schain plotting operations
152 152
153 153 This class should never be use directtly you must subclass a new operation,
154 154 children classes must be defined as follow:
155 155
156 156 ExamplePlot(Plot):
157 157
158 158 CODE = 'code'
159 159 colormap = 'jet'
160 160 plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
161 161
162 162 def setup(self):
163 163 pass
164 164
165 165 def plot(self):
166 166 pass
167 167
168 168 """
169 169
170 170 CODE = 'Figure'
171 171 colormap = 'jet'
172 172 bgcolor = 'white'
173 173 buffering = True
174 174 __missing = 1E30
175 175
176 176 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
177 177 'showprofile']
178 178
179 179 def __init__(self):
180 180
181 181 Operation.__init__(self)
182 182 self.isConfig = False
183 183 self.isPlotConfig = False
184 184 self.save_time = 0
185 185 self.sender_time = 0
186 186 self.data = None
187 187 self.firsttime = True
188 188 self.sender_queue = deque(maxlen=10)
189 189 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
190 190
191 191 def __fmtTime(self, x, pos):
192 192 '''
193 193 '''
194 194
195 195 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
196 196
197 197 def __setup(self, **kwargs):
198 198 '''
199 199 Initialize variables
200 200 '''
201 201
202 202 self.figures = []
203 203 self.axes = []
204 204 self.cb_axes = []
205 205 self.localtime = kwargs.pop('localtime', True)
206 206 self.show = kwargs.get('show', True)
207 207 self.save = kwargs.get('save', False)
208 208 self.save_period = kwargs.get('save_period', 0)
209 209 self.colormap = kwargs.get('colormap', self.colormap)
210 210 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
211 211 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
212 212 self.colormaps = kwargs.get('colormaps', None)
213 213 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
214 214 self.showprofile = kwargs.get('showprofile', False)
215 215 self.title = kwargs.get('wintitle', self.CODE.upper())
216 216 self.cb_label = kwargs.get('cb_label', None)
217 217 self.cb_labels = kwargs.get('cb_labels', None)
218 218 self.labels = kwargs.get('labels', None)
219 219 self.xaxis = kwargs.get('xaxis', 'frequency')
220 220 self.zmin = kwargs.get('zmin', None)
221 221 self.zmax = kwargs.get('zmax', None)
222 222 self.zlimits = kwargs.get('zlimits', None)
223 223 self.xlimits = kwargs.get('xlimits', None)
224 224 self.xstep_given = kwargs.get('xstep_given', None)
225 225 self.ystep_given = kwargs.get('ystep_given', None)
226 226 self.autoxticks = kwargs.get('autoxticks', True)
227 227 self.xmin = kwargs.get('xmin', None)
228 228 self.xmax = kwargs.get('xmax', None)
229 229 self.xrange = kwargs.get('xrange', 12)
230 230 self.xscale = kwargs.get('xscale', None)
231 231 self.ymin = kwargs.get('ymin', None)
232 232 self.ymax = kwargs.get('ymax', None)
233 233 self.yscale = kwargs.get('yscale', None)
234 234 self.xlabel = kwargs.get('xlabel', None)
235 235 self.attr_time = kwargs.get('attr_time', 'utctime')
236 236 self.attr_data = kwargs.get('attr_data', 'data_param')
237 237 self.decimation = kwargs.get('decimation', None)
238 238 self.oneFigure = kwargs.get('oneFigure', True)
239 239 self.width = kwargs.get('width', None)
240 240 self.height = kwargs.get('height', None)
241 241 self.colorbar = kwargs.get('colorbar', True)
242 242 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
243 243 self.channels = kwargs.get('channels', None)
244 244 self.titles = kwargs.get('titles', [])
245 245 self.polar = False
246 246 self.type = kwargs.get('type', 'iq')
247 247 self.grid = kwargs.get('grid', False)
248 248 self.pause = kwargs.get('pause', False)
249 249 self.save_code = kwargs.get('save_code', self.CODE)
250 250 self.throttle = kwargs.get('throttle', 0)
251 251 self.exp_code = kwargs.get('exp_code', None)
252 252 self.server = kwargs.get('server', False)
253 253 self.sender_period = kwargs.get('sender_period', 60)
254 254 self.tag = kwargs.get('tag', '')
255 255 self.height_index = kwargs.get('height_index', None)
256 256 self.__throttle_plot = apply_throttle(self.throttle)
257 257 code = self.attr_data if self.attr_data else self.CODE
258 258 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
259 #self.EEJtype = kwargs.get('EEJtype', 2)
259 260
260 261 if self.server:
261 262 if not self.server.startswith('tcp://'):
262 263 self.server = 'tcp://{}'.format(self.server)
263 264 log.success(
264 265 'Sending to server: {}'.format(self.server),
265 266 self.name
266 267 )
267 268
268 269 if isinstance(self.attr_data, str):
269 270 self.attr_data = [self.attr_data]
270 271
271 272 def __setup_plot(self):
272 273 '''
273 274 Common setup for all figures, here figures and axes are created
274 275 '''
275 276
276 277 self.setup()
277 278
278 279 self.time_label = 'LT' if self.localtime else 'UTC'
279 280
280 281 if self.width is None:
281 282 self.width = 8
282 283
283 284 self.figures = []
284 285 self.axes = []
285 286 self.cb_axes = []
286 287 self.pf_axes = []
287 288 self.cmaps = []
288 289
289 290 size = '15%' if self.ncols == 1 else '30%'
290 291 pad = '4%' if self.ncols == 1 else '8%'
291 292
292 293 if self.oneFigure:
293 294 if self.height is None:
294 295 self.height = 1.4 * self.nrows + 1
295 296 fig = plt.figure(figsize=(self.width, self.height),
296 297 edgecolor='k',
297 298 facecolor='w')
298 299 self.figures.append(fig)
299 300 for n in range(self.nplots):
300 301 ax = fig.add_subplot(self.nrows, self.ncols,
301 302 n + 1, polar=self.polar)
302 303 ax.tick_params(labelsize=8)
303 304 ax.firsttime = True
304 305 ax.index = 0
305 306 ax.press = None
306 307 self.axes.append(ax)
307 308 if self.showprofile:
308 309 cax = self.__add_axes(ax, size=size, pad=pad)
309 310 cax.tick_params(labelsize=8)
310 311 self.pf_axes.append(cax)
311 312 else:
312 313 if self.height is None:
313 314 self.height = 3
314 315 for n in range(self.nplots):
315 316 fig = plt.figure(figsize=(self.width, self.height),
316 317 edgecolor='k',
317 318 facecolor='w')
318 319 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
319 320 ax.tick_params(labelsize=8)
320 321 ax.firsttime = True
321 322 ax.index = 0
322 323 ax.press = None
323 324 self.figures.append(fig)
324 325 self.axes.append(ax)
325 326 if self.showprofile:
326 327 cax = self.__add_axes(ax, size=size, pad=pad)
327 328 cax.tick_params(labelsize=8)
328 329 self.pf_axes.append(cax)
329 330
330 331 for n in range(self.nrows):
331 332 if self.colormaps is not None:
332 333 cmap = plt.get_cmap(self.colormaps[n])
333 334 else:
334 335 cmap = plt.get_cmap(self.colormap)
335 336 cmap.set_bad(self.bgcolor, 1.)
336 337 self.cmaps.append(cmap)
337 338
338 339 def __add_axes(self, ax, size='30%', pad='8%'):
339 340 '''
340 341 Add new axes to the given figure
341 342 '''
342 343 divider = make_axes_locatable(ax)
343 344 nax = divider.new_horizontal(size=size, pad=pad)
344 345 ax.figure.add_axes(nax)
345 346 return nax
346 347
347 348 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
348 349 '''
349 350 Create a masked array for missing data
350 351 '''
351 352 if x_buffer.shape[0] < 2:
352 353 return x_buffer, y_buffer, z_buffer
353 354
354 355 deltas = x_buffer[1:] - x_buffer[0:-1]
355 356 x_median = numpy.median(deltas)
356 357
357 358 index = numpy.where(deltas > 5 * x_median)
358 359
359 360 if len(index[0]) != 0:
360 361 z_buffer[::, index[0], ::] = self.__missing
361 362 z_buffer = numpy.ma.masked_inside(z_buffer,
362 363 0.99 * self.__missing,
363 364 1.01 * self.__missing)
364 365
365 366 return x_buffer, y_buffer, z_buffer
366 367
367 368 def decimate(self):
368 369
369 370 # dx = int(len(self.x)/self.__MAXNUMX) + 1
370 371 dy = int(len(self.y) / self.decimation) + 1
371 372
372 373 # x = self.x[::dx]
373 374 x = self.x
374 375 y = self.y[::dy]
375 376 z = self.z[::, ::, ::dy]
376 377
377 378 return x, y, z
378 379
379 380 def format(self):
380 381 '''
381 382 Set min and max values, labels, ticks and titles
382 383 '''
383 384
384 385 for n, ax in enumerate(self.axes):
385 386 if ax.firsttime:
386 387 if self.xaxis != 'time':
387 388 xmin = self.xmin
388 389 xmax = self.xmax
389 390 else:
390 391 xmin = self.tmin
391 392 xmax = self.tmin + self.xrange*60*60
392 393 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
393 394 ax.xaxis.set_major_locator(LinearLocator(9))
394 395 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
395 396 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
396 397 ax.set_facecolor(self.bgcolor)
397 398 if self.xscale:
398 399 ax.xaxis.set_major_formatter(FuncFormatter(
399 400 lambda x, pos: '{0:g}'.format(x*self.xscale)))
400 401 if self.yscale:
401 402 ax.yaxis.set_major_formatter(FuncFormatter(
402 403 lambda x, pos: '{0:g}'.format(x*self.yscale)))
403 404 if self.xlabel is not None:
404 405 ax.set_xlabel(self.xlabel)
405 406 if self.ylabel is not None:
406 407 ax.set_ylabel(self.ylabel)
407 408 if self.showprofile:
408 409 self.pf_axes[n].set_ylim(ymin, ymax)
409 410 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
410 411 self.pf_axes[n].set_xlabel('dB')
411 412 self.pf_axes[n].grid(b=True, axis='x')
412 413 [tick.set_visible(False)
413 414 for tick in self.pf_axes[n].get_yticklabels()]
414 415 if self.colorbar:
415 416 ax.cbar = plt.colorbar(
416 417 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
417 418 ax.cbar.ax.tick_params(labelsize=8)
418 419 ax.cbar.ax.press = None
419 420 if self.cb_label:
420 421 ax.cbar.set_label(self.cb_label, size=8)
421 422 elif self.cb_labels:
422 423 ax.cbar.set_label(self.cb_labels[n], size=8)
423 424 else:
424 425 ax.cbar = None
425 426 ax.set_xlim(xmin, xmax)
426 427 ax.set_ylim(ymin, ymax)
427 428 ax.firsttime = False
428 429 if self.grid:
429 430 ax.grid(True)
430 431 if not self.polar:
431 432 ax.set_title('{} {} {}'.format(
432 433 self.titles[n],
433 434 self.getDateTime(self.data.max_time).strftime(
434 435 '%Y-%m-%d %H:%M:%S'),
435 436 self.time_label),
436 437 size=8)
437 438 else:
438 439 ax.set_title('{}'.format(self.titles[n]), size=8)
439 440 ax.set_ylim(0, 90)
440 441 ax.set_yticks(numpy.arange(0, 90, 20))
441 442 ax.yaxis.labelpad = 40
442 443
443 444 if self.firsttime:
444 445 for n, fig in enumerate(self.figures):
445 446 fig.subplots_adjust(**self.plots_adjust)
446 447 self.firsttime = False
447 448
448 449 def clear_figures(self):
449 450 '''
450 451 Reset axes for redraw plots
451 452 '''
452 453
453 454 for ax in self.axes+self.pf_axes+self.cb_axes:
454 455 ax.clear()
455 456 ax.firsttime = True
456 457 if hasattr(ax, 'cbar') and ax.cbar:
457 458 ax.cbar.remove()
458 459
459 460 def __plot(self):
460 461 '''
461 462 Main function to plot, format and save figures
462 463 '''
463 464
464 465 self.plot()
465 466 self.format()
466 467
467 468 for n, fig in enumerate(self.figures):
468 469 if self.nrows == 0 or self.nplots == 0:
469 470 log.warning('No data', self.name)
470 471 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
471 472 fig.canvas.manager.set_window_title(self.CODE)
472 473 continue
473 474
474 475 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
475 476 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
476 477 fig.canvas.draw()
477 478 if self.show:
478 479 fig.show()
479 480 figpause(0.01)
480 481
481 482 if self.save:
482 483 self.save_figure(n)
483 484
484 485 if self.server:
485 486 self.send_to_server()
486 487
487 488 def __update(self, dataOut, timestamp):
488 489 '''
489 490 '''
490 491
491 492 metadata = {
492 493 'yrange': dataOut.heightList,
493 494 'interval': dataOut.timeInterval,
494 495 'channels': dataOut.channelList
495 496 }
496 497
497 498 data, meta = self.update(dataOut)
498 499 metadata.update(meta)
499 500 self.data.update(data, timestamp, metadata)
500 501
501 502 def save_figure(self, n):
502 503 '''
503 504 '''
504 505
505 506 if (self.data.max_time - self.save_time) <= self.save_period:
506 507 return
507 508
508 509 self.save_time = self.data.max_time
509 510
510 511 fig = self.figures[n]
511 512
512 513 if self.throttle == 0:
513 514 figname = os.path.join(
514 515 self.save,
515 516 self.save_code,
516 517 '{}_{}.png'.format(
517 518 self.save_code,
518 519 self.getDateTime(self.data.max_time).strftime(
519 520 '%Y%m%d_%H%M%S'
520 521 ),
521 522 )
522 523 )
523 524 log.log('Saving figure: {}'.format(figname), self.name)
524 525 if not os.path.isdir(os.path.dirname(figname)):
525 526 os.makedirs(os.path.dirname(figname))
526 527 fig.savefig(figname)
527 528
528 529 figname = os.path.join(
529 530 self.save,
530 531 #self.save_code,
531 532 '{}_{}.png'.format(
532 533 self.save_code,
533 534 self.getDateTime(self.data.min_time).strftime(
534 535 '%Y%m%d'
535 536 ),
536 537 )
537 538 )
538 539 log.log('Saving figure: {}'.format(figname), self.name)
539 540 if not os.path.isdir(os.path.dirname(figname)):
540 541 os.makedirs(os.path.dirname(figname))
541 542 fig.savefig(figname)
542 543
543 544 def send_to_server(self):
544 545 '''
545 546 '''
546 547
547 548 if self.exp_code == None:
548 549 log.warning('Missing `exp_code` skipping sending to server...')
549 550
550 551 last_time = self.data.max_time
551 552 interval = last_time - self.sender_time
552 553 if interval < self.sender_period:
553 554 return
554 555
555 556 self.sender_time = last_time
556 557
557 558 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
558 559 for attr in attrs:
559 560 value = getattr(self, attr)
560 561 if value:
561 562 if isinstance(value, (numpy.float32, numpy.float64)):
562 563 value = round(float(value), 2)
563 564 self.data.meta[attr] = value
564 565 if self.colormap == 'jet':
565 566 self.data.meta['colormap'] = 'Jet'
566 567 elif 'RdBu' in self.colormap:
567 568 self.data.meta['colormap'] = 'RdBu'
568 569 else:
569 570 self.data.meta['colormap'] = 'Viridis'
570 571 self.data.meta['interval'] = int(interval)
571 572
572 573 self.sender_queue.append(last_time)
573 574
574 575 while True:
575 576 try:
576 577 tm = self.sender_queue.popleft()
577 578 except IndexError:
578 579 break
579 580 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
580 581 self.socket.send_string(msg)
581 582 socks = dict(self.poll.poll(2000))
582 583 if socks.get(self.socket) == zmq.POLLIN:
583 584 reply = self.socket.recv_string()
584 585 if reply == 'ok':
585 586 log.log("Response from server ok", self.name)
586 587 time.sleep(0.1)
587 588 continue
588 589 else:
589 590 log.warning(
590 591 "Malformed reply from server: {}".format(reply), self.name)
591 592 else:
592 593 log.warning(
593 594 "No response from server, retrying...", self.name)
594 595 self.sender_queue.appendleft(tm)
595 596 self.socket.setsockopt(zmq.LINGER, 0)
596 597 self.socket.close()
597 598 self.poll.unregister(self.socket)
598 599 self.socket = self.context.socket(zmq.REQ)
599 600 self.socket.connect(self.server)
600 601 self.poll.register(self.socket, zmq.POLLIN)
601 602 break
602 603
603 604 def setup(self):
604 605 '''
605 606 This method should be implemented in the child class, the following
606 607 attributes should be set:
607 608
608 609 self.nrows: number of rows
609 610 self.ncols: number of cols
610 611 self.nplots: number of plots (channels or pairs)
611 612 self.ylabel: label for Y axes
612 613 self.titles: list of axes title
613 614
614 615 '''
615 616 raise NotImplementedError
616 617
617 618 def plot(self):
618 619 '''
619 620 Must be defined in the child class, the actual plotting method
620 621 '''
621 622 raise NotImplementedError
622 623
623 624 def update(self, dataOut):
624 625 '''
625 626 Must be defined in the child class, update self.data with new data
626 627 '''
627 628
628 629 data = {
629 630 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
630 631 }
631 632 meta = {}
632 633
633 634 return data, meta
634 635
635 636 def run(self, dataOut, **kwargs):
636 637 '''
637 638 Main plotting routine
638 639 '''
639 640
640 641 if self.isConfig is False:
641 642 self.__setup(**kwargs)
642 643
643 644 if self.localtime:
644 645 self.getDateTime = datetime.datetime.fromtimestamp
645 646 else:
646 647 self.getDateTime = datetime.datetime.utcfromtimestamp
647 648
648 649 self.data.setup()
649 650 self.isConfig = True
650 651 if self.server:
651 652 self.context = zmq.Context()
652 653 self.socket = self.context.socket(zmq.REQ)
653 654 self.socket.connect(self.server)
654 655 self.poll = zmq.Poller()
655 656 self.poll.register(self.socket, zmq.POLLIN)
656 657
657 658 tm = getattr(dataOut, self.attr_time)
658 659
659 660 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
660 661 self.save_time = tm
661 662 self.__plot()
662 663 self.tmin += self.xrange*60*60
663 664 self.data.setup()
664 665 self.clear_figures()
665 666
666 667 self.__update(dataOut, tm)
667 668
668 669 if self.isPlotConfig is False:
669 670 self.__setup_plot()
670 671 self.isPlotConfig = True
671 672 if self.xaxis == 'time':
672 673 dt = self.getDateTime(tm)
673 674 if self.xmin is None:
674 675 self.tmin = tm
675 676 self.xmin = dt.hour
676 677 minutes = (self.xmin-int(self.xmin)) * 60
677 678 seconds = (minutes - int(minutes)) * 60
678 679 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
679 680 datetime.datetime(1970, 1, 1)).total_seconds()
680 681 if self.localtime:
681 682 self.tmin += time.timezone
682 683
683 684 if self.xmin is not None and self.xmax is not None:
684 685 self.xrange = self.xmax - self.xmin
685 686
686 687 if self.throttle == 0:
687 688 self.__plot()
688 689 else:
689 690 self.__throttle_plot(self.__plot)#, coerce=coerce)
690 691
691 692 def close(self):
692 693
693 694 if self.data and not self.data.flagNoData:
694 695 self.save_time = 0
695 696 self.__plot()
696 697 if self.data and not self.data.flagNoData and self.pause:
697 698 figpause(10)
Show file before | Show file after | Hide comments
@@ -1,381 +1,494
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
7 7 from schainpy.utils import log
8 8
9 9 EARTH_RADIUS = 6.3710e3
10 10
11 11
12 12 def ll2xy(lat1, lon1, lat2, lon2):
13 13
14 14 p = 0.017453292519943295
15 15 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
16 16 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
17 17 r = 12742 * numpy.arcsin(numpy.sqrt(a))
18 18 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
19 19 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
20 20 theta = -theta + numpy.pi/2
21 21 return r*numpy.cos(theta), r*numpy.sin(theta)
22 22
23 23
24 24 def km2deg(km):
25 25 '''
26 26 Convert distance in km to degrees
27 27 '''
28 28
29 29 return numpy.rad2deg(km/EARTH_RADIUS)
30 30
31 31
32 32
33 33 class SpectralMomentsPlot(SpectraPlot):
34 34 '''
35 35 Plot for Spectral Moments
36 36 '''
37 37 CODE = 'spc_moments'
38 38 # colormap = 'jet'
39 39 # plot_type = 'pcolor'
40 40
41 41 class DobleGaussianPlot(SpectraPlot):
42 42 '''
43 43 Plot for Double Gaussian Plot
44 44 '''
45 45 CODE = 'gaussian_fit'
46 46 # colormap = 'jet'
47 47 # plot_type = 'pcolor'
48 48
49 49
50 50 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
51 51 '''
52 52 Plot SpectraCut with Double Gaussian Fit
53 53 '''
54 54 CODE = 'cut_gaussian_fit'
55 55
56 56
57 57 class SpectralFitObliquePlot(SpectraPlot):
58 58 '''
59 59 Plot for Spectral Oblique
60 60 '''
61 61 CODE = 'spc_moments'
62 62 colormap = 'jet'
63 63 plot_type = 'pcolor'
64 64
65 65
66 66
67 67 class SnrPlot(RTIPlot):
68 68 '''
69 69 Plot for SNR Data
70 70 '''
71 71
72 72 CODE = 'snr'
73 73 colormap = 'jet'
74 74
75 75 def update(self, dataOut):
76 76
77 77 data = {
78 'snr': 10*numpy.log10(dataOut.data_snr)
78 'snr': 10*numpy.log10(dataOut.data_snr)
79 79 }
80 80
81 81 return data, {}
82 82
83 83 class DopplerPlot(RTIPlot):
84 84 '''
85 85 Plot for DOPPLER Data (1st moment)
86 86 '''
87 87
88 88 CODE = 'dop'
89 89 colormap = 'jet'
90 90
91 91 def update(self, dataOut):
92 92
93 93 data = {
94 'dop': 10*numpy.log10(dataOut.data_dop)
94 'dop': 10*numpy.log10(dataOut.data_dop)
95 }
96
97 return data, {}
98
99 class DopplerEEJPlot_V0(RTIPlot):
100 '''
101 Written by R. Flores
102 '''
103 '''
104 Plot for EEJ
105 '''
106
107 CODE = 'dop'
108 colormap = 'RdBu_r'
109 colormap = 'jet'
110
111 def setup(self):
112
113 self.xaxis = 'time'
114 self.ncols = 1
115 self.nrows = len(self.data.channels)
116 self.nplots = len(self.data.channels)
117 self.ylabel = 'Range [km]'
118 self.xlabel = 'Time'
119 self.cb_label = '(m/s)'
120 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
121 self.titles = ['{} Channel {}'.format(
122 self.CODE.upper(), x) for x in range(self.nrows)]
123
124 def update(self, dataOut):
125 #print(self.EEJtype)
126
127 if self.EEJtype == 1:
128 data = {
129 'dop': dataOut.Oblique_params[:,-2,:]
130 }
131 elif self.EEJtype == 2:
132 data = {
133 'dop': dataOut.Oblique_params[:,-1,:]
134 }
135
136 return data, {}
137
138 class DopplerEEJPlot(RTIPlot):
139 '''
140 Written by R. Flores
141 '''
142 '''
143 Plot for Doppler Shift EEJ
144 '''
145
146 CODE = 'dop'
147 colormap = 'RdBu_r'
148 #colormap = 'jet'
149
150 def setup(self):
151
152 self.xaxis = 'time'
153 self.ncols = 1
154 self.nrows = 2
155 self.nplots = 2
156 self.ylabel = 'Range [km]'
157 self.xlabel = 'Time'
158 self.cb_label = '(m/s)'
159 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
160 self.titles = ['{} EJJ Type {} /'.format(
161 self.CODE.upper(), x) for x in range(1,1+self.nrows)]
162
163 def update(self, dataOut):
164
165 if dataOut.mode == 11: #Double Gaussian
166 doppler = numpy.append(dataOut.Oblique_params[:,1,:],dataOut.Oblique_params[:,4,:],axis=0)
167 elif dataOut.mode == 9: #Double Skew Gaussian
168 doppler = numpy.append(dataOut.Oblique_params[:,-2,:],dataOut.Oblique_params[:,-1,:],axis=0)
169 data = {
170 'dop': doppler
171 }
172
173 return data, {}
174
175 class SpcWidthEEJPlot(RTIPlot):
176 '''
177 Written by R. Flores
178 '''
179 '''
180 Plot for EEJ Spectral Width
181 '''
182
183 CODE = 'width'
184 colormap = 'RdBu_r'
185 colormap = 'jet'
186
187 def setup(self):
188
189 self.xaxis = 'time'
190 self.ncols = 1
191 self.nrows = 2
192 self.nplots = 2
193 self.ylabel = 'Range [km]'
194 self.xlabel = 'Time'
195 self.cb_label = '(m/s)'
196 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
197 self.titles = ['{} EJJ Type {} /'.format(
198 self.CODE.upper(), x) for x in range(1,1+self.nrows)]
199
200 def update(self, dataOut):
201
202 if dataOut.mode == 11: #Double Gaussian
203 width = numpy.append(dataOut.Oblique_params[:,2,:],dataOut.Oblique_params[:,5,:],axis=0)
204 elif dataOut.mode == 9: #Double Skew Gaussian
205 width = numpy.append(dataOut.Oblique_params[:,2,:],dataOut.Oblique_params[:,6,:],axis=0)
206 data = {
207 'width': width
95 208 }
96 209
97 210 return data, {}
98 211
99 212 class PowerPlot(RTIPlot):
100 213 '''
101 214 Plot for Power Data (0 moment)
102 215 '''
103 216
104 217 CODE = 'pow'
105 218 colormap = 'jet'
106 219
107 220 def update(self, dataOut):
108 221
109 222 data = {
110 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
223 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
111 224 }
112 225
113 226 return data, {}
114 227
115 228 class SpectralWidthPlot(RTIPlot):
116 229 '''
117 230 Plot for Spectral Width Data (2nd moment)
118 231 '''
119 232
120 233 CODE = 'width'
121 234 colormap = 'jet'
122 235
123 236 def update(self, dataOut):
124 237
125 238 data = {
126 239 'width': dataOut.data_width
127 240 }
128 241
129 242 return data, {}
130 243
131 244 class SkyMapPlot(Plot):
132 245 '''
133 246 Plot for meteors detection data
134 247 '''
135 248
136 249 CODE = 'param'
137 250
138 251 def setup(self):
139 252
140 253 self.ncols = 1
141 254 self.nrows = 1
142 255 self.width = 7.2
143 256 self.height = 7.2
144 257 self.nplots = 1
145 258 self.xlabel = 'Zonal Zenith Angle (deg)'
146 259 self.ylabel = 'Meridional Zenith Angle (deg)'
147 260 self.polar = True
148 261 self.ymin = -180
149 262 self.ymax = 180
150 263 self.colorbar = False
151 264
152 265 def plot(self):
153 266
154 267 arrayParameters = numpy.concatenate(self.data['param'])
155 268 error = arrayParameters[:, -1]
156 269 indValid = numpy.where(error == 0)[0]
157 270 finalMeteor = arrayParameters[indValid, :]
158 271 finalAzimuth = finalMeteor[:, 3]
159 272 finalZenith = finalMeteor[:, 4]
160 273
161 274 x = finalAzimuth * numpy.pi / 180
162 275 y = finalZenith
163 276
164 277 ax = self.axes[0]
165 278
166 279 if ax.firsttime:
167 280 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
168 281 else:
169 282 ax.plot.set_data(x, y)
170 283
171 284 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
172 285 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
173 286 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
174 287 dt2,
175 288 len(x))
176 289 self.titles[0] = title
177 290
178 291
179 292 class GenericRTIPlot(Plot):
180 293 '''
181 294 Plot for data_xxxx object
182 295 '''
183 296
184 297 CODE = 'param'
185 298 colormap = 'viridis'
186 299 plot_type = 'pcolorbuffer'
187 300
188 301 def setup(self):
189 302 self.xaxis = 'time'
190 303 self.ncols = 1
191 304 self.nrows = self.data.shape('param')[0]
192 305 self.nplots = self.nrows
193 306 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
194 307
195 308 if not self.xlabel:
196 309 self.xlabel = 'Time'
197 310
198 311 self.ylabel = 'Range [km]'
199 312 if not self.titles:
200 313 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
201 314
202 315 def update(self, dataOut):
203 316
204 317 data = {
205 318 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
206 319 }
207 320
208 321 meta = {}
209 322
210 323 return data, meta
211
324
212 325 def plot(self):
213 326 # self.data.normalize_heights()
214 327 self.x = self.data.times
215 328 self.y = self.data.yrange
216 329 self.z = self.data['param']
217 330
218 331 self.z = numpy.ma.masked_invalid(self.z)
219 332
220 333 if self.decimation is None:
221 334 x, y, z = self.fill_gaps(self.x, self.y, self.z)
222 335 else:
223 336 x, y, z = self.fill_gaps(*self.decimate())
224 337
225 338 for n, ax in enumerate(self.axes):
226 339
227 340 self.zmax = self.zmax if self.zmax is not None else numpy.max(
228 341 self.z[n])
229 342 self.zmin = self.zmin if self.zmin is not None else numpy.min(
230 343 self.z[n])
231 344
232 345 if ax.firsttime:
233 346 if self.zlimits is not None:
234 347 self.zmin, self.zmax = self.zlimits[n]
235 348
236 349 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
237 350 vmin=self.zmin,
238 351 vmax=self.zmax,
239 352 cmap=self.cmaps[n]
240 353 )
241 354 else:
242 355 if self.zlimits is not None:
243 356 self.zmin, self.zmax = self.zlimits[n]
244 357 ax.collections.remove(ax.collections[0])
245 358 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
246 359 vmin=self.zmin,
247 360 vmax=self.zmax,
248 361 cmap=self.cmaps[n]
249 362 )
250 363
251 364
252 365 class PolarMapPlot(Plot):
253 366 '''
254 367 Plot for weather radar
255 368 '''
256 369
257 370 CODE = 'param'
258 371 colormap = 'seismic'
259 372
260 373 def setup(self):
261 374 self.ncols = 1
262 375 self.nrows = 1
263 376 self.width = 9
264 377 self.height = 8
265 378 self.mode = self.data.meta['mode']
266 379 if self.channels is not None:
267 380 self.nplots = len(self.channels)
268 381 self.nrows = len(self.channels)
269 382 else:
270 383 self.nplots = self.data.shape(self.CODE)[0]
271 384 self.nrows = self.nplots
272 385 self.channels = list(range(self.nplots))
273 386 if self.mode == 'E':
274 387 self.xlabel = 'Longitude'
275 388 self.ylabel = 'Latitude'
276 389 else:
277 390 self.xlabel = 'Range (km)'
278 391 self.ylabel = 'Height (km)'
279 392 self.bgcolor = 'white'
280 393 self.cb_labels = self.data.meta['units']
281 394 self.lat = self.data.meta['latitude']
282 395 self.lon = self.data.meta['longitude']
283 396 self.xmin, self.xmax = float(
284 397 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
285 398 self.ymin, self.ymax = float(
286 399 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
287 400 # self.polar = True
288 401
289 402 def plot(self):
290 403
291 404 for n, ax in enumerate(self.axes):
292 405 data = self.data['param'][self.channels[n]]
293 406
294 407 zeniths = numpy.linspace(
295 408 0, self.data.meta['max_range'], data.shape[1])
296 409 if self.mode == 'E':
297 410 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
298 411 r, theta = numpy.meshgrid(zeniths, azimuths)
299 412 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
300 413 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
301 414 x = km2deg(x) + self.lon
302 415 y = km2deg(y) + self.lat
303 416 else:
304 417 azimuths = numpy.radians(self.data.yrange)
305 418 r, theta = numpy.meshgrid(zeniths, azimuths)
306 419 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
307 420 self.y = zeniths
308 421
309 422 if ax.firsttime:
310 423 if self.zlimits is not None:
311 424 self.zmin, self.zmax = self.zlimits[n]
312 425 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
313 426 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
314 427 vmin=self.zmin,
315 428 vmax=self.zmax,
316 429 cmap=self.cmaps[n])
317 430 else:
318 431 if self.zlimits is not None:
319 432 self.zmin, self.zmax = self.zlimits[n]
320 433 ax.collections.remove(ax.collections[0])
321 434 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
322 435 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
323 436 vmin=self.zmin,
324 437 vmax=self.zmax,
325 438 cmap=self.cmaps[n])
326 439
327 440 if self.mode == 'A':
328 441 continue
329 442
330 443 # plot district names
331 444 f = open('/data/workspace/schain_scripts/distrito.csv')
332 445 for line in f:
333 446 label, lon, lat = [s.strip() for s in line.split(',') if s]
334 447 lat = float(lat)
335 448 lon = float(lon)
336 449 # ax.plot(lon, lat, '.b', ms=2)
337 450 ax.text(lon, lat, label.decode('utf8'), ha='center',
338 451 va='bottom', size='8', color='black')
339 452
340 453 # plot limites
341 454 limites = []
342 455 tmp = []
343 456 for line in open('/data/workspace/schain_scripts/lima.csv'):
344 457 if '#' in line:
345 458 if tmp:
346 459 limites.append(tmp)
347 460 tmp = []
348 461 continue
349 462 values = line.strip().split(',')
350 463 tmp.append((float(values[0]), float(values[1])))
351 464 for points in limites:
352 465 ax.add_patch(
353 466 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
354 467
355 468 # plot Cuencas
356 469 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
357 470 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
358 471 values = [line.strip().split(',') for line in f]
359 472 points = [(float(s[0]), float(s[1])) for s in values]
360 473 ax.add_patch(Polygon(points, ec='b', fc='none'))
361 474
362 475 # plot grid
363 476 for r in (15, 30, 45, 60):
364 477 ax.add_artist(plt.Circle((self.lon, self.lat),
365 478 km2deg(r), color='0.6', fill=False, lw=0.2))
366 479 ax.text(
367 480 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
368 481 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
369 482 '{}km'.format(r),
370 483 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
371 484
372 485 if self.mode == 'E':
373 486 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
374 487 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
375 488 else:
376 489 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
377 490 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
378 491
379 492 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
380 493 self.titles = ['{} {}'.format(
381 494 self.data.parameters[x], title) for x in self.channels]
Show file before | Show file after | Hide comments
@@ -1,1290 +1,1315
1 1 # Copyright (c) 2012-2021 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13
14 14 class SpectraPlot(Plot):
15 15 '''
16 16 Plot for Spectra data
17 17 '''
18 18
19 19 CODE = 'spc'
20 20 colormap = 'jet'
21 21 plot_type = 'pcolor'
22 22 buffering = False
23 23
24 24 def setup(self):
25 25
26 26 self.nplots = len(self.data.channels)
27 27 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
28 28 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
29 29 self.height = 2.6 * self.nrows
30 30 self.cb_label = 'dB'
31 31 if self.showprofile:
32 32 self.width = 4 * self.ncols
33 33 else:
34 34 self.width = 3.5 * self.ncols
35 35 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
36 36 self.ylabel = 'Range [km]'
37 37
38 38 def update(self, dataOut):
39 39
40 40 data = {}
41 41 meta = {}
42
42 43 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
44 #print("Spc: ",spc[0])
45 #exit(1)
43 46 data['spc'] = spc
44 47 data['rti'] = dataOut.getPower()
48 #print(data['rti'][0])
49 #exit(1)
45 50 #print("NormFactor: ",dataOut.normFactor)
46 51 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
47 52 if hasattr(dataOut, 'LagPlot'): #Double Pulse
48 53 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
49 54 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=46,ymax_index=max_hei_id)/dataOut.normFactor)
50 55 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=40,ymax_index=max_hei_id)/dataOut.normFactor)
51 56 data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)/dataOut.normFactor)
52 57 data['noise'][0] = 10*numpy.log10(dataOut.getNoise(ymin_index=53)[0]/dataOut.normFactor)
53 58 #data['noise'][1] = 22.035507
54 59 else:
55 60 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
56 61 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=26,ymax_index=44)/dataOut.normFactor)
57 62 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
58 63
59 64 if self.CODE == 'spc_moments':
60 65 data['moments'] = dataOut.moments
61 66 if self.CODE == 'gaussian_fit':
62 67 data['gaussfit'] = dataOut.DGauFitParams
63 68
64 69 return data, meta
65 70
66 71 def plot(self):
67 72
68 73 if self.xaxis == "frequency":
69 74 x = self.data.xrange[0]
70 75 self.xlabel = "Frequency (kHz)"
71 76 elif self.xaxis == "time":
72 77 x = self.data.xrange[1]
73 78 self.xlabel = "Time (ms)"
74 79 else:
75 80 x = self.data.xrange[2]
76 81 self.xlabel = "Velocity (m/s)"
77 82
78 83 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
79 84 x = self.data.xrange[2]
80 85 self.xlabel = "Velocity (m/s)"
81 86
82 87 self.titles = []
83 88
84 89 y = self.data.yrange
85 90 self.y = y
86 91
87 92 data = self.data[-1]
88 93 z = data['spc']
89 94
90 95 self.CODE2 = 'spc_oblique'
91 96
92 97
93 98 for n, ax in enumerate(self.axes):
94 99 noise = data['noise'][n]
95 100 if self.CODE == 'spc_moments':
96 101 mean = data['moments'][n, 1]
97 102 if self.CODE == 'gaussian_fit':
98 103 gau0 = data['gaussfit'][n][2,:,0]
99 104 gau1 = data['gaussfit'][n][2,:,1]
100 105 if ax.firsttime:
101 106 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
102 107 self.xmin = self.xmin if self.xmin else numpy.nanmin(x)#-self.xmax
103 108 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
104 109 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
105 110
106 111 ax.plt = ax.pcolormesh(x, y, z[n].T,
107 112 vmin=self.zmin,
108 113 vmax=self.zmax,
109 114 cmap=plt.get_cmap(self.colormap),
110 115 )
111 116
112 117 if self.showprofile:
113 118 ax.plt_profile = self.pf_axes[n].plot(
114 119 data['rti'][n], y)[0]
115 120 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
116 121 color="k", linestyle="dashed", lw=1)[0]
117 122 if self.CODE == 'spc_moments':
118 123 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
119 124 if self.CODE == 'gaussian_fit':
120 125 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
121 126 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
122 127 else:
123 128 ax.plt.set_array(z[n].T.ravel())
124 129 if self.showprofile:
125 130 ax.plt_profile.set_data(data['rti'][n], y)
126 131 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
127 132 if self.CODE == 'spc_moments':
128 133 ax.plt_mean.set_data(mean, y)
129 134 if self.CODE == 'gaussian_fit':
130 135 ax.plt_gau0.set_data(gau0, y)
131 136 ax.plt_gau1.set_data(gau1, y)
132 137 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
133 138
134 139 class SpectraObliquePlot(Plot):
135 140 '''
136 141 Plot for Spectra data
137 142 '''
138 143
139 144 CODE = 'spc_oblique'
140 145 colormap = 'jet'
141 146 plot_type = 'pcolor'
142 147
143 148 def setup(self):
144 149 self.xaxis = "oblique"
145 150 self.nplots = len(self.data.channels)
146 151 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
147 152 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
148 153 self.height = 2.6 * self.nrows
149 154 self.cb_label = 'dB'
150 155 if self.showprofile:
151 156 self.width = 4 * self.ncols
152 157 else:
153 158 self.width = 3.5 * self.ncols
154 159 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
155 160 self.ylabel = 'Range [km]'
156 161
157 162 def update(self, dataOut):
158 163
159 164 data = {}
160 165 meta = {}
161 166 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
162 167 data['spc'] = spc
163 168 data['rti'] = dataOut.getPower()
164 169 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
165 170 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
166
167 data['shift1'] = dataOut.Oblique_params[0][1]
168 data['shift2'] = dataOut.Oblique_params[0][4]
169 data['shift1_error'] = dataOut.Oblique_param_errors[0][1]
170 data['shift2_error'] = dataOut.Oblique_param_errors[0][4]
171 '''
172 data['shift1'] = dataOut.Oblique_params[0,-2,:]
173 data['shift2'] = dataOut.Oblique_params[0,-1,:]
174 data['shift1_error'] = dataOut.Oblique_param_errors[0,-2,:]
175 data['shift2_error'] = dataOut.Oblique_param_errors[0,-1,:]
176 '''
177 '''
178 data['shift1'] = dataOut.Oblique_params[0,1,:]
179 data['shift2'] = dataOut.Oblique_params[0,4,:]
180 data['shift1_error'] = dataOut.Oblique_param_errors[0,1,:]
181 data['shift2_error'] = dataOut.Oblique_param_errors[0,4,:]
182 '''
183 data['shift1'] = dataOut.Dop_EEJ_T1[0]
184 data['shift2'] = dataOut.Dop_EEJ_T2[0]
185 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
186 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
171 187
172 188 return data, meta
173 189
174 190 def plot(self):
175 191
176 192 if self.xaxis == "frequency":
177 193 x = self.data.xrange[0]
178 194 self.xlabel = "Frequency (kHz)"
179 195 elif self.xaxis == "time":
180 196 x = self.data.xrange[1]
181 197 self.xlabel = "Time (ms)"
182 198 else:
183 199 x = self.data.xrange[2]
184 200 self.xlabel = "Velocity (m/s)"
185 201
186 202 self.titles = []
187 203
188 204 y = self.data.yrange
189 205 self.y = y
190 z = self.data['spc']
206
207 data = self.data[-1]
208 z = data['spc']
191 209
192 210 for n, ax in enumerate(self.axes):
193 211 noise = self.data['noise'][n][-1]
194 shift1 = self.data['shift1']
195 shift2 = self.data['shift2']
196 err1 = self.data['shift1_error']
197 err2 = self.data['shift2_error']
212 shift1 = data['shift1']
213 #print(shift1)
214 shift2 = data['shift2']
215 err1 = data['shift1_error']
216 err2 = data['shift2_error']
198 217 if ax.firsttime:
218
199 219 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
200 220 self.xmin = self.xmin if self.xmin else -self.xmax
201 221 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
202 222 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
203 223 ax.plt = ax.pcolormesh(x, y, z[n].T,
204 224 vmin=self.zmin,
205 225 vmax=self.zmax,
206 226 cmap=plt.get_cmap(self.colormap)
207 227 )
208 228
209 229 if self.showprofile:
210 230 ax.plt_profile = self.pf_axes[n].plot(
211 231 self.data['rti'][n][-1], y)[0]
212 232 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
213 233 color="k", linestyle="dashed", lw=1)[0]
214 234
215 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=0.2, marker='x', linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
216 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
235 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
236 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
237 #print("plotter1: ", self.ploterr1,shift1)
238
217 239 else:
240 #print("else plotter1: ", self.ploterr1,shift1)
218 241 self.ploterr1.remove()
219 242 self.ploterr2.remove()
220 243 ax.plt.set_array(z[n].T.ravel())
221 244 if self.showprofile:
222 245 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
223 246 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
224 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
225 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
247 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
248 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
226 249
227 250 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
228 251
229 252
230 253 class CrossSpectraPlot(Plot):
231 254
232 255 CODE = 'cspc'
233 256 colormap = 'jet'
234 257 plot_type = 'pcolor'
235 258 zmin_coh = None
236 259 zmax_coh = None
237 260 zmin_phase = None
238 261 zmax_phase = None
239 262
240 263 def setup(self):
241 264
242 265 self.ncols = 4
243 266 self.nplots = len(self.data.pairs) * 2
244 267 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
245 268 self.width = 3.1 * self.ncols
246 269 self.height = 5 * self.nrows
247 270 self.ylabel = 'Range [km]'
248 271 self.showprofile = False
249 272 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
250 273
251 274 def update(self, dataOut):
252 275
253 276 data = {}
254 277 meta = {}
255 278
256 279 spc = dataOut.data_spc
257 280 cspc = dataOut.data_cspc
258 281 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
259 282 meta['pairs'] = dataOut.pairsList
260 283
261 284 tmp = []
262 285
263 286 for n, pair in enumerate(meta['pairs']):
264 287 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
265 288 coh = numpy.abs(out)
266 289 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
267 290 tmp.append(coh)
268 291 tmp.append(phase)
269 292
270 293 data['cspc'] = numpy.array(tmp)
271 294
272 295 return data, meta
273 296
274 297 def plot(self):
275 298
276 299 if self.xaxis == "frequency":
277 300 x = self.data.xrange[0]
278 301 self.xlabel = "Frequency (kHz)"
279 302 elif self.xaxis == "time":
280 303 x = self.data.xrange[1]
281 304 self.xlabel = "Time (ms)"
282 305 else:
283 306 x = self.data.xrange[2]
284 307 self.xlabel = "Velocity (m/s)"
285 308
286 309 self.titles = []
287 310
288 311 y = self.data.yrange
289 312 self.y = y
290 313
291 314 data = self.data[-1]
292 315 cspc = data['cspc']
293 316
294 317 for n in range(len(self.data.pairs)):
295 318 pair = self.data.pairs[n]
296 319 coh = cspc[n*2]
297 320 phase = cspc[n*2+1]
298 321 ax = self.axes[2 * n]
299 322 if ax.firsttime:
300 323 ax.plt = ax.pcolormesh(x, y, coh.T,
301 324 vmin=0,
302 325 vmax=1,
303 326 cmap=plt.get_cmap(self.colormap_coh)
304 327 )
305 328 else:
306 329 ax.plt.set_array(coh.T.ravel())
307 330 self.titles.append(
308 331 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
309 332
310 333 ax = self.axes[2 * n + 1]
311 334 if ax.firsttime:
312 335 ax.plt = ax.pcolormesh(x, y, phase.T,
313 336 vmin=-180,
314 337 vmax=180,
315 338 cmap=plt.get_cmap(self.colormap_phase)
316 339 )
317 340 else:
318 341 ax.plt.set_array(phase.T.ravel())
319 342 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
320 343
321 344
322 345 class CrossSpectra4Plot(Plot):
323 346
324 347 CODE = 'cspc'
325 348 colormap = 'jet'
326 349 plot_type = 'pcolor'
327 350 zmin_coh = None
328 351 zmax_coh = None
329 352 zmin_phase = None
330 353 zmax_phase = None
331 354
332 355 def setup(self):
333 356
334 357 self.ncols = 4
335 358 self.nrows = len(self.data.pairs)
336 359 self.nplots = self.nrows * 4
337 360 self.width = 3.1 * self.ncols
338 361 self.height = 5 * self.nrows
339 362 self.ylabel = 'Range [km]'
340 363 self.showprofile = False
341 364 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
342 365
343 366 def plot(self):
344 367
345 368 if self.xaxis == "frequency":
346 369 x = self.data.xrange[0]
347 370 self.xlabel = "Frequency (kHz)"
348 371 elif self.xaxis == "time":
349 372 x = self.data.xrange[1]
350 373 self.xlabel = "Time (ms)"
351 374 else:
352 375 x = self.data.xrange[2]
353 376 self.xlabel = "Velocity (m/s)"
354 377
355 378 self.titles = []
356 379
357 380
358 381 y = self.data.heights
359 382 self.y = y
360 383 nspc = self.data['spc']
361 384 #print(numpy.shape(self.data['spc']))
362 385 spc = self.data['cspc'][0]
363 386 #print(numpy.shape(nspc))
364 387 #exit()
365 388 #nspc[1,:,:] = numpy.flip(nspc[1,:,:],axis=0)
366 389 #print(numpy.shape(spc))
367 390 #exit()
368 391 cspc = self.data['cspc'][1]
369 392
370 393 #xflip=numpy.flip(x)
371 394 #print(numpy.shape(cspc))
372 395 #exit()
373 396
374 397 for n in range(self.nrows):
375 398 noise = self.data['noise'][:,-1]
376 399 pair = self.data.pairs[n]
377 400 #print(pair)
378 401 #exit()
379 402 ax = self.axes[4 * n]
380 403 if ax.firsttime:
381 404 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
382 405 self.xmin = self.xmin if self.xmin else -self.xmax
383 406 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
384 407 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
385 408 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
386 409 vmin=self.zmin,
387 410 vmax=self.zmax,
388 411 cmap=plt.get_cmap(self.colormap)
389 412 )
390 413 else:
391 414 #print(numpy.shape(nspc[pair[0]].T))
392 415 #exit()
393 416 ax.plt.set_array(nspc[pair[0]].T.ravel())
394 417 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
395 418
396 419 ax = self.axes[4 * n + 1]
397 420
398 421 if ax.firsttime:
399 422 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
400 423 vmin=self.zmin,
401 424 vmax=self.zmax,
402 425 cmap=plt.get_cmap(self.colormap)
403 426 )
404 427 else:
405 428
406 429 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
407 430 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
408 431
409 432 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
410 433 coh = numpy.abs(out)
411 434 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
412 435
413 436 ax = self.axes[4 * n + 2]
414 437 if ax.firsttime:
415 438 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
416 439 vmin=0,
417 440 vmax=1,
418 441 cmap=plt.get_cmap(self.colormap_coh)
419 442 )
420 443 else:
421 444 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
422 445 self.titles.append(
423 446 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
424 447
425 448 ax = self.axes[4 * n + 3]
426 449 if ax.firsttime:
427 450 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
428 451 vmin=-180,
429 452 vmax=180,
430 453 cmap=plt.get_cmap(self.colormap_phase)
431 454 )
432 455 else:
433 456 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
434 457 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
435 458
436 459
437 460 class CrossSpectra2Plot(Plot):
438 461
439 462 CODE = 'cspc'
440 463 colormap = 'jet'
441 464 plot_type = 'pcolor'
442 465 zmin_coh = None
443 466 zmax_coh = None
444 467 zmin_phase = None
445 468 zmax_phase = None
446 469
447 470 def setup(self):
448 471
449 472 self.ncols = 1
450 473 self.nrows = len(self.data.pairs)
451 474 self.nplots = self.nrows * 1
452 475 self.width = 3.1 * self.ncols
453 476 self.height = 5 * self.nrows
454 477 self.ylabel = 'Range [km]'
455 478 self.showprofile = False
456 479 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
457 480
458 481 def plot(self):
459 482
460 483 if self.xaxis == "frequency":
461 484 x = self.data.xrange[0]
462 485 self.xlabel = "Frequency (kHz)"
463 486 elif self.xaxis == "time":
464 487 x = self.data.xrange[1]
465 488 self.xlabel = "Time (ms)"
466 489 else:
467 490 x = self.data.xrange[2]
468 491 self.xlabel = "Velocity (m/s)"
469 492
470 493 self.titles = []
471 494
472 495
473 496 y = self.data.heights
474 497 self.y = y
475 498 #nspc = self.data['spc']
476 499 #print(numpy.shape(self.data['spc']))
477 500 #spc = self.data['cspc'][0]
478 501 #print(numpy.shape(spc))
479 502 #exit()
480 503 cspc = self.data['cspc'][1]
481 504 #print(numpy.shape(cspc))
482 505 #exit()
483 506
484 507 for n in range(self.nrows):
485 508 noise = self.data['noise'][:,-1]
486 509 pair = self.data.pairs[n]
487 510 #print(pair) #exit()
488 511
489 512
490 513
491 514 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
492 515
493 516 #print(out[:,53])
494 517 #exit()
495 518 cross = numpy.abs(out)
496 519 z = cross/self.data.nFactor
497 520 #print("here")
498 521 #print(dataOut.data_spc[0,0,0])
499 522 #exit()
500 523
501 524 cross = 10*numpy.log10(z)
502 525 #print(numpy.shape(cross))
503 526 #print(cross[0,:])
504 527 #print(self.data.nFactor)
505 528 #exit()
506 529 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
507 530
508 531 ax = self.axes[1 * n]
509 532 if ax.firsttime:
510 533 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
511 534 self.xmin = self.xmin if self.xmin else -self.xmax
512 535 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
513 536 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
514 537 ax.plt = ax.pcolormesh(x, y, cross.T,
515 538 vmin=self.zmin,
516 539 vmax=self.zmax,
517 540 cmap=plt.get_cmap(self.colormap)
518 541 )
519 542 else:
520 543 ax.plt.set_array(cross.T.ravel())
521 544 self.titles.append(
522 545 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
523 546
524 547
525 548 class CrossSpectra3Plot(Plot):
526 549
527 550 CODE = 'cspc'
528 551 colormap = 'jet'
529 552 plot_type = 'pcolor'
530 553 zmin_coh = None
531 554 zmax_coh = None
532 555 zmin_phase = None
533 556 zmax_phase = None
534 557
535 558 def setup(self):
536 559
537 560 self.ncols = 3
538 561 self.nrows = len(self.data.pairs)
539 562 self.nplots = self.nrows * 3
540 563 self.width = 3.1 * self.ncols
541 564 self.height = 5 * self.nrows
542 565 self.ylabel = 'Range [km]'
543 566 self.showprofile = False
544 567 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
545 568
546 569 def plot(self):
547 570
548 571 if self.xaxis == "frequency":
549 572 x = self.data.xrange[0]
550 573 self.xlabel = "Frequency (kHz)"
551 574 elif self.xaxis == "time":
552 575 x = self.data.xrange[1]
553 576 self.xlabel = "Time (ms)"
554 577 else:
555 578 x = self.data.xrange[2]
556 579 self.xlabel = "Velocity (m/s)"
557 580
558 581 self.titles = []
559 582
560 583
561 584 y = self.data.heights
562 585 self.y = y
563 586 #nspc = self.data['spc']
564 587 #print(numpy.shape(self.data['spc']))
565 588 #spc = self.data['cspc'][0]
566 589 #print(numpy.shape(spc))
567 590 #exit()
568 591 cspc = self.data['cspc'][1]
569 592 #print(numpy.shape(cspc))
570 593 #exit()
571 594
572 595 for n in range(self.nrows):
573 596 noise = self.data['noise'][:,-1]
574 597 pair = self.data.pairs[n]
575 598 #print(pair) #exit()
576 599
577 600
578 601
579 602 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
580 603
581 604 #print(out[:,53])
582 605 #exit()
583 606 cross = numpy.abs(out)
584 607 z = cross/self.data.nFactor
585 608 cross = 10*numpy.log10(z)
586 609
587 610 out_r= out.real/self.data.nFactor
588 611 #out_r = 10*numpy.log10(out_r)
589 612
590 613 out_i= out.imag/self.data.nFactor
591 614 #out_i = 10*numpy.log10(out_i)
592 615 #print(numpy.shape(cross))
593 616 #print(cross[0,:])
594 617 #print(self.data.nFactor)
595 618 #exit()
596 619 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
597 620
598 621 ax = self.axes[3 * n]
599 622 if ax.firsttime:
600 623 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
601 624 self.xmin = self.xmin if self.xmin else -self.xmax
602 625 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
603 626 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
604 627 ax.plt = ax.pcolormesh(x, y, cross.T,
605 628 vmin=self.zmin,
606 629 vmax=self.zmax,
607 630 cmap=plt.get_cmap(self.colormap)
608 631 )
609 632 else:
610 633 ax.plt.set_array(cross.T.ravel())
611 634 self.titles.append(
612 635 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
613 636
614 637 ax = self.axes[3 * n + 1]
615 638 if ax.firsttime:
616 639 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
617 640 self.xmin = self.xmin if self.xmin else -self.xmax
618 641 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
619 642 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
620 643 ax.plt = ax.pcolormesh(x, y, out_r.T,
621 644 vmin=-1.e6,
622 645 vmax=0,
623 646 cmap=plt.get_cmap(self.colormap)
624 647 )
625 648 else:
626 649 ax.plt.set_array(out_r.T.ravel())
627 650 self.titles.append(
628 651 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
629 652
630 653 ax = self.axes[3 * n + 2]
631 654
632 655
633 656 if ax.firsttime:
634 657 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
635 658 self.xmin = self.xmin if self.xmin else -self.xmax
636 659 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
637 660 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
638 661 ax.plt = ax.pcolormesh(x, y, out_i.T,
639 662 vmin=-1.e6,
640 663 vmax=1.e6,
641 664 cmap=plt.get_cmap(self.colormap)
642 665 )
643 666 else:
644 667 ax.plt.set_array(out_i.T.ravel())
645 668 self.titles.append(
646 669 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
647 670
648 671 class RTIPlot(Plot):
649 672 '''
650 673 Plot for RTI data
651 674 '''
652 675
653 676 CODE = 'rti'
654 677 colormap = 'jet'
655 678 plot_type = 'pcolorbuffer'
656 679
657 680 def setup(self):
658 681 self.xaxis = 'time'
659 682 self.ncols = 1
660 683 self.nrows = len(self.data.channels)
661 684 self.nplots = len(self.data.channels)
662 685 self.ylabel = 'Range [km]'
663 686 self.xlabel = 'Time'
664 687 self.cb_label = 'dB'
665 688 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
666 689 self.titles = ['{} Channel {}'.format(
667 690 self.CODE.upper(), x) for x in range(self.nrows)]
668 691
669 692 def update(self, dataOut):
670 693
671 694 data = {}
672 695 meta = {}
673 696 data['rti'] = dataOut.getPower()
697 #print(numpy.shape(data['rti']))
674 698
675 699 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
676 700
677 701 return data, meta
678 702
679 703 def plot(self):
680 704
681 705 self.x = self.data.times
682 706 self.y = self.data.yrange
683 707 self.z = self.data[self.CODE]
684 708
685 709 self.z = numpy.ma.masked_invalid(self.z)
686 710
687 711 if self.decimation is None:
688 712 x, y, z = self.fill_gaps(self.x, self.y, self.z)
689 713 else:
690 714 x, y, z = self.fill_gaps(*self.decimate())
691 715
692 716 for n, ax in enumerate(self.axes):
693 717 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
694 718 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
719
695 720 if ax.firsttime:
696 721 ax.plt = ax.pcolormesh(x, y, z[n].T,
697 722 vmin=self.zmin,
698 723 vmax=self.zmax,
699 724 cmap=plt.get_cmap(self.colormap)
700 725 )
701 726 if self.showprofile:
702 727 ax.plot_profile = self.pf_axes[n].plot(
703 728 self.data['rti'][n][-1], self.y)[0]
704 729 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
705 730 color="k", linestyle="dashed", lw=1)[0]
706 731 else:
707 732 ax.collections.remove(ax.collections[0])
708 733 ax.plt = ax.pcolormesh(x, y, z[n].T,
709 734 vmin=self.zmin,
710 735 vmax=self.zmax,
711 736 cmap=plt.get_cmap(self.colormap)
712 737 )
713 738 if self.showprofile:
714 739 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
715 740 ax.plot_noise.set_data(numpy.repeat(
716 741 self.data['noise'][n][-1], len(self.y)), self.y)
717 742
718 743
719 744 class SpectrogramPlot(Plot):
720 745 '''
721 746 Plot for Spectrogram data
722 747 '''
723 748
724 749 CODE = 'Spectrogram_Profile'
725 750 colormap = 'binary'
726 751 plot_type = 'pcolorbuffer'
727 752
728 753 def setup(self):
729 754 self.xaxis = 'time'
730 755 self.ncols = 1
731 756 self.nrows = len(self.data.channels)
732 757 self.nplots = len(self.data.channels)
733 758 self.xlabel = 'Time'
734 759 #self.cb_label = 'dB'
735 760 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
736 761 self.titles = []
737 762
738 763 #self.titles = ['{} Channel {} \n H = {} km ({} - {})'.format(
739 764 #self.CODE.upper(), x, self.data.heightList[self.data.hei], self.data.heightList[self.data.hei],self.data.heightList[self.data.hei]+(self.data.DH*self.data.nProfiles)) for x in range(self.nrows)]
740 765
741 766 self.titles = ['{} Channel {}'.format(
742 767 self.CODE.upper(), x) for x in range(self.nrows)]
743 768
744 769
745 770 def update(self, dataOut):
746 771 data = {}
747 772 meta = {}
748 773
749 774 maxHei = 1620#+12000
750 775 indb = numpy.where(dataOut.heightList <= maxHei)
751 776 hei = indb[0][-1]
752 777 #print(dataOut.heightList)
753 778
754 779 factor = dataOut.nIncohInt
755 780 z = dataOut.data_spc[:,:,hei] / factor
756 781 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
757 782 #buffer = 10 * numpy.log10(z)
758 783
759 784 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
760 785
761 786
762 787 #self.hei = hei
763 788 #self.heightList = dataOut.heightList
764 789 #self.DH = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
765 790 #self.nProfiles = dataOut.nProfiles
766 791
767 792 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
768 793
769 794 data['hei'] = hei
770 795 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
771 796 data['nProfiles'] = dataOut.nProfiles
772 797 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
773 798 '''
774 799 import matplotlib.pyplot as plt
775 800 plt.plot(10 * numpy.log10(z[0,:]))
776 801 plt.show()
777 802
778 803 from time import sleep
779 804 sleep(10)
780 805 '''
781 806 return data, meta
782 807
783 808 def plot(self):
784 809
785 810 self.x = self.data.times
786 811 self.z = self.data[self.CODE]
787 812 self.y = self.data.xrange[0]
788 813
789 814 hei = self.data['hei'][-1]
790 815 DH = self.data['DH'][-1]
791 816 nProfiles = self.data['nProfiles'][-1]
792 817
793 818 self.ylabel = "Frequency (kHz)"
794 819
795 820 self.z = numpy.ma.masked_invalid(self.z)
796 821
797 822 if self.decimation is None:
798 823 x, y, z = self.fill_gaps(self.x, self.y, self.z)
799 824 else:
800 825 x, y, z = self.fill_gaps(*self.decimate())
801 826
802 827 for n, ax in enumerate(self.axes):
803 828 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
804 829 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
805 830 data = self.data[-1]
806 831 if ax.firsttime:
807 832 ax.plt = ax.pcolormesh(x, y, z[n].T,
808 833 vmin=self.zmin,
809 834 vmax=self.zmax,
810 835 cmap=plt.get_cmap(self.colormap)
811 836 )
812 837 else:
813 838 ax.collections.remove(ax.collections[0])
814 839 ax.plt = ax.pcolormesh(x, y, z[n].T,
815 840 vmin=self.zmin,
816 841 vmax=self.zmax,
817 842 cmap=plt.get_cmap(self.colormap)
818 843 )
819 844
820 845 #self.titles.append('Spectrogram')
821 846
822 847 #self.titles.append('{} Channel {} \n H = {} km ({} - {})'.format(
823 848 #self.CODE.upper(), x, y[hei], y[hei],y[hei]+(DH*nProfiles)))
824 849
825 850
826 851
827 852
828 853 class CoherencePlot(RTIPlot):
829 854 '''
830 855 Plot for Coherence data
831 856 '''
832 857
833 858 CODE = 'coh'
834 859
835 860 def setup(self):
836 861 self.xaxis = 'time'
837 862 self.ncols = 1
838 863 self.nrows = len(self.data.pairs)
839 864 self.nplots = len(self.data.pairs)
840 865 self.ylabel = 'Range [km]'
841 866 self.xlabel = 'Time'
842 867 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
843 868 if self.CODE == 'coh':
844 869 self.cb_label = ''
845 870 self.titles = [
846 871 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
847 872 else:
848 873 self.cb_label = 'Degrees'
849 874 self.titles = [
850 875 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
851 876
852 877 def update(self, dataOut):
853 878
854 879 data = {}
855 880 meta = {}
856 881 data['coh'] = dataOut.getCoherence()
857 882 meta['pairs'] = dataOut.pairsList
858 883
859 884 return data, meta
860 885
861 886 class PhasePlot(CoherencePlot):
862 887 '''
863 888 Plot for Phase map data
864 889 '''
865 890
866 891 CODE = 'phase'
867 892 colormap = 'seismic'
868 893
869 894 def update(self, dataOut):
870 895
871 896 data = {}
872 897 meta = {}
873 898 data['phase'] = dataOut.getCoherence(phase=True)
874 899 meta['pairs'] = dataOut.pairsList
875 900
876 901 return data, meta
877 902
878 903 class NoisePlot(Plot):
879 904 '''
880 905 Plot for noise
881 906 '''
882 907
883 908 CODE = 'noise'
884 909 plot_type = 'scatterbuffer'
885 910
886 911 def setup(self):
887 912 self.xaxis = 'time'
888 913 self.ncols = 1
889 914 self.nrows = 1
890 915 self.nplots = 1
891 916 self.ylabel = 'Intensity [dB]'
892 917 self.xlabel = 'Time'
893 918 self.titles = ['Noise']
894 919 self.colorbar = False
895 920 self.plots_adjust.update({'right': 0.85 })
896 921
897 922 def update(self, dataOut):
898 923
899 924 data = {}
900 925 meta = {}
901 926 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
902 927 meta['yrange'] = numpy.array([])
903 928
904 929 return data, meta
905 930
906 931 def plot(self):
907 932
908 933 x = self.data.times
909 934 xmin = self.data.min_time
910 935 xmax = xmin + self.xrange * 60 * 60
911 936 Y = self.data['noise']
912 937
913 938 if self.axes[0].firsttime:
914 939 self.ymin = numpy.nanmin(Y) - 5
915 940 self.ymax = numpy.nanmax(Y) + 5
916 941 for ch in self.data.channels:
917 942 y = Y[ch]
918 943 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
919 944 plt.legend(bbox_to_anchor=(1.18, 1.0))
920 945 else:
921 946 for ch in self.data.channels:
922 947 y = Y[ch]
923 948 self.axes[0].lines[ch].set_data(x, y)
924 949
925 950 self.ymin = numpy.nanmin(Y) - 5
926 951 self.ymax = numpy.nanmax(Y) + 10
927 952
928 953
929 954 class PowerProfilePlot(Plot):
930 955
931 956 CODE = 'pow_profile'
932 957 plot_type = 'scatter'
933 958
934 959 def setup(self):
935 960
936 961 self.ncols = 1
937 962 self.nrows = 1
938 963 self.nplots = 1
939 964 self.height = 4
940 965 self.width = 3
941 966 self.ylabel = 'Range [km]'
942 967 self.xlabel = 'Intensity [dB]'
943 968 self.titles = ['Power Profile']
944 969 self.colorbar = False
945 970
946 971 def update(self, dataOut):
947 972
948 973 data = {}
949 974 meta = {}
950 975 data[self.CODE] = dataOut.getPower()
951 976
952 977 return data, meta
953 978
954 979 def plot(self):
955 980
956 981 y = self.data.yrange
957 982 self.y = y
958 983
959 984 x = self.data[-1][self.CODE]
960 985
961 986 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
962 987 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
963 988
964 989 if self.axes[0].firsttime:
965 990 for ch in self.data.channels:
966 991 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
967 992 plt.legend()
968 993 else:
969 994 for ch in self.data.channels:
970 995 self.axes[0].lines[ch].set_data(x[ch], y)
971 996
972 997
973 998 class SpectraCutPlot(Plot):
974 999
975 1000 CODE = 'spc_cut'
976 1001 plot_type = 'scatter'
977 1002 buffering = False
978 1003
979 1004 def setup(self):
980 1005
981 1006 self.nplots = len(self.data.channels)
982 1007 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
983 1008 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
984 1009 self.width = 3.4 * self.ncols + 1.5
985 1010 self.height = 3 * self.nrows
986 1011 self.ylabel = 'Power [dB]'
987 1012 self.colorbar = False
988 1013 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
989 1014
990 1015 def update(self, dataOut):
991 1016
992 1017 data = {}
993 1018 meta = {}
994 1019 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
995 1020 data['spc'] = spc
996 1021 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
997 1022 if self.CODE == 'cut_gaussian_fit':
998 1023 data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
999 1024 data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
1000 1025 return data, meta
1001 1026
1002 1027 def plot(self):
1003 1028 if self.xaxis == "frequency":
1004 1029 x = self.data.xrange[0][1:]
1005 1030 self.xlabel = "Frequency (kHz)"
1006 1031 elif self.xaxis == "time":
1007 1032 x = self.data.xrange[1]
1008 1033 self.xlabel = "Time (ms)"
1009 1034 else:
1010 1035 x = self.data.xrange[2][:-1]
1011 1036 self.xlabel = "Velocity (m/s)"
1012 1037
1013 1038 if self.CODE == 'cut_gaussian_fit':
1014 1039 x = self.data.xrange[2][:-1]
1015 1040 self.xlabel = "Velocity (m/s)"
1016 1041
1017 1042 self.titles = []
1018 1043
1019 1044 y = self.data.yrange
1020 1045 data = self.data[-1]
1021 1046 z = data['spc']
1022 1047
1023 1048 if self.height_index:
1024 1049 index = numpy.array(self.height_index)
1025 1050 else:
1026 1051 index = numpy.arange(0, len(y), int((len(y))/9))
1027 1052
1028 1053 for n, ax in enumerate(self.axes):
1029 1054 if self.CODE == 'cut_gaussian_fit':
1030 1055 gau0 = data['gauss_fit0']
1031 1056 gau1 = data['gauss_fit1']
1032 1057 if ax.firsttime:
1033 1058 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1034 1059 self.xmin = self.xmin if self.xmin else -self.xmax
1035 1060 self.ymin = self.ymin if self.ymin else numpy.nanmin(z[:,:,index])
1036 1061 self.ymax = self.ymax if self.ymax else numpy.nanmax(z[:,:,index])
1037 1062 #print(self.ymax)
1038 1063 #print(z[n, :, index])
1039 1064 ax.plt = ax.plot(x, z[n, :, index].T, lw=0.25)
1040 1065 if self.CODE == 'cut_gaussian_fit':
1041 1066 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
1042 1067 for i, line in enumerate(ax.plt_gau0):
1043 1068 line.set_color(ax.plt[i].get_color())
1044 1069 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
1045 1070 for i, line in enumerate(ax.plt_gau1):
1046 1071 line.set_color(ax.plt[i].get_color())
1047 1072 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1048 1073 self.figures[0].legend(ax.plt, labels, loc='center right')
1049 1074 else:
1050 1075 for i, line in enumerate(ax.plt):
1051 1076 line.set_data(x, z[n, :, index[i]].T)
1052 1077 for i, line in enumerate(ax.plt_gau0):
1053 1078 line.set_data(x, gau0[n, :, index[i]].T)
1054 1079 line.set_color(ax.plt[i].get_color())
1055 1080 for i, line in enumerate(ax.plt_gau1):
1056 1081 line.set_data(x, gau1[n, :, index[i]].T)
1057 1082 line.set_color(ax.plt[i].get_color())
1058 1083 self.titles.append('CH {}'.format(n))
1059 1084
1060 1085
1061 1086 class BeaconPhase(Plot):
1062 1087
1063 1088 __isConfig = None
1064 1089 __nsubplots = None
1065 1090
1066 1091 PREFIX = 'beacon_phase'
1067 1092
1068 1093 def __init__(self):
1069 1094 Plot.__init__(self)
1070 1095 self.timerange = 24*60*60
1071 1096 self.isConfig = False
1072 1097 self.__nsubplots = 1
1073 1098 self.counter_imagwr = 0
1074 1099 self.WIDTH = 800
1075 1100 self.HEIGHT = 400
1076 1101 self.WIDTHPROF = 120
1077 1102 self.HEIGHTPROF = 0
1078 1103 self.xdata = None
1079 1104 self.ydata = None
1080 1105
1081 1106 self.PLOT_CODE = BEACON_CODE
1082 1107
1083 1108 self.FTP_WEI = None
1084 1109 self.EXP_CODE = None
1085 1110 self.SUB_EXP_CODE = None
1086 1111 self.PLOT_POS = None
1087 1112
1088 1113 self.filename_phase = None
1089 1114
1090 1115 self.figfile = None
1091 1116
1092 1117 self.xmin = None
1093 1118 self.xmax = None
1094 1119
1095 1120 def getSubplots(self):
1096 1121
1097 1122 ncol = 1
1098 1123 nrow = 1
1099 1124
1100 1125 return nrow, ncol
1101 1126
1102 1127 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1103 1128
1104 1129 self.__showprofile = showprofile
1105 1130 self.nplots = nplots
1106 1131
1107 1132 ncolspan = 7
1108 1133 colspan = 6
1109 1134 self.__nsubplots = 2
1110 1135
1111 1136 self.createFigure(id = id,
1112 1137 wintitle = wintitle,
1113 1138 widthplot = self.WIDTH+self.WIDTHPROF,
1114 1139 heightplot = self.HEIGHT+self.HEIGHTPROF,
1115 1140 show=show)
1116 1141
1117 1142 nrow, ncol = self.getSubplots()
1118 1143
1119 1144 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1120 1145
1121 1146 def save_phase(self, filename_phase):
1122 1147 f = open(filename_phase,'w+')
1123 1148 f.write('\n\n')
1124 1149 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1125 1150 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1126 1151 f.close()
1127 1152
1128 1153 def save_data(self, filename_phase, data, data_datetime):
1129 1154 f=open(filename_phase,'a')
1130 1155 timetuple_data = data_datetime.timetuple()
1131 1156 day = str(timetuple_data.tm_mday)
1132 1157 month = str(timetuple_data.tm_mon)
1133 1158 year = str(timetuple_data.tm_year)
1134 1159 hour = str(timetuple_data.tm_hour)
1135 1160 minute = str(timetuple_data.tm_min)
1136 1161 second = str(timetuple_data.tm_sec)
1137 1162 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1138 1163 f.close()
1139 1164
1140 1165 def plot(self):
1141 1166 log.warning('TODO: Not yet implemented...')
1142 1167
1143 1168 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1144 1169 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1145 1170 timerange=None,
1146 1171 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1147 1172 server=None, folder=None, username=None, password=None,
1148 1173 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1149 1174
1150 1175 if dataOut.flagNoData:
1151 1176 return dataOut
1152 1177
1153 1178 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1154 1179 return
1155 1180
1156 1181 if pairsList == None:
1157 1182 pairsIndexList = dataOut.pairsIndexList[:10]
1158 1183 else:
1159 1184 pairsIndexList = []
1160 1185 for pair in pairsList:
1161 1186 if pair not in dataOut.pairsList:
1162 1187 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1163 1188 pairsIndexList.append(dataOut.pairsList.index(pair))
1164 1189
1165 1190 if pairsIndexList == []:
1166 1191 return
1167 1192
1168 1193 # if len(pairsIndexList) > 4:
1169 1194 # pairsIndexList = pairsIndexList[0:4]
1170 1195
1171 1196 hmin_index = None
1172 1197 hmax_index = None
1173 1198
1174 1199 if hmin != None and hmax != None:
1175 1200 indexes = numpy.arange(dataOut.nHeights)
1176 1201 hmin_list = indexes[dataOut.heightList >= hmin]
1177 1202 hmax_list = indexes[dataOut.heightList <= hmax]
1178 1203
1179 1204 if hmin_list.any():
1180 1205 hmin_index = hmin_list[0]
1181 1206
1182 1207 if hmax_list.any():
1183 1208 hmax_index = hmax_list[-1]+1
1184 1209
1185 1210 x = dataOut.getTimeRange()
1186 1211
1187 1212 thisDatetime = dataOut.datatime
1188 1213
1189 1214 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1190 1215 xlabel = "Local Time"
1191 1216 ylabel = "Phase (degrees)"
1192 1217
1193 1218 update_figfile = False
1194 1219
1195 1220 nplots = len(pairsIndexList)
1196 1221 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1197 1222 phase_beacon = numpy.zeros(len(pairsIndexList))
1198 1223 for i in range(nplots):
1199 1224 pair = dataOut.pairsList[pairsIndexList[i]]
1200 1225 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1201 1226 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1202 1227 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1203 1228 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1204 1229 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1205 1230
1206 1231 if dataOut.beacon_heiIndexList:
1207 1232 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1208 1233 else:
1209 1234 phase_beacon[i] = numpy.average(phase)
1210 1235
1211 1236 if not self.isConfig:
1212 1237
1213 1238 nplots = len(pairsIndexList)
1214 1239
1215 1240 self.setup(id=id,
1216 1241 nplots=nplots,
1217 1242 wintitle=wintitle,
1218 1243 showprofile=showprofile,
1219 1244 show=show)
1220 1245
1221 1246 if timerange != None:
1222 1247 self.timerange = timerange
1223 1248
1224 1249 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1225 1250
1226 1251 if ymin == None: ymin = 0
1227 1252 if ymax == None: ymax = 360
1228 1253
1229 1254 self.FTP_WEI = ftp_wei
1230 1255 self.EXP_CODE = exp_code
1231 1256 self.SUB_EXP_CODE = sub_exp_code
1232 1257 self.PLOT_POS = plot_pos
1233 1258
1234 1259 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1235 1260 self.isConfig = True
1236 1261 self.figfile = figfile
1237 1262 self.xdata = numpy.array([])
1238 1263 self.ydata = numpy.array([])
1239 1264
1240 1265 update_figfile = True
1241 1266
1242 1267 #open file beacon phase
1243 1268 path = '%s%03d' %(self.PREFIX, self.id)
1244 1269 beacon_file = os.path.join(path,'%s.txt'%self.name)
1245 1270 self.filename_phase = os.path.join(figpath,beacon_file)
1246 1271 #self.save_phase(self.filename_phase)
1247 1272
1248 1273
1249 1274 #store data beacon phase
1250 1275 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1251 1276
1252 1277 self.setWinTitle(title)
1253 1278
1254 1279
1255 1280 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1256 1281
1257 1282 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1258 1283
1259 1284 axes = self.axesList[0]
1260 1285
1261 1286 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1262 1287
1263 1288 if len(self.ydata)==0:
1264 1289 self.ydata = phase_beacon.reshape(-1,1)
1265 1290 else:
1266 1291 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1267 1292
1268 1293
1269 1294 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1270 1295 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1271 1296 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1272 1297 XAxisAsTime=True, grid='both'
1273 1298 )
1274 1299
1275 1300 self.draw()
1276 1301
1277 1302 if dataOut.ltctime >= self.xmax:
1278 1303 self.counter_imagwr = wr_period
1279 1304 self.isConfig = False
1280 1305 update_figfile = True
1281 1306
1282 1307 self.save(figpath=figpath,
1283 1308 figfile=figfile,
1284 1309 save=save,
1285 1310 ftp=ftp,
1286 1311 wr_period=wr_period,
1287 1312 thisDatetime=thisDatetime,
1288 1313 update_figfile=update_figfile)
1289 1314
1290 1315 return dataOut
Show file before | Show file after | Hide comments
@@ -1,1285 +1,1285
1 1
2 2 import os
3 3 import time
4 4 import math
5 5 import datetime
6 6 import numpy
7 7 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator #YONG
8 8
9 9 from .jroplot_spectra import RTIPlot, NoisePlot
10 10
11 11 from schainpy.utils import log
12 12 from .plotting_codes import *
13 13
14 14 from schainpy.model.graphics.jroplot_base import Plot, plt
15 15
16 16 import matplotlib.pyplot as plt
17 17 import matplotlib.colors as colors
18 18 from matplotlib.ticker import MultipleLocator
19 19
20 20
21 21 class RTIDPPlot(RTIPlot):
22 22
23 23 '''Plot for RTI Double Pulse Experiment
24 24 '''
25 25
26 26 CODE = 'RTIDP'
27 27 colormap = 'jet'
28 28 plot_name = 'RTI'
29 29 plot_type = 'pcolorbuffer'
30 30
31 31 def setup(self):
32 32 self.xaxis = 'time'
33 33 self.ncols = 1
34 34 self.nrows = 3
35 35 self.nplots = self.nrows
36 36
37 37 self.ylabel = 'Range [km]'
38 38 self.xlabel = 'Time (LT)'
39 39
40 40 self.cb_label = 'Intensity (dB)'
41 41
42 42 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
43 43
44 44 self.titles = ['{} Channel {}'.format(
45 45 self.plot_name.upper(), '0x1'),'{} Channel {}'.format(
46 46 self.plot_name.upper(), '0'),'{} Channel {}'.format(
47 47 self.plot_name.upper(), '1')]
48 48
49 49 def update(self, dataOut):
50 50
51 51 data = {}
52 52 meta = {}
53 53 data['rti'] = dataOut.data_for_RTI_DP
54 54 data['NDP'] = dataOut.NDP
55 55
56 56 return data, meta
57 57
58 58 def plot(self):
59 59
60 60 NDP = self.data['NDP'][-1]
61 61 self.x = self.data.times
62 62 self.y = self.data.yrange[0:NDP]
63 63 self.z = self.data['rti']
64 64 self.z = numpy.ma.masked_invalid(self.z)
65 65
66 66 if self.decimation is None:
67 67 x, y, z = self.fill_gaps(self.x, self.y, self.z)
68 68 else:
69 69 x, y, z = self.fill_gaps(*self.decimate())
70 70
71 71 for n, ax in enumerate(self.axes):
72 72
73 73 self.zmax = self.zmax if self.zmax is not None else numpy.max(
74 74 self.z[1][0,12:40])
75 75 self.zmin = self.zmin if self.zmin is not None else numpy.min(
76 76 self.z[1][0,12:40])
77 77
78 78 if ax.firsttime:
79 79
80 80 if self.zlimits is not None:
81 81 self.zmin, self.zmax = self.zlimits[n]
82 82
83 83 ax.plt = ax.pcolormesh(x, y, z[n].T,
84 84 vmin=self.zmin,
85 85 vmax=self.zmax,
86 86 cmap=plt.get_cmap(self.colormap)
87 87 )
88 88 else:
89 89 #if self.zlimits is not None:
90 90 #self.zmin, self.zmax = self.zlimits[n]
91 91 ax.collections.remove(ax.collections[0])
92 92 ax.plt = ax.pcolormesh(x, y, z[n].T,
93 93 vmin=self.zmin,
94 94 vmax=self.zmax,
95 95 cmap=plt.get_cmap(self.colormap)
96 96 )
97 97
98 98
99 99 class RTILPPlot(RTIPlot):
100 100
101 101 '''
102 102 Plot for RTI Long Pulse
103 103 '''
104 104
105 105 CODE = 'RTILP'
106 106 colormap = 'jet'
107 107 plot_name = 'RTI LP'
108 108 plot_type = 'pcolorbuffer'
109 109
110 110 def setup(self):
111 111 self.xaxis = 'time'
112 112 self.ncols = 1
113 113 self.nrows = 4
114 114 self.nplots = self.nrows
115 115
116 116 self.ylabel = 'Range [km]'
117 117 self.xlabel = 'Time (LT)'
118 118
119 119 self.cb_label = 'Intensity (dB)'
120 120
121 121 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
122 122
123 123
124 124 self.titles = ['{} Channel {}'.format(
125 125 self.plot_name.upper(), '0'),'{} Channel {}'.format(
126 126 self.plot_name.upper(), '1'),'{} Channel {}'.format(
127 127 self.plot_name.upper(), '2'),'{} Channel {}'.format(
128 128 self.plot_name.upper(), '3')]
129 129
130 130
131 131 def update(self, dataOut):
132 132
133 133 data = {}
134 134 meta = {}
135 135 data['rti'] = dataOut.data_for_RTI_LP
136 136 data['NRANGE'] = dataOut.NRANGE
137 137
138 138 return data, meta
139 139
140 140 def plot(self):
141 141
142 142 NRANGE = self.data['NRANGE'][-1]
143 143 self.x = self.data.times
144 144 self.y = self.data.yrange[0:NRANGE]
145 145
146 146 self.z = self.data['rti']
147 147
148 148 self.z = numpy.ma.masked_invalid(self.z)
149 149
150 150 if self.decimation is None:
151 151 x, y, z = self.fill_gaps(self.x, self.y, self.z)
152 152 else:
153 153 x, y, z = self.fill_gaps(*self.decimate())
154 154
155 155 for n, ax in enumerate(self.axes):
156 156
157 157 self.zmax = self.zmax if self.zmax is not None else numpy.max(
158 158 self.z[1][0,12:40])
159 159 self.zmin = self.zmin if self.zmin is not None else numpy.min(
160 160 self.z[1][0,12:40])
161 161
162 162 if ax.firsttime:
163 163
164 164 if self.zlimits is not None:
165 165 self.zmin, self.zmax = self.zlimits[n]
166 166
167 167
168 168 ax.plt = ax.pcolormesh(x, y, z[n].T,
169 169 vmin=self.zmin,
170 170 vmax=self.zmax,
171 171 cmap=plt.get_cmap(self.colormap)
172 172 )
173 173
174 174 else:
175 175 #if self.zlimits is not None:
176 176 #self.zmin, self.zmax = self.zlimits[n]
177 177 ax.collections.remove(ax.collections[0])
178 178 ax.plt = ax.pcolormesh(x, y, z[n].T,
179 179 vmin=self.zmin,
180 180 vmax=self.zmax,
181 181 cmap=plt.get_cmap(self.colormap)
182 182 )
183 183
184 184
185 185 class DenRTIPlot(RTIPlot):
186 186
187 187 '''
188 188 Plot for Den
189 189 '''
190 190
191 191 CODE = 'denrti'
192 192 colormap = 'jet'
193 193
194 194 def setup(self):
195 195 self.xaxis = 'time'
196 196 self.ncols = 1
197 197 self.nrows = self.data.shape(self.CODE)[0]
198 198 self.nplots = self.nrows
199 199
200 200 self.ylabel = 'Range [km]'
201 201 self.xlabel = 'Time (LT)'
202 202
203 203 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
204 204
205 205 if self.CODE == 'denrti':
206 206 self.cb_label = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
207 207
208 208
209 209 self.titles = ['Electron Density RTI']
210 210
211 211 def update(self, dataOut):
212 212
213 213 data = {}
214 214 meta = {}
215 215
216 data['denrti'] = dataOut.DensityFinal
216 data['denrti'] = dataOut.DensityFinal*1.e-6 #To Plot in cm^-3
217 217
218 218 return data, meta
219 219
220 220 def plot(self):
221 221
222 222 self.x = self.data.times
223 223 self.y = self.data.yrange
224 224
225 225 self.z = self.data[self.CODE]
226 226
227 227 self.z = numpy.ma.masked_invalid(self.z)
228 228
229 229 if self.decimation is None:
230 230 x, y, z = self.fill_gaps(self.x, self.y, self.z)
231 231 else:
232 232 x, y, z = self.fill_gaps(*self.decimate())
233 233
234 234 for n, ax in enumerate(self.axes):
235 235
236 236 self.zmax = self.zmax if self.zmax is not None else numpy.max(
237 237 self.z[n])
238 238 self.zmin = self.zmin if self.zmin is not None else numpy.min(
239 239 self.z[n])
240 240
241 241 if ax.firsttime:
242 242
243 243 if self.zlimits is not None:
244 244 self.zmin, self.zmax = self.zlimits[n]
245 245 if numpy.log10(self.zmin)<0:
246 246 self.zmin=1
247 247 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
248 248 vmin=self.zmin,
249 249 vmax=self.zmax,
250 250 cmap=self.cmaps[n],
251 251 norm=colors.LogNorm()
252 252 )
253 253
254 254 else:
255 255 if self.zlimits is not None:
256 256 self.zmin, self.zmax = self.zlimits[n]
257 257 ax.collections.remove(ax.collections[0])
258 258 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
259 259 vmin=self.zmin,
260 260 vmax=self.zmax,
261 261 cmap=self.cmaps[n],
262 262 norm=colors.LogNorm()
263 263 )
264 264
265 265
266 266 class ETempRTIPlot(RTIPlot):
267 267
268 268 '''
269 269 Plot for Electron Temperature
270 270 '''
271 271
272 272 CODE = 'ETemp'
273 273 colormap = 'jet'
274 274
275 275 def setup(self):
276 276 self.xaxis = 'time'
277 277 self.ncols = 1
278 278 self.nrows = self.data.shape(self.CODE)[0]
279 279 self.nplots = self.nrows
280 280
281 281 self.ylabel = 'Range [km]'
282 282 self.xlabel = 'Time (LT)'
283 283 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
284 284 if self.CODE == 'ETemp':
285 285 self.cb_label = 'Electron Temperature (K)'
286 286 self.titles = ['Electron Temperature RTI']
287 287 if self.CODE == 'ITemp':
288 288 self.cb_label = 'Ion Temperature (K)'
289 289 self.titles = ['Ion Temperature RTI']
290 290 if self.CODE == 'HeFracLP':
291 291 self.cb_label='He+ Fraction'
292 292 self.titles = ['He+ Fraction RTI']
293 293 self.zmax=0.16
294 294 if self.CODE== 'HFracLP':
295 295 self.cb_label='H+ Fraction'
296 296 self.titles = ['H+ Fraction RTI']
297 297
298 298 def update(self, dataOut):
299 299
300 300 data = {}
301 301 meta = {}
302 302
303 303 data['ETemp'] = dataOut.ElecTempFinal
304 304
305 305 return data, meta
306 306
307 307 def plot(self):
308 308
309 309 self.x = self.data.times
310 310 self.y = self.data.yrange
311 311
312 312
313 313 self.z = self.data[self.CODE]
314 314
315 315 self.z = numpy.ma.masked_invalid(self.z)
316 316
317 317 if self.decimation is None:
318 318 x, y, z = self.fill_gaps(self.x, self.y, self.z)
319 319 else:
320 320 x, y, z = self.fill_gaps(*self.decimate())
321 321
322 322 for n, ax in enumerate(self.axes):
323 323
324 324 self.zmax = self.zmax if self.zmax is not None else numpy.max(
325 325 self.z[n])
326 326 self.zmin = self.zmin if self.zmin is not None else numpy.min(
327 327 self.z[n])
328 328
329 329 if ax.firsttime:
330 330
331 331 if self.zlimits is not None:
332 332 self.zmin, self.zmax = self.zlimits[n]
333 333
334 334 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
335 335 vmin=self.zmin,
336 336 vmax=self.zmax,
337 337 cmap=self.cmaps[n]
338 338 )
339 339 #plt.tight_layout()
340 340
341 341 else:
342 342 if self.zlimits is not None:
343 343 self.zmin, self.zmax = self.zlimits[n]
344 344 ax.collections.remove(ax.collections[0])
345 345 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
346 346 vmin=self.zmin,
347 347 vmax=self.zmax,
348 348 cmap=self.cmaps[n]
349 349 )
350 350
351 351
352 352 class ITempRTIPlot(ETempRTIPlot):
353 353
354 354 '''
355 355 Plot for Ion Temperature
356 356 '''
357 357
358 358 CODE = 'ITemp'
359 359 colormap = 'jet'
360 360 plot_name = 'Ion Temperature'
361 361
362 362 def update(self, dataOut):
363 363
364 364 data = {}
365 365 meta = {}
366 366
367 367 data['ITemp'] = dataOut.IonTempFinal
368 368
369 369 return data, meta
370 370
371 371
372 372 class HFracRTIPlot(ETempRTIPlot):
373 373
374 374 '''
375 375 Plot for H+ LP
376 376 '''
377 377
378 378 CODE = 'HFracLP'
379 379 colormap = 'jet'
380 380 plot_name = 'H+ Frac'
381 381
382 382 def update(self, dataOut):
383 383
384 384 data = {}
385 385 meta = {}
386 386 data['HFracLP'] = dataOut.PhyFinal
387 387
388 388 return data, meta
389 389
390 390
391 391 class HeFracRTIPlot(ETempRTIPlot):
392 392
393 393 '''
394 394 Plot for He+ LP
395 395 '''
396 396
397 397 CODE = 'HeFracLP'
398 398 colormap = 'jet'
399 399 plot_name = 'He+ Frac'
400 400
401 401 def update(self, dataOut):
402 402
403 403 data = {}
404 404 meta = {}
405 405 data['HeFracLP'] = dataOut.PheFinal
406 406
407 407 return data, meta
408 408
409 409
410 410 class TempsDPPlot(Plot):
411 411 '''
412 412 Plot for Electron - Ion Temperatures
413 413 '''
414 414
415 415 CODE = 'tempsDP'
416 416 #plot_name = 'Temperatures'
417 417 plot_type = 'scatterbuffer'
418 418
419 419 def setup(self):
420 420
421 421 self.ncols = 1
422 422 self.nrows = 1
423 423 self.nplots = 1
424 424 self.ylabel = 'Range [km]'
425 425 self.xlabel = 'Temperature (K)'
426 426 self.titles = ['Electron/Ion Temperatures']
427 427 self.width = 3.5
428 428 self.height = 5.5
429 429 self.colorbar = False
430 430 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
431 431
432 432 def update(self, dataOut):
433 433 data = {}
434 434 meta = {}
435 435
436 436 data['Te'] = dataOut.te2
437 437 data['Ti'] = dataOut.ti2
438 438 data['Te_error'] = dataOut.ete2
439 439 data['Ti_error'] = dataOut.eti2
440 440
441 441 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
442 442
443 443 return data, meta
444 444
445 445 def plot(self):
446 446
447 447 y = self.data.yrange
448 448
449 449 self.xmin = -100
450 450 self.xmax = 5000
451 451
452 452 ax = self.axes[0]
453 453
454 454 data = self.data[-1]
455 455
456 456 Te = data['Te']
457 457 Ti = data['Ti']
458 458 errTe = data['Te_error']
459 459 errTi = data['Ti_error']
460 460
461 461 if ax.firsttime:
462 462 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
463 463 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
464 464 plt.legend(loc='lower right')
465 465 self.ystep_given = 50
466 466 ax.yaxis.set_minor_locator(MultipleLocator(15))
467 467 ax.grid(which='minor')
468 468
469 469 else:
470 470 self.clear_figures()
471 471 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
472 472 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
473 473 plt.legend(loc='lower right')
474 474 ax.yaxis.set_minor_locator(MultipleLocator(15))
475 475
476 476
477 477 class TempsHPPlot(Plot):
478 478 '''
479 479 Plot for Temperatures Hybrid Experiment
480 480 '''
481 481
482 482 CODE = 'temps_LP'
483 483 #plot_name = 'Temperatures'
484 484 plot_type = 'scatterbuffer'
485 485
486 486
487 487 def setup(self):
488 488
489 489 self.ncols = 1
490 490 self.nrows = 1
491 491 self.nplots = 1
492 492 self.ylabel = 'Range [km]'
493 493 self.xlabel = 'Temperature (K)'
494 494 self.titles = ['Electron/Ion Temperatures']
495 495 self.width = 3.5
496 496 self.height = 6.5
497 497 self.colorbar = False
498 498 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
499 499
500 500 def update(self, dataOut):
501 501 data = {}
502 502 meta = {}
503 503
504 504
505 505 data['Te'] = numpy.concatenate((dataOut.te2[:dataOut.cut],dataOut.te[dataOut.cut:]))
506 506 data['Ti'] = numpy.concatenate((dataOut.ti2[:dataOut.cut],dataOut.ti[dataOut.cut:]))
507 507 data['Te_error'] = numpy.concatenate((dataOut.ete2[:dataOut.cut],dataOut.ete[dataOut.cut:]))
508 508 data['Ti_error'] = numpy.concatenate((dataOut.eti2[:dataOut.cut],dataOut.eti[dataOut.cut:]))
509 509
510 510 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
511 511
512 512 return data, meta
513 513
514 514 def plot(self):
515 515
516 516
517 517 self.y = self.data.yrange
518 518 self.xmin = -100
519 519 self.xmax = 4500
520 520 ax = self.axes[0]
521 521
522 522 data = self.data[-1]
523 523
524 524 Te = data['Te']
525 525 Ti = data['Ti']
526 526 errTe = data['Te_error']
527 527 errTi = data['Ti_error']
528 528
529 529 if ax.firsttime:
530 530
531 531 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
532 532 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
533 533 plt.legend(loc='lower right')
534 534 self.ystep_given = 200
535 535 ax.yaxis.set_minor_locator(MultipleLocator(15))
536 536 ax.grid(which='minor')
537 537
538 538 else:
539 539 self.clear_figures()
540 540 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
541 541 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
542 542 plt.legend(loc='lower right')
543 543 ax.yaxis.set_minor_locator(MultipleLocator(15))
544 544 ax.grid(which='minor')
545 545
546 546
547 547 class FracsHPPlot(Plot):
548 548 '''
549 549 Plot for Composition LP
550 550 '''
551 551
552 552 CODE = 'fracs_LP'
553 553 plot_type = 'scatterbuffer'
554 554
555 555
556 556 def setup(self):
557 557
558 558 self.ncols = 1
559 559 self.nrows = 1
560 560 self.nplots = 1
561 561 self.ylabel = 'Range [km]'
562 562 self.xlabel = 'Frac'
563 563 self.titles = ['Composition']
564 564 self.width = 3.5
565 565 self.height = 6.5
566 566 self.colorbar = False
567 567 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
568 568
569 569 def update(self, dataOut):
570 570 data = {}
571 571 meta = {}
572 572
573 573 #aux_nan=numpy.zeros(dataOut.cut,'float32')
574 574 #aux_nan[:]=numpy.nan
575 575 #data['ph'] = numpy.concatenate((aux_nan,dataOut.ph[dataOut.cut:]))
576 576 #data['eph'] = numpy.concatenate((aux_nan,dataOut.eph[dataOut.cut:]))
577 577
578 578 data['ph'] = dataOut.ph[dataOut.cut:]
579 579 data['eph'] = dataOut.eph[dataOut.cut:]
580 580 data['phe'] = dataOut.phe[dataOut.cut:]
581 581 data['ephe'] = dataOut.ephe[dataOut.cut:]
582 582
583 583 data['cut'] = dataOut.cut
584 584
585 585 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
586 586
587 587
588 588 return data, meta
589 589
590 590 def plot(self):
591 591
592 592 data = self.data[-1]
593 593
594 594 ph = data['ph']
595 595 eph = data['eph']
596 596 phe = data['phe']
597 597 ephe = data['ephe']
598 598 cut = data['cut']
599 599 self.y = self.data.yrange
600 600
601 601 self.xmin = 0
602 602 self.xmax = 1
603 603 ax = self.axes[0]
604 604
605 605 if ax.firsttime:
606 606
607 607 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
608 608 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
609 609 plt.legend(loc='lower right')
610 610 self.xstep_given = 0.2
611 611 self.ystep_given = 200
612 612 ax.yaxis.set_minor_locator(MultipleLocator(15))
613 613 ax.grid(which='minor')
614 614
615 615 else:
616 616 self.clear_figures()
617 617 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
618 618 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
619 619 plt.legend(loc='lower right')
620 620 ax.yaxis.set_minor_locator(MultipleLocator(15))
621 621
622 622 class EDensityPlot(Plot):
623 623 '''
624 624 Plot for electron density
625 625 '''
626 626
627 627 CODE = 'den'
628 628 #plot_name = 'Electron Density'
629 629 plot_type = 'scatterbuffer'
630 630
631 631 def setup(self):
632 632
633 633 self.ncols = 1
634 634 self.nrows = 1
635 635 self.nplots = 1
636 636 self.ylabel = 'Range [km]'
637 637 self.xlabel = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
638 638 self.titles = ['Electron Density']
639 639 self.width = 3.5
640 640 self.height = 5.5
641 641 self.colorbar = False
642 642 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
643 643
644 644 def update(self, dataOut):
645 645 data = {}
646 646 meta = {}
647 647
648 648 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
649 649 data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS]
650 650 data['den_error'] = dataOut.sdp2[:dataOut.NSHTS]
651 651 #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS]
652 652
653 653 data['NSHTS'] = dataOut.NSHTS
654 654
655 655 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
656 656
657 657 return data, meta
658 658
659 659 def plot(self):
660 660
661 661 y = self.data.yrange
662 662
663 663 self.xmin = 1e3
664 664 self.xmax = 1e7
665 665
666 666 ax = self.axes[0]
667 667
668 668 data = self.data[-1]
669 669
670 670 DenPow = data['den_power']
671 671 DenFar = data['den_Faraday']
672 672 errDenPow = data['den_error']
673 673 #errFaraday = data['err_Faraday']
674 674
675 675 NSHTS = data['NSHTS']
676 676
677 677 if self.CODE == 'denLP':
678 678 DenPowLP = data['den_LP']
679 679 errDenPowLP = data['den_LP_error']
680 680 cut = data['cut']
681 681
682 682 if ax.firsttime:
683 683 self.autoxticks=False
684 684 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
685 685 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2)
686 686 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
687 687 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2)
688 688
689 689 if self.CODE=='denLP':
690 690 ax.errorbar(DenPowLP[cut:], y[cut:], xerr=errDenPowLP[cut:], fmt='r^-',elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
691 691
692 692 plt.legend(loc='upper left',fontsize=8.5)
693 693 #plt.legend(loc='lower left',fontsize=8.5)
694 694 ax.set_xscale("log", nonposx='clip')
695 695 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
696 696 self.ystep_given=100
697 697 if self.CODE=='denLP':
698 698 self.ystep_given=200
699 699 ax.set_yticks(grid_y_ticks,minor=True)
700 700 ax.grid(which='minor')
701 701
702 702 else:
703 703 dataBefore = self.data[-2]
704 704 DenPowBefore = dataBefore['den_power']
705 705 self.clear_figures()
706 706 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
707 707 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2)
708 708 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
709 709 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2)
710 710 ax.errorbar(DenPowBefore, y[:NSHTS], elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed")
711 711
712 712 if self.CODE=='denLP':
713 713 ax.errorbar(DenPowLP[cut:], y[cut:], fmt='r^-', xerr=errDenPowLP[cut:],elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
714 714
715 715 ax.set_xscale("log", nonposx='clip')
716 716 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
717 717 ax.set_yticks(grid_y_ticks,minor=True)
718 718 ax.grid(which='minor')
719 719 plt.legend(loc='upper left',fontsize=8.5)
720 720 #plt.legend(loc='lower left',fontsize=8.5)
721 721
722 722 class FaradayAnglePlot(Plot):
723 723 '''
724 724 Plot for electron density
725 725 '''
726 726
727 727 CODE = 'angle'
728 728 plot_name = 'Faraday Angle'
729 729 plot_type = 'scatterbuffer'
730 730
731 731 def setup(self):
732 732
733 733 self.ncols = 1
734 734 self.nrows = 1
735 735 self.nplots = 1
736 736 self.ylabel = 'Range [km]'
737 737 self.xlabel = 'Faraday Angle (º)'
738 738 self.titles = ['Electron Density']
739 739 self.width = 3.5
740 740 self.height = 5.5
741 741 self.colorbar = False
742 742 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
743 743
744 744 def update(self, dataOut):
745 745 data = {}
746 746 meta = {}
747 747
748 748 data['angle'] = numpy.degrees(dataOut.phi)
749 749 #'''
750 750 print(dataOut.phi_uwrp)
751 751 print(data['angle'])
752 752 exit(1)
753 753 #'''
754 754 data['dphi'] = dataOut.dphi_uc*10
755 755 #print(dataOut.dphi)
756 756
757 757 #data['NSHTS'] = dataOut.NSHTS
758 758
759 759 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
760 760
761 761 return data, meta
762 762
763 763 def plot(self):
764 764
765 765 data = self.data[-1]
766 766 self.x = data[self.CODE]
767 767 dphi = data['dphi']
768 768 self.y = self.data.yrange
769 769 self.xmin = -360#-180
770 770 self.xmax = 360#180
771 771 ax = self.axes[0]
772 772
773 773 if ax.firsttime:
774 774 self.autoxticks=False
775 775 #if self.CODE=='den':
776 776 ax.plot(self.x, self.y,marker='o',color='g',linewidth=1.0,markersize=2)
777 777 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
778 778
779 779 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
780 780 self.ystep_given=100
781 781 if self.CODE=='denLP':
782 782 self.ystep_given=200
783 783 ax.set_yticks(grid_y_ticks,minor=True)
784 784 ax.grid(which='minor')
785 785 #plt.tight_layout()
786 786 else:
787 787
788 788 self.clear_figures()
789 789 #if self.CODE=='den':
790 790 #print(numpy.shape(self.x))
791 791 ax.plot(self.x, self.y, marker='o',color='g',linewidth=1.0, markersize=2)
792 792 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
793 793
794 794 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
795 795 ax.set_yticks(grid_y_ticks,minor=True)
796 796 ax.grid(which='minor')
797 797
798 798 class EDensityHPPlot(EDensityPlot):
799 799
800 800 '''
801 801 Plot for Electron Density Hybrid Experiment
802 802 '''
803 803
804 804 CODE = 'denLP'
805 805 plot_name = 'Electron Density'
806 806 plot_type = 'scatterbuffer'
807 807
808 808 def update(self, dataOut):
809 809 data = {}
810 810 meta = {}
811 811
812 812 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
813 813 data['den_Faraday']=dataOut.dphi[:dataOut.NSHTS]
814 814 data['den_error']=dataOut.sdp2[:dataOut.NSHTS]
815 815 data['den_LP']=dataOut.ne[:dataOut.NACF]
816 816 data['den_LP_error']=dataOut.ene[:dataOut.NACF]*dataOut.ne[:dataOut.NACF]*0.434
817 817 #self.ene=10**dataOut.ene[:dataOut.NACF]
818 818 data['NSHTS']=dataOut.NSHTS
819 819 data['cut']=dataOut.cut
820 820
821 821 return data, meta
822 822
823 823
824 824 class ACFsPlot(Plot):
825 825 '''
826 826 Plot for ACFs Double Pulse Experiment
827 827 '''
828 828
829 829 CODE = 'acfs'
830 830 #plot_name = 'ACF'
831 831 plot_type = 'scatterbuffer'
832 832
833 833
834 834 def setup(self):
835 835 self.ncols = 1
836 836 self.nrows = 1
837 837 self.nplots = 1
838 838 self.ylabel = 'Range [km]'
839 839 self.xlabel = 'Lag (ms)'
840 840 self.titles = ['ACFs']
841 841 self.width = 3.5
842 842 self.height = 5.5
843 843 self.colorbar = False
844 844 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
845 845
846 846 def update(self, dataOut):
847 847 data = {}
848 848 meta = {}
849 849
850 850 data['ACFs'] = dataOut.acfs_to_plot
851 851 data['ACFs_error'] = dataOut.acfs_error_to_plot
852 852 data['lags'] = dataOut.lags_to_plot
853 853 data['Lag_contaminated_1'] = dataOut.x_igcej_to_plot
854 854 data['Lag_contaminated_2'] = dataOut.x_ibad_to_plot
855 855 data['Height_contaminated_1'] = dataOut.y_igcej_to_plot
856 856 data['Height_contaminated_2'] = dataOut.y_ibad_to_plot
857 857
858 858 meta['yrange'] = numpy.array([])
859 859 #meta['NSHTS'] = dataOut.NSHTS
860 860 #meta['DPL'] = dataOut.DPL
861 861 data['NSHTS'] = dataOut.NSHTS #This is metadata
862 862 data['DPL'] = dataOut.DPL #This is metadata
863 863
864 864 return data, meta
865 865
866 866 def plot(self):
867 867
868 868 data = self.data[-1]
869 869 #NSHTS = self.meta['NSHTS']
870 870 #DPL = self.meta['DPL']
871 871 NSHTS = data['NSHTS'] #This is metadata
872 872 DPL = data['DPL'] #This is metadata
873 873
874 874 lags = data['lags']
875 875 ACFs = data['ACFs']
876 876 errACFs = data['ACFs_error']
877 877 BadLag1 = data['Lag_contaminated_1']
878 878 BadLag2 = data['Lag_contaminated_2']
879 879 BadHei1 = data['Height_contaminated_1']
880 880 BadHei2 = data['Height_contaminated_2']
881 881
882 882 self.xmin = 0.0
883 883 self.xmax = 2.0
884 884 self.y = ACFs
885 885
886 886 ax = self.axes[0]
887 887
888 888 if ax.firsttime:
889 889
890 890 for i in range(NSHTS):
891 891 x_aux = numpy.isfinite(lags[i,:])
892 892 y_aux = numpy.isfinite(ACFs[i,:])
893 893 yerr_aux = numpy.isfinite(errACFs[i,:])
894 894 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
895 895 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
896 896 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
897 897 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
898 898 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
899 899 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',marker='o',linewidth=1.0,markersize=2)
900 900 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
901 901 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
902 902
903 903 self.xstep_given = (self.xmax-self.xmin)/(DPL-1)
904 904 self.ystep_given = 50
905 905 ax.yaxis.set_minor_locator(MultipleLocator(15))
906 906 ax.grid(which='minor')
907 907
908 908 else:
909 909 self.clear_figures()
910 910 for i in range(NSHTS):
911 911 x_aux = numpy.isfinite(lags[i,:])
912 912 y_aux = numpy.isfinite(ACFs[i,:])
913 913 yerr_aux = numpy.isfinite(errACFs[i,:])
914 914 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
915 915 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
916 916 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
917 917 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
918 918 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
919 919 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],linewidth=1.0,markersize=2,color='b',marker='o')
920 920 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
921 921 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
922 922 ax.yaxis.set_minor_locator(MultipleLocator(15))
923 923
924 924 class ACFsLPPlot(Plot):
925 925 '''
926 926 Plot for ACFs Double Pulse Experiment
927 927 '''
928 928
929 929 CODE = 'acfs_LP'
930 930 #plot_name = 'ACF'
931 931 plot_type = 'scatterbuffer'
932 932
933 933
934 934 def setup(self):
935 935 self.ncols = 1
936 936 self.nrows = 1
937 937 self.nplots = 1
938 938 self.ylabel = 'Range [km]'
939 939 self.xlabel = 'Lag (ms)'
940 940 self.titles = ['ACFs']
941 941 self.width = 3.5
942 942 self.height = 5.5
943 943 self.colorbar = False
944 944 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
945 945
946 946 def update(self, dataOut):
947 947 data = {}
948 948 meta = {}
949 949
950 950 aux=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
951 951 errors=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
952 952 lags_LP_to_plot=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
953 953
954 954 for i in range(dataOut.NACF):
955 955 for j in range(dataOut.IBITS):
956 956 if numpy.abs(dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0])<1.0:
957 957 aux[i,j]=dataOut.output_LP_integrated.real[j,i,0]/dataOut.output_LP_integrated.real[0,i,0]
958 958 aux[i,j]=max(min(aux[i,j],1.0),-1.0)*dataOut.DH+dataOut.heightList[i]
959 959 lags_LP_to_plot[i,j]=dataOut.lags_LP[j]
960 960 errors[i,j]=dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0]*dataOut.DH
961 961 else:
962 962 aux[i,j]=numpy.nan
963 963 lags_LP_to_plot[i,j]=numpy.nan
964 964 errors[i,j]=numpy.nan
965 965
966 966 data['ACFs'] = aux
967 967 data['ACFs_error'] = errors
968 968 data['lags'] = lags_LP_to_plot
969 969
970 970 meta['yrange'] = numpy.array([])
971 971 #meta['NACF'] = dataOut.NACF
972 972 #meta['NLAG'] = dataOut.NLAG
973 973 data['NACF'] = dataOut.NACF #This is metadata
974 974 data['NLAG'] = dataOut.NLAG #This is metadata
975 975
976 976 return data, meta
977 977
978 978 def plot(self):
979 979
980 980 data = self.data[-1]
981 981 #NACF = self.meta['NACF']
982 982 #NLAG = self.meta['NLAG']
983 983 NACF = data['NACF'] #This is metadata
984 984 NLAG = data['NLAG'] #This is metadata
985 985
986 986 lags = data['lags']
987 987 ACFs = data['ACFs']
988 988 errACFs = data['ACFs_error']
989 989
990 990 self.xmin = 0.0
991 991 self.xmax = 1.5
992 992
993 993 self.y = ACFs
994 994
995 995 ax = self.axes[0]
996 996
997 997 if ax.firsttime:
998 998
999 999 for i in range(NACF):
1000 1000 x_aux = numpy.isfinite(lags[i,:])
1001 1001 y_aux = numpy.isfinite(ACFs[i,:])
1002 1002 yerr_aux = numpy.isfinite(errACFs[i,:])
1003 1003
1004 1004 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1005 1005 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1006 1006
1007 1007 #self.xstep_given = (self.xmax-self.xmin)/(self.data.NLAG-1)
1008 1008 self.xstep_given=0.3
1009 1009 self.ystep_given = 200
1010 1010 ax.yaxis.set_minor_locator(MultipleLocator(15))
1011 1011 ax.grid(which='minor')
1012 1012
1013 1013 else:
1014 1014 self.clear_figures()
1015 1015
1016 1016 for i in range(NACF):
1017 1017 x_aux = numpy.isfinite(lags[i,:])
1018 1018 y_aux = numpy.isfinite(ACFs[i,:])
1019 1019 yerr_aux = numpy.isfinite(errACFs[i,:])
1020 1020
1021 1021 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1022 1022 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1023 1023
1024 1024 ax.yaxis.set_minor_locator(MultipleLocator(15))
1025 1025
1026 1026
1027 1027 class CrossProductsPlot(Plot):
1028 1028 '''
1029 1029 Plot for cross products
1030 1030 '''
1031 1031
1032 1032 CODE = 'crossprod'
1033 1033 plot_name = 'Cross Products'
1034 1034 plot_type = 'scatterbuffer'
1035 1035
1036 1036 def setup(self):
1037 1037
1038 1038 self.ncols = 3
1039 1039 self.nrows = 1
1040 1040 self.nplots = 3
1041 1041 self.ylabel = 'Range [km]'
1042 1042 self.titles = []
1043 1043 self.width = 3.5*self.nplots
1044 1044 self.height = 5.5
1045 1045 self.colorbar = False
1046 1046 self.plots_adjust.update({'wspace':.3, 'left': 0.12, 'right': 0.92, 'bottom': 0.1})
1047 1047
1048 1048
1049 1049 def update(self, dataOut):
1050 1050
1051 1051 data = {}
1052 1052 meta = {}
1053 1053
1054 1054 data['crossprod'] = dataOut.crossprods
1055 1055 data['NDP'] = dataOut.NDP
1056 1056
1057 1057 return data, meta
1058 1058
1059 1059 def plot(self):
1060 1060
1061 1061 NDP = self.data['NDP'][-1]
1062 1062 x = self.data['crossprod'][:,-1,:,:,:,:]
1063 1063 y = self.data.yrange[0:NDP]
1064 1064
1065 1065 for n, ax in enumerate(self.axes):
1066 1066
1067 1067 self.xmin=numpy.min(numpy.concatenate((x[n][0,20:30,0,0],x[n][1,20:30,0,0],x[n][2,20:30,0,0],x[n][3,20:30,0,0])))
1068 1068 self.xmax=numpy.max(numpy.concatenate((x[n][0,20:30,0,0],x[n][1,20:30,0,0],x[n][2,20:30,0,0],x[n][3,20:30,0,0])))
1069 1069
1070 1070 if ax.firsttime:
1071 1071
1072 1072 self.autoxticks=False
1073 1073 if n==0:
1074 1074 label1='kax'
1075 1075 label2='kay'
1076 1076 label3='kbx'
1077 1077 label4='kby'
1078 1078 self.xlimits=[(self.xmin,self.xmax)]
1079 1079 elif n==1:
1080 1080 label1='kax2'
1081 1081 label2='kay2'
1082 1082 label3='kbx2'
1083 1083 label4='kby2'
1084 1084 self.xlimits.append((self.xmin,self.xmax))
1085 1085 elif n==2:
1086 1086 label1='kaxay'
1087 1087 label2='kbxby'
1088 1088 label3='kaxbx'
1089 1089 label4='kaxby'
1090 1090 self.xlimits.append((self.xmin,self.xmax))
1091 1091
1092 1092 ax.plotline1 = ax.plot(x[n][0,:,0,0], y, color='r',linewidth=2.0, label=label1)
1093 1093 ax.plotline2 = ax.plot(x[n][1,:,0,0], y, color='k',linewidth=2.0, label=label2)
1094 1094 ax.plotline3 = ax.plot(x[n][2,:,0,0], y, color='b',linewidth=2.0, label=label3)
1095 1095 ax.plotline4 = ax.plot(x[n][3,:,0,0], y, color='m',linewidth=2.0, label=label4)
1096 1096 ax.legend(loc='upper right')
1097 1097 ax.set_xlim(self.xmin, self.xmax)
1098 1098 self.titles.append('{}'.format(self.plot_name.upper()))
1099 1099
1100 1100 else:
1101 1101
1102 1102 if n==0:
1103 1103 self.xlimits=[(self.xmin,self.xmax)]
1104 1104 else:
1105 1105 self.xlimits.append((self.xmin,self.xmax))
1106 1106
1107 1107 ax.set_xlim(self.xmin, self.xmax)
1108 1108
1109 1109 ax.plotline1[0].set_data(x[n][0,:,0,0],y)
1110 1110 ax.plotline2[0].set_data(x[n][1,:,0,0],y)
1111 1111 ax.plotline3[0].set_data(x[n][2,:,0,0],y)
1112 1112 ax.plotline4[0].set_data(x[n][3,:,0,0],y)
1113 1113 self.titles.append('{}'.format(self.plot_name.upper()))
1114 1114
1115 1115
1116 1116 class CrossProductsLPPlot(Plot):
1117 1117 '''
1118 1118 Plot for cross products LP
1119 1119 '''
1120 1120
1121 1121 CODE = 'crossprodslp'
1122 1122 plot_name = 'Cross Products LP'
1123 1123 plot_type = 'scatterbuffer'
1124 1124
1125 1125
1126 1126 def setup(self):
1127 1127
1128 1128 self.ncols = 2
1129 1129 self.nrows = 1
1130 1130 self.nplots = 2
1131 1131 self.ylabel = 'Range [km]'
1132 1132 self.xlabel = 'dB'
1133 1133 self.width = 3.5*self.nplots
1134 1134 self.height = 5.5
1135 1135 self.colorbar = False
1136 1136 self.titles = []
1137 1137 self.plots_adjust.update({'wspace': .8 ,'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1138 1138
1139 1139 def update(self, dataOut):
1140 1140 data = {}
1141 1141 meta = {}
1142 1142
1143 1143 data['crossprodslp'] = 10*numpy.log10(numpy.abs(dataOut.output_LP))
1144 1144
1145 1145 data['NRANGE'] = dataOut.NRANGE #This is metadata
1146 1146 data['NLAG'] = dataOut.NLAG #This is metadata
1147 1147
1148 1148 return data, meta
1149 1149
1150 1150 def plot(self):
1151 1151
1152 1152 NRANGE = self.data['NRANGE'][-1]
1153 1153 NLAG = self.data['NLAG'][-1]
1154 1154
1155 1155 x = self.data[self.CODE][:,-1,:,:]
1156 1156 self.y = self.data.yrange[0:NRANGE]
1157 1157
1158 1158 label_array=numpy.array(['lag '+ str(x) for x in range(NLAG)])
1159 1159 color_array=['r','k','g','b','c','m','y','orange','steelblue','purple','peru','darksalmon','grey','limegreen','olive','midnightblue']
1160 1160
1161 1161
1162 1162 for n, ax in enumerate(self.axes):
1163 1163
1164 1164 self.xmin=28#30
1165 1165 self.xmax=70#70
1166 1166 #self.xmin=numpy.min(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1167 1167 #self.xmax=numpy.max(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1168 1168
1169 1169 if ax.firsttime:
1170 1170
1171 1171 self.autoxticks=False
1172 1172 if n == 0:
1173 1173 self.plotline_array=numpy.zeros((2,NLAG),dtype=object)
1174 1174
1175 1175 for i in range(NLAG):
1176 1176 self.plotline_array[n,i], = ax.plot(x[i,:,n], self.y, color=color_array[i],linewidth=1.0, label=label_array[i])
1177 1177
1178 1178 ax.legend(loc='upper right')
1179 1179 ax.set_xlim(self.xmin, self.xmax)
1180 1180 if n==0:
1181 1181 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1182 1182 if n==1:
1183 1183 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1184 1184 else:
1185 1185 for i in range(NLAG):
1186 1186 self.plotline_array[n,i].set_data(x[i,:,n],self.y)
1187 1187
1188 1188 if n==0:
1189 1189 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1190 1190 if n==1:
1191 1191 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1192 1192
1193 1193
1194 1194 class NoiseDPPlot(NoisePlot):
1195 1195 '''
1196 1196 Plot for noise Double Pulse
1197 1197 '''
1198 1198
1199 1199 CODE = 'noise'
1200 1200 #plot_name = 'Noise'
1201 1201 #plot_type = 'scatterbuffer'
1202 1202
1203 1203 def update(self, dataOut):
1204 1204
1205 1205 data = {}
1206 1206 meta = {}
1207 1207 data['noise'] = 10*numpy.log10(dataOut.noise_final)
1208 1208
1209 1209 return data, meta
1210 1210
1211 1211
1212 1212 class XmitWaveformPlot(Plot):
1213 1213 '''
1214 1214 Plot for xmit waveform
1215 1215 '''
1216 1216
1217 1217 CODE = 'xmit'
1218 1218 plot_name = 'Xmit Waveform'
1219 1219 plot_type = 'scatterbuffer'
1220 1220
1221 1221
1222 1222 def setup(self):
1223 1223
1224 1224 self.ncols = 1
1225 1225 self.nrows = 1
1226 1226 self.nplots = 1
1227 1227 self.ylabel = ''
1228 1228 self.xlabel = 'Number of Lag'
1229 1229 self.width = 5.5
1230 1230 self.height = 3.5
1231 1231 self.colorbar = False
1232 1232 self.plots_adjust.update({'right': 0.85 })
1233 1233 self.titles = [self.plot_name]
1234 1234 #self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1235 1235
1236 1236 #if not self.titles:
1237 1237 #self.titles = self.data.parameters \
1238 1238 #if self.data.parameters else ['{}'.format(self.plot_name.upper())]
1239 1239
1240 1240 def update(self, dataOut):
1241 1241
1242 1242 data = {}
1243 1243 meta = {}
1244 1244
1245 1245 y_1=numpy.arctan2(dataOut.output_LP[:,0,2].imag,dataOut.output_LP[:,0,2].real)* 180 / (numpy.pi*10)
1246 1246 y_2=numpy.abs(dataOut.output_LP[:,0,2])
1247 1247 norm=numpy.max(y_2)
1248 1248 norm=max(norm,0.1)
1249 1249 y_2=y_2/norm
1250 1250
1251 1251 meta['yrange'] = numpy.array([])
1252 1252
1253 1253 data['xmit'] = numpy.vstack((y_1,y_2))
1254 1254 data['NLAG'] = dataOut.NLAG
1255 1255
1256 1256 return data, meta
1257 1257
1258 1258 def plot(self):
1259 1259
1260 1260 data = self.data[-1]
1261 1261 NLAG = data['NLAG']
1262 1262 x = numpy.arange(0,NLAG,1,'float32')
1263 1263 y = data['xmit']
1264 1264
1265 1265 self.xmin = 0
1266 1266 self.xmax = NLAG-1
1267 1267 self.ymin = -1.0
1268 1268 self.ymax = 1.0
1269 1269 ax = self.axes[0]
1270 1270
1271 1271 if ax.firsttime:
1272 1272 ax.plotline0=ax.plot(x,y[0,:],color='blue')
1273 1273 ax.plotline1=ax.plot(x,y[1,:],color='red')
1274 1274 secax=ax.secondary_xaxis(location=0.5)
1275 1275 secax.xaxis.tick_bottom()
1276 1276 secax.tick_params( labelleft=False, labeltop=False,
1277 1277 labelright=False, labelbottom=False)
1278 1278
1279 1279 self.xstep_given = 3
1280 1280 self.ystep_given = .25
1281 1281 secax.set_xticks(numpy.linspace(self.xmin, self.xmax, 6)) #only works on matplotlib.version>3.2
1282 1282
1283 1283 else:
1284 1284 ax.plotline0[0].set_data(x,y[0,:])
1285 1285 ax.plotline1[0].set_data(x,y[1,:])
Show file before | Show file after | Hide comments
@@ -1,251 +1,252
1 1 '''
2 2 Base clases to create Processing units and operations, the MPDecorator
3 3 must be used in plotting and writing operations to allow to run as an
4 4 external process.
5 5 '''
6 6
7 7 import os
8 8 import inspect
9 9 import zmq
10 10 import time
11 11 import pickle
12 12 import traceback
13 13 from threading import Thread
14 14 from multiprocessing import Process, Queue
15 15 from schainpy.utils import log
16 16
17 17 import copy
18 18
19 19 QUEUE_SIZE = int(os.environ.get('QUEUE_MAX_SIZE', '10'))
20 20
21 21 class ProcessingUnit(object):
22 22 '''
23 23 Base class to create Signal Chain Units
24 24 '''
25 25
26 26 proc_type = 'processing'
27 27
28 28 def __init__(self):
29 29
30 30 self.dataIn = None
31 31 self.dataOut = None
32 32 self.isConfig = False
33 33 self.operations = []
34 34 self.name = 'Test'
35 35 self.inputs = []
36 36
37 37 def setInput(self, unit):
38 38
39 39 attr = 'dataIn'
40 40 for i, u in enumerate(unit):
41 41 if i==0:
42 42 #print(u.dataOut.flagNoData)
43 43 #exit(1)
44 44 self.dataIn = u.dataOut#.copy()
45 45 self.inputs.append('dataIn')
46 46 else:
47 47 setattr(self, 'dataIn{}'.format(i), u.dataOut)#.copy())
48 48 self.inputs.append('dataIn{}'.format(i))
49 49
50 50
51 51 def getAllowedArgs(self):
52 52 if hasattr(self, '__attrs__'):
53 53 return self.__attrs__
54 54 else:
55 55 return inspect.getargspec(self.run).args
56 56
57 57 def addOperation(self, conf, operation):
58 58 '''
59 59 '''
60 60
61 61 self.operations.append((operation, conf.type, conf.getKwargs()))
62 62
63 63 def getOperationObj(self, objId):
64 64
65 65 if objId not in list(self.operations.keys()):
66 66 return None
67 67
68 68 return self.operations[objId]
69 69
70 70 def call(self, **kwargs):
71 71 '''
72 72 '''
73 #print("call")
73
74 74 try:
75 75 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
76 76 #if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error and not self.dataIn.runNextUnit:
77 77 if self.dataIn.runNextUnit:
78 78 #print("SUCCESSSSSSS")
79 79 #exit(1)
80 80 return not self.dataIn.isReady()
81 81 else:
82 82 return self.dataIn.isReady()
83 83 elif self.dataIn is None or not self.dataIn.error:
84 84 #print([getattr(self, at) for at in self.inputs])
85 85 #print("Elif 1")
86 86 self.run(**kwargs)
87 87 elif self.dataIn.error:
88 88 #print("Elif 2")
89 89 self.dataOut.error = self.dataIn.error
90 90 self.dataOut.flagNoData = True
91 91 except:
92 92 #print("Except")
93 93 err = traceback.format_exc()
94 94 if 'SchainWarning' in err:
95 95 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), self.name)
96 96 elif 'SchainError' in err:
97 97 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), self.name)
98 98 else:
99 99 log.error(err, self.name)
100 100 self.dataOut.error = True
101 101 #print("before op")
102 102 for op, optype, opkwargs in self.operations:
103 103 aux = self.dataOut.copy()
104 104 #aux = copy.deepcopy(self.dataOut)
105 105 #print("**********************Before",op)
106 106 if optype == 'other' and not self.dataOut.flagNoData:
107 107 #print("**********************Other",op)
108 108 #print(self.dataOut.flagNoData)
109 109 self.dataOut = op.run(self.dataOut, **opkwargs)
110 110 elif optype == 'external' and not self.dataOut.flagNoData:
111 111 op.queue.put(aux)
112 112 elif optype == 'external' and self.dataOut.error:
113 113 op.queue.put(aux)
114 114 #elif optype == 'external' and self.dataOut.isReady():
115 115 #op.queue.put(copy.deepcopy(self.dataOut))
116 116 #print(not self.dataOut.isReady())
117 117
118 118 try:
119 119 if self.dataOut.runNextUnit:
120 120 runNextUnit = self.dataOut.runNextUnit
121 121 #print(self.operations)
122 122 #print("Tru")
123 123
124 124 else:
125 125 runNextUnit = self.dataOut.isReady()
126 126 except:
127 127 runNextUnit = self.dataOut.isReady()
128 #exit(1)
128 129 #if not self.dataOut.isReady():
129 130 #return 'Error' if self.dataOut.error else input()
130 131 #print("NexT",runNextUnit)
131 132 #print("error: ",self.dataOut.error)
132 133 return 'Error' if self.dataOut.error else runNextUnit# self.dataOut.isReady()
133 134
134 135 def setup(self):
135 136
136 137 raise NotImplementedError
137 138
138 139 def run(self):
139 140
140 141 raise NotImplementedError
141 142
142 143 def close(self):
143 144
144 145 return
145 146
146 147
147 148 class Operation(object):
148 149
149 150 '''
150 151 '''
151 152
152 153 proc_type = 'operation'
153 154
154 155 def __init__(self):
155 156
156 157 self.id = None
157 158 self.isConfig = False
158 159
159 160 if not hasattr(self, 'name'):
160 161 self.name = self.__class__.__name__
161 162
162 163 def getAllowedArgs(self):
163 164 if hasattr(self, '__attrs__'):
164 165 return self.__attrs__
165 166 else:
166 167 return inspect.getargspec(self.run).args
167 168
168 169 def setup(self):
169 170
170 171 self.isConfig = True
171 172
172 173 raise NotImplementedError
173 174
174 175 def run(self, dataIn, **kwargs):
175 176 """
176 177 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
177 178 atributos del objeto dataIn.
178 179
179 180 Input:
180 181
181 182 dataIn : objeto del tipo JROData
182 183
183 184 Return:
184 185
185 186 None
186 187
187 188 Affected:
188 189 __buffer : buffer de recepcion de datos.
189 190
190 191 """
191 192 if not self.isConfig:
192 193 self.setup(**kwargs)
193 194
194 195 raise NotImplementedError
195 196
196 197 def close(self):
197 198
198 199 return
199 200
200 201
201 202 def MPDecorator(BaseClass):
202 203 """
203 204 Multiprocessing class decorator
204 205
205 206 This function add multiprocessing features to a BaseClass.
206 207 """
207 208
208 209 class MPClass(BaseClass, Process):
209 210
210 211 def __init__(self, *args, **kwargs):
211 212 super(MPClass, self).__init__()
212 213 Process.__init__(self)
213 214
214 215 self.args = args
215 216 self.kwargs = kwargs
216 217 self.t = time.time()
217 218 self.op_type = 'external'
218 219 self.name = BaseClass.__name__
219 220 self.__doc__ = BaseClass.__doc__
220 221
221 222 if 'plot' in self.name.lower() and not self.name.endswith('_'):
222 223 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
223 224
224 225 self.start_time = time.time()
225 226 self.err_queue = args[3]
226 227 self.queue = Queue(maxsize=QUEUE_SIZE)
227 228 self.myrun = BaseClass.run
228 229
229 230 def run(self):
230 231
231 232 while True:
232 233
233 234 dataOut = self.queue.get()
234 235
235 236 if not dataOut.error:
236 237 try:
237 238 BaseClass.run(self, dataOut, **self.kwargs)
238 239 except:
239 240 err = traceback.format_exc()
240 241 log.error(err, self.name)
241 242 else:
242 243 break
243 244
244 245 self.close()
245 246
246 247 def close(self):
247 248
248 249 BaseClass.close(self)
249 250 log.success('Done...(Time:{:4.2f} secs)'.format(time.time() - self.start_time), self.name)
250 251
251 252 return MPClass
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,910 +1,957
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15
16 16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 17 from schainpy.model.data.jrodata import Spectra
18 18 from schainpy.model.data.jrodata import hildebrand_sekhon
19 19 from schainpy.utils import log
20 20
21 21
22 22 class SpectraProc(ProcessingUnit):
23 23
24 24 def __init__(self):
25 25
26 26 ProcessingUnit.__init__(self)
27 27
28 28 self.buffer = None
29 29 self.firstdatatime = None
30 30 self.profIndex = 0
31 31 self.dataOut = Spectra()
32 32 self.id_min = None
33 33 self.id_max = None
34 34 self.setupReq = False #Agregar a todas las unidades de proc
35 35
36 36 def __updateSpecFromVoltage(self):
37 37
38 38 self.dataOut.timeZone = self.dataIn.timeZone
39 39 self.dataOut.dstFlag = self.dataIn.dstFlag
40 40 self.dataOut.errorCount = self.dataIn.errorCount
41 41 self.dataOut.useLocalTime = self.dataIn.useLocalTime
42 42 try:
43 43 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
44 44 except:
45 45 pass
46 46 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
47 47 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
48 48 self.dataOut.channelList = self.dataIn.channelList
49 49 self.dataOut.heightList = self.dataIn.heightList
50 50 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
51 51 self.dataOut.nProfiles = self.dataOut.nFFTPoints
52 52 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
53 53 self.dataOut.utctime = self.firstdatatime
54 54 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
55 55 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
56 56 self.dataOut.flagShiftFFT = False
57 57 self.dataOut.nCohInt = self.dataIn.nCohInt
58 58 self.dataOut.nIncohInt = 1
59 59 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
60 60 self.dataOut.frequency = self.dataIn.frequency
61 61 self.dataOut.realtime = self.dataIn.realtime
62 62 self.dataOut.azimuth = self.dataIn.azimuth
63 63 self.dataOut.zenith = self.dataIn.zenith
64 64 self.dataOut.beam.codeList = self.dataIn.beam.codeList
65 65 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
66 66 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
67 67 self.dataOut.runNextUnit = self.dataIn.runNextUnit
68 68 try:
69 69 self.dataOut.step = self.dataIn.step
70 70 except:
71 71 pass
72 72
73 73 def __getFft(self):
74 74 """
75 75 Convierte valores de Voltaje a Spectra
76 76
77 77 Affected:
78 78 self.dataOut.data_spc
79 79 self.dataOut.data_cspc
80 80 self.dataOut.data_dc
81 81 self.dataOut.heightList
82 82 self.profIndex
83 83 self.buffer
84 84 self.dataOut.flagNoData
85 85 """
86 86 fft_volt = numpy.fft.fft(
87 87 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
88 88 fft_volt = fft_volt.astype(numpy.dtype('complex'))
89 89 dc = fft_volt[:, 0, :]
90 90
91 91 # calculo de self-spectra
92 92 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
93 93 spc = fft_volt * numpy.conjugate(fft_volt)
94 94 spc = spc.real
95 95
96 96 blocksize = 0
97 97 blocksize += dc.size
98 98 blocksize += spc.size
99 99
100 100 cspc = None
101 101 pairIndex = 0
102 102 if self.dataOut.pairsList != None:
103 103 # calculo de cross-spectra
104 104 cspc = numpy.zeros(
105 105 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
106 106 for pair in self.dataOut.pairsList:
107 107 if pair[0] not in self.dataOut.channelList:
108 108 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
109 109 str(pair), str(self.dataOut.channelList)))
110 110 if pair[1] not in self.dataOut.channelList:
111 111 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
112 112 str(pair), str(self.dataOut.channelList)))
113 113
114 114 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
115 115 numpy.conjugate(fft_volt[pair[1], :, :])
116 116 pairIndex += 1
117 117 blocksize += cspc.size
118 118
119 119 self.dataOut.data_spc = spc
120 120 self.dataOut.data_cspc = cspc
121 121 self.dataOut.data_dc = dc
122 122 self.dataOut.blockSize = blocksize
123 123 self.dataOut.flagShiftFFT = False
124 124
125 125 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
126
126
127 127 self.dataIn.runNextUnit = runNextUnit
128 128 if self.dataIn.type == "Spectra":
129
129 130 self.dataOut.copy(self.dataIn)
130 131 if shift_fft:
131 132 #desplaza a la derecha en el eje 2 determinadas posiciones
132 133 shift = int(self.dataOut.nFFTPoints/2)
133 134 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
134 135
135 136 if self.dataOut.data_cspc is not None:
136 137 #desplaza a la derecha en el eje 2 determinadas posiciones
137 138 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
138 139 if pairsList:
139 140 self.__selectPairs(pairsList)
140 141
141 142 elif self.dataIn.type == "Voltage":
142 143
143 144 self.dataOut.flagNoData = True
144 145
145 146 if nFFTPoints == None:
146 147 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
147 148
148 149 if nProfiles == None:
149 150 nProfiles = nFFTPoints
150
151 #print(self.dataOut.ipp)
152 #exit(1)
151 153 if ippFactor == None:
152 154 self.dataOut.ippFactor = 1
155 #if ippFactor is not None:
156 #self.dataOut.ippFactor = ippFactor
157 #print(ippFactor)
158 #print(self.dataOut.ippFactor)
159 #exit(1)
153 160
154 161 self.dataOut.nFFTPoints = nFFTPoints
155 162
156 163 if self.buffer is None:
157 164 self.buffer = numpy.zeros((self.dataIn.nChannels,
158 165 nProfiles,
159 166 self.dataIn.nHeights),
160 167 dtype='complex')
161 168
162 169 if self.dataIn.flagDataAsBlock:
163 170 nVoltProfiles = self.dataIn.data.shape[1]
164 171
165 172 if nVoltProfiles == nProfiles:
166 173 self.buffer = self.dataIn.data.copy()
167 174 self.profIndex = nVoltProfiles
168 175
169 176 elif nVoltProfiles < nProfiles:
170 177
171 178 if self.profIndex == 0:
172 179 self.id_min = 0
173 180 self.id_max = nVoltProfiles
174 181 #print(self.id_min)
175 182 #print(self.id_max)
176 183 #print(numpy.shape(self.buffer))
177 184 self.buffer[:, self.id_min:self.id_max,
178 185 :] = self.dataIn.data
179 186 self.profIndex += nVoltProfiles
180 187 self.id_min += nVoltProfiles
181 188 self.id_max += nVoltProfiles
182 189 else:
183 190 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
184 191 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
185 192 self.dataOut.flagNoData = True
186 193 else:
187 194 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
188 195 self.profIndex += 1
189 196
190 197 if self.firstdatatime == None:
191 198 self.firstdatatime = self.dataIn.utctime
192 199
193 200 if self.profIndex == nProfiles:
194 201 self.__updateSpecFromVoltage()
195 202 if pairsList == None:
196 203 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
197 204 else:
198 205 self.dataOut.pairsList = pairsList
199 206 self.__getFft()
200 207 self.dataOut.flagNoData = False
201 208 self.firstdatatime = None
202 209 self.profIndex = 0
203 210 else:
204 211 raise ValueError("The type of input object '%s' is not valid".format(
205 212 self.dataIn.type))
206 213
207 214
208 215 def __selectPairs(self, pairsList):
209 216
210 217 if not pairsList:
211 218 return
212 219
213 220 pairs = []
214 221 pairsIndex = []
215 222
216 223 for pair in pairsList:
217 224 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
218 225 continue
219 226 pairs.append(pair)
220 227 pairsIndex.append(pairs.index(pair))
221 228
222 229 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
223 230 self.dataOut.pairsList = pairs
224 231
225 232 return
226 233
227 234 def selectFFTs(self, minFFT, maxFFT ):
228 235 """
229 236 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
230 237 minFFT<= FFT <= maxFFT
231 238 """
232 239
233 240 if (minFFT > maxFFT):
234 241 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
235 242
236 243 if (minFFT < self.dataOut.getFreqRange()[0]):
237 244 minFFT = self.dataOut.getFreqRange()[0]
238 245
239 246 if (maxFFT > self.dataOut.getFreqRange()[-1]):
240 247 maxFFT = self.dataOut.getFreqRange()[-1]
241 248
242 249 minIndex = 0
243 250 maxIndex = 0
244 251 FFTs = self.dataOut.getFreqRange()
245 252
246 253 inda = numpy.where(FFTs >= minFFT)
247 254 indb = numpy.where(FFTs <= maxFFT)
248 255
249 256 try:
250 257 minIndex = inda[0][0]
251 258 except:
252 259 minIndex = 0
253 260
254 261 try:
255 262 maxIndex = indb[0][-1]
256 263 except:
257 264 maxIndex = len(FFTs)
258 265
259 266 self.selectFFTsByIndex(minIndex, maxIndex)
260 267
261 268 return 1
262 269
263 270 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
264 271 newheis = numpy.where(
265 272 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
266 273
267 274 if hei_ref != None:
268 275 newheis = numpy.where(self.dataOut.heightList > hei_ref)
269 276
270 277 minIndex = min(newheis[0])
271 278 maxIndex = max(newheis[0])
272 279 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
273 280 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
274 281
275 282 # determina indices
276 283 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
277 284 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
278 285 avg_dB = 10 * \
279 286 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
280 287 beacon_dB = numpy.sort(avg_dB)[-nheis:]
281 288 beacon_heiIndexList = []
282 289 for val in avg_dB.tolist():
283 290 if val >= beacon_dB[0]:
284 291 beacon_heiIndexList.append(avg_dB.tolist().index(val))
285 292
286 293 #data_spc = data_spc[:,:,beacon_heiIndexList]
287 294 data_cspc = None
288 295 if self.dataOut.data_cspc is not None:
289 296 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
290 297 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
291 298
292 299 data_dc = None
293 300 if self.dataOut.data_dc is not None:
294 301 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
295 302 #data_dc = data_dc[:,beacon_heiIndexList]
296 303
297 304 self.dataOut.data_spc = data_spc
298 305 self.dataOut.data_cspc = data_cspc
299 306 self.dataOut.data_dc = data_dc
300 307 self.dataOut.heightList = heightList
301 308 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
302 309
303 310 return 1
304 311
305 312 def selectFFTsByIndex(self, minIndex, maxIndex):
306 313 """
307 314
308 315 """
309 316
310 317 if (minIndex < 0) or (minIndex > maxIndex):
311 318 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
312 319
313 320 if (maxIndex >= self.dataOut.nProfiles):
314 321 maxIndex = self.dataOut.nProfiles-1
315 322
316 323 #Spectra
317 324 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
318 325
319 326 data_cspc = None
320 327 if self.dataOut.data_cspc is not None:
321 328 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
322 329
323 330 data_dc = None
324 331 if self.dataOut.data_dc is not None:
325 332 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
326 333
327 334 self.dataOut.data_spc = data_spc
328 335 self.dataOut.data_cspc = data_cspc
329 336 self.dataOut.data_dc = data_dc
330 337
331 338 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
332 339 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
333 340 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
334 341
335 342 return 1
336 343
337 344 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
338 345 # validacion de rango
339 346 print("NOISeeee")
340 347 if minHei == None:
341 348 minHei = self.dataOut.heightList[0]
342 349
343 350 if maxHei == None:
344 351 maxHei = self.dataOut.heightList[-1]
345 352
346 353 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
347 354 print('minHei: %.2f is out of the heights range' % (minHei))
348 355 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
349 356 minHei = self.dataOut.heightList[0]
350 357
351 358 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
352 359 print('maxHei: %.2f is out of the heights range' % (maxHei))
353 360 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
354 361 maxHei = self.dataOut.heightList[-1]
355 362
356 363 # validacion de velocidades
357 364 velrange = self.dataOut.getVelRange(1)
358 365
359 366 if minVel == None:
360 367 minVel = velrange[0]
361 368
362 369 if maxVel == None:
363 370 maxVel = velrange[-1]
364 371
365 372 if (minVel < velrange[0]) or (minVel > maxVel):
366 373 print('minVel: %.2f is out of the velocity range' % (minVel))
367 374 print('minVel is setting to %.2f' % (velrange[0]))
368 375 minVel = velrange[0]
369 376
370 377 if (maxVel > velrange[-1]) or (maxVel < minVel):
371 378 print('maxVel: %.2f is out of the velocity range' % (maxVel))
372 379 print('maxVel is setting to %.2f' % (velrange[-1]))
373 380 maxVel = velrange[-1]
374 381
375 382 # seleccion de indices para rango
376 383 minIndex = 0
377 384 maxIndex = 0
378 385 heights = self.dataOut.heightList
379 386
380 387 inda = numpy.where(heights >= minHei)
381 388 indb = numpy.where(heights <= maxHei)
382 389
383 390 try:
384 391 minIndex = inda[0][0]
385 392 except:
386 393 minIndex = 0
387 394
388 395 try:
389 396 maxIndex = indb[0][-1]
390 397 except:
391 398 maxIndex = len(heights)
392 399
393 400 if (minIndex < 0) or (minIndex > maxIndex):
394 401 raise ValueError("some value in (%d,%d) is not valid" % (
395 402 minIndex, maxIndex))
396 403
397 404 if (maxIndex >= self.dataOut.nHeights):
398 405 maxIndex = self.dataOut.nHeights - 1
399 406
400 407 # seleccion de indices para velocidades
401 408 indminvel = numpy.where(velrange >= minVel)
402 409 indmaxvel = numpy.where(velrange <= maxVel)
403 410 try:
404 411 minIndexVel = indminvel[0][0]
405 412 except:
406 413 minIndexVel = 0
407 414
408 415 try:
409 416 maxIndexVel = indmaxvel[0][-1]
410 417 except:
411 418 maxIndexVel = len(velrange)
412 419
413 420 # seleccion del espectro
414 421 data_spc = self.dataOut.data_spc[:,
415 422 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
416 423 # estimacion de ruido
417 424 noise = numpy.zeros(self.dataOut.nChannels)
418 425
419 426 for channel in range(self.dataOut.nChannels):
420 427 daux = data_spc[channel, :, :]
421 428 sortdata = numpy.sort(daux, axis=None)
422 429 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
423 430
424 431 self.dataOut.noise_estimation = noise.copy()
425 432
426 433 return 1
427 434
435 class GetSNR(Operation):
436 '''
437 Written by R. Flores
438 '''
439 """Operation to get SNR.
440
441 Parameters:
442 -----------
443
444 Example
445 --------
446
447 op = proc_unit.addOperation(name='GetSNR', optype='other')
448
449 """
450
451 def __init__(self, **kwargs):
452
453 Operation.__init__(self, **kwargs)
454
455
456 def run(self,dataOut):
457
458 noise = dataOut.getNoise()
459 maxdB = 16
460
461 normFactor = 24
462
463 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.normFactor)
464 dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints)
465
466 dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr)
467 #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr)
468 #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0)
469 #print(dataOut.data_snr.shape)
470 #exit(1)
471
472
473 return dataOut
474
428 475 class removeDC(Operation):
429 476
430 477 def run(self, dataOut, mode=2):
431 478 self.dataOut = dataOut
432 479 jspectra = self.dataOut.data_spc
433 480 jcspectra = self.dataOut.data_cspc
434 481
435 482 num_chan = jspectra.shape[0]
436 483 num_hei = jspectra.shape[2]
437 484
438 485 if jcspectra is not None:
439 486 jcspectraExist = True
440 487 num_pairs = jcspectra.shape[0]
441 488 else:
442 489 jcspectraExist = False
443 490
444 491 freq_dc = int(jspectra.shape[1] / 2)
445 492 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
446 493 ind_vel = ind_vel.astype(int)
447 494
448 495 if ind_vel[0] < 0:
449 496 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
450 497
451 498 if mode == 1:
452 499 jspectra[:, freq_dc, :] = (
453 500 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
454 501
455 502 if jcspectraExist:
456 503 jcspectra[:, freq_dc, :] = (
457 504 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
458 505
459 506 if mode == 2:
460 507
461 508 vel = numpy.array([-2, -1, 1, 2])
462 509 xx = numpy.zeros([4, 4])
463 510
464 511 for fil in range(4):
465 512 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
466 513
467 514 xx_inv = numpy.linalg.inv(xx)
468 515 xx_aux = xx_inv[0, :]
469 516
470 517 for ich in range(num_chan):
471 518 yy = jspectra[ich, ind_vel, :]
472 519 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
473 520
474 521 junkid = jspectra[ich, freq_dc, :] <= 0
475 522 cjunkid = sum(junkid)
476 523
477 524 if cjunkid.any():
478 525 jspectra[ich, freq_dc, junkid.nonzero()] = (
479 526 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
480 527
481 528 if jcspectraExist:
482 529 for ip in range(num_pairs):
483 530 yy = jcspectra[ip, ind_vel, :]
484 531 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
485 532
486 533 self.dataOut.data_spc = jspectra
487 534 self.dataOut.data_cspc = jcspectra
488 535
489 536 return self.dataOut
490 537
491 538 class removeInterference(Operation):
492 539
493 540 def removeInterference2(self):
494 541
495 542 cspc = self.dataOut.data_cspc
496 543 spc = self.dataOut.data_spc
497 544 Heights = numpy.arange(cspc.shape[2])
498 545 realCspc = numpy.abs(cspc)
499 546
500 547 for i in range(cspc.shape[0]):
501 548 LinePower= numpy.sum(realCspc[i], axis=0)
502 549 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
503 550 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
504 551 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
505 552 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
506 553 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
507 554
508 555
509 556 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
510 557 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
511 558 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
512 559 cspc[i,InterferenceRange,:] = numpy.NaN
513 560
514 561 self.dataOut.data_cspc = cspc
515 562
516 563 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
517 564
518 565 jspectra = self.dataOut.data_spc
519 566 jcspectra = self.dataOut.data_cspc
520 567 jnoise = self.dataOut.getNoise()
521 568 num_incoh = self.dataOut.nIncohInt
522 569
523 570 num_channel = jspectra.shape[0]
524 571 num_prof = jspectra.shape[1]
525 572 num_hei = jspectra.shape[2]
526 573
527 574 # hei_interf
528 575 if hei_interf is None:
529 576 count_hei = int(num_hei / 2)
530 577 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
531 578 hei_interf = numpy.asarray(hei_interf)[0]
532 579 # nhei_interf
533 580 if (nhei_interf == None):
534 581 nhei_interf = 5
535 582 if (nhei_interf < 1):
536 583 nhei_interf = 1
537 584 if (nhei_interf > count_hei):
538 585 nhei_interf = count_hei
539 586 if (offhei_interf == None):
540 587 offhei_interf = 0
541 588
542 589 ind_hei = list(range(num_hei))
543 590 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
544 591 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
545 592 mask_prof = numpy.asarray(list(range(num_prof)))
546 593 num_mask_prof = mask_prof.size
547 594 comp_mask_prof = [0, num_prof / 2]
548 595
549 596 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
550 597 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
551 598 jnoise = numpy.nan
552 599 noise_exist = jnoise[0] < numpy.Inf
553 600
554 601 # Subrutina de Remocion de la Interferencia
555 602 for ich in range(num_channel):
556 603 # Se ordena los espectros segun su potencia (menor a mayor)
557 604 power = jspectra[ich, mask_prof, :]
558 605 power = power[:, hei_interf]
559 606 power = power.sum(axis=0)
560 607 psort = power.ravel().argsort()
561 608
562 609 # Se estima la interferencia promedio en los Espectros de Potencia empleando
563 610 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
564 611 offhei_interf, nhei_interf + offhei_interf))]]]
565 612
566 613 if noise_exist:
567 614 # tmp_noise = jnoise[ich] / num_prof
568 615 tmp_noise = jnoise[ich]
569 616 junkspc_interf = junkspc_interf - tmp_noise
570 617 #junkspc_interf[:,comp_mask_prof] = 0
571 618
572 619 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
573 620 jspc_interf = jspc_interf.transpose()
574 621 # Calculando el espectro de interferencia promedio
575 622 noiseid = numpy.where(
576 623 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
577 624 noiseid = noiseid[0]
578 625 cnoiseid = noiseid.size
579 626 interfid = numpy.where(
580 627 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
581 628 interfid = interfid[0]
582 629 cinterfid = interfid.size
583 630
584 631 if (cnoiseid > 0):
585 632 jspc_interf[noiseid] = 0
586 633
587 634 # Expandiendo los perfiles a limpiar
588 635 if (cinterfid > 0):
589 636 new_interfid = (
590 637 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
591 638 new_interfid = numpy.asarray(new_interfid)
592 639 new_interfid = {x for x in new_interfid}
593 640 new_interfid = numpy.array(list(new_interfid))
594 641 new_cinterfid = new_interfid.size
595 642 else:
596 643 new_cinterfid = 0
597 644
598 645 for ip in range(new_cinterfid):
599 646 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
600 647 jspc_interf[new_interfid[ip]
601 648 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
602 649
603 650 jspectra[ich, :, ind_hei] = jspectra[ich, :,
604 651 ind_hei] - jspc_interf # Corregir indices
605 652
606 653 # Removiendo la interferencia del punto de mayor interferencia
607 654 ListAux = jspc_interf[mask_prof].tolist()
608 655 maxid = ListAux.index(max(ListAux))
609 656
610 657 if cinterfid > 0:
611 658 for ip in range(cinterfid * (interf == 2) - 1):
612 659 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
613 660 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
614 661 cind = len(ind)
615 662
616 663 if (cind > 0):
617 664 jspectra[ich, interfid[ip], ind] = tmp_noise * \
618 665 (1 + (numpy.random.uniform(cind) - 0.5) /
619 666 numpy.sqrt(num_incoh))
620 667
621 668 ind = numpy.array([-2, -1, 1, 2])
622 669 xx = numpy.zeros([4, 4])
623 670
624 671 for id1 in range(4):
625 672 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
626 673
627 674 xx_inv = numpy.linalg.inv(xx)
628 675 xx = xx_inv[:, 0]
629 676 ind = (ind + maxid + num_mask_prof) % num_mask_prof
630 677 yy = jspectra[ich, mask_prof[ind], :]
631 678 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
632 679 yy.transpose(), xx)
633 680
634 681 indAux = (jspectra[ich, :, :] < tmp_noise *
635 682 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
636 683 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
637 684 (1 - 1 / numpy.sqrt(num_incoh))
638 685
639 686 # Remocion de Interferencia en el Cross Spectra
640 687 if jcspectra is None:
641 688 return jspectra, jcspectra
642 689 num_pairs = int(jcspectra.size / (num_prof * num_hei))
643 690 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
644 691
645 692 for ip in range(num_pairs):
646 693
647 694 #-------------------------------------------
648 695
649 696 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
650 697 cspower = cspower[:, hei_interf]
651 698 cspower = cspower.sum(axis=0)
652 699
653 700 cspsort = cspower.ravel().argsort()
654 701 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
655 702 offhei_interf, nhei_interf + offhei_interf))]]]
656 703 junkcspc_interf = junkcspc_interf.transpose()
657 704 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
658 705
659 706 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
660 707
661 708 median_real = int(numpy.median(numpy.real(
662 709 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
663 710 median_imag = int(numpy.median(numpy.imag(
664 711 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
665 712 comp_mask_prof = [int(e) for e in comp_mask_prof]
666 713 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
667 714 median_real, median_imag)
668 715
669 716 for iprof in range(num_prof):
670 717 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
671 718 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
672 719
673 720 # Removiendo la Interferencia
674 721 jcspectra[ip, :, ind_hei] = jcspectra[ip,
675 722 :, ind_hei] - jcspc_interf
676 723
677 724 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
678 725 maxid = ListAux.index(max(ListAux))
679 726
680 727 ind = numpy.array([-2, -1, 1, 2])
681 728 xx = numpy.zeros([4, 4])
682 729
683 730 for id1 in range(4):
684 731 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
685 732
686 733 xx_inv = numpy.linalg.inv(xx)
687 734 xx = xx_inv[:, 0]
688 735
689 736 ind = (ind + maxid + num_mask_prof) % num_mask_prof
690 737 yy = jcspectra[ip, mask_prof[ind], :]
691 738 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
692 739
693 740 # Guardar Resultados
694 741 self.dataOut.data_spc = jspectra
695 742 self.dataOut.data_cspc = jcspectra
696 743
697 744 return 1
698 745
699 746 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
700 747
701 748 self.dataOut = dataOut
702 749
703 750 if mode == 1:
704 751 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
705 752 elif mode == 2:
706 753 self.removeInterference2()
707 754
708 755 return self.dataOut
709 756
710 757
711 758 class IncohInt(Operation):
712 759
713 760 __profIndex = 0
714 761 __withOverapping = False
715 762
716 763 __byTime = False
717 764 __initime = None
718 765 __lastdatatime = None
719 766 __integrationtime = None
720 767
721 768 __buffer_spc = None
722 769 __buffer_cspc = None
723 770 __buffer_dc = None
724 771
725 772 __dataReady = False
726 773
727 774 __timeInterval = None
728 775
729 776 n = None
730 777
731 778 def __init__(self):
732 779
733 780 Operation.__init__(self)
734 781
735 782 def setup(self, n=None, timeInterval=None, overlapping=False):
736 783 """
737 784 Set the parameters of the integration class.
738 785
739 786 Inputs:
740 787
741 788 n : Number of coherent integrations
742 789 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
743 790 overlapping :
744 791
745 792 """
746 793
747 794 self.__initime = None
748 795 self.__lastdatatime = 0
749 796
750 797 self.__buffer_spc = 0
751 798 self.__buffer_cspc = 0
752 799 self.__buffer_dc = 0
753 800
754 801 self.__profIndex = 0
755 802 self.__dataReady = False
756 803 self.__byTime = False
757 804
758 805 if n is None and timeInterval is None:
759 806 raise ValueError("n or timeInterval should be specified ...")
760 807
761 808 if n is not None:
762 809 self.n = int(n)
763 810 else:
764 811
765 812 self.__integrationtime = int(timeInterval)
766 813 self.n = None
767 814 self.__byTime = True
768 815
769 816 def putData(self, data_spc, data_cspc, data_dc):
770 817 """
771 818 Add a profile to the __buffer_spc and increase in one the __profileIndex
772 819
773 820 """
774 821
775 822 self.__buffer_spc += data_spc
776 823
777 824 if data_cspc is None:
778 825 self.__buffer_cspc = None
779 826 else:
780 827 self.__buffer_cspc += data_cspc
781 828
782 829 if data_dc is None:
783 830 self.__buffer_dc = None
784 831 else:
785 832 self.__buffer_dc += data_dc
786 833
787 834 self.__profIndex += 1
788 835
789 836 return
790 837
791 838 def pushData(self):
792 839 """
793 840 Return the sum of the last profiles and the profiles used in the sum.
794 841
795 842 Affected:
796 843
797 844 self.__profileIndex
798 845
799 846 """
800 847
801 848 data_spc = self.__buffer_spc
802 849 data_cspc = self.__buffer_cspc
803 850 data_dc = self.__buffer_dc
804 851 n = self.__profIndex
805 852
806 853 self.__buffer_spc = 0
807 854 self.__buffer_cspc = 0
808 855 self.__buffer_dc = 0
809 856 self.__profIndex = 0
810 857
811 858 return data_spc, data_cspc, data_dc, n
812 859
813 860 def byProfiles(self, *args):
814 861
815 862 self.__dataReady = False
816 863 avgdata_spc = None
817 864 avgdata_cspc = None
818 865 avgdata_dc = None
819 866
820 867 self.putData(*args)
821 868
822 869 if self.__profIndex == self.n:
823 870
824 871 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
825 872 self.n = n
826 873 self.__dataReady = True
827 874
828 875 return avgdata_spc, avgdata_cspc, avgdata_dc
829 876
830 877 def byTime(self, datatime, *args):
831 878
832 879 self.__dataReady = False
833 880 avgdata_spc = None
834 881 avgdata_cspc = None
835 882 avgdata_dc = None
836 883
837 884 self.putData(*args)
838 885
839 886 if (datatime - self.__initime) >= self.__integrationtime:
840 887 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
841 888 self.n = n
842 889 self.__dataReady = True
843 890
844 891 return avgdata_spc, avgdata_cspc, avgdata_dc
845 892
846 893 def integrate(self, datatime, *args):
847 894
848 895 if self.__profIndex == 0:
849 896 self.__initime = datatime
850 897
851 898 if self.__byTime:
852 899 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
853 900 datatime, *args)
854 901 else:
855 902 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
856 903
857 904 if not self.__dataReady:
858 905 return None, None, None, None
859 906
860 907 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
861 908
862 909 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
863 910 if n == 1:
864 911 return dataOut
865 912 print("JERE")
866 913 dataOut.flagNoData = True
867 914
868 915 if not self.isConfig:
869 916 self.setup(n, timeInterval, overlapping)
870 917 self.isConfig = True
871 918
872 919 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
873 920 dataOut.data_spc,
874 921 dataOut.data_cspc,
875 922 dataOut.data_dc)
876 923
877 924 if self.__dataReady:
878 925
879 926 dataOut.data_spc = avgdata_spc
880 927 print(numpy.sum(dataOut.data_spc))
881 928 exit(1)
882 929 dataOut.data_cspc = avgdata_cspc
883 930 dataOut.data_dc = avgdata_dc
884 931 dataOut.nIncohInt *= self.n
885 932 dataOut.utctime = avgdatatime
886 933 dataOut.flagNoData = False
887 934
888 935 return dataOut
889 936
890 937 class dopplerFlip(Operation):
891 938
892 939 def run(self, dataOut):
893 940 # arreglo 1: (num_chan, num_profiles, num_heights)
894 941 self.dataOut = dataOut
895 942 # JULIA-oblicua, indice 2
896 943 # arreglo 2: (num_profiles, num_heights)
897 944 jspectra = self.dataOut.data_spc[2]
898 945 jspectra_tmp = numpy.zeros(jspectra.shape)
899 946 num_profiles = jspectra.shape[0]
900 947 freq_dc = int(num_profiles / 2)
901 948 # Flip con for
902 949 for j in range(num_profiles):
903 950 jspectra_tmp[num_profiles-j-1]= jspectra[j]
904 951 # Intercambio perfil de DC con perfil inmediato anterior
905 952 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
906 953 jspectra_tmp[freq_dc]= jspectra[freq_dc]
907 954 # canal modificado es re-escrito en el arreglo de canales
908 955 self.dataOut.data_spc[2] = jspectra_tmp
909 956
910 957 return self.dataOut
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now