##// END OF EJS Templates
schain
RSS
last update
rflores -
r1600:2c146cb976d3
parent child
Show More
Side by Side Unified
Context file:
r1600:2c146cb976d3 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,1082
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 #pwcode = 64
552 #print("pwcode: ", pwcode)
553 #exit(1)
551 554 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
552 555 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
553 556
554 557 return normFactor
555 558
556 559 @property
557 560 def flag_cspc(self):
558 561
559 562 if self.data_cspc is None:
560 563 return True
561 564
562 565 return False
563 566
564 567 @property
565 568 def flag_dc(self):
566 569
567 570 if self.data_dc is None:
568 571 return True
569 572
570 573 return False
571 574
572 575 @property
573 576 def timeInterval(self):
574 577
575 578 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
576 579 if self.nmodes:
577 580 return self.nmodes*timeInterval
578 581 else:
579 582 return timeInterval
580 583
581 584 def getPower(self):
582 585
583 586 factor = self.normFactor
584 587 z = self.data_spc / factor
585 588 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
586 589 avg = numpy.average(z, axis=1)
587 590
588 591 return 10 * numpy.log10(avg)
589 592
590 593 def getCoherence(self, pairsList=None, phase=False):
591 594
592 595 z = []
593 596 if pairsList is None:
594 597 pairsIndexList = self.pairsIndexList
595 598 else:
596 599 pairsIndexList = []
597 600 for pair in pairsList:
598 601 if pair not in self.pairsList:
599 602 raise ValueError("Pair %s is not in dataOut.pairsList" % (
600 603 pair))
601 604 pairsIndexList.append(self.pairsList.index(pair))
602 605 for i in range(len(pairsIndexList)):
603 606 pair = self.pairsList[pairsIndexList[i]]
604 607 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
605 608 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
606 609 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
607 610 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
608 611 if phase:
609 612 data = numpy.arctan2(avgcoherenceComplex.imag,
610 613 avgcoherenceComplex.real) * 180 / numpy.pi
611 614 else:
612 615 data = numpy.abs(avgcoherenceComplex)
613 616
614 617 z.append(data)
615 618
616 619 return numpy.array(z)
617 620
618 621 def setValue(self, value):
619 622
620 623 print("This property should not be initialized")
621 624
622 625 return
623 626
624 627 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
625 628
626 629
627 630 class SpectraHeis(Spectra):
628 631
629 632 def __init__(self):
630 633
631 634 self.radarControllerHeaderObj = RadarControllerHeader()
632 635 self.systemHeaderObj = SystemHeader()
633 636 self.type = "SpectraHeis"
634 637 self.nProfiles = None
635 638 self.heightList = None
636 639 self.channelList = None
637 640 self.flagNoData = True
638 641 self.flagDiscontinuousBlock = False
639 642 self.utctime = None
640 643 self.blocksize = None
641 644 self.profileIndex = 0
642 645 self.nCohInt = 1
643 646 self.nIncohInt = 1
644 647
645 648 @property
646 649 def normFactor(self):
647 650 pwcode = 1
648 651 if self.flagDecodeData:
649 652 pwcode = numpy.sum(self.code[0]**2)
650 653
651 654 normFactor = self.nIncohInt * self.nCohInt * pwcode
652 655
653 656 return normFactor
654 657
655 658 @property
656 659 def timeInterval(self):
657 660
658 661 return self.ippSeconds * self.nCohInt * self.nIncohInt
659 662
660 663
661 664 class Fits(JROData):
662 665
663 666 def __init__(self):
664 667
665 668 self.type = "Fits"
666 669 self.nProfiles = None
667 670 self.heightList = None
668 671 self.channelList = None
669 672 self.flagNoData = True
670 673 self.utctime = None
671 674 self.nCohInt = 1
672 675 self.nIncohInt = 1
673 676 self.useLocalTime = True
674 677 self.profileIndex = 0
675 678 self.timeZone = 0
676 679
677 680 def getTimeRange(self):
678 681
679 682 datatime = []
680 683
681 684 datatime.append(self.ltctime)
682 685 datatime.append(self.ltctime + self.timeInterval)
683 686
684 687 datatime = numpy.array(datatime)
685 688
686 689 return datatime
687 690
688 691 def getChannelIndexList(self):
689 692
690 693 return list(range(self.nChannels))
691 694
692 695 def getNoise(self, type=1):
693 696
694 697
695 698 if type == 1:
696 699 noise = self.getNoisebyHildebrand()
697 700
698 701 if type == 2:
699 702 noise = self.getNoisebySort()
700 703
701 704 if type == 3:
702 705 noise = self.getNoisebyWindow()
703 706
704 707 return noise
705 708
706 709 @property
707 710 def timeInterval(self):
708 711
709 712 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
710 713
711 714 return timeInterval
712 715
713 716 @property
714 717 def ippSeconds(self):
715 718 '''
716 719 '''
717 720 return self.ipp_sec
718 721
719 722 noise = property(getNoise, "I'm the 'nHeights' property.")
720 723
721 724
722 725 class Correlation(JROData):
723 726
724 727 def __init__(self):
725 728 '''
726 729 Constructor
727 730 '''
728 731 self.radarControllerHeaderObj = RadarControllerHeader()
729 732 self.systemHeaderObj = SystemHeader()
730 733 self.type = "Correlation"
731 734 self.data = None
732 735 self.dtype = None
733 736 self.nProfiles = None
734 737 self.heightList = None
735 738 self.channelList = None
736 739 self.flagNoData = True
737 740 self.flagDiscontinuousBlock = False
738 741 self.utctime = None
739 742 self.timeZone = 0
740 743 self.dstFlag = None
741 744 self.errorCount = None
742 745 self.blocksize = None
743 746 self.flagDecodeData = False # asumo q la data no esta decodificada
744 747 self.flagDeflipData = False # asumo q la data no esta sin flip
745 748 self.pairsList = None
746 749 self.nPoints = None
747 750
748 751 def getPairsList(self):
749 752
750 753 return self.pairsList
751 754
752 755 def getNoise(self, mode=2):
753 756
754 757 indR = numpy.where(self.lagR == 0)[0][0]
755 758 indT = numpy.where(self.lagT == 0)[0][0]
756 759
757 760 jspectra0 = self.data_corr[:, :, indR, :]
758 761 jspectra = copy.copy(jspectra0)
759 762
760 763 num_chan = jspectra.shape[0]
761 764 num_hei = jspectra.shape[2]
762 765
763 766 freq_dc = jspectra.shape[1] / 2
764 767 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
765 768
766 769 if ind_vel[0] < 0:
767 770 ind_vel[list(range(0, 1))] = ind_vel[list(
768 771 range(0, 1))] + self.num_prof
769 772
770 773 if mode == 1:
771 774 jspectra[:, freq_dc, :] = (
772 775 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
773 776
774 777 if mode == 2:
775 778
776 779 vel = numpy.array([-2, -1, 1, 2])
777 780 xx = numpy.zeros([4, 4])
778 781
779 782 for fil in range(4):
780 783 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
781 784
782 785 xx_inv = numpy.linalg.inv(xx)
783 786 xx_aux = xx_inv[0, :]
784 787
785 788 for ich in range(num_chan):
786 789 yy = jspectra[ich, ind_vel, :]
787 790 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
788 791
789 792 junkid = jspectra[ich, freq_dc, :] <= 0
790 793 cjunkid = sum(junkid)
791 794
792 795 if cjunkid.any():
793 796 jspectra[ich, freq_dc, junkid.nonzero()] = (
794 797 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
795 798
796 799 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
797 800
798 801 return noise
799 802
800 803 @property
801 804 def timeInterval(self):
802 805
803 806 return self.ippSeconds * self.nCohInt * self.nProfiles
804 807
805 808 def splitFunctions(self):
806 809
807 810 pairsList = self.pairsList
808 811 ccf_pairs = []
809 812 acf_pairs = []
810 813 ccf_ind = []
811 814 acf_ind = []
812 815 for l in range(len(pairsList)):
813 816 chan0 = pairsList[l][0]
814 817 chan1 = pairsList[l][1]
815 818
816 819 # Obteniendo pares de Autocorrelacion
817 820 if chan0 == chan1:
818 821 acf_pairs.append(chan0)
819 822 acf_ind.append(l)
820 823 else:
821 824 ccf_pairs.append(pairsList[l])
822 825 ccf_ind.append(l)
823 826
824 827 data_acf = self.data_cf[acf_ind]
825 828 data_ccf = self.data_cf[ccf_ind]
826 829
827 830 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
828 831
829 832 @property
830 833 def normFactor(self):
831 834 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
832 835 acf_pairs = numpy.array(acf_pairs)
833 836 normFactor = numpy.zeros((self.nPairs, self.nHeights))
834 837
835 838 for p in range(self.nPairs):
836 839 pair = self.pairsList[p]
837 840
838 841 ch0 = pair[0]
839 842 ch1 = pair[1]
840 843
841 844 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
842 845 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
843 846 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
844 847
845 848 return normFactor
846 849
847 850
848 851 class Parameters(Spectra):
849 852
850 853 groupList = None # List of Pairs, Groups, etc
851 854 data_param = None # Parameters obtained
852 855 data_pre = None # Data Pre Parametrization
853 856 data_SNR = None # Signal to Noise Ratio
854 857 abscissaList = None # Abscissa, can be velocities, lags or time
855 858 utctimeInit = None # Initial UTC time
856 859 paramInterval = None # Time interval to calculate Parameters in seconds
857 860 useLocalTime = True
858 861 # Fitting
859 862 data_error = None # Error of the estimation
860 863 constants = None
861 864 library = None
862 865 # Output signal
863 866 outputInterval = None # Time interval to calculate output signal in seconds
864 867 data_output = None # Out signal
865 868 nAvg = None
866 869 noise_estimation = None
867 870 GauSPC = None # Fit gaussian SPC
868 871
869 872 def __init__(self):
870 873 '''
871 874 Constructor
872 875 '''
873 876 self.radarControllerHeaderObj = RadarControllerHeader()
874 877 self.systemHeaderObj = SystemHeader()
875 878 self.type = "Parameters"
876 879 self.timeZone = 0
877 880
878 881 def getTimeRange1(self, interval):
879 882
880 883 datatime = []
881 884
882 885 if self.useLocalTime:
883 886 time1 = self.utctimeInit - self.timeZone * 60
884 887 else:
885 888 time1 = self.utctimeInit
886 889
887 890 datatime.append(time1)
888 891 datatime.append(time1 + interval)
889 892 datatime = numpy.array(datatime)
890 893
891 894 return datatime
892 895
893 896 @property
894 897 def timeInterval(self):
895 898
896 899 if hasattr(self, 'timeInterval1'):
897 900 return self.timeInterval1
898 901 else:
899 902 return self.paramInterval
900 903
901 904
902 905 def setValue(self, value):
903 906
904 907 print("This property should not be initialized")
905 908
906 909 return
907 910
908 911 def getNoise(self):
909 912
910 913 return self.spc_noise
911 914
912 915 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
913 916
914 917
915 918 class PlotterData(object):
916 919 '''
917 920 Object to hold data to be plotted
918 921 '''
919 922
920 923 MAXNUMX = 200
921 924 MAXNUMY = 200
922 925
923 926 def __init__(self, code, exp_code, localtime=True):
924 927
925 928 self.key = code
926 929 self.exp_code = exp_code
927 930 self.ready = False
928 931 self.flagNoData = False
929 932 self.localtime = localtime
930 933 self.data = {}
931 934 self.meta = {}
932 935 self.__heights = []
933 936
934 937 def __str__(self):
935 938 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
936 939 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
937 940
938 941 def __len__(self):
939 942 return len(self.data)
940 943
941 944 def __getitem__(self, key):
942 945 if isinstance(key, int):
943 946 return self.data[self.times[key]]
944 947 elif isinstance(key, str):
945 948 ret = numpy.array([self.data[x][key] for x in self.times])
946 949 if ret.ndim > 1:
947 950 ret = numpy.swapaxes(ret, 0, 1)
948 951 return ret
949 952
950 953 def __contains__(self, key):
951 954 return key in self.data[self.min_time]
952 955
953 956 def setup(self):
954 957 '''
955 958 Configure object
956 959 '''
957 960 self.type = ''
958 961 self.ready = False
959 962 del self.data
960 963 self.data = {}
961 964 self.__heights = []
962 965 self.__all_heights = set()
963 966
964 967 def shape(self, key):
965 968 '''
966 969 Get the shape of the one-element data for the given key
967 970 '''
968 971
969 972 if len(self.data[self.min_time][key]):
970 973 return self.data[self.min_time][key].shape
971 974 return (0,)
972 975
973 976 def update(self, data, tm, meta={}):
974 977 '''
975 978 Update data object with new dataOut
976 979 '''
977 980
978 981 self.data[tm] = data
979 982
980 983 for key, value in meta.items():
981 984 setattr(self, key, value)
982 985
983 986 def normalize_heights(self):
984 987 '''
985 988 Ensure same-dimension of the data for different heighList
986 989 '''
987 990
988 991 H = numpy.array(list(self.__all_heights))
989 992 H.sort()
990 993 for key in self.data:
991 994 shape = self.shape(key)[:-1] + H.shape
992 995 for tm, obj in list(self.data[key].items()):
993 996 h = self.__heights[self.times.tolist().index(tm)]
994 997 if H.size == h.size:
995 998 continue
996 999 index = numpy.where(numpy.in1d(H, h))[0]
997 1000 dummy = numpy.zeros(shape) + numpy.nan
998 1001 if len(shape) == 2:
999 1002 dummy[:, index] = obj
1000 1003 else:
1001 1004 dummy[index] = obj
1002 1005 self.data[key][tm] = dummy
1003 1006
1004 1007 self.__heights = [H for tm in self.times]
1005 1008
1006 1009 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1007 1010 '''
1008 1011 Convert data to json
1009 1012 '''
1010 1013
1011 1014 meta = {}
1012 1015 meta['xrange'] = []
1013 1016 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1014 1017 tmp = self.data[tm][self.key]
1015 1018 shape = tmp.shape
1016 1019 if len(shape) == 2:
1017 1020 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1018 1021 elif len(shape) == 3:
1019 1022 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1020 1023 data = self.roundFloats(
1021 1024 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1022 1025 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1023 1026 else:
1024 1027 data = self.roundFloats(self.data[tm][self.key].tolist())
1025 1028
1026 1029 ret = {
1027 1030 'plot': plot_name,
1028 1031 'code': self.exp_code,
1029 1032 'time': float(tm),
1030 1033 'data': data,
1031 1034 }
1032 1035 meta['type'] = plot_type
1033 1036 meta['interval'] = float(self.interval)
1034 1037 meta['localtime'] = self.localtime
1035 1038 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1036 1039 meta.update(self.meta)
1037 1040 ret['metadata'] = meta
1038 1041 return json.dumps(ret)
1039 1042
1040 1043 @property
1041 1044 def times(self):
1042 1045 '''
1043 1046 Return the list of times of the current data
1044 1047 '''
1045 1048
1046 1049 ret = [t for t in self.data]
1047 1050 ret.sort()
1048 1051 return numpy.array(ret)
1049 1052
1050 1053 @property
1051 1054 def min_time(self):
1052 1055 '''
1053 1056 Return the minimun time value
1054 1057 '''
1055 1058
1056 1059 return self.times[0]
1057 1060
1058 1061 @property
1059 1062 def max_time(self):
1060 1063 '''
1061 1064 Return the maximun time value
1062 1065 '''
1063 1066
1064 1067 return self.times[-1]
1065 1068
1066 1069 # @property
1067 1070 # def heights(self):
1068 1071 # '''
1069 1072 # Return the list of heights of the current data
1070 1073 # '''
1071 1074
1072 1075 # return numpy.array(self.__heights[-1])
1073 1076
1074 1077 @staticmethod
1075 1078 def roundFloats(obj):
1076 1079 if isinstance(obj, list):
1077 1080 return list(map(PlotterData.roundFloats, obj))
1078 1081 elif isinstance(obj, float):
1079 1082 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,1344 +1,1348
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 #import collections.abc
12 12
13 13 from schainpy.model.graphics.jroplot_base import Plot, plt, log
14 14
15 15 class SpectraPlot(Plot):
16 16 '''
17 17 Plot for Spectra data
18 18 '''
19 19
20 20 CODE = 'spc'
21 21 colormap = 'jet'
22 22 plot_type = 'pcolor'
23 23 buffering = False
24 24
25 25 def setup(self):
26 26
27 27 self.nplots = len(self.data.channels)
28 28 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
29 29 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
30 30 self.height = 2.6 * self.nrows
31 31 self.cb_label = 'dB'
32 32 if self.showprofile:
33 33 self.width = 4 * self.ncols
34 34 else:
35 35 self.width = 3.5 * self.ncols
36 36 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
37 37 self.ylabel = 'Range [km]'
38 38
39 39 def update(self, dataOut):
40 40
41 41 data = {}
42 42 meta = {}
43 43
44 44 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
45 #print("dataOut.normFactor: ", dataOut.normFactor)
46 #print("spc: ", dataOut.data_spc[0,0,0])
47 #spc = 10*numpy.log10(dataOut.data_spc)
45 48 #print("Spc: ",spc[0])
46 49 #exit(1)
47 50 data['spc'] = spc
48 51 data['rti'] = dataOut.getPower()
49 52 #print(data['rti'][0])
50 53 #exit(1)
51 54 #print("NormFactor: ",dataOut.normFactor)
52 55 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
53 56 if hasattr(dataOut, 'LagPlot'): #Double Pulse
54 57 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
55 58 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=46,ymax_index=max_hei_id)/dataOut.normFactor)
56 59 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=40,ymax_index=max_hei_id)/dataOut.normFactor)
57 60 data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)/dataOut.normFactor)
58 61 data['noise'][0] = 10*numpy.log10(dataOut.getNoise(ymin_index=53)[0]/dataOut.normFactor)
59 62 #data['noise'][1] = 22.035507
60 63 else:
61 64 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
62 65 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=26,ymax_index=44)/dataOut.normFactor)
63 66 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
64 67
65 68 if self.CODE == 'spc_moments':
66 69 data['moments'] = dataOut.moments
67 70 if self.CODE == 'gaussian_fit':
68 71 data['gaussfit'] = dataOut.DGauFitParams
69 72
70 73 return data, meta
71 74
72 75 def plot(self):
73 76
74 77 if self.xaxis == "frequency":
75 78 x = self.data.xrange[0]
76 79 self.xlabel = "Frequency (kHz)"
77 80 elif self.xaxis == "time":
78 81 x = self.data.xrange[1]
79 82 self.xlabel = "Time (ms)"
80 83 else:
81 84 x = self.data.xrange[2]
82 85 self.xlabel = "Velocity (m/s)"
83 86
84 87 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
85 88 x = self.data.xrange[2]
86 89 self.xlabel = "Velocity (m/s)"
87 90
88 91 self.titles = []
89 92
90 93 y = self.data.yrange
91 94 self.y = y
92 95
93 96 data = self.data[-1]
94 97 z = data['spc']
95 98
96 99 self.CODE2 = 'spc_oblique'
97 100
98 101 for n, ax in enumerate(self.axes):
99 102 noise = data['noise'][n]
100 103 if self.CODE == 'spc_moments':
101 104 mean = data['moments'][n, 1]
102 105 if self.CODE == 'gaussian_fit':
103 106 gau0 = data['gaussfit'][n][2,:,0]
104 107 gau1 = data['gaussfit'][n][2,:,1]
105 108 if ax.firsttime:
106 109 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
107 110 self.xmin = self.xmin if self.xmin else numpy.nanmin(x)#-self.xmax
108 111 #self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
109 112 #self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
110 113 if self.zlimits is not None:
111 114 self.zmin, self.zmax = self.zlimits[n]
112 115
113 116 ax.plt = ax.pcolormesh(x, y, z[n].T,
114 117 vmin=self.zmin,
115 118 vmax=self.zmax,
116 119 cmap=plt.get_cmap(self.colormap),
117 120 )
118 121
119 122 if self.showprofile:
120 123 ax.plt_profile = self.pf_axes[n].plot(
121 124 data['rti'][n], y)[0]
122 125 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
123 126 color="k", linestyle="dashed", lw=1)[0]
124 127 if self.CODE == 'spc_moments':
125 128 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
126 129 if self.CODE == 'gaussian_fit':
127 130 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
128 131 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
129 132 else:
130 133 if self.zlimits is not None:
131 134 self.zmin, self.zmax = self.zlimits[n]
132 135 ax.plt.set_array(z[n].T.ravel())
133 136 if self.showprofile:
134 137 ax.plt_profile.set_data(data['rti'][n], y)
135 138 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
136 139 if self.CODE == 'spc_moments':
137 140 ax.plt_mean.set_data(mean, y)
138 141 if self.CODE == 'gaussian_fit':
139 142 ax.plt_gau0.set_data(gau0, y)
140 143 ax.plt_gau1.set_data(gau1, y)
141 144 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
142 145
143 146 class SpectraObliquePlot(Plot):
144 147 '''
145 148 Plot for Spectra data
146 149 '''
147 150
148 151 CODE = 'spc_oblique'
149 152 colormap = 'jet'
150 153 plot_type = 'pcolor'
151 154
152 155 def setup(self):
153 156 self.xaxis = "oblique"
154 157 self.nplots = len(self.data.channels)
155 158 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
156 159 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
157 160 self.height = 2.6 * self.nrows
158 161 self.cb_label = 'dB'
159 162 if self.showprofile:
160 163 self.width = 4 * self.ncols
161 164 else:
162 165 self.width = 3.5 * self.ncols
163 166 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
164 167 self.ylabel = 'Range [km]'
165 168
166 169 def update(self, dataOut):
167 170
168 171 data = {}
169 172 meta = {}
170 173
171 174 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
172 175 data['spc'] = spc
173 176 data['rti'] = dataOut.getPower()
174 177 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
175 178 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
176 179 '''
177 180 data['shift1'] = dataOut.Oblique_params[0,-2,:]
178 181 data['shift2'] = dataOut.Oblique_params[0,-1,:]
179 182 data['shift1_error'] = dataOut.Oblique_param_errors[0,-2,:]
180 183 data['shift2_error'] = dataOut.Oblique_param_errors[0,-1,:]
181 184 '''
182 185 '''
183 186 data['shift1'] = dataOut.Oblique_params[0,1,:]
184 187 data['shift2'] = dataOut.Oblique_params[0,4,:]
185 188 data['shift1_error'] = dataOut.Oblique_param_errors[0,1,:]
186 189 data['shift2_error'] = dataOut.Oblique_param_errors[0,4,:]
187 190 '''
188 191 data['shift1'] = dataOut.Dop_EEJ_T1[0]
189 192 data['shift2'] = dataOut.Dop_EEJ_T2[0]
190 193 data['max_val_2'] = dataOut.Oblique_params[0,-1,:]
191 194 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
192 195 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
193 196
194 197 return data, meta
195 198
196 199 def plot(self):
197 200
198 201 if self.xaxis == "frequency":
199 202 x = self.data.xrange[0]
200 203 self.xlabel = "Frequency (kHz)"
201 204 elif self.xaxis == "time":
202 205 x = self.data.xrange[1]
203 206 self.xlabel = "Time (ms)"
204 207 else:
205 208 x = self.data.xrange[2]
206 209 self.xlabel = "Velocity (m/s)"
207 210
208 211 self.titles = []
209 212
210 213 y = self.data.yrange
211 214 self.y = y
212 215
213 216 data = self.data[-1]
214 217 z = data['spc']
215 218
216 219 for n, ax in enumerate(self.axes):
217 220 noise = self.data['noise'][n][-1]
218 221 shift1 = data['shift1']
219 222 #print(shift1)
220 223 shift2 = data['shift2']
221 224 max_val_2 = data['max_val_2']
222 225 err1 = data['shift1_error']
223 226 err2 = data['shift2_error']
224 227 if ax.firsttime:
225 228
226 229 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
227 230 self.xmin = self.xmin if self.xmin else -self.xmax
228 231 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
229 232 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
230 233 ax.plt = ax.pcolormesh(x, y, z[n].T,
231 234 vmin=self.zmin,
232 235 vmax=self.zmax,
233 236 cmap=plt.get_cmap(self.colormap)
234 237 )
235 238
236 239 if self.showprofile:
237 240 ax.plt_profile = self.pf_axes[n].plot(
238 241 self.data['rti'][n][-1], y)[0]
239 242 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
240 243 color="k", linestyle="dashed", lw=1)[0]
241 244
242 245 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)
243 246 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)
244 247 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
245 248
246 249 #print("plotter1: ", self.ploterr1,shift1)
247 250
248 251 else:
249 252 #print("else plotter1: ", self.ploterr1,shift1)
250 253 self.ploterr1.remove()
251 254 self.ploterr2.remove()
252 255 self.ploterr3.remove()
253 256 ax.plt.set_array(z[n].T.ravel())
254 257 if self.showprofile:
255 258 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
256 259 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
257 260 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)
258 261 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)
259 262 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
260 263
261 264 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
262 265
263 266
264 267 class CrossSpectraPlot(Plot):
265 268
266 269 CODE = 'cspc'
267 270 colormap = 'jet'
268 271 plot_type = 'pcolor'
269 272 zmin_coh = None
270 273 zmax_coh = None
271 274 zmin_phase = None
272 275 zmax_phase = None
273 276
274 277 def setup(self):
275 278
276 279 self.ncols = 4
277 280 self.nplots = len(self.data.pairs) * 2
278 281 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
279 282 self.width = 3.1 * self.ncols
280 283 self.height = 5 * self.nrows
281 284 self.ylabel = 'Range [km]'
282 285 self.showprofile = False
283 286 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
284 287
285 288 def update(self, dataOut):
286 289
287 290 data = {}
288 291 meta = {}
289 292
290 293 spc = dataOut.data_spc
291 294 cspc = dataOut.data_cspc
292 295 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
293 296 meta['pairs'] = dataOut.pairsList
294 297
295 298 tmp = []
296 299
297 300 for n, pair in enumerate(meta['pairs']):
298 301 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
299 302 coh = numpy.abs(out)
300 303 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
301 304 tmp.append(coh)
302 305 tmp.append(phase)
303 306
304 307 data['cspc'] = numpy.array(tmp)
305 308
306 309 return data, meta
307 310
308 311 def plot(self):
309 312
310 313 if self.xaxis == "frequency":
311 314 x = self.data.xrange[0]
312 315 self.xlabel = "Frequency (kHz)"
313 316 elif self.xaxis == "time":
314 317 x = self.data.xrange[1]
315 318 self.xlabel = "Time (ms)"
316 319 else:
317 320 x = self.data.xrange[2]
318 321 self.xlabel = "Velocity (m/s)"
319 322
320 323 self.titles = []
321 324
322 325 y = self.data.yrange
323 326 self.y = y
324 327
325 328 data = self.data[-1]
326 329 cspc = data['cspc']
327 330
328 331 for n in range(len(self.data.pairs)):
329 332 pair = self.data.pairs[n]
330 333 coh = cspc[n*2]
331 334 phase = cspc[n*2+1]
332 335 ax = self.axes[2 * n]
333 336 if ax.firsttime:
334 337 ax.plt = ax.pcolormesh(x, y, coh.T,
335 338 vmin=0,
336 339 vmax=1,
337 340 cmap=plt.get_cmap(self.colormap_coh)
338 341 )
339 342 else:
340 343 ax.plt.set_array(coh.T.ravel())
341 344 self.titles.append(
342 345 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
343 346
344 347 ax = self.axes[2 * n + 1]
345 348 if ax.firsttime:
346 349 ax.plt = ax.pcolormesh(x, y, phase.T,
347 350 vmin=-180,
348 351 vmax=180,
349 352 cmap=plt.get_cmap(self.colormap_phase)
350 353 )
351 354 else:
352 355 ax.plt.set_array(phase.T.ravel())
353 356 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
354 357
355 358
356 359 class CrossSpectra4Plot(Plot):
357 360
358 361 CODE = 'cspc'
359 362 colormap = 'jet'
360 363 plot_type = 'pcolor'
361 364 zmin_coh = None
362 365 zmax_coh = None
363 366 zmin_phase = None
364 367 zmax_phase = None
365 368
366 369 def setup(self):
367 370
368 371 self.ncols = 4
369 372 self.nrows = len(self.data.pairs)
370 373 self.nplots = self.nrows * 4
371 374 self.width = 3.1 * self.ncols
372 375 self.height = 5 * self.nrows
373 376 self.ylabel = 'Range [km]'
374 377 self.showprofile = False
375 378 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
376 379
377 380 def plot(self):
378 381
379 382 if self.xaxis == "frequency":
380 383 x = self.data.xrange[0]
381 384 self.xlabel = "Frequency (kHz)"
382 385 elif self.xaxis == "time":
383 386 x = self.data.xrange[1]
384 387 self.xlabel = "Time (ms)"
385 388 else:
386 389 x = self.data.xrange[2]
387 390 self.xlabel = "Velocity (m/s)"
388 391
389 392 self.titles = []
390 393
391 394
392 395 y = self.data.heights
393 396 self.y = y
394 397 nspc = self.data['spc']
395 398 #print(numpy.shape(self.data['spc']))
396 399 spc = self.data['cspc'][0]
397 400 #print(numpy.shape(nspc))
398 401 #exit()
399 402 #nspc[1,:,:] = numpy.flip(nspc[1,:,:],axis=0)
400 403 #print(numpy.shape(spc))
401 404 #exit()
402 405 cspc = self.data['cspc'][1]
403 406
404 407 #xflip=numpy.flip(x)
405 408 #print(numpy.shape(cspc))
406 409 #exit()
407 410
408 411 for n in range(self.nrows):
409 412 noise = self.data['noise'][:,-1]
410 413 pair = self.data.pairs[n]
411 414 #print(pair)
412 415 #exit()
413 416 ax = self.axes[4 * n]
414 417 if ax.firsttime:
415 418 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
416 419 self.xmin = self.xmin if self.xmin else -self.xmax
417 420 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
418 421 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
419 422 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
420 423 vmin=self.zmin,
421 424 vmax=self.zmax,
422 425 cmap=plt.get_cmap(self.colormap)
423 426 )
424 427 else:
425 428 #print(numpy.shape(nspc[pair[0]].T))
426 429 #exit()
427 430 ax.plt.set_array(nspc[pair[0]].T.ravel())
428 431 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
429 432
430 433 ax = self.axes[4 * n + 1]
431 434
432 435 if ax.firsttime:
433 436 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
434 437 vmin=self.zmin,
435 438 vmax=self.zmax,
436 439 cmap=plt.get_cmap(self.colormap)
437 440 )
438 441 else:
439 442
440 443 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
441 444 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
442 445
443 446 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
444 447 coh = numpy.abs(out)
445 448 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
446 449
447 450 ax = self.axes[4 * n + 2]
448 451 if ax.firsttime:
449 452 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
450 453 vmin=0,
451 454 vmax=1,
452 455 cmap=plt.get_cmap(self.colormap_coh)
453 456 )
454 457 else:
455 458 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
456 459 self.titles.append(
457 460 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
458 461
459 462 ax = self.axes[4 * n + 3]
460 463 if ax.firsttime:
461 464 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
462 465 vmin=-180,
463 466 vmax=180,
464 467 cmap=plt.get_cmap(self.colormap_phase)
465 468 )
466 469 else:
467 470 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
468 471 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
469 472
470 473
471 474 class CrossSpectra2Plot(Plot):
472 475
473 476 CODE = 'cspc'
474 477 colormap = 'jet'
475 478 plot_type = 'pcolor'
476 479 zmin_coh = None
477 480 zmax_coh = None
478 481 zmin_phase = None
479 482 zmax_phase = None
480 483
481 484 def setup(self):
482 485
483 486 self.ncols = 1
484 487 self.nrows = len(self.data.pairs)
485 488 self.nplots = self.nrows * 1
486 489 self.width = 3.1 * self.ncols
487 490 self.height = 5 * self.nrows
488 491 self.ylabel = 'Range [km]'
489 492 self.showprofile = False
490 493 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
491 494
492 495 def plot(self):
493 496
494 497 if self.xaxis == "frequency":
495 498 x = self.data.xrange[0]
496 499 self.xlabel = "Frequency (kHz)"
497 500 elif self.xaxis == "time":
498 501 x = self.data.xrange[1]
499 502 self.xlabel = "Time (ms)"
500 503 else:
501 504 x = self.data.xrange[2]
502 505 self.xlabel = "Velocity (m/s)"
503 506
504 507 self.titles = []
505 508
506 509
507 510 y = self.data.heights
508 511 self.y = y
509 512 #nspc = self.data['spc']
510 513 #print(numpy.shape(self.data['spc']))
511 514 #spc = self.data['cspc'][0]
512 515 #print(numpy.shape(spc))
513 516 #exit()
514 517 cspc = self.data['cspc'][1]
515 518 #print(numpy.shape(cspc))
516 519 #exit()
517 520
518 521 for n in range(self.nrows):
519 522 noise = self.data['noise'][:,-1]
520 523 pair = self.data.pairs[n]
521 524 #print(pair) #exit()
522 525
523 526
524 527
525 528 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
526 529
527 530 #print(out[:,53])
528 531 #exit()
529 532 cross = numpy.abs(out)
530 533 z = cross/self.data.nFactor
531 534 #print("here")
532 535 #print(dataOut.data_spc[0,0,0])
533 536 #exit()
534 537
535 538 cross = 10*numpy.log10(z)
536 539 #print(numpy.shape(cross))
537 540 #print(cross[0,:])
538 541 #print(self.data.nFactor)
539 542 #exit()
540 543 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
541 544
542 545 ax = self.axes[1 * n]
543 546 if ax.firsttime:
544 547 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
545 548 self.xmin = self.xmin if self.xmin else -self.xmax
546 549 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
547 550 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
548 551 ax.plt = ax.pcolormesh(x, y, cross.T,
549 552 vmin=self.zmin,
550 553 vmax=self.zmax,
551 554 cmap=plt.get_cmap(self.colormap)
552 555 )
553 556 else:
554 557 ax.plt.set_array(cross.T.ravel())
555 558 self.titles.append(
556 559 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
557 560
558 561
559 562 class CrossSpectra3Plot(Plot):
560 563
561 564 CODE = 'cspc'
562 565 colormap = 'jet'
563 566 plot_type = 'pcolor'
564 567 zmin_coh = None
565 568 zmax_coh = None
566 569 zmin_phase = None
567 570 zmax_phase = None
568 571
569 572 def setup(self):
570 573
571 574 self.ncols = 3
572 575 self.nrows = len(self.data.pairs)
573 576 self.nplots = self.nrows * 3
574 577 self.width = 3.1 * self.ncols
575 578 self.height = 5 * self.nrows
576 579 self.ylabel = 'Range [km]'
577 580 self.showprofile = False
578 581 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
579 582
580 583 def plot(self):
581 584
582 585 if self.xaxis == "frequency":
583 586 x = self.data.xrange[0]
584 587 self.xlabel = "Frequency (kHz)"
585 588 elif self.xaxis == "time":
586 589 x = self.data.xrange[1]
587 590 self.xlabel = "Time (ms)"
588 591 else:
589 592 x = self.data.xrange[2]
590 593 self.xlabel = "Velocity (m/s)"
591 594
592 595 self.titles = []
593 596
594 597
595 598 y = self.data.heights
596 599 self.y = y
597 600 #nspc = self.data['spc']
598 601 #print(numpy.shape(self.data['spc']))
599 602 #spc = self.data['cspc'][0]
600 603 #print(numpy.shape(spc))
601 604 #exit()
602 605 cspc = self.data['cspc'][1]
603 606 #print(numpy.shape(cspc))
604 607 #exit()
605 608
606 609 for n in range(self.nrows):
607 610 noise = self.data['noise'][:,-1]
608 611 pair = self.data.pairs[n]
609 612 #print(pair) #exit()
610 613
611 614
612 615
613 616 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
614 617
615 618 #print(out[:,53])
616 619 #exit()
617 620 cross = numpy.abs(out)
618 621 z = cross/self.data.nFactor
619 622 cross = 10*numpy.log10(z)
620 623
621 624 out_r= out.real/self.data.nFactor
622 625 #out_r = 10*numpy.log10(out_r)
623 626
624 627 out_i= out.imag/self.data.nFactor
625 628 #out_i = 10*numpy.log10(out_i)
626 629 #print(numpy.shape(cross))
627 630 #print(cross[0,:])
628 631 #print(self.data.nFactor)
629 632 #exit()
630 633 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
631 634
632 635 ax = self.axes[3 * n]
633 636 if ax.firsttime:
634 637 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
635 638 self.xmin = self.xmin if self.xmin else -self.xmax
636 639 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
637 640 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
638 641 ax.plt = ax.pcolormesh(x, y, cross.T,
639 642 vmin=self.zmin,
640 643 vmax=self.zmax,
641 644 cmap=plt.get_cmap(self.colormap)
642 645 )
643 646 else:
644 647 ax.plt.set_array(cross.T.ravel())
645 648 self.titles.append(
646 649 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
647 650
648 651 ax = self.axes[3 * n + 1]
649 652 if ax.firsttime:
650 653 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
651 654 self.xmin = self.xmin if self.xmin else -self.xmax
652 655 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
653 656 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
654 657 ax.plt = ax.pcolormesh(x, y, out_r.T,
655 658 vmin=-1.e6,
656 659 vmax=0,
657 660 cmap=plt.get_cmap(self.colormap)
658 661 )
659 662 else:
660 663 ax.plt.set_array(out_r.T.ravel())
661 664 self.titles.append(
662 665 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
663 666
664 667 ax = self.axes[3 * n + 2]
665 668
666 669
667 670 if ax.firsttime:
668 671 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
669 672 self.xmin = self.xmin if self.xmin else -self.xmax
670 673 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
671 674 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
672 675 ax.plt = ax.pcolormesh(x, y, out_i.T,
673 676 vmin=-1.e6,
674 677 vmax=1.e6,
675 678 cmap=plt.get_cmap(self.colormap)
676 679 )
677 680 else:
678 681 ax.plt.set_array(out_i.T.ravel())
679 682 self.titles.append(
680 683 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
681 684
682 685 class RTIPlot(Plot):
683 686 '''
684 687 Plot for RTI data
685 688 '''
686 689
687 690 CODE = 'rti'
688 691 colormap = 'jet'
689 692 plot_type = 'pcolorbuffer'
690 693
691 694 def setup(self):
692 695 self.xaxis = 'time'
693 696 self.ncols = 1
694 697 self.nrows = len(self.data.channels)
695 698 self.nplots = len(self.data.channels)
696 699 self.ylabel = 'Range [km]'
697 700 self.xlabel = 'Time'
698 701 self.cb_label = 'dB'
699 702 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
700 703 self.titles = ['{} Channel {}'.format(
701 704 self.CODE.upper(), x) for x in range(self.nrows)]
702 705
703 706 def update(self, dataOut):
704 707
705 708 data = {}
706 709 meta = {}
707 710 data['rti'] = dataOut.getPower()
708 711 #print(numpy.shape(data['rti']))
709 712
710 713 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
711 714
712 715 return data, meta
713 716
714 717 def plot(self):
715 718
716 719 self.x = self.data.times
717 720 self.y = self.data.yrange
718 721 self.z = self.data[self.CODE]
719
722 #print("Inside RTI: ", self.z)
720 723 self.z = numpy.ma.masked_invalid(self.z)
721 724
722 725 if self.decimation is None:
723 726 x, y, z = self.fill_gaps(self.x, self.y, self.z)
724 727 else:
725 728 x, y, z = self.fill_gaps(*self.decimate())
726
729 #print("self.z: ", self.z)
730 #exit(1)
727 731 '''
728 732 if not isinstance(self.zmin, collections.abc.Sequence):
729 733 if not self.zmin:
730 734 self.zmin = [numpy.min(self.z)]*len(self.axes)
731 735 else:
732 736 self.zmin = [self.zmin]*len(self.axes)
733 737
734 738 if not isinstance(self.zmax, collections.abc.Sequence):
735 739 if not self.zmax:
736 740 self.zmax = [numpy.max(self.z)]*len(self.axes)
737 741 else:
738 742 self.zmax = [self.zmax]*len(self.axes)
739 743 '''
740 744 for n, ax in enumerate(self.axes):
741 745
742 746 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
743 747 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
744 748
745 749 if ax.firsttime:
746 750 if self.zlimits is not None:
747 751 self.zmin, self.zmax = self.zlimits[n]
748 752 ax.plt = ax.pcolormesh(x, y, z[n].T,
749 753 vmin=self.zmin,
750 754 vmax=self.zmax,
751 755 cmap=plt.get_cmap(self.colormap)
752 756 )
753 757 if self.showprofile:
754 758 ax.plot_profile = self.pf_axes[n].plot(
755 759 self.data['rti'][n][-1], self.y)[0]
756 760 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
757 761 color="k", linestyle="dashed", lw=1)[0]
758 762 else:
759 763 if self.zlimits is not None:
760 764 self.zmin, self.zmax = self.zlimits[n]
761 765 ax.collections.remove(ax.collections[0])
762 766 ax.plt = ax.pcolormesh(x, y, z[n].T,
763 767 vmin=self.zmin,
764 768 vmax=self.zmax,
765 769 cmap=plt.get_cmap(self.colormap)
766 770 )
767 771 if self.showprofile:
768 772 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
769 773 ax.plot_noise.set_data(numpy.repeat(
770 774 self.data['noise'][n][-1], len(self.y)), self.y)
771 775
772 776
773 777 class SpectrogramPlot(Plot):
774 778 '''
775 779 Plot for Spectrogram data
776 780 '''
777 781
778 782 CODE = 'Spectrogram_Profile'
779 783 colormap = 'binary'
780 784 plot_type = 'pcolorbuffer'
781 785
782 786 def setup(self):
783 787 self.xaxis = 'time'
784 788 self.ncols = 1
785 789 self.nrows = len(self.data.channels)
786 790 self.nplots = len(self.data.channels)
787 791 self.xlabel = 'Time'
788 792 #self.cb_label = 'dB'
789 793 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
790 794 self.titles = []
791 795
792 796 #self.titles = ['{} Channel {} \n H = {} km ({} - {})'.format(
793 797 #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)]
794 798
795 799 self.titles = ['{} Channel {}'.format(
796 800 self.CODE.upper(), x) for x in range(self.nrows)]
797 801
798 802
799 803 def update(self, dataOut):
800 804 data = {}
801 805 meta = {}
802 806
803 807 maxHei = 1620#+12000
804 808 indb = numpy.where(dataOut.heightList <= maxHei)
805 809 hei = indb[0][-1]
806 810 #print(dataOut.heightList)
807 811
808 812 factor = dataOut.nIncohInt
809 813 z = dataOut.data_spc[:,:,hei] / factor
810 814 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
811 815 #buffer = 10 * numpy.log10(z)
812 816
813 817 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
814 818
815 819
816 820 #self.hei = hei
817 821 #self.heightList = dataOut.heightList
818 822 #self.DH = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
819 823 #self.nProfiles = dataOut.nProfiles
820 824
821 825 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
822 826
823 827 data['hei'] = hei
824 828 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
825 829 data['nProfiles'] = dataOut.nProfiles
826 830 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
827 831 '''
828 832 import matplotlib.pyplot as plt
829 833 plt.plot(10 * numpy.log10(z[0,:]))
830 834 plt.show()
831 835
832 836 from time import sleep
833 837 sleep(10)
834 838 '''
835 839 return data, meta
836 840
837 841 def plot(self):
838 842
839 843 self.x = self.data.times
840 844 self.z = self.data[self.CODE]
841 845 self.y = self.data.xrange[0]
842 846
843 847 hei = self.data['hei'][-1]
844 848 DH = self.data['DH'][-1]
845 849 nProfiles = self.data['nProfiles'][-1]
846 850
847 851 self.ylabel = "Frequency (kHz)"
848 852
849 853 self.z = numpy.ma.masked_invalid(self.z)
850 854
851 855 if self.decimation is None:
852 856 x, y, z = self.fill_gaps(self.x, self.y, self.z)
853 857 else:
854 858 x, y, z = self.fill_gaps(*self.decimate())
855 859
856 860 for n, ax in enumerate(self.axes):
857 861 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
858 862 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
859 863 data = self.data[-1]
860 864 if ax.firsttime:
861 865 ax.plt = ax.pcolormesh(x, y, z[n].T,
862 866 vmin=self.zmin,
863 867 vmax=self.zmax,
864 868 cmap=plt.get_cmap(self.colormap)
865 869 )
866 870 else:
867 871 ax.collections.remove(ax.collections[0])
868 872 ax.plt = ax.pcolormesh(x, y, z[n].T,
869 873 vmin=self.zmin,
870 874 vmax=self.zmax,
871 875 cmap=plt.get_cmap(self.colormap)
872 876 )
873 877
874 878 #self.titles.append('Spectrogram')
875 879
876 880 #self.titles.append('{} Channel {} \n H = {} km ({} - {})'.format(
877 881 #self.CODE.upper(), x, y[hei], y[hei],y[hei]+(DH*nProfiles)))
878 882
879 883
880 884
881 885
882 886 class CoherencePlot(RTIPlot):
883 887 '''
884 888 Plot for Coherence data
885 889 '''
886 890
887 891 CODE = 'coh'
888 892
889 893 def setup(self):
890 894 self.xaxis = 'time'
891 895 self.ncols = 1
892 896 self.nrows = len(self.data.pairs)
893 897 self.nplots = len(self.data.pairs)
894 898 self.ylabel = 'Range [km]'
895 899 self.xlabel = 'Time'
896 900 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
897 901 if self.CODE == 'coh':
898 902 self.cb_label = ''
899 903 self.titles = [
900 904 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
901 905 else:
902 906 self.cb_label = 'Degrees'
903 907 self.titles = [
904 908 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
905 909
906 910 def update(self, dataOut):
907 911
908 912 data = {}
909 913 meta = {}
910 914 data['coh'] = dataOut.getCoherence()
911 915 meta['pairs'] = dataOut.pairsList
912 916
913 917 return data, meta
914 918
915 919 class PhasePlot(CoherencePlot):
916 920 '''
917 921 Plot for Phase map data
918 922 '''
919 923
920 924 CODE = 'phase'
921 925 colormap = 'seismic'
922 926
923 927 def update(self, dataOut):
924 928
925 929 data = {}
926 930 meta = {}
927 931 data['phase'] = dataOut.getCoherence(phase=True)
928 932 meta['pairs'] = dataOut.pairsList
929 933
930 934 return data, meta
931 935
932 936 class NoisePlot(Plot):
933 937 '''
934 938 Plot for noise
935 939 '''
936 940
937 941 CODE = 'noise'
938 942 plot_type = 'scatterbuffer'
939 943
940 944 def setup(self):
941 945 self.xaxis = 'time'
942 946 self.ncols = 1
943 947 self.nrows = 1
944 948 self.nplots = 1
945 949 self.ylabel = 'Intensity [dB]'
946 950 self.xlabel = 'Time'
947 951 self.titles = ['Noise']
948 952 self.colorbar = False
949 953 self.plots_adjust.update({'right': 0.85 })
950 954
951 955 def update(self, dataOut):
952 956
953 957 data = {}
954 958 meta = {}
955 959 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
956 960 meta['yrange'] = numpy.array([])
957 961
958 962 return data, meta
959 963
960 964 def plot(self):
961 965
962 966 x = self.data.times
963 967 xmin = self.data.min_time
964 968 xmax = xmin + self.xrange * 60 * 60
965 969 Y = self.data['noise']
966 970
967 971 if self.axes[0].firsttime:
968 972 self.ymin = numpy.nanmin(Y) - 5
969 973 self.ymax = numpy.nanmax(Y) + 5
970 974 for ch in self.data.channels:
971 975 y = Y[ch]
972 976 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
973 977 plt.legend(bbox_to_anchor=(1.18, 1.0))
974 978 else:
975 979 for ch in self.data.channels:
976 980 y = Y[ch]
977 981 self.axes[0].lines[ch].set_data(x, y)
978 982
979 983 self.ymin = numpy.nanmin(Y) - 5
980 984 self.ymax = numpy.nanmax(Y) + 10
981 985
982 986
983 987 class PowerProfilePlot(Plot):
984 988
985 989 CODE = 'pow_profile'
986 990 plot_type = 'scatter'
987 991
988 992 def setup(self):
989 993
990 994 self.ncols = 1
991 995 self.nrows = 1
992 996 self.nplots = 1
993 997 self.height = 4
994 998 self.width = 3
995 999 self.ylabel = 'Range [km]'
996 1000 self.xlabel = 'Intensity [dB]'
997 1001 self.titles = ['Power Profile']
998 1002 self.colorbar = False
999 1003
1000 1004 def update(self, dataOut):
1001 1005
1002 1006 data = {}
1003 1007 meta = {}
1004 1008 data[self.CODE] = dataOut.getPower()
1005 1009
1006 1010 return data, meta
1007 1011
1008 1012 def plot(self):
1009 1013
1010 1014 y = self.data.yrange
1011 1015 self.y = y
1012 1016
1013 1017 x = self.data[-1][self.CODE]
1014 1018
1015 1019 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
1016 1020 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
1017 1021
1018 1022 if self.axes[0].firsttime:
1019 1023 for ch in self.data.channels:
1020 1024 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
1021 1025 plt.legend()
1022 1026 else:
1023 1027 for ch in self.data.channels:
1024 1028 self.axes[0].lines[ch].set_data(x[ch], y)
1025 1029
1026 1030
1027 1031 class SpectraCutPlot(Plot):
1028 1032
1029 1033 CODE = 'spc_cut'
1030 1034 plot_type = 'scatter'
1031 1035 buffering = False
1032 1036
1033 1037 def setup(self):
1034 1038
1035 1039 self.nplots = len(self.data.channels)
1036 1040 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
1037 1041 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
1038 1042 self.width = 3.4 * self.ncols + 1.5
1039 1043 self.height = 3 * self.nrows
1040 1044 self.ylabel = 'Power [dB]'
1041 1045 self.colorbar = False
1042 1046 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
1043 1047
1044 1048 def update(self, dataOut):
1045 1049
1046 1050 data = {}
1047 1051 meta = {}
1048 1052 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
1049 1053 data['spc'] = spc
1050 1054 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
1051 1055 if self.CODE == 'cut_gaussian_fit':
1052 1056 data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
1053 1057 data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
1054 1058 return data, meta
1055 1059
1056 1060 def plot(self):
1057 1061 if self.xaxis == "frequency":
1058 1062 x = self.data.xrange[0][1:]
1059 1063 self.xlabel = "Frequency (kHz)"
1060 1064 elif self.xaxis == "time":
1061 1065 x = self.data.xrange[1]
1062 1066 self.xlabel = "Time (ms)"
1063 1067 else:
1064 1068 x = self.data.xrange[2][:-1]
1065 1069 self.xlabel = "Velocity (m/s)"
1066 1070
1067 1071 if self.CODE == 'cut_gaussian_fit':
1068 1072 x = self.data.xrange[2][:-1]
1069 1073 self.xlabel = "Velocity (m/s)"
1070 1074
1071 1075 self.titles = []
1072 1076
1073 1077 y = self.data.yrange
1074 1078 data = self.data[-1]
1075 1079 z = data['spc']
1076 1080
1077 1081 if self.height_index:
1078 1082 index = numpy.array(self.height_index)
1079 1083 else:
1080 1084 index = numpy.arange(0, len(y), int((len(y))/9))
1081 1085
1082 1086 for n, ax in enumerate(self.axes):
1083 1087 if self.CODE == 'cut_gaussian_fit':
1084 1088 gau0 = data['gauss_fit0']
1085 1089 gau1 = data['gauss_fit1']
1086 1090 if ax.firsttime:
1087 1091 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1088 1092 self.xmin = self.xmin if self.xmin else -self.xmax
1089 1093 self.ymin = self.ymin if self.ymin else numpy.nanmin(z[:,:,index])
1090 1094 self.ymax = self.ymax if self.ymax else numpy.nanmax(z[:,:,index])
1091 1095 #print(self.ymax)
1092 1096 #print(z[n, :, index])
1093 1097 ax.plt = ax.plot(x, z[n, :, index].T, lw=0.25)
1094 1098 if self.CODE == 'cut_gaussian_fit':
1095 1099 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
1096 1100 for i, line in enumerate(ax.plt_gau0):
1097 1101 line.set_color(ax.plt[i].get_color())
1098 1102 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
1099 1103 for i, line in enumerate(ax.plt_gau1):
1100 1104 line.set_color(ax.plt[i].get_color())
1101 1105 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1102 1106 self.figures[0].legend(ax.plt, labels, loc='center right')
1103 1107 else:
1104 1108 for i, line in enumerate(ax.plt):
1105 1109 line.set_data(x, z[n, :, index[i]].T)
1106 1110 for i, line in enumerate(ax.plt_gau0):
1107 1111 line.set_data(x, gau0[n, :, index[i]].T)
1108 1112 line.set_color(ax.plt[i].get_color())
1109 1113 for i, line in enumerate(ax.plt_gau1):
1110 1114 line.set_data(x, gau1[n, :, index[i]].T)
1111 1115 line.set_color(ax.plt[i].get_color())
1112 1116 self.titles.append('CH {}'.format(n))
1113 1117
1114 1118
1115 1119 class BeaconPhase(Plot):
1116 1120
1117 1121 __isConfig = None
1118 1122 __nsubplots = None
1119 1123
1120 1124 PREFIX = 'beacon_phase'
1121 1125
1122 1126 def __init__(self):
1123 1127 Plot.__init__(self)
1124 1128 self.timerange = 24*60*60
1125 1129 self.isConfig = False
1126 1130 self.__nsubplots = 1
1127 1131 self.counter_imagwr = 0
1128 1132 self.WIDTH = 800
1129 1133 self.HEIGHT = 400
1130 1134 self.WIDTHPROF = 120
1131 1135 self.HEIGHTPROF = 0
1132 1136 self.xdata = None
1133 1137 self.ydata = None
1134 1138
1135 1139 self.PLOT_CODE = BEACON_CODE
1136 1140
1137 1141 self.FTP_WEI = None
1138 1142 self.EXP_CODE = None
1139 1143 self.SUB_EXP_CODE = None
1140 1144 self.PLOT_POS = None
1141 1145
1142 1146 self.filename_phase = None
1143 1147
1144 1148 self.figfile = None
1145 1149
1146 1150 self.xmin = None
1147 1151 self.xmax = None
1148 1152
1149 1153 def getSubplots(self):
1150 1154
1151 1155 ncol = 1
1152 1156 nrow = 1
1153 1157
1154 1158 return nrow, ncol
1155 1159
1156 1160 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1157 1161
1158 1162 self.__showprofile = showprofile
1159 1163 self.nplots = nplots
1160 1164
1161 1165 ncolspan = 7
1162 1166 colspan = 6
1163 1167 self.__nsubplots = 2
1164 1168
1165 1169 self.createFigure(id = id,
1166 1170 wintitle = wintitle,
1167 1171 widthplot = self.WIDTH+self.WIDTHPROF,
1168 1172 heightplot = self.HEIGHT+self.HEIGHTPROF,
1169 1173 show=show)
1170 1174
1171 1175 nrow, ncol = self.getSubplots()
1172 1176
1173 1177 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1174 1178
1175 1179 def save_phase(self, filename_phase):
1176 1180 f = open(filename_phase,'w+')
1177 1181 f.write('\n\n')
1178 1182 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1179 1183 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1180 1184 f.close()
1181 1185
1182 1186 def save_data(self, filename_phase, data, data_datetime):
1183 1187 f=open(filename_phase,'a')
1184 1188 timetuple_data = data_datetime.timetuple()
1185 1189 day = str(timetuple_data.tm_mday)
1186 1190 month = str(timetuple_data.tm_mon)
1187 1191 year = str(timetuple_data.tm_year)
1188 1192 hour = str(timetuple_data.tm_hour)
1189 1193 minute = str(timetuple_data.tm_min)
1190 1194 second = str(timetuple_data.tm_sec)
1191 1195 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1192 1196 f.close()
1193 1197
1194 1198 def plot(self):
1195 1199 log.warning('TODO: Not yet implemented...')
1196 1200
1197 1201 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1198 1202 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1199 1203 timerange=None,
1200 1204 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1201 1205 server=None, folder=None, username=None, password=None,
1202 1206 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1203 1207
1204 1208 if dataOut.flagNoData:
1205 1209 return dataOut
1206 1210
1207 1211 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1208 1212 return
1209 1213
1210 1214 if pairsList == None:
1211 1215 pairsIndexList = dataOut.pairsIndexList[:10]
1212 1216 else:
1213 1217 pairsIndexList = []
1214 1218 for pair in pairsList:
1215 1219 if pair not in dataOut.pairsList:
1216 1220 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1217 1221 pairsIndexList.append(dataOut.pairsList.index(pair))
1218 1222
1219 1223 if pairsIndexList == []:
1220 1224 return
1221 1225
1222 1226 # if len(pairsIndexList) > 4:
1223 1227 # pairsIndexList = pairsIndexList[0:4]
1224 1228
1225 1229 hmin_index = None
1226 1230 hmax_index = None
1227 1231
1228 1232 if hmin != None and hmax != None:
1229 1233 indexes = numpy.arange(dataOut.nHeights)
1230 1234 hmin_list = indexes[dataOut.heightList >= hmin]
1231 1235 hmax_list = indexes[dataOut.heightList <= hmax]
1232 1236
1233 1237 if hmin_list.any():
1234 1238 hmin_index = hmin_list[0]
1235 1239
1236 1240 if hmax_list.any():
1237 1241 hmax_index = hmax_list[-1]+1
1238 1242
1239 1243 x = dataOut.getTimeRange()
1240 1244
1241 1245 thisDatetime = dataOut.datatime
1242 1246
1243 1247 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1244 1248 xlabel = "Local Time"
1245 1249 ylabel = "Phase (degrees)"
1246 1250
1247 1251 update_figfile = False
1248 1252
1249 1253 nplots = len(pairsIndexList)
1250 1254 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1251 1255 phase_beacon = numpy.zeros(len(pairsIndexList))
1252 1256 for i in range(nplots):
1253 1257 pair = dataOut.pairsList[pairsIndexList[i]]
1254 1258 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1255 1259 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1256 1260 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1257 1261 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1258 1262 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1259 1263
1260 1264 if dataOut.beacon_heiIndexList:
1261 1265 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1262 1266 else:
1263 1267 phase_beacon[i] = numpy.average(phase)
1264 1268
1265 1269 if not self.isConfig:
1266 1270
1267 1271 nplots = len(pairsIndexList)
1268 1272
1269 1273 self.setup(id=id,
1270 1274 nplots=nplots,
1271 1275 wintitle=wintitle,
1272 1276 showprofile=showprofile,
1273 1277 show=show)
1274 1278
1275 1279 if timerange != None:
1276 1280 self.timerange = timerange
1277 1281
1278 1282 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1279 1283
1280 1284 if ymin == None: ymin = 0
1281 1285 if ymax == None: ymax = 360
1282 1286
1283 1287 self.FTP_WEI = ftp_wei
1284 1288 self.EXP_CODE = exp_code
1285 1289 self.SUB_EXP_CODE = sub_exp_code
1286 1290 self.PLOT_POS = plot_pos
1287 1291
1288 1292 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1289 1293 self.isConfig = True
1290 1294 self.figfile = figfile
1291 1295 self.xdata = numpy.array([])
1292 1296 self.ydata = numpy.array([])
1293 1297
1294 1298 update_figfile = True
1295 1299
1296 1300 #open file beacon phase
1297 1301 path = '%s%03d' %(self.PREFIX, self.id)
1298 1302 beacon_file = os.path.join(path,'%s.txt'%self.name)
1299 1303 self.filename_phase = os.path.join(figpath,beacon_file)
1300 1304 #self.save_phase(self.filename_phase)
1301 1305
1302 1306
1303 1307 #store data beacon phase
1304 1308 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1305 1309
1306 1310 self.setWinTitle(title)
1307 1311
1308 1312
1309 1313 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1310 1314
1311 1315 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1312 1316
1313 1317 axes = self.axesList[0]
1314 1318
1315 1319 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1316 1320
1317 1321 if len(self.ydata)==0:
1318 1322 self.ydata = phase_beacon.reshape(-1,1)
1319 1323 else:
1320 1324 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1321 1325
1322 1326
1323 1327 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1324 1328 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1325 1329 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1326 1330 XAxisAsTime=True, grid='both'
1327 1331 )
1328 1332
1329 1333 self.draw()
1330 1334
1331 1335 if dataOut.ltctime >= self.xmax:
1332 1336 self.counter_imagwr = wr_period
1333 1337 self.isConfig = False
1334 1338 update_figfile = True
1335 1339
1336 1340 self.save(figpath=figpath,
1337 1341 figfile=figfile,
1338 1342 save=save,
1339 1343 ftp=ftp,
1340 1344 wr_period=wr_period,
1341 1345 thisDatetime=thisDatetime,
1342 1346 update_figfile=update_figfile)
1343 1347
1344 1348 return dataOut
Show file before | Show file after | Hide comments
@@ -1,1345 +1,1345
1 1
2 2 import os
3 3 import time
4 4 import math
5 5 import datetime
6 6 import numpy
7 7
8 8 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator #YONG
9 9
10 10 from .jroplot_spectra import RTIPlot, NoisePlot
11 11
12 12 from schainpy.utils import log
13 13 from .plotting_codes import *
14 14
15 15 from schainpy.model.graphics.jroplot_base import Plot, plt
16 16
17 17 import matplotlib.pyplot as plt
18 18 import matplotlib.colors as colors
19 19 from matplotlib.ticker import MultipleLocator, LogLocator, NullFormatter
20 20
21 21 class RTIDPPlot(RTIPlot):
22 22 '''
23 23 Written by R. Flores
24 24 '''
25 25 '''Plot for RTI Double Pulse Experiment Using Cross Products Analysis
26 26 '''
27 27
28 28 CODE = 'RTIDP'
29 29 colormap = 'jet'
30 30 plot_name = 'RTI'
31 31 plot_type = 'pcolorbuffer'
32 32
33 33 def setup(self):
34 34 self.xaxis = 'time'
35 35 self.ncols = 1
36 36 self.nrows = 3
37 37 self.nplots = self.nrows
38 38
39 39 self.ylabel = 'Range [km]'
40 40 self.xlabel = 'Time (LT)'
41 41
42 42 self.cb_label = 'Intensity (dB)'
43 43
44 44 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
45 45
46 46 self.titles = ['{} Channel {}'.format(
47 47 self.plot_name.upper(), '0x1'),'{} Channel {}'.format(
48 48 self.plot_name.upper(), '0'),'{} Channel {}'.format(
49 49 self.plot_name.upper(), '1')]
50 50
51 51 def update(self, dataOut):
52 52
53 53 data = {}
54 54 meta = {}
55 55 data['rti'] = dataOut.data_for_RTI_DP
56 56 data['NDP'] = dataOut.NDP
57 57
58 58 return data, meta
59 59
60 60 def plot(self):
61 61
62 62 NDP = self.data['NDP'][-1]
63 63 self.x = self.data.times
64 64 self.y = self.data.yrange[0:NDP]
65 65 self.z = self.data['rti']
66 66 self.z = numpy.ma.masked_invalid(self.z)
67 67
68 68 if self.decimation is None:
69 69 x, y, z = self.fill_gaps(self.x, self.y, self.z)
70 70 else:
71 71 x, y, z = self.fill_gaps(*self.decimate())
72 72
73 73 for n, ax in enumerate(self.axes):
74 74
75 75 self.zmax = self.zmax if self.zmax is not None else numpy.max(
76 76 self.z[1][0,12:40])
77 77 self.zmin = self.zmin if self.zmin is not None else numpy.min(
78 78 self.z[1][0,12:40])
79 79
80 80 if ax.firsttime:
81 81
82 82 if self.zlimits is not None:
83 83 self.zmin, self.zmax = self.zlimits[n]
84 84
85 85 ax.plt = ax.pcolormesh(x, y, z[n].T,
86 86 vmin=self.zmin,
87 87 vmax=self.zmax,
88 88 cmap=plt.get_cmap(self.colormap)
89 89 )
90 90 else:
91 91 #if self.zlimits is not None:
92 92 #self.zmin, self.zmax = self.zlimits[n]
93 93 ax.collections.remove(ax.collections[0])
94 94 ax.plt = ax.pcolormesh(x, y, z[n].T,
95 95 vmin=self.zmin,
96 96 vmax=self.zmax,
97 97 cmap=plt.get_cmap(self.colormap)
98 98 )
99 99
100 100
101 101 class RTILPPlot(RTIPlot):
102 102 '''
103 103 Written by R. Flores
104 104 '''
105 105 '''
106 106 Plot for RTI Long Pulse Using Cross Products Analysis
107 107 '''
108 108
109 109 CODE = 'RTILP'
110 110 colormap = 'jet'
111 111 plot_name = 'RTI LP'
112 112 plot_type = 'pcolorbuffer'
113 113
114 114 def setup(self):
115 115 self.xaxis = 'time'
116 116 self.ncols = 1
117 117 self.nrows = 2
118 118 self.nplots = self.nrows
119 119
120 120 self.ylabel = 'Range [km]'
121 121 self.xlabel = 'Time (LT)'
122 122
123 123 self.cb_label = 'Intensity (dB)'
124 124
125 125 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
126 126
127 127
128 128 self.titles = ['{} Channel {}'.format(
129 129 self.plot_name.upper(), '0'),'{} Channel {}'.format(
130 130 self.plot_name.upper(), '1'),'{} Channel {}'.format(
131 131 self.plot_name.upper(), '2'),'{} Channel {}'.format(
132 132 self.plot_name.upper(), '3')]
133 133
134 134
135 135 def update(self, dataOut):
136 136
137 137 data = {}
138 138 meta = {}
139 139 data['rti'] = dataOut.data_for_RTI_LP
140 140 data['NRANGE'] = dataOut.NRANGE
141 141
142 142 return data, meta
143 143
144 144 def plot(self):
145 145
146 146 NRANGE = self.data['NRANGE'][-1]
147 147 self.x = self.data.times
148 148 self.y = self.data.yrange[0:NRANGE]
149 149
150 150 self.z = self.data['rti']
151 151
152 152 self.z = numpy.ma.masked_invalid(self.z)
153 153
154 154 if self.decimation is None:
155 155 x, y, z = self.fill_gaps(self.x, self.y, self.z)
156 156 else:
157 157 x, y, z = self.fill_gaps(*self.decimate())
158 158
159 159 for n, ax in enumerate(self.axes):
160 160
161 161 self.zmax = self.zmax if self.zmax is not None else numpy.max(
162 162 self.z[1][0,12:40])
163 163 self.zmin = self.zmin if self.zmin is not None else numpy.min(
164 164 self.z[1][0,12:40])
165 165
166 166 if ax.firsttime:
167 167
168 168 if self.zlimits is not None:
169 169 self.zmin, self.zmax = self.zlimits[n]
170 170
171 171
172 172 ax.plt = ax.pcolormesh(x, y, z[n].T,
173 173 vmin=self.zmin,
174 174 vmax=self.zmax,
175 175 cmap=plt.get_cmap(self.colormap)
176 176 )
177 177
178 178 else:
179 179 if self.zlimits is not None:
180 180 self.zmin, self.zmax = self.zlimits[n]
181 181 ax.collections.remove(ax.collections[0])
182 182 ax.plt = ax.pcolormesh(x, y, z[n].T,
183 183 vmin=self.zmin,
184 184 vmax=self.zmax,
185 185 cmap=plt.get_cmap(self.colormap)
186 186 )
187 187
188 188
189 189 class DenRTIPlot(RTIPlot):
190 190 '''
191 191 Written by R. Flores
192 192 '''
193 193 '''
194 194 Plot for Den
195 195 '''
196 196
197 197 CODE = 'denrti'
198 198 colormap = 'jet'
199 199
200 200 def setup(self):
201 201 self.xaxis = 'time'
202 202 self.ncols = 1
203 203 self.nrows = self.data.shape(self.CODE)[0]
204 204 self.nplots = self.nrows
205 205
206 206 self.ylabel = 'Range [km]'
207 207 self.xlabel = 'Time (LT)'
208 208
209 209 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
210 210
211 211 if self.CODE == 'denrti':
212 212 self.cb_label = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
213 213
214 214 self.titles = ['Electron Density RTI']
215 215
216 216 def update(self, dataOut):
217 217
218 218 data = {}
219 219 meta = {}
220 220
221 221 data['denrti'] = dataOut.DensityFinal*1.e-6 #To Plot in cm^-3
222 222
223 223 return data, meta
224 224
225 225 def plot(self):
226 226
227 227 self.x = self.data.times
228 228 self.y = self.data.yrange
229 229
230 230 self.z = self.data[self.CODE]
231 231
232 232 self.z = numpy.ma.masked_invalid(self.z)
233 233
234 234 if self.decimation is None:
235 235 x, y, z = self.fill_gaps(self.x, self.y, self.z)
236 236 else:
237 237 x, y, z = self.fill_gaps(*self.decimate())
238 238
239 239 for n, ax in enumerate(self.axes):
240 240
241 241 self.zmax = self.zmax if self.zmax is not None else numpy.max(
242 242 self.z[n])
243 243 self.zmin = self.zmin if self.zmin is not None else numpy.min(
244 244 self.z[n])
245 245
246 246 if ax.firsttime:
247 247
248 248 if self.zlimits is not None:
249 249 self.zmin, self.zmax = self.zlimits[n]
250 250 if numpy.log10(self.zmin)<0:
251 251 self.zmin=1
252 252 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
253 253 vmin=self.zmin,
254 254 vmax=self.zmax,
255 255 cmap=self.cmaps[n],
256 256 norm=colors.LogNorm()
257 257 )
258 258
259 259 else:
260 260 if self.zlimits is not None:
261 261 self.zmin, self.zmax = self.zlimits[n]
262 262 ax.collections.remove(ax.collections[0])
263 263 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
264 264 vmin=self.zmin,
265 265 vmax=self.zmax,
266 266 cmap=self.cmaps[n],
267 267 norm=colors.LogNorm()
268 268 )
269 269
270 270
271 271 class ETempRTIPlot(RTIPlot):
272 272 '''
273 273 Written by R. Flores
274 274 '''
275 275 '''
276 276 Plot for Electron Temperature
277 277 '''
278 278
279 279 CODE = 'ETemp'
280 280 colormap = 'jet'
281 281
282 282 def setup(self):
283 283 self.xaxis = 'time'
284 284 self.ncols = 1
285 285 self.nrows = self.data.shape(self.CODE)[0]
286 286 self.nplots = self.nrows
287 287
288 288 self.ylabel = 'Range [km]'
289 289 self.xlabel = 'Time (LT)'
290 290 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
291 291 if self.CODE == 'ETemp':
292 292 self.cb_label = 'Electron Temperature (K)'
293 293 self.titles = ['Electron Temperature RTI']
294 294 if self.CODE == 'ITemp':
295 295 self.cb_label = 'Ion Temperature (K)'
296 296 self.titles = ['Ion Temperature RTI']
297 297 if self.CODE == 'HeFracLP':
298 298 self.cb_label ='He+ Fraction'
299 299 self.titles = ['He+ Fraction RTI']
300 300 self.zmax=0.16
301 301 if self.CODE == 'HFracLP':
302 302 self.cb_label ='H+ Fraction'
303 303 self.titles = ['H+ Fraction RTI']
304 304
305 305 def update(self, dataOut):
306 306
307 307 data = {}
308 308 meta = {}
309 309
310 310 data['ETemp'] = dataOut.ElecTempFinal
311 311
312 312 return data, meta
313 313
314 314 def plot(self):
315 315
316 316 self.x = self.data.times
317 317 self.y = self.data.yrange
318 318 self.z = self.data[self.CODE]
319 319
320 320 self.z = numpy.ma.masked_invalid(self.z)
321 321
322 322 if self.decimation is None:
323 323 x, y, z = self.fill_gaps(self.x, self.y, self.z)
324 324 else:
325 325 x, y, z = self.fill_gaps(*self.decimate())
326 326
327 327 for n, ax in enumerate(self.axes):
328 328
329 329 self.zmax = self.zmax if self.zmax is not None else numpy.max(
330 330 self.z[n])
331 331 self.zmin = self.zmin if self.zmin is not None else numpy.min(
332 332 self.z[n])
333 333
334 334 if ax.firsttime:
335 335
336 336 if self.zlimits is not None:
337 337 self.zmin, self.zmax = self.zlimits[n]
338 338
339 339 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
340 340 vmin=self.zmin,
341 341 vmax=self.zmax,
342 342 cmap=self.cmaps[n]
343 343 )
344 344 #plt.tight_layout()
345 345
346 346 else:
347 347 if self.zlimits is not None:
348 348 self.zmin, self.zmax = self.zlimits[n]
349 349 ax.collections.remove(ax.collections[0])
350 350 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
351 351 vmin=self.zmin,
352 352 vmax=self.zmax,
353 353 cmap=self.cmaps[n]
354 354 )
355 355
356 356
357 357 class ITempRTIPlot(ETempRTIPlot):
358 358 '''
359 359 Written by R. Flores
360 360 '''
361 361 '''
362 362 Plot for Ion Temperature
363 363 '''
364 364
365 365 CODE = 'ITemp'
366 366 colormap = 'jet'
367 367 plot_name = 'Ion Temperature'
368 368
369 369 def update(self, dataOut):
370 370
371 371 data = {}
372 372 meta = {}
373 373
374 374 data['ITemp'] = dataOut.IonTempFinal
375 375
376 376 return data, meta
377 377
378 378
379 379 class HFracRTIPlot(ETempRTIPlot):
380 380 '''
381 381 Written by R. Flores
382 382 '''
383 383 '''
384 384 Plot for H+ LP
385 385 '''
386 386
387 387 CODE = 'HFracLP'
388 388 colormap = 'jet'
389 389 plot_name = 'H+ Frac'
390 390
391 391 def update(self, dataOut):
392 392
393 393 data = {}
394 394 meta = {}
395 395 data['HFracLP'] = dataOut.PhyFinal
396 396
397 397 return data, meta
398 398
399 399
400 400 class HeFracRTIPlot(ETempRTIPlot):
401 401 '''
402 402 Written by R. Flores
403 403 '''
404 404 '''
405 405 Plot for He+ LP
406 406 '''
407 407
408 408 CODE = 'HeFracLP'
409 409 colormap = 'jet'
410 410 plot_name = 'He+ Frac'
411 411
412 412 def update(self, dataOut):
413 413
414 414 data = {}
415 415 meta = {}
416 416 data['HeFracLP'] = dataOut.PheFinal
417 417
418 418 return data, meta
419 419
420 420
421 421 class TempsDPPlot(Plot):
422 422 '''
423 423 Written by R. Flores
424 424 '''
425 425 '''
426 426 Plot for Electron - Ion Temperatures
427 427 '''
428 428
429 429 CODE = 'tempsDP'
430 430 #plot_name = 'Temperatures'
431 431 plot_type = 'scatterbuffer'
432 432
433 433 def setup(self):
434 434
435 435 self.ncols = 1
436 436 self.nrows = 1
437 437 self.nplots = 1
438 438 self.ylabel = 'Range [km]'
439 439 self.xlabel = 'Temperature (K)'
440 440 self.titles = ['Electron/Ion Temperatures']
441 441 self.width = 3.5
442 442 self.height = 5.5
443 443 self.colorbar = False
444 444 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
445 445
446 446 def update(self, dataOut):
447 447 data = {}
448 448 meta = {}
449 449
450 450 data['Te'] = dataOut.te2
451 451 data['Ti'] = dataOut.ti2
452 452 data['Te_error'] = dataOut.ete2
453 453 data['Ti_error'] = dataOut.eti2
454 454
455 455 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
456 456
457 457 return data, meta
458 458
459 459 def plot(self):
460 460
461 461 y = self.data.yrange
462 462
463 463 self.xmin = -100
464 464 self.xmax = 5000
465 465
466 466 ax = self.axes[0]
467 467
468 468 data = self.data[-1]
469 469
470 470 Te = data['Te']
471 471 Ti = data['Ti']
472 472 errTe = data['Te_error']
473 473 errTi = data['Ti_error']
474 474
475 475 if ax.firsttime:
476 476 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
477 477 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
478 478 plt.legend(loc='lower right')
479 479 self.ystep_given = 50
480 480 ax.yaxis.set_minor_locator(MultipleLocator(15))
481 481 ax.grid(which='minor')
482 482
483 483 else:
484 484 self.clear_figures()
485 485 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
486 486 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
487 487 plt.legend(loc='lower right')
488 488 ax.yaxis.set_minor_locator(MultipleLocator(15))
489 489
490 490
491 491 class TempsHPPlot(Plot):
492 492 '''
493 493 Written by R. Flores
494 494 '''
495 495 '''
496 496 Plot for Temperatures Hybrid Experiment
497 497 '''
498 498
499 499 CODE = 'temps_LP'
500 500 #plot_name = 'Temperatures'
501 501 plot_type = 'scatterbuffer'
502 502
503 503
504 504 def setup(self):
505 505
506 506 self.ncols = 1
507 507 self.nrows = 1
508 508 self.nplots = 1
509 509 self.ylabel = 'Range [km]'
510 510 self.xlabel = 'Temperature (K)'
511 511 self.titles = ['Electron/Ion Temperatures']
512 512 self.width = 3.5
513 513 self.height = 6.5
514 514 self.colorbar = False
515 515 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
516 516
517 517 def update(self, dataOut):
518 518 data = {}
519 519 meta = {}
520 520
521 521
522 522 data['Te'] = numpy.concatenate((dataOut.te2[:dataOut.cut],dataOut.te[dataOut.cut:]))
523 523 data['Ti'] = numpy.concatenate((dataOut.ti2[:dataOut.cut],dataOut.ti[dataOut.cut:]))
524 524 data['Te_error'] = numpy.concatenate((dataOut.ete2[:dataOut.cut],dataOut.ete[dataOut.cut:]))
525 525 data['Ti_error'] = numpy.concatenate((dataOut.eti2[:dataOut.cut],dataOut.eti[dataOut.cut:]))
526 526
527 527 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
528 528
529 529 return data, meta
530 530
531 531 def plot(self):
532 532
533 533
534 534 self.y = self.data.yrange
535 535 self.xmin = -100
536 536 self.xmax = 4500
537 537 ax = self.axes[0]
538 538
539 539 data = self.data[-1]
540 540
541 541 Te = data['Te']
542 542 Ti = data['Ti']
543 543 errTe = data['Te_error']
544 544 errTi = data['Ti_error']
545 545
546 546 if ax.firsttime:
547 547
548 548 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
549 549 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
550 550 plt.legend(loc='lower right')
551 551 self.ystep_given = 200
552 552 ax.yaxis.set_minor_locator(MultipleLocator(15))
553 553 ax.grid(which='minor')
554 554
555 555 else:
556 556 self.clear_figures()
557 557 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
558 558 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
559 559 plt.legend(loc='lower right')
560 560 ax.yaxis.set_minor_locator(MultipleLocator(15))
561 561 ax.grid(which='minor')
562 562
563 563
564 564 class FracsHPPlot(Plot):
565 565 '''
566 566 Written by R. Flores
567 567 '''
568 568 '''
569 569 Plot for Composition LP
570 570 '''
571 571
572 572 CODE = 'fracs_LP'
573 573 plot_type = 'scatterbuffer'
574 574
575 575
576 576 def setup(self):
577 577
578 578 self.ncols = 1
579 579 self.nrows = 1
580 580 self.nplots = 1
581 581 self.ylabel = 'Range [km]'
582 582 self.xlabel = 'Frac'
583 583 self.titles = ['Composition']
584 584 self.width = 3.5
585 585 self.height = 6.5
586 586 self.colorbar = False
587 587 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
588 588
589 589 def update(self, dataOut):
590 590 data = {}
591 591 meta = {}
592 592
593 593 #aux_nan=numpy.zeros(dataOut.cut,'float32')
594 594 #aux_nan[:]=numpy.nan
595 595 #data['ph'] = numpy.concatenate((aux_nan,dataOut.ph[dataOut.cut:]))
596 596 #data['eph'] = numpy.concatenate((aux_nan,dataOut.eph[dataOut.cut:]))
597 597
598 598 data['ph'] = dataOut.ph[dataOut.cut:]
599 599 data['eph'] = dataOut.eph[dataOut.cut:]
600 600 data['phe'] = dataOut.phe[dataOut.cut:]
601 601 data['ephe'] = dataOut.ephe[dataOut.cut:]
602 602
603 603 data['cut'] = dataOut.cut
604 604
605 605 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
606 606
607 607
608 608 return data, meta
609 609
610 610 def plot(self):
611 611
612 612 data = self.data[-1]
613 613
614 614 ph = data['ph']
615 615 eph = data['eph']
616 616 phe = data['phe']
617 617 ephe = data['ephe']
618 618 cut = data['cut']
619 619 self.y = self.data.yrange
620 620
621 621 self.xmin = 0
622 622 self.xmax = 1
623 623 ax = self.axes[0]
624 624
625 625 if ax.firsttime:
626 626
627 627 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
628 628 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
629 629 plt.legend(loc='lower right')
630 630 self.xstep_given = 0.2
631 631 self.ystep_given = 200
632 632 ax.yaxis.set_minor_locator(MultipleLocator(15))
633 633 ax.grid(which='minor')
634 634
635 635 else:
636 636 self.clear_figures()
637 637 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
638 638 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
639 639 plt.legend(loc='lower right')
640 640 ax.yaxis.set_minor_locator(MultipleLocator(15))
641 641 ax.grid(which='minor')
642 642
643 643 class EDensityPlot(Plot):
644 644 '''
645 645 Written by R. Flores
646 646 '''
647 647 '''
648 648 Plot for electron density
649 649 '''
650 650
651 651 CODE = 'den'
652 652 #plot_name = 'Electron Density'
653 653 plot_type = 'scatterbuffer'
654 654
655 655 def setup(self):
656 656
657 657 self.ncols = 1
658 658 self.nrows = 1
659 659 self.nplots = 1
660 660 self.ylabel = 'Range [km]'
661 661 self.xlabel = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
662 662 self.titles = ['Electron Density']
663 663 self.width = 3.5
664 664 self.height = 5.5
665 665 self.colorbar = False
666 666 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
667 667
668 668 def update(self, dataOut):
669 669 data = {}
670 670 meta = {}
671 671
672 672 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
673 673 data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS]
674 674 data['den_error'] = dataOut.sdp2[:dataOut.NSHTS]
675 675 #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS]
676 print(numpy.shape(data['den_power']))
677 print(numpy.shape(data['den_Faraday']))
678 print(numpy.shape(data['den_error']))
676 #print(numpy.shape(data['den_power']))
677 #print(numpy.shape(data['den_Faraday']))
678 #print(numpy.shape(data['den_error']))
679 679
680 680 data['NSHTS'] = dataOut.NSHTS
681 681
682 682 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
683 683
684 684 return data, meta
685 685
686 686 def plot(self):
687 687
688 688 y = self.data.yrange
689 689
690 690 #self.xmin = 1e3
691 691 #self.xmax = 1e7
692 692
693 693 ax = self.axes[0]
694 694
695 695 data = self.data[-1]
696 696
697 697 DenPow = data['den_power']
698 698 DenFar = data['den_Faraday']
699 699 errDenPow = data['den_error']
700 700 #errFaraday = data['err_Faraday']
701 701
702 702 NSHTS = data['NSHTS']
703 703
704 704 if self.CODE == 'denLP':
705 705 DenPowLP = data['den_LP']
706 706 errDenPowLP = data['den_LP_error']
707 707 cut = data['cut']
708 708
709 709 if ax.firsttime:
710 710 self.autoxticks=False
711 711 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
712 712 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2,linestyle='-')
713 713 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
714 714 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2,linestyle='-')
715 715
716 716 if self.CODE=='denLP':
717 717 ax.errorbar(DenPowLP[cut:], y[cut:], xerr=errDenPowLP[cut:], fmt='r^-',elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
718 718
719 719 plt.legend(loc='upper left',fontsize=8.5)
720 720 #plt.legend(loc='lower left',fontsize=8.5)
721 721 ax.set_xscale("log", nonposx='clip')
722 722 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
723 723 self.ystep_given=100
724 724 if self.CODE=='denLP':
725 725 self.ystep_given=200
726 726 ax.set_yticks(grid_y_ticks,minor=True)
727 727 locmaj = LogLocator(base=10,numticks=12)
728 728 ax.xaxis.set_major_locator(locmaj)
729 729 locmin = LogLocator(base=10.0,subs=(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12)
730 730 ax.xaxis.set_minor_locator(locmin)
731 731 ax.xaxis.set_minor_formatter(NullFormatter())
732 732 ax.grid(which='minor')
733 733
734 734 else:
735 735 dataBefore = self.data[-2]
736 736 DenPowBefore = dataBefore['den_power']
737 737 self.clear_figures()
738 738 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
739 739 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2,linestyle='-')
740 740 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
741 741 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2,linestyle='-')
742 742 ax.errorbar(DenPowBefore, y[:NSHTS], elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed")
743 743
744 744 if self.CODE=='denLP':
745 745 ax.errorbar(DenPowLP[cut:], y[cut:], fmt='r^-', xerr=errDenPowLP[cut:],elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
746 746
747 747 ax.set_xscale("log", nonposx='clip')
748 748 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
749 749 ax.set_yticks(grid_y_ticks,minor=True)
750 750 locmaj = LogLocator(base=10,numticks=12)
751 751 ax.xaxis.set_major_locator(locmaj)
752 752 locmin = LogLocator(base=10.0,subs=(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12)
753 753 ax.xaxis.set_minor_locator(locmin)
754 754 ax.xaxis.set_minor_formatter(NullFormatter())
755 755 ax.grid(which='minor')
756 756 plt.legend(loc='upper left',fontsize=8.5)
757 757 #plt.legend(loc='lower left',fontsize=8.5)
758 758
759 759 class FaradayAnglePlot(Plot):
760 760 '''
761 761 Written by R. Flores
762 762 '''
763 763 '''
764 764 Plot for electron density
765 765 '''
766 766
767 767 CODE = 'angle'
768 768 plot_name = 'Faraday Angle'
769 769 plot_type = 'scatterbuffer'
770 770
771 771 def setup(self):
772 772
773 773 self.ncols = 1
774 774 self.nrows = 1
775 775 self.nplots = 1
776 776 self.ylabel = 'Range [km]'
777 777 self.xlabel = 'Faraday Angle (º)'
778 778 self.titles = ['Electron Density']
779 779 self.width = 3.5
780 780 self.height = 5.5
781 781 self.colorbar = False
782 782 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
783 783
784 784 def update(self, dataOut):
785 785 data = {}
786 786 meta = {}
787 787
788 788 data['angle'] = numpy.degrees(dataOut.phi)
789 789 #'''
790 790 #print(dataOut.phi_uwrp)
791 791 #print(data['angle'])
792 792 #exit(1)
793 793 #'''
794 794 data['dphi'] = dataOut.dphi_uc*10
795 795 #print(dataOut.dphi)
796 796
797 797 #data['NSHTS'] = dataOut.NSHTS
798 798
799 799 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
800 800
801 801 return data, meta
802 802
803 803 def plot(self):
804 804
805 805 data = self.data[-1]
806 806 self.x = data[self.CODE]
807 807 dphi = data['dphi']
808 808 self.y = self.data.yrange
809 809 self.xmin = -360#-180
810 810 self.xmax = 360#180
811 811 ax = self.axes[0]
812 812
813 813 if ax.firsttime:
814 814 self.autoxticks=False
815 815 #if self.CODE=='den':
816 816 ax.plot(self.x, self.y,marker='o',color='g',linewidth=1.0,markersize=2)
817 817 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
818 818
819 819 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
820 820 self.ystep_given=100
821 821 if self.CODE=='denLP':
822 822 self.ystep_given=200
823 823 ax.set_yticks(grid_y_ticks,minor=True)
824 824 ax.grid(which='minor')
825 825 #plt.tight_layout()
826 826 else:
827 827
828 828 self.clear_figures()
829 829 #if self.CODE=='den':
830 830 #print(numpy.shape(self.x))
831 831 ax.plot(self.x, self.y, marker='o',color='g',linewidth=1.0, markersize=2)
832 832 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
833 833
834 834 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
835 835 ax.set_yticks(grid_y_ticks,minor=True)
836 836 ax.grid(which='minor')
837 837
838 838 class EDensityHPPlot(EDensityPlot):
839 839 '''
840 840 Written by R. Flores
841 841 '''
842 842 '''
843 843 Plot for Electron Density Hybrid Experiment
844 844 '''
845 845
846 846 CODE = 'denLP'
847 847 plot_name = 'Electron Density'
848 848 plot_type = 'scatterbuffer'
849 849
850 850 def update(self, dataOut):
851 851 data = {}
852 852 meta = {}
853 853
854 854 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
855 855 data['den_Faraday']=dataOut.dphi[:dataOut.NSHTS]
856 856 data['den_error']=dataOut.sdp2[:dataOut.NSHTS]
857 857 data['den_LP']=dataOut.ne[:dataOut.NACF]
858 858 data['den_LP_error']=dataOut.ene[:dataOut.NACF]*dataOut.ne[:dataOut.NACF]*0.434
859 859 #self.ene=10**dataOut.ene[:dataOut.NACF]
860 860 data['NSHTS']=dataOut.NSHTS
861 861 data['cut']=dataOut.cut
862 862
863 863 return data, meta
864 864
865 865
866 866 class ACFsPlot(Plot):
867 867 '''
868 868 Written by R. Flores
869 869 '''
870 870 '''
871 871 Plot for ACFs Double Pulse Experiment
872 872 '''
873 873
874 874 CODE = 'acfs'
875 875 #plot_name = 'ACF'
876 876 plot_type = 'scatterbuffer'
877 877
878 878
879 879 def setup(self):
880 880 self.ncols = 1
881 881 self.nrows = 1
882 882 self.nplots = 1
883 883 self.ylabel = 'Range [km]'
884 884 self.xlabel = 'Lag (ms)'
885 885 self.titles = ['ACFs']
886 886 self.width = 3.5
887 887 self.height = 5.5
888 888 self.colorbar = False
889 889 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
890 890
891 891 def update(self, dataOut):
892 892 data = {}
893 893 meta = {}
894 894
895 895 data['ACFs'] = dataOut.acfs_to_plot
896 896 data['ACFs_error'] = dataOut.acfs_error_to_plot
897 897 data['lags'] = dataOut.lags_to_plot
898 898 data['Lag_contaminated_1'] = dataOut.x_igcej_to_plot
899 899 data['Lag_contaminated_2'] = dataOut.x_ibad_to_plot
900 900 data['Height_contaminated_1'] = dataOut.y_igcej_to_plot
901 901 data['Height_contaminated_2'] = dataOut.y_ibad_to_plot
902 902
903 903 meta['yrange'] = numpy.array([])
904 904 #meta['NSHTS'] = dataOut.NSHTS
905 905 #meta['DPL'] = dataOut.DPL
906 906 data['NSHTS'] = dataOut.NSHTS #This is metadata
907 907 data['DPL'] = dataOut.DPL #This is metadata
908 908
909 909 return data, meta
910 910
911 911 def plot(self):
912 912
913 913 data = self.data[-1]
914 914 #NSHTS = self.meta['NSHTS']
915 915 #DPL = self.meta['DPL']
916 916 NSHTS = data['NSHTS'] #This is metadata
917 917 DPL = data['DPL'] #This is metadata
918 918
919 919 lags = data['lags']
920 920 ACFs = data['ACFs']
921 921 errACFs = data['ACFs_error']
922 922 BadLag1 = data['Lag_contaminated_1']
923 923 BadLag2 = data['Lag_contaminated_2']
924 924 BadHei1 = data['Height_contaminated_1']
925 925 BadHei2 = data['Height_contaminated_2']
926 926
927 927 self.xmin = 0.0
928 928 self.xmax = 2.0
929 929 self.y = ACFs
930 930
931 931 ax = self.axes[0]
932 932
933 933 if ax.firsttime:
934 934
935 935 for i in range(NSHTS):
936 936 x_aux = numpy.isfinite(lags[i,:])
937 937 y_aux = numpy.isfinite(ACFs[i,:])
938 938 yerr_aux = numpy.isfinite(errACFs[i,:])
939 939 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
940 940 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
941 941 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
942 942 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
943 943 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
944 944 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',marker='o',linewidth=1.0,markersize=2)
945 945 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
946 946 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
947 947
948 948 self.xstep_given = (self.xmax-self.xmin)/(DPL-1)
949 949 self.ystep_given = 50
950 950 ax.yaxis.set_minor_locator(MultipleLocator(15))
951 951 ax.grid(which='minor')
952 952
953 953 else:
954 954 self.clear_figures()
955 955 for i in range(NSHTS):
956 956 x_aux = numpy.isfinite(lags[i,:])
957 957 y_aux = numpy.isfinite(ACFs[i,:])
958 958 yerr_aux = numpy.isfinite(errACFs[i,:])
959 959 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
960 960 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
961 961 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
962 962 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
963 963 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
964 964 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],linewidth=1.0,markersize=2,color='b',marker='o')
965 965 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
966 966 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
967 967 ax.yaxis.set_minor_locator(MultipleLocator(15))
968 968
969 969 class ACFsLPPlot(Plot):
970 970 '''
971 971 Written by R. Flores
972 972 '''
973 973 '''
974 974 Plot for ACFs Double Pulse Experiment
975 975 '''
976 976
977 977 CODE = 'acfs_LP'
978 978 #plot_name = 'ACF'
979 979 plot_type = 'scatterbuffer'
980 980
981 981
982 982 def setup(self):
983 983 self.ncols = 1
984 984 self.nrows = 1
985 985 self.nplots = 1
986 986 self.ylabel = 'Range [km]'
987 987 self.xlabel = 'Lag (ms)'
988 988 self.titles = ['ACFs']
989 989 self.width = 3.5
990 990 self.height = 5.5
991 991 self.colorbar = False
992 992 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
993 993
994 994 def update(self, dataOut):
995 995 data = {}
996 996 meta = {}
997 997
998 998 aux=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
999 999 errors=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
1000 1000 lags_LP_to_plot=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
1001 1001
1002 1002 for i in range(dataOut.NACF):
1003 1003 for j in range(dataOut.IBITS):
1004 1004 if numpy.abs(dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0])<1.0:
1005 1005 aux[i,j]=dataOut.output_LP_integrated.real[j,i,0]/dataOut.output_LP_integrated.real[0,i,0]
1006 1006 aux[i,j]=max(min(aux[i,j],1.0),-1.0)*dataOut.DH+dataOut.heightList[i]
1007 1007 lags_LP_to_plot[i,j]=dataOut.lags_LP[j]
1008 1008 errors[i,j]=dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0]*dataOut.DH
1009 1009 else:
1010 1010 aux[i,j]=numpy.nan
1011 1011 lags_LP_to_plot[i,j]=numpy.nan
1012 1012 errors[i,j]=numpy.nan
1013 1013
1014 1014 data['ACFs'] = aux
1015 1015 data['ACFs_error'] = errors
1016 1016 data['lags'] = lags_LP_to_plot
1017 1017
1018 1018 meta['yrange'] = numpy.array([])
1019 1019 #meta['NACF'] = dataOut.NACF
1020 1020 #meta['NLAG'] = dataOut.NLAG
1021 1021 data['NACF'] = dataOut.NACF #This is metadata
1022 1022 data['NLAG'] = dataOut.NLAG #This is metadata
1023 1023
1024 1024 return data, meta
1025 1025
1026 1026 def plot(self):
1027 1027
1028 1028 data = self.data[-1]
1029 1029 #NACF = self.meta['NACF']
1030 1030 #NLAG = self.meta['NLAG']
1031 1031 NACF = data['NACF'] #This is metadata
1032 1032 NLAG = data['NLAG'] #This is metadata
1033 1033
1034 1034 lags = data['lags']
1035 1035 ACFs = data['ACFs']
1036 1036 errACFs = data['ACFs_error']
1037 1037
1038 1038 self.xmin = 0.0
1039 1039 self.xmax = 1.5
1040 1040
1041 1041 self.y = ACFs
1042 1042
1043 1043 ax = self.axes[0]
1044 1044
1045 1045 if ax.firsttime:
1046 1046
1047 1047 for i in range(NACF):
1048 1048 x_aux = numpy.isfinite(lags[i,:])
1049 1049 y_aux = numpy.isfinite(ACFs[i,:])
1050 1050 yerr_aux = numpy.isfinite(errACFs[i,:])
1051 1051
1052 1052 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1053 1053 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1054 1054
1055 1055 #self.xstep_given = (self.xmax-self.xmin)/(self.data.NLAG-1)
1056 1056 self.xstep_given=0.3
1057 1057 self.ystep_given = 200
1058 1058 ax.yaxis.set_minor_locator(MultipleLocator(15))
1059 1059 ax.grid(which='minor')
1060 1060
1061 1061 else:
1062 1062 self.clear_figures()
1063 1063
1064 1064 for i in range(NACF):
1065 1065 x_aux = numpy.isfinite(lags[i,:])
1066 1066 y_aux = numpy.isfinite(ACFs[i,:])
1067 1067 yerr_aux = numpy.isfinite(errACFs[i,:])
1068 1068
1069 1069 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1070 1070 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1071 1071
1072 1072 ax.yaxis.set_minor_locator(MultipleLocator(15))
1073 1073
1074 1074
1075 1075 class CrossProductsPlot(Plot):
1076 1076 '''
1077 1077 Written by R. Flores
1078 1078 '''
1079 1079 '''
1080 1080 Plot for cross products
1081 1081 '''
1082 1082
1083 1083 CODE = 'crossprod'
1084 1084 plot_name = 'Cross Products'
1085 1085 plot_type = 'scatterbuffer'
1086 1086
1087 1087 def setup(self):
1088 1088
1089 1089 self.ncols = 3
1090 1090 self.nrows = 1
1091 1091 self.nplots = 3
1092 1092 self.ylabel = 'Range [km]'
1093 1093 self.titles = []
1094 1094 self.width = 3.5*self.nplots
1095 1095 self.height = 5.5
1096 1096 self.colorbar = False
1097 1097 self.plots_adjust.update({'wspace':.3, 'left': 0.12, 'right': 0.92, 'bottom': 0.1})
1098 1098
1099 1099
1100 1100 def update(self, dataOut):
1101 1101
1102 1102 data = {}
1103 1103 meta = {}
1104 1104
1105 1105 data['crossprod'] = dataOut.crossprods
1106 1106 data['NDP'] = dataOut.NDP
1107 1107
1108 1108 return data, meta
1109 1109
1110 1110 def plot(self):
1111 1111
1112 1112 NDP = self.data['NDP'][-1]
1113 1113 x = self.data['crossprod'][:,-1,:,:,:,:]
1114 1114 y = self.data.yrange[0:NDP]
1115 1115
1116 1116 for n, ax in enumerate(self.axes):
1117 1117
1118 1118 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])))
1119 1119 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])))
1120 1120
1121 1121 if ax.firsttime:
1122 1122
1123 1123 self.autoxticks=False
1124 1124 if n==0:
1125 1125 label1='kax'
1126 1126 label2='kay'
1127 1127 label3='kbx'
1128 1128 label4='kby'
1129 1129 self.xlimits=[(self.xmin,self.xmax)]
1130 1130 elif n==1:
1131 1131 label1='kax2'
1132 1132 label2='kay2'
1133 1133 label3='kbx2'
1134 1134 label4='kby2'
1135 1135 self.xlimits.append((self.xmin,self.xmax))
1136 1136 elif n==2:
1137 1137 label1='kaxay'
1138 1138 label2='kbxby'
1139 1139 label3='kaxbx'
1140 1140 label4='kaxby'
1141 1141 self.xlimits.append((self.xmin,self.xmax))
1142 1142
1143 1143 ax.plotline1 = ax.plot(x[n][0,:,0,0], y, color='r',linewidth=2.0, label=label1)
1144 1144 ax.plotline2 = ax.plot(x[n][1,:,0,0], y, color='k',linewidth=2.0, label=label2)
1145 1145 ax.plotline3 = ax.plot(x[n][2,:,0,0], y, color='b',linewidth=2.0, label=label3)
1146 1146 ax.plotline4 = ax.plot(x[n][3,:,0,0], y, color='m',linewidth=2.0, label=label4)
1147 1147 ax.legend(loc='upper right')
1148 1148 ax.set_xlim(self.xmin, self.xmax)
1149 1149 self.titles.append('{}'.format(self.plot_name.upper()))
1150 1150
1151 1151 else:
1152 1152
1153 1153 if n==0:
1154 1154 self.xlimits=[(self.xmin,self.xmax)]
1155 1155 else:
1156 1156 self.xlimits.append((self.xmin,self.xmax))
1157 1157
1158 1158 ax.set_xlim(self.xmin, self.xmax)
1159 1159
1160 1160 ax.plotline1[0].set_data(x[n][0,:,0,0],y)
1161 1161 ax.plotline2[0].set_data(x[n][1,:,0,0],y)
1162 1162 ax.plotline3[0].set_data(x[n][2,:,0,0],y)
1163 1163 ax.plotline4[0].set_data(x[n][3,:,0,0],y)
1164 1164 self.titles.append('{}'.format(self.plot_name.upper()))
1165 1165
1166 1166
1167 1167 class CrossProductsLPPlot(Plot):
1168 1168 '''
1169 1169 Written by R. Flores
1170 1170 '''
1171 1171 '''
1172 1172 Plot for cross products LP
1173 1173 '''
1174 1174
1175 1175 CODE = 'crossprodslp'
1176 1176 plot_name = 'Cross Products LP'
1177 1177 plot_type = 'scatterbuffer'
1178 1178
1179 1179
1180 1180 def setup(self):
1181 1181
1182 1182 self.ncols = 2
1183 1183 self.nrows = 1
1184 1184 self.nplots = 2
1185 1185 self.ylabel = 'Range [km]'
1186 1186 self.xlabel = 'dB'
1187 1187 self.width = 3.5*self.nplots
1188 1188 self.height = 5.5
1189 1189 self.colorbar = False
1190 1190 self.titles = []
1191 1191 self.plots_adjust.update({'wspace': .8 ,'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1192 1192
1193 1193 def update(self, dataOut):
1194 1194 data = {}
1195 1195 meta = {}
1196 1196
1197 1197 data['crossprodslp'] = 10*numpy.log10(numpy.abs(dataOut.output_LP))
1198 1198
1199 1199 data['NRANGE'] = dataOut.NRANGE #This is metadata
1200 1200 data['NLAG'] = dataOut.NLAG #This is metadata
1201 1201
1202 1202 return data, meta
1203 1203
1204 1204 def plot(self):
1205 1205
1206 1206 NRANGE = self.data['NRANGE'][-1]
1207 1207 NLAG = self.data['NLAG'][-1]
1208 1208
1209 1209 x = self.data[self.CODE][:,-1,:,:]
1210 1210 self.y = self.data.yrange[0:NRANGE]
1211 1211
1212 1212 label_array=numpy.array(['lag '+ str(x) for x in range(NLAG)])
1213 1213 color_array=['r','k','g','b','c','m','y','orange','steelblue','purple','peru','darksalmon','grey','limegreen','olive','midnightblue']
1214 1214
1215 1215
1216 1216 for n, ax in enumerate(self.axes):
1217 1217
1218 1218 self.xmin=28#30
1219 1219 self.xmax=70#70
1220 1220 #self.xmin=numpy.min(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1221 1221 #self.xmax=numpy.max(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1222 1222
1223 1223 if ax.firsttime:
1224 1224
1225 1225 self.autoxticks=False
1226 1226 if n == 0:
1227 1227 self.plotline_array=numpy.zeros((2,NLAG),dtype=object)
1228 1228
1229 1229 for i in range(NLAG):
1230 1230 self.plotline_array[n,i], = ax.plot(x[i,:,n], self.y, color=color_array[i],linewidth=1.0, label=label_array[i])
1231 1231
1232 1232 ax.legend(loc='upper right')
1233 1233 ax.set_xlim(self.xmin, self.xmax)
1234 1234 if n==0:
1235 1235 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1236 1236 if n==1:
1237 1237 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1238 1238 else:
1239 1239 for i in range(NLAG):
1240 1240 self.plotline_array[n,i].set_data(x[i,:,n],self.y)
1241 1241
1242 1242 if n==0:
1243 1243 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1244 1244 if n==1:
1245 1245 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1246 1246
1247 1247
1248 1248 class NoiseDPPlot(NoisePlot):
1249 1249 '''
1250 1250 Written by R. Flores
1251 1251 '''
1252 1252 '''
1253 1253 Plot for noise Double Pulse
1254 1254 '''
1255 1255
1256 1256 CODE = 'noise'
1257 1257 #plot_name = 'Noise'
1258 1258 #plot_type = 'scatterbuffer'
1259 1259
1260 1260 def update(self, dataOut):
1261 1261
1262 1262 data = {}
1263 1263 meta = {}
1264 1264 data['noise'] = 10*numpy.log10(dataOut.noise_final)
1265 1265
1266 1266 return data, meta
1267 1267
1268 1268
1269 1269 class XmitWaveformPlot(Plot):
1270 1270 '''
1271 1271 Written by R. Flores
1272 1272 '''
1273 1273 '''
1274 1274 Plot for xmit waveform
1275 1275 '''
1276 1276
1277 1277 CODE = 'xmit'
1278 1278 plot_name = 'Xmit Waveform'
1279 1279 plot_type = 'scatterbuffer'
1280 1280
1281 1281
1282 1282 def setup(self):
1283 1283
1284 1284 self.ncols = 1
1285 1285 self.nrows = 1
1286 1286 self.nplots = 1
1287 1287 self.ylabel = ''
1288 1288 self.xlabel = 'Number of Lag'
1289 1289 self.width = 5.5
1290 1290 self.height = 3.5
1291 1291 self.colorbar = False
1292 1292 self.plots_adjust.update({'right': 0.85 })
1293 1293 self.titles = [self.plot_name]
1294 1294 #self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1295 1295
1296 1296 #if not self.titles:
1297 1297 #self.titles = self.data.parameters \
1298 1298 #if self.data.parameters else ['{}'.format(self.plot_name.upper())]
1299 1299
1300 1300 def update(self, dataOut):
1301 1301
1302 1302 data = {}
1303 1303 meta = {}
1304 1304
1305 1305 y_1=numpy.arctan2(dataOut.output_LP[:,0,2].imag,dataOut.output_LP[:,0,2].real)* 180 / (numpy.pi*10)
1306 1306 y_2=numpy.abs(dataOut.output_LP[:,0,2])
1307 1307 norm=numpy.max(y_2)
1308 1308 norm=max(norm,0.1)
1309 1309 y_2=y_2/norm
1310 1310
1311 1311 meta['yrange'] = numpy.array([])
1312 1312
1313 1313 data['xmit'] = numpy.vstack((y_1,y_2))
1314 1314 data['NLAG'] = dataOut.NLAG
1315 1315
1316 1316 return data, meta
1317 1317
1318 1318 def plot(self):
1319 1319
1320 1320 data = self.data[-1]
1321 1321 NLAG = data['NLAG']
1322 1322 x = numpy.arange(0,NLAG,1,'float32')
1323 1323 y = data['xmit']
1324 1324
1325 1325 self.xmin = 0
1326 1326 self.xmax = NLAG-1
1327 1327 self.ymin = -1.0
1328 1328 self.ymax = 1.0
1329 1329 ax = self.axes[0]
1330 1330
1331 1331 if ax.firsttime:
1332 1332 ax.plotline0=ax.plot(x,y[0,:],color='blue')
1333 1333 ax.plotline1=ax.plot(x,y[1,:],color='red')
1334 1334 secax=ax.secondary_xaxis(location=0.5)
1335 1335 secax.xaxis.tick_bottom()
1336 1336 secax.tick_params( labelleft=False, labeltop=False,
1337 1337 labelright=False, labelbottom=False)
1338 1338
1339 1339 self.xstep_given = 3
1340 1340 self.ystep_given = .25
1341 1341 secax.set_xticks(numpy.linspace(self.xmin, self.xmax, 6)) #only works on matplotlib.version>3.2
1342 1342
1343 1343 else:
1344 1344 ax.plotline0[0].set_data(x,y[0,:])
1345 1345 ax.plotline1[0].set_data(x,y[1,:])
Show file before | Show file after | Hide comments
@@ -1,648 +1,649
1 1 '''
2 2 Created on Aug 1, 2017
3 3
4 4 @author: Juan C. Espinoza
5 5 '''
6 6
7 7 import os
8 8 import sys
9 9 import time
10 10 import json
11 11 import glob
12 12 import datetime
13 13
14 14 import numpy
15 15 import h5py
16 16
17 17 import schainpy.admin
18 18 from schainpy.model.io.jroIO_base import LOCALTIME, Reader
19 19 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
20 20 from schainpy.model.data.jrodata import Parameters
21 21 from schainpy.utils import log
22 22
23 23 try:
24 24 import madrigal.cedar
25 25 except:
26 26 pass
27 27
28 28 try:
29 29 basestring
30 30 except:
31 31 basestring = str
32 32
33 33 DEF_CATALOG = {
34 34 'principleInvestigator': 'Marco Milla',
35 35 'expPurpose': '',
36 36 'cycleTime': '',
37 37 'correlativeExp': '',
38 38 'sciRemarks': '',
39 39 'instRemarks': ''
40 40 }
41 41
42 42 DEF_HEADER = {
43 43 'kindatDesc': '',
44 44 'analyst': 'Jicamarca User',
45 45 'comments': '',
46 46 'history': ''
47 47 }
48 48
49 49 MNEMONICS = {
50 50 10: 'jro',
51 51 11: 'jbr',
52 14: 'jmp', #Added by R. Flores
52 53 840: 'jul',
53 54 13: 'jas',
54 55 1000: 'pbr',
55 56 1001: 'hbr',
56 57 1002: 'obr',
57 58 400: 'clr'
58 59
59 60 }
60 61
61 62 UT1970 = datetime.datetime(1970, 1, 1) - datetime.timedelta(seconds=time.timezone)
62 63
63 64 def load_json(obj):
64 65 '''
65 66 Parse json as string instead of unicode
66 67 '''
67 68
68 69 if isinstance(obj, str):
69 70 iterable = json.loads(obj)
70 71 else:
71 72 iterable = obj
72 73
73 74 if isinstance(iterable, dict):
74 75 return {str(k): load_json(v) if isinstance(v, dict) else str(v) if isinstance(v, basestring) else v
75 76 for k, v in list(iterable.items())}
76 77 elif isinstance(iterable, (list, tuple)):
77 78 return [str(v) if isinstance(v, basestring) else v for v in iterable]
78 79
79 80 return iterable
80 81
81 82
82 83 class MADReader(Reader, ProcessingUnit):
83 84
84 85 def __init__(self):
85 86
86 87 ProcessingUnit.__init__(self)
87 88
88 89 self.dataOut = Parameters()
89 90 self.counter_records = 0
90 91 self.nrecords = None
91 92 self.flagNoMoreFiles = 0
92 93 self.filename = None
93 94 self.intervals = set()
94 95 self.datatime = datetime.datetime(1900,1,1)
95 96 self.format = None
96 97 self.filefmt = "***%Y%m%d*******"
97 98
98 99 def setup(self, **kwargs):
99 100
100 101 self.set_kwargs(**kwargs)
101 102 self.oneDDict = load_json(self.oneDDict)
102 103 self.twoDDict = load_json(self.twoDDict)
103 104 self.ind2DList = load_json(self.ind2DList)
104 105 self.independentParam = self.ind2DList[0]
105 106
106 107 if self.path is None:
107 108 raise ValueError('The path is not valid')
108 109
109 110 self.open_file = open
110 111 self.open_mode = 'rb'
111 112
112 113 if self.format is None:
113 114 raise ValueError('The format is not valid choose simple or hdf5')
114 115 elif self.format.lower() in ('simple', 'txt'):
115 116 self.ext = '.txt'
116 117 elif self.format.lower() in ('cedar',):
117 118 self.ext = '.001'
118 119 else:
119 120 self.ext = '.hdf5'
120 121 self.open_file = h5py.File
121 122 self.open_mode = 'r'
122 123
123 124 if self.online:
124 125 log.log("Searching files in online mode...", self.name)
125 126
126 127 for nTries in range(self.nTries):
127 128 fullpath = self.searchFilesOnLine(self.path, self.startDate,
128 129 self.endDate, self.expLabel, self.ext, self.walk,
129 130 self.filefmt, self.folderfmt)
130 131
131 132 try:
132 133 fullpath = next(fullpath)
133 134 except:
134 135 fullpath = None
135 136
136 137 if fullpath:
137 138 break
138 139
139 140 log.warning(
140 141 'Waiting {} sec for a valid file in {}: try {} ...'.format(
141 142 self.delay, self.path, nTries + 1),
142 143 self.name)
143 144 time.sleep(self.delay)
144 145
145 146 if not(fullpath):
146 147 raise schainpy.admin.SchainError(
147 148 'There isn\'t any valid file in {}'.format(self.path))
148 149
149 150 else:
150 151 log.log("Searching files in {}".format(self.path), self.name)
151 152 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
152 153 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
153 154
154 155 self.setNextFile()
155 156
156 157 def readFirstHeader(self):
157 158 '''Read header and data'''
158 159
159 160 self.parseHeader()
160 161 self.parseData()
161 162 self.blockIndex = 0
162 163
163 164 return
164 165
165 166 def parseHeader(self):
166 167 '''
167 168 '''
168 169
169 170 self.output = {}
170 171 self.version = '2'
171 172 s_parameters = None
172 173 if self.ext == '.txt':
173 174 self.parameters = [s.strip().lower() for s in self.fp.readline().decode().strip().split(' ') if s]
174 175 elif self.ext == '.hdf5':
175 176 self.metadata = self.fp['Metadata']
176 177 if '_record_layout' in self.metadata:
177 178 s_parameters = [s[0].lower().decode() for s in self.metadata['Independent Spatial Parameters']]
178 179 self.version = '3'
179 180 self.parameters = [s[0].lower().decode() for s in self.metadata['Data Parameters']]
180 181
181 182 log.success('Parameters found: {}'.format(self.parameters),
182 183 'MADReader')
183 184 if s_parameters:
184 185 log.success('Spatial parameters found: {}'.format(s_parameters),
185 186 'MADReader')
186 187
187 188 for param in list(self.oneDDict.keys()):
188 189 if param.lower() not in self.parameters:
189 190 log.warning(
190 191 'Parameter {} not found will be ignored'.format(
191 192 param),
192 193 'MADReader')
193 194 self.oneDDict.pop(param, None)
194 195
195 196 for param, value in list(self.twoDDict.items()):
196 197 if param.lower() not in self.parameters:
197 198 log.warning(
198 199 'Parameter {} not found, it will be ignored'.format(
199 200 param),
200 201 'MADReader')
201 202 self.twoDDict.pop(param, None)
202 203 continue
203 204 if isinstance(value, list):
204 205 if value[0] not in self.output:
205 206 self.output[value[0]] = []
206 207 self.output[value[0]].append([])
207 208
208 209 def parseData(self):
209 210 '''
210 211 '''
211 212
212 213 if self.ext == '.txt':
213 214 self.data = numpy.genfromtxt(self.fp, missing_values=('missing'))
214 215 self.nrecords = self.data.shape[0]
215 216 self.ranges = numpy.unique(self.data[:,self.parameters.index(self.independentParam.lower())])
216 217 self.counter_records = 0
217 218 elif self.ext == '.hdf5':
218 219 self.data = self.fp['Data']
219 220 self.ranges = numpy.unique(self.data['Table Layout'][self.independentParam.lower()])
220 221 self.times = numpy.unique(self.data['Table Layout']['ut1_unix'])
221 222 self.counter_records = int(self.data['Table Layout']['recno'][0])
222 223 self.nrecords = int(self.data['Table Layout']['recno'][-1])
223 224
224 225 def readNextBlock(self):
225 226
226 227 while True:
227 228 self.flagDiscontinuousBlock = 0
228 229 if self.counter_records == self.nrecords:
229 230 self.setNextFile()
230 231
231 232 self.readBlock()
232 233
233 234 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
234 235 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
235 236 log.warning(
236 237 'Reading Record No. {}/{} -> {} [Skipping]'.format(
237 238 self.counter_records,
238 239 self.nrecords,
239 240 self.datatime.ctime()),
240 241 'MADReader')
241 242 continue
242 243 break
243 244
244 245 log.log(
245 246 'Reading Record No. {}/{} -> {}'.format(
246 247 self.counter_records,
247 248 self.nrecords,
248 249 self.datatime.ctime()),
249 250 'MADReader')
250 251
251 252 return 1
252 253
253 254 def readBlock(self):
254 255 '''
255 256 '''
256 257 dum = []
257 258 if self.ext == '.txt':
258 259 dt = self.data[self.counter_records][:6].astype(int)
259 260 if datetime.datetime(dt[0], dt[1], dt[2], dt[3], dt[4], dt[5]).date() > self.datatime.date():
260 261 self.flagDiscontinuousBlock = 1
261 262 self.datatime = datetime.datetime(dt[0], dt[1], dt[2], dt[3], dt[4], dt[5])
262 263 while True:
263 264 dt = self.data[self.counter_records][:6].astype(int)
264 265 datatime = datetime.datetime(dt[0], dt[1], dt[2], dt[3], dt[4], dt[5])
265 266 if datatime == self.datatime:
266 267 dum.append(self.data[self.counter_records])
267 268 self.counter_records += 1
268 269 if self.counter_records == self.nrecords:
269 270 break
270 271 continue
271 272 self.intervals.add((datatime-self.datatime).seconds)
272 273 break
273 274 elif self.ext == '.hdf5':
274 275 datatime = datetime.datetime.utcfromtimestamp(
275 276 self.times[self.counter_records])
276 277 dum = self.data['Table Layout'][self.data['Table Layout']['recno']==self.counter_records]
277 278 self.intervals.add((datatime-self.datatime).seconds)
278 279 if datatime.date()>self.datatime.date():
279 280 self.flagDiscontinuousBlock = 1
280 281 self.datatime = datatime
281 282 self.counter_records += 1
282 283
283 284 self.buffer = numpy.array(dum)
284 285 return
285 286
286 287 def set_output(self):
287 288 '''
288 289 Storing data from buffer to dataOut object
289 290 '''
290 291
291 292 parameters = [None for __ in self.parameters]
292 293
293 294 for param, attr in list(self.oneDDict.items()):
294 295 x = self.parameters.index(param.lower())
295 296 setattr(self.dataOut, attr, self.buffer[0][x])
296 297
297 298 for param, value in list(self.twoDDict.items()):
298 299 dummy = numpy.zeros(self.ranges.shape) + numpy.nan
299 300 if self.ext == '.txt':
300 301 x = self.parameters.index(param.lower())
301 302 y = self.parameters.index(self.independentParam.lower())
302 303 ranges = self.buffer[:,y]
303 304 #if self.ranges.size == ranges.size:
304 305 # continue
305 306 index = numpy.where(numpy.in1d(self.ranges, ranges))[0]
306 307 dummy[index] = self.buffer[:,x]
307 308 else:
308 309 ranges = self.buffer[self.independentParam.lower()]
309 310 index = numpy.where(numpy.in1d(self.ranges, ranges))[0]
310 311 dummy[index] = self.buffer[param.lower()]
311 312
312 313 if isinstance(value, str):
313 314 if value not in self.independentParam:
314 315 setattr(self.dataOut, value, dummy.reshape(1,-1))
315 316 elif isinstance(value, list):
316 317 self.output[value[0]][value[1]] = dummy
317 318 parameters[value[1]] = param
318 319 for key, value in list(self.output.items()):
319 320 setattr(self.dataOut, key, numpy.array(value))
320 321
321 322 self.dataOut.parameters = [s for s in parameters if s]
322 323 self.dataOut.heightList = self.ranges
323 324 self.dataOut.utctime = (self.datatime - datetime.datetime(1970, 1, 1)).total_seconds()
324 325 self.dataOut.utctimeInit = self.dataOut.utctime
325 326 self.dataOut.paramInterval = min(self.intervals)
326 327 self.dataOut.useLocalTime = False
327 328 self.dataOut.flagNoData = False
328 329 self.dataOut.nrecords = self.nrecords
329 330 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
330 331
331 332 def getData(self):
332 333 '''
333 334 Storing data from databuffer to dataOut object
334 335 '''
335 336
336 337 if not self.readNextBlock():
337 338 self.dataOut.flagNoData = True
338 339 return 0
339 340
340 341 self.set_output()
341 342
342 343 return 1
343 344
344 345 def run(self, **kwargs):
345 346
346 347 if not(self.isConfig):
347 348 self.setup(**kwargs)
348 349 self.isConfig = True
349 350
350 351 self.getData()
351 352
352 353 return
353 354
354 355 @MPDecorator
355 356 class MADWriter(Operation):
356 357 '''Writing module for Madrigal files
357 358
358 359 type: external
359 360
360 361 Inputs:
361 362 path path where files will be created
362 363 oneDDict json of one-dimensional parameters in record where keys
363 364 are Madrigal codes (integers or mnemonics) and values the corresponding
364 365 dataOut attribute e.g: {
365 366 'gdlatr': 'lat',
366 367 'gdlonr': 'lon',
367 368 'gdlat2':'lat',
368 369 'glon2':'lon'}
369 370 ind2DList list of independent spatial two-dimensional parameters e.g:
370 371 ['heigthList']
371 372 twoDDict json of two-dimensional parameters in record where keys
372 373 are Madrigal codes (integers or mnemonics) and values the corresponding
373 374 dataOut attribute if multidimensional array specify as tupple
374 375 ('attr', pos) e.g: {
375 376 'gdalt': 'heightList',
376 377 'vn1p2': ('data_output', 0),
377 378 'vn2p2': ('data_output', 1),
378 379 'vn3': ('data_output', 2),
379 380 'snl': ('data_SNR', 'db')
380 381 }
381 382 metadata json of madrigal metadata (kinst, kindat, catalog and header)
382 383 format hdf5, cedar
383 384 blocks number of blocks per file'''
384 385
385 386 __attrs__ = ['path', 'oneDDict', 'ind2DList', 'twoDDict','metadata', 'format', 'blocks']
386 387 missing = -32767
387 388 currentDay = None
388 389
389 390 def __init__(self):
390 391
391 392 Operation.__init__(self)
392 393 self.dataOut = Parameters()
393 394 self.counter = 0
394 395 self.path = None
395 396 self.fp = None
396 397
397 398 def run(self, dataOut, path, oneDDict, ind2DList='[]', twoDDict='{}',
398 399 metadata='{}', format='cedar', **kwargs):
399 400
400 401
401 402 #if dataOut.AUX==1: #Modified
402 403
403 404 if not self.isConfig:
404 405 self.setup(path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs)
405 406 self.isConfig = True
406 407
407 408 self.dataOut = dataOut
408 409 self.putData()
409 410
410 411 return 1
411 412
412 413 def setup(self, path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs):
413 414 '''
414 415 Configure Operation
415 416 '''
416 417
417 418 self.path = path
418 419 self.blocks = kwargs.get('blocks', None)
419 420 self.counter = 0
420 421 self.oneDDict = load_json(oneDDict)
421 422 self.twoDDict = load_json(twoDDict)
422 423 self.ind2DList = load_json(ind2DList)
423 424 meta = load_json(metadata)
424 425 self.kinst = meta.get('kinst')
425 426 self.kindat = meta.get('kindat')
426 427 self.catalog = meta.get('catalog', DEF_CATALOG)
427 428 self.header = meta.get('header', DEF_HEADER)
428 429 if format == 'cedar':
429 430 self.ext = '.dat'
430 431 self.extra_args = {}
431 432 elif format == 'hdf5':
432 433 self.ext = '.hdf5'
433 434 self.extra_args = {'ind2DList': self.ind2DList}
434 435
435 436 self.keys = [k.lower() for k in self.twoDDict]
436 437 if 'range' in self.keys:
437 438 self.keys.remove('range')
438 439 if 'gdalt' in self.keys:
439 440 self.keys.remove('gdalt')
440 441
441 442 def setFile(self):
442 443 '''
443 444 Create new cedar file object
444 445 '''
445 446
446 447 self.mnemonic = MNEMONICS[self.kinst] #TODO get mnemonic from madrigal
447 448 date = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
448 449 #if self.dataOut.input_dat_type:
449 450 #date=datetime.datetime.fromtimestamp(self.dataOut.TimeBlockSeconds_for_dp_power)
450 451 #print("date",date)
451 452
452 453 filename = '{}{}{}'.format(self.mnemonic,
453 454 date.strftime('%Y%m%d_%H%M%S'),
454 455 self.ext)
455 456
456 457 self.fullname = os.path.join(self.path, filename)
457 458
458 459 if os.path.isfile(self.fullname) :
459 460 log.warning(
460 461 'Destination file {} already exists, previous file deleted.'.format(
461 462 self.fullname),
462 463 'MADWriter')
463 464 os.remove(self.fullname)
464 465
465 466 try:
466 467 log.success(
467 468 'Creating file: {}'.format(self.fullname),
468 469 'MADWriter')
469 470 if not os.path.exists(self.path):
470 471 os.makedirs(self.path)
471 472 self.fp = madrigal.cedar.MadrigalCedarFile(self.fullname, True)
472 473
473 474
474 475 except ValueError as e:
475 476 log.error(
476 477 'Impossible to create a cedar object with "madrigal.cedar.MadrigalCedarFile"',
477 478 'MADWriter')
478 479 return
479 480
480 481 return 1
481 482
482 483 def writeBlock(self):
483 484 '''
484 485 Add data records to cedar file taking data from oneDDict and twoDDict
485 486 attributes.
486 487 Allowed parameters in: parcodes.tab
487 488 '''
488 489 #self.dataOut.paramInterval=2
489 490 startTime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
490 491
491 492 endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
492 493
493 494 #if self.dataOut.input_dat_type:
494 495 #if self.dataOut.experiment=="DP":
495 496 #startTime=datetime.datetime.fromtimestamp(self.dataOut.TimeBlockSeconds_for_dp_power)
496 497 #endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
497 498
498 499
499 500 #print("2: ",startTime)
500 501 #print(endTime)
501 502 heights = self.dataOut.heightList
502 503 #print(heights)
503 504 #exit(1)
504 505 #print(self.blocks)
505 506 #print(startTime)
506 507 #print(endTime)
507 508 #print(heights)
508 509 #input()
509 510 if self.ext == '.dat':
510 511 for key, value in list(self.twoDDict.items()):
511 512 if isinstance(value, str):
512 513 data = getattr(self.dataOut, value)
513 514 invalid = numpy.isnan(data)
514 515 data[invalid] = self.missing
515 516 elif isinstance(value, (tuple, list)):
516 517 attr, key = value
517 518 data = getattr(self.dataOut, attr)
518 519 invalid = numpy.isnan(data)
519 520 data[invalid] = self.missing
520 521
521 522 out = {}
522 523 for key, value in list(self.twoDDict.items()):
523 524 key = key.lower()
524 525 if isinstance(value, str):
525 526 if 'db' in value.lower():
526 527 tmp = getattr(self.dataOut, value.replace('_db', ''))
527 528 SNRavg = numpy.average(tmp, axis=0)
528 529 tmp = 10*numpy.log10(SNRavg)
529 530 else:
530 531 tmp = getattr(self.dataOut, value)
531 532 out[key] = tmp.flatten()[:len(heights)]
532 533 elif isinstance(value, (tuple, list)):
533 534 attr, x = value
534 535 data = getattr(self.dataOut, attr)
535 536 #print(x)
536 537 #print(len(heights))
537 538 #print(data[int(x)][:len(heights)])
538 539 #print(numpy.shape(out))
539 540 #print(numpy.shape(data))
540 541
541 542 out[key] = data[int(x)][:len(heights)]
542 543
543 544 a = numpy.array([out[k] for k in self.keys])
544 545 #print(a)
545 546 nrows = numpy.array([numpy.isnan(a[:, x]).all() for x in range(len(heights))])
546 547 index = numpy.where(nrows == False)[0]
547 548
548 549 #print(startTime.minute)
549 550 rec = madrigal.cedar.MadrigalDataRecord(
550 551 self.kinst,
551 552 self.kindat,
552 553 startTime.year,
553 554 startTime.month,
554 555 startTime.day,
555 556 startTime.hour,
556 557 startTime.minute,
557 558 startTime.second,
558 559 startTime.microsecond/10000,
559 560 endTime.year,
560 561 endTime.month,
561 562 endTime.day,
562 563 endTime.hour,
563 564 endTime.minute,
564 565 endTime.second,
565 566 endTime.microsecond/10000,
566 567 list(self.oneDDict.keys()),
567 568 list(self.twoDDict.keys()),
568 569 len(index),
569 570 **self.extra_args
570 571 )
571 572 #print("rec",rec)
572 573 # Setting 1d values
573 574 for key in self.oneDDict:
574 575 rec.set1D(key, getattr(self.dataOut, self.oneDDict[key]))
575 576
576 577 # Setting 2d values
577 578 nrec = 0
578 579 for n in index:
579 580 for key in out:
580 581 rec.set2D(key, nrec, out[key][n])
581 582 nrec += 1
582 583
583 584 self.fp.append(rec)
584 585 if self.ext == '.hdf5' and self.counter %2 == 0 and self.counter > 0:
585 586 #print("here")
586 587 self.fp.dump()
587 588 if self.counter % 20 == 0 and self.counter > 0:
588 589 #self.fp.write()
589 590 log.log(
590 591 'Writing {} records'.format(
591 592 self.counter),
592 593 'MADWriter')
593 594
594 595 def setHeader(self):
595 596 '''
596 597 Create an add catalog and header to cedar file
597 598 '''
598 599
599 600 log.success('Closing file {}'.format(self.fullname), 'MADWriter')
600 601
601 602 if self.ext == '.dat':
602 603 self.fp.write()
603 604 else:
604 605 self.fp.dump()
605 606 self.fp.close()
606 607
607 608 header = madrigal.cedar.CatalogHeaderCreator(self.fullname)
608 609 header.createCatalog(**self.catalog)
609 610 header.createHeader(**self.header)
610 611 header.write()
611 612
612 613 def timeFlag(self):
613 614 currentTime = self.dataOut.utctime
614 615 timeTuple = time.localtime(currentTime)
615 616 dataDay = timeTuple.tm_yday
616 617
617 618 if self.currentDay is None:
618 619 self.currentDay = dataDay
619 620 return False
620 621
621 622 #Si el dia es diferente
622 623 if dataDay != self.currentDay:
623 624 self.currentDay = dataDay
624 625 return True
625 626
626 627 else:
627 628 return False
628 629
629 630 def putData(self):
630 631
631 632 if self.dataOut.flagNoData:
632 633 return 0
633 634
634 635 if self.dataOut.flagDiscontinuousBlock or self.counter == self.blocks or self.timeFlag():
635 636 if self.counter > 0:
636 637 self.setHeader()
637 638 self.counter = 0
638 639
639 640 if self.counter == 0:
640 641 self.setFile()
641 642
642 643 self.writeBlock()
643 644 self.counter += 1
644 645
645 646 def close(self):
646 647
647 648 if self.counter > 0:
648 649 self.setHeader()
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
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
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
General Comments 0
You need to be logged in to leave comments. Login now