##// END OF EJS Templates
schain
RSS
Drifts tested
joabAM -
r1740:7f5b085e2124
parent child
Show More
Side by Side Unified
Context file:
r1740:7f5b085e2124 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,1156 +1,1159
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,ProcessingHeader
21 21 from schainpy.model.data import _noise
22 22 SPEED_OF_LIGHT = 3e8
23 23
24 24
25 25 def getNumpyDtype(dataTypeCode):
26 26
27 27 if dataTypeCode == 0:
28 28 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
29 29 elif dataTypeCode == 1:
30 30 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
31 31 elif dataTypeCode == 2:
32 32 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
33 33 elif dataTypeCode == 3:
34 34 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
35 35 elif dataTypeCode == 4:
36 36 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
37 37 elif dataTypeCode == 5:
38 38 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
39 39 else:
40 40 raise ValueError('dataTypeCode was not defined')
41 41
42 42 return numpyDtype
43 43
44 44
45 45 def getDataTypeCode(numpyDtype):
46 46
47 47 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
48 48 datatype = 0
49 49 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
50 50 datatype = 1
51 51 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
52 52 datatype = 2
53 53 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
54 54 datatype = 3
55 55 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
56 56 datatype = 4
57 57 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
58 58 datatype = 5
59 59 else:
60 60 datatype = None
61 61
62 62 return datatype
63 63
64 64
65 65 def hildebrand_sekhon(data, navg):
66 66 """
67 67 This method is for the objective determination of the noise level in Doppler spectra. This
68 68 implementation technique is based on the fact that the standard deviation of the spectral
69 69 densities is equal to the mean spectral density for white Gaussian noise
70 70
71 71 Inputs:
72 72 Data : heights
73 73 navg : numbers of averages
74 74
75 75 Return:
76 76 mean : noise's level
77 77 """
78 78
79 79 sortdata = numpy.sort(data, axis=None)
80 80 '''
81 81 lenOfData = len(sortdata)
82 82 nums_min = lenOfData*0.2
83 83
84 84 if nums_min <= 5:
85 85
86 86 nums_min = 5
87 87
88 88 sump = 0.
89 89 sumq = 0.
90 90
91 91 j = 0
92 92 cont = 1
93 93
94 94 while((cont == 1)and(j < lenOfData)):
95 95
96 96 sump += sortdata[j]
97 97 sumq += sortdata[j]**2
98 98
99 99 if j > nums_min:
100 100 rtest = float(j)/(j-1) + 1.0/navg
101 101 if ((sumq*j) > (rtest*sump**2)):
102 102 j = j - 1
103 103 sump = sump - sortdata[j]
104 104 sumq = sumq - sortdata[j]**2
105 105 cont = 0
106 106
107 107 j += 1
108 108
109 109 lnoise = sump / j
110 110 '''
111 111 return _noise.hildebrand_sekhon(sortdata, navg)
112 112
113 113
114 114 class Beam:
115 115
116 116 def __init__(self):
117 117 self.codeList = []
118 118 self.azimuthList = []
119 119 self.zenithList = []
120 120
121 121
122 122 class GenericData(object):
123 123
124 124 flagNoData = True
125 125
126 126 def copy(self, inputObj=None):
127 127
128 128 if inputObj == None:
129 129 return copy.deepcopy(self)
130 130
131 131 for key in list(inputObj.__dict__.keys()):
132 132
133 133 attribute = inputObj.__dict__[key]
134 134
135 135 # If this attribute is a tuple or list
136 136 if type(inputObj.__dict__[key]) in (tuple, list):
137 137 self.__dict__[key] = attribute[:]
138 138 continue
139 139
140 140 # If this attribute is another object or instance
141 141 if hasattr(attribute, '__dict__'):
142 142 self.__dict__[key] = attribute.copy()
143 143 continue
144 144
145 145 self.__dict__[key] = inputObj.__dict__[key]
146 146
147 147 def deepcopy(self):
148 148
149 149 return copy.deepcopy(self)
150 150
151 151 def isEmpty(self):
152 152
153 153 return self.flagNoData
154 154
155 155 def isReady(self):
156 156
157 157 return not self.flagNoData
158 158
159 159
160 160 class JROData(GenericData):
161 161
162 162 systemHeaderObj = SystemHeader()
163 163 radarControllerHeaderObj = RadarControllerHeader()
164 164 type = None
165 165 datatype = None # dtype but in string
166 166 nProfiles = None
167 167 heightList = None
168 168 channelList = None
169 169 flagDiscontinuousBlock = False
170 170 useLocalTime = False
171 171 utctime = None
172 172 timeZone = None
173 173 dstFlag = None
174 174 errorCount = None
175 175 blocksize = None
176 176 flagDecodeData = False # asumo q la data no esta decodificada
177 177 flagDeflipData = False # asumo q la data no esta sin flip
178 178 flagShiftFFT = False
179 179 nCohInt = None
180 180 windowOfFilter = 1
181 181 C = 3e8
182 182 frequency = 49.92e6
183 183 realtime = False
184 184 beacon_heiIndexList = None
185 185 last_block = None
186 186 blocknow = None
187 187 azimuth = None
188 188 zenith = None
189 189 beam = Beam()
190 190 profileIndex = None
191 191 error = None
192 192 data = None
193 193 nmodes = None
194 194 metadata_list = ['heightList', 'timeZone', 'type']
195 195
196 196 ippFactor = 1 #Added to correct the freq and vel range for AMISR data
197 197 useInputBuffer = False
198 198 buffer_empty = True
199 199 codeList = []
200 200 azimuthList = []
201 201 elevationList = []
202 202 last_noise = None
203 203 __ipp = None
204 204 __ippSeconds = None
205 205 sampled_heightsFFT = None
206 206 pulseLength_TxA = None
207 207 deltaHeight = None
208 208 __code = None
209 209 __nCode = None
210 210 __nBaud = None
211 211 unitsDescription = "The units of the parameters are according to the International System of units (Seconds, Meter, Hertz, ...), except \
212 212 the parameters related to distances such as heightList, or heightResolution wich are in Km"
213 213
214 214
215 215
216 216 def __str__(self):
217 217
218 218 return '{} - {}'.format(self.type, self.datatime())
219 219
220 220 def getNoise(self):
221 221
222 222 raise NotImplementedError
223 223
224 224 @property
225 225 def nChannels(self):
226 226
227 227 return len(self.channelList)
228 228
229 229 @property
230 230 def channelIndexList(self):
231 231
232 232 return list(range(self.nChannels))
233 233
234 234 @property
235 235 def nHeights(self):
236 236
237 237 return len(self.heightList)
238 238
239 239 def getDeltaH(self):
240 240
241 241 return self.heightList[1] - self.heightList[0]
242 242
243 243 @property
244 244 def ltctime(self):
245 245
246 246 if self.useLocalTime:
247 247 return self.utctime - self.timeZone * 60
248 248
249 249 return self.utctime
250 250
251 251 @property
252 252 def datatime(self):
253 253
254 254 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
255 255 return datatimeValue
256 256
257 257 def getTimeRange(self):
258 258
259 259 datatime = []
260 260
261 261 datatime.append(self.ltctime)
262 262 datatime.append(self.ltctime + self.timeInterval + 1)
263 263
264 264 datatime = numpy.array(datatime)
265 265
266 266 return datatime
267 267
268 268 def getFmaxTimeResponse(self):
269 269
270 270 period = (10**-6) * self.getDeltaH() / (0.15)
271 271
272 272 PRF = 1. / (period * self.nCohInt)
273 273
274 274 fmax = PRF
275 275
276 276 return fmax
277 277
278 278 def getFmax(self):
279 279 PRF = 1. / (self.__ippSeconds * self.nCohInt)
280 280
281 281 fmax = PRF
282 282 return fmax
283 283
284 284 def getVmax(self):
285 285
286 286 _lambda = self.C / self.frequency
287 287
288 288 vmax = self.getFmax() * _lambda / 2
289 289
290 290 return vmax
291 291
292 292 ## Radar Controller Header must be immutable
293 293 @property
294 294 def ippSeconds(self):
295 295 '''
296 296 '''
297 297 #return self.radarControllerHeaderObj.ippSeconds
298 298 return self.__ippSeconds
299 299
300 300 @ippSeconds.setter
301 301 def ippSeconds(self, ippSeconds):
302 302 '''
303 303 '''
304 304 #self.radarControllerHeaderObj.ippSeconds = ippSeconds
305 305 self.__ippSeconds = ippSeconds
306 306 self.__ipp = ippSeconds*SPEED_OF_LIGHT/2000.0
307 307
308 308 @property
309 309 def code(self):
310 310 '''
311 311 '''
312 312 # return self.radarControllerHeaderObj.code
313 313 return self.__code
314 314
315 315 @code.setter
316 316 def code(self, code):
317 317 '''
318 318 '''
319 319 # self.radarControllerHeaderObj.code = code
320 320 self.__code = code
321 321
322 322 @property
323 323 def nCode(self):
324 324 '''
325 325 '''
326 326 # return self.radarControllerHeaderObj.nCode
327 327 return self.__nCode
328 328
329 329 @nCode.setter
330 330 def nCode(self, ncode):
331 331 '''
332 332 '''
333 333 # self.radarControllerHeaderObj.nCode = ncode
334 334 self.__nCode = ncode
335 335
336 336 @property
337 337 def nBaud(self):
338 338 '''
339 339 '''
340 340 # return self.radarControllerHeaderObj.nBaud
341 341 return self.__nBaud
342 342
343 343 @nBaud.setter
344 344 def nBaud(self, nbaud):
345 345 '''
346 346 '''
347 347 # self.radarControllerHeaderObj.nBaud = nbaud
348 348 self.__nBaud = nbaud
349 349
350 350 @property
351 351 def ipp(self):
352 352 '''
353 353 '''
354 354 # return self.radarControllerHeaderObj.ipp
355 355 return self.__ipp
356 356
357 357 @ipp.setter
358 358 def ipp(self, ipp):
359 359 '''
360 360 '''
361 361 # self.radarControllerHeaderObj.ipp = ipp
362 362 self.__ipp = ipp
363 363
364 364 @property
365 365 def metadata(self):
366 366 '''
367 367 '''
368 368
369 369 return {attr: getattr(self, attr) for attr in self.metadata_list}
370 370
371 371
372 372 class Voltage(JROData):
373 373
374 374 dataPP_POW = None
375 375 dataPP_DOP = None
376 376 dataPP_WIDTH = None
377 377 dataPP_SNR = None
378 378
379 379 # To use oper
380 380 flagProfilesByRange = False
381 381 nProfilesByRange = None
382 382 max_nIncohInt = 1
383 383
384 384 def __init__(self):
385 385 '''
386 386 Constructor
387 387 '''
388 388
389 389 self.useLocalTime = True
390 390 self.radarControllerHeaderObj = RadarControllerHeader()
391 391 self.systemHeaderObj = SystemHeader()
392 392 self.processingHeaderObj = ProcessingHeader()
393 393 self.type = "Voltage"
394 394 self.data = None
395 395 self.nProfiles = None
396 396 self.heightList = None
397 397 self.channelList = None
398 398 self.flagNoData = True
399 399 self.flagDiscontinuousBlock = False
400 400 self.utctime = None
401 401 self.timeZone = 0
402 402 self.dstFlag = None
403 403 self.errorCount = None
404 404 self.nCohInt = None
405 405 self.blocksize = None
406 406 self.flagCohInt = False
407 407 self.flagDecodeData = False # asumo q la data no esta decodificada
408 408 self.flagDeflipData = False # asumo q la data no esta sin flip
409 409 self.flagShiftFFT = False
410 410 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
411 411 self.profileIndex = 0
412 412 self.ippFactor=1
413 413 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
414 414 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
415 415
416 416 def getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None):
417 417 """
418 418 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
419 419
420 420 Return:
421 421 noiselevel
422 422 """
423 423
424 424 if channel != None:
425 425 data = self.data[channel,ymin_index:ymax_index]
426 426 nChannels = 1
427 427 else:
428 428 data = self.data[:,ymin_index:ymax_index]
429 429 nChannels = self.nChannels
430 430
431 431 noise = numpy.zeros(nChannels)
432 432 power = data * numpy.conjugate(data)
433 433
434 434 for thisChannel in range(nChannels):
435 435 if nChannels == 1:
436 436 daux = power[:].real
437 437 else:
438 438 daux = power[thisChannel, :].real
439 439 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
440 440
441 441 return noise
442 442
443 443 def getNoise(self, type=1, channel=None,ymin_index=None, ymax_index=None):
444 444
445 445 if type == 1:
446 446 noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index)
447 447
448 448 return noise
449 449
450 450 def getPower(self, channel=None):
451 451
452 452 if channel != None:
453 453 data = self.data[channel]
454 454 else:
455 455 data = self.data
456 456
457 457 power = data * numpy.conjugate(data)
458 458 powerdB = 10 * numpy.log10(power.real)
459 459 powerdB = numpy.squeeze(powerdB)
460 460
461 461 return powerdB
462 @property
463 def data_pow(self):
464 return self.getPower()
462 465
463 466 @property
464 467 def timeInterval(self):
465 468
466 469 return self.ippSeconds * self.nCohInt
467 470
468 471 noise = property(getNoise, "I'm the 'nHeights' property.")
469 472
470 473
471 474 class Spectra(JROData):
472 475
473 476 data_outlier = None
474 477 flagProfilesByRange = False
475 478 nProfilesByRange = None
476 479
477 480 def __init__(self):
478 481 '''
479 482 Constructor
480 483 '''
481 484
482 485 self.data_dc = None
483 486 self.data_spc = None
484 487 self.data_cspc = None
485 488 self.useLocalTime = True
486 489 self.radarControllerHeaderObj = RadarControllerHeader()
487 490 self.systemHeaderObj = SystemHeader()
488 491 self.processingHeaderObj = ProcessingHeader()
489 492 self.type = "Spectra"
490 493 self.timeZone = 0
491 494 self.nProfiles = None
492 495 self.heightList = None
493 496 self.channelList = None
494 497 self.pairsList = None
495 498 self.flagNoData = True
496 499 self.flagDiscontinuousBlock = False
497 500 self.utctime = None
498 501 self.nCohInt = None
499 502 self.nIncohInt = None
500 503 self.blocksize = None
501 504 self.nFFTPoints = None
502 505 self.wavelength = None
503 506 self.flagDecodeData = False # asumo q la data no esta decodificada
504 507 self.flagDeflipData = False # asumo q la data no esta sin flip
505 508 self.flagShiftFFT = False
506 509 self.ippFactor = 1
507 510 self.beacon_heiIndexList = []
508 511 self.noise_estimation = None
509 512 self.codeList = []
510 513 self.azimuthList = []
511 514 self.elevationList = []
512 515 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
513 516 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
514 517
515 518 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
516 519 """
517 520 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
518 521
519 522 Return:
520 523 noiselevel
521 524 """
522 525
523 526 noise = numpy.zeros(self.nChannels)
524 527
525 528 for channel in range(self.nChannels):
526 529 daux = self.data_spc[channel,
527 530 xmin_index:xmax_index, ymin_index:ymax_index]
528 531 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
529 532 noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt[channel])
530 533
531 534 return noise
532 535
533 536 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
534 537
535 538 if self.noise_estimation is not None:
536 539 # this was estimated by getNoise Operation defined in jroproc_spectra.py
537 540 return self.noise_estimation
538 541 else:
539 542 noise = self.getNoisebyHildebrand(
540 543 xmin_index, xmax_index, ymin_index, ymax_index)
541 544 return noise
542 545
543 546 def getFreqRangeTimeResponse(self, extrapoints=0):
544 547
545 548 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
546 549 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
547 550
548 551 return freqrange
549 552
550 553 def getAcfRange(self, extrapoints=0):
551 554
552 555 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
553 556 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
554 557
555 558 return freqrange
556 559
557 560 def getFreqRange(self, extrapoints=0):
558 561
559 562 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
560 563 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
561 564
562 565 return freqrange
563 566
564 567 def getVelRange(self, extrapoints=0):
565 568
566 569 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
567 570 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
568 571
569 572 if self.nmodes:
570 573 return velrange/self.nmodes
571 574 else:
572 575 return velrange
573 576
574 577 @property
575 578 def nPairs(self):
576 579
577 580 return len(self.pairsList)
578 581
579 582 @property
580 583 def pairsIndexList(self):
581 584
582 585 return list(range(self.nPairs))
583 586
584 587 @property
585 588 def normFactor(self):
586 589
587 590 pwcode = 1
588 591 if self.flagDecodeData:
589 592 try:
590 593 pwcode = numpy.sum(self.code[0]**2)
591 594 except Exception as e:
592 595 log.warning("Failed pwcode read, setting to 1")
593 596 pwcode = 1
594 597 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
595 598 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
596 599 if self.flagProfilesByRange:
597 600 normFactor *= (self.nProfilesByRange/self.nProfilesByRange.max())
598 601 return normFactor
599 602
600 603 @property
601 604 def flag_cspc(self):
602 605
603 606 if self.data_cspc is None:
604 607 return True
605 608
606 609 return False
607 610
608 611 @property
609 612 def flag_dc(self):
610 613
611 614 if self.data_dc is None:
612 615 return True
613 616
614 617 return False
615 618
616 619 @property
617 620 def timeInterval(self):
618 621
619 622 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
620 623 if self.nmodes:
621 624 return self.nmodes*timeInterval
622 625 else:
623 626 return timeInterval
624 627
625 628 def getPower(self):
626 629
627 630 factor = self.normFactor
628 631 power = numpy.zeros( (self.nChannels,self.nHeights) )
629 632 for ch in range(self.nChannels):
630 633 z = None
631 634 if hasattr(factor,'shape'):
632 635 if factor.ndim > 1:
633 636 z = self.data_spc[ch]/factor[ch]
634 637 else:
635 638 z = self.data_spc[ch]/factor
636 639 else:
637 640 z = self.data_spc[ch]/factor
638 641 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
639 642 avg = numpy.average(z, axis=0)
640 643 power[ch] = 10 * numpy.log10(avg)
641 644 return power
642 645
643 646 @property
644 647 def max_nIncohInt(self):
645 648
646 649 ints = numpy.zeros(self.nChannels)
647 650 for ch in range(self.nChannels):
648 651 if hasattr(self.nIncohInt,'shape'):
649 652 if self.nIncohInt.ndim > 1:
650 653 ints[ch,] = self.nIncohInt[ch].max()
651 654 else:
652 655 ints[ch,] = self.nIncohInt
653 656 self.nIncohInt = int(self.nIncohInt)
654 657 else:
655 658 ints[ch,] = self.nIncohInt
656 659
657 660 return ints
658 661
659 662 def getCoherence(self, pairsList=None, phase=False):
660 663
661 664 z = []
662 665 if pairsList is None:
663 666 pairsIndexList = self.pairsIndexList
664 667 else:
665 668 pairsIndexList = []
666 669 for pair in pairsList:
667 670 if pair not in self.pairsList:
668 671 raise ValueError("Pair %s is not in dataOut.pairsList" % (
669 672 pair))
670 673 pairsIndexList.append(self.pairsList.index(pair))
671 674 for i in range(len(pairsIndexList)):
672 675 pair = self.pairsList[pairsIndexList[i]]
673 676 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
674 677 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
675 678 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
676 679 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
677 680 if phase:
678 681 data = numpy.arctan2(avgcoherenceComplex.imag,
679 682 avgcoherenceComplex.real) * 180 / numpy.pi
680 683 else:
681 684 data = numpy.abs(avgcoherenceComplex)
682 685
683 686 z.append(data)
684 687
685 688 return numpy.array(z)
686 689
687 690 def setValue(self, value):
688 691
689 692 print("This property should not be initialized", value)
690 693
691 694 return
692 695
693 696 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
694 697
695 698
696 699 class SpectraHeis(Spectra):
697 700
698 701 def __init__(self):
699 702
700 703 self.radarControllerHeaderObj = RadarControllerHeader()
701 704 self.systemHeaderObj = SystemHeader()
702 705 self.type = "SpectraHeis"
703 706 self.nProfiles = None
704 707 self.heightList = None
705 708 self.channelList = None
706 709 self.flagNoData = True
707 710 self.flagDiscontinuousBlock = False
708 711 self.utctime = None
709 712 self.blocksize = None
710 713 self.profileIndex = 0
711 714 self.nCohInt = 1
712 715 self.nIncohInt = 1
713 716
714 717 @property
715 718 def normFactor(self):
716 719 pwcode = 1
717 720 if self.flagDecodeData:
718 721 pwcode = numpy.sum(self.code[0]**2)
719 722
720 723 normFactor = self.nIncohInt * self.nCohInt * pwcode
721 724
722 725 return normFactor
723 726
724 727 @property
725 728 def timeInterval(self):
726 729
727 730 return self.ippSeconds * self.nCohInt * self.nIncohInt
728 731
729 732
730 733 class Fits(JROData):
731 734
732 735 def __init__(self):
733 736
734 737 self.type = "Fits"
735 738 self.nProfiles = None
736 739 self.heightList = None
737 740 self.channelList = None
738 741 self.flagNoData = True
739 742 self.utctime = None
740 743 self.nCohInt = 1
741 744 self.nIncohInt = 1
742 745 self.useLocalTime = True
743 746 self.profileIndex = 0
744 747 self.timeZone = 0
745 748
746 749 def getTimeRange(self):
747 750
748 751 datatime = []
749 752
750 753 datatime.append(self.ltctime)
751 754 datatime.append(self.ltctime + self.timeInterval)
752 755
753 756 datatime = numpy.array(datatime)
754 757
755 758 return datatime
756 759
757 760 def getChannelIndexList(self):
758 761
759 762 return list(range(self.nChannels))
760 763
761 764 def getNoise(self, type=1):
762 765
763 766
764 767 if type == 1:
765 768 noise = self.getNoisebyHildebrand()
766 769
767 770 if type == 2:
768 771 noise = self.getNoisebySort()
769 772
770 773 if type == 3:
771 774 noise = self.getNoisebyWindow()
772 775
773 776 return noise
774 777
775 778 @property
776 779 def timeInterval(self):
777 780
778 781 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
779 782
780 783 return timeInterval
781 784
782 785 @property
783 786 def ippSeconds(self):
784 787 '''
785 788 '''
786 789 return self.ipp_sec
787 790
788 791 noise = property(getNoise, "I'm the 'nHeights' property.")
789 792
790 793
791 794 class Correlation(JROData):
792 795
793 796 def __init__(self):
794 797 '''
795 798 Constructor
796 799 '''
797 800 self.radarControllerHeaderObj = RadarControllerHeader()
798 801 self.systemHeaderObj = SystemHeader()
799 802 self.type = "Correlation"
800 803 self.data = None
801 804 self.dtype = None
802 805 self.nProfiles = None
803 806 self.heightList = None
804 807 self.channelList = None
805 808 self.flagNoData = True
806 809 self.flagDiscontinuousBlock = False
807 810 self.utctime = None
808 811 self.timeZone = 0
809 812 self.dstFlag = None
810 813 self.errorCount = None
811 814 self.blocksize = None
812 815 self.flagDecodeData = False # asumo q la data no esta decodificada
813 816 self.flagDeflipData = False # asumo q la data no esta sin flip
814 817 self.pairsList = None
815 818 self.nPoints = None
816 819
817 820 def getPairsList(self):
818 821
819 822 return self.pairsList
820 823
821 824 def getNoise(self, mode=2):
822 825
823 826 indR = numpy.where(self.lagR == 0)[0][0]
824 827 indT = numpy.where(self.lagT == 0)[0][0]
825 828
826 829 jspectra0 = self.data_corr[:, :, indR, :]
827 830 jspectra = copy.copy(jspectra0)
828 831
829 832 num_chan = jspectra.shape[0]
830 833 num_hei = jspectra.shape[2]
831 834
832 835 freq_dc = jspectra.shape[1] / 2
833 836 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
834 837
835 838 if ind_vel[0] < 0:
836 839 ind_vel[list(range(0, 1))] = ind_vel[list(
837 840 range(0, 1))] + self.num_prof
838 841
839 842 if mode == 1:
840 843 jspectra[:, freq_dc, :] = (
841 844 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
842 845
843 846 if mode == 2:
844 847
845 848 vel = numpy.array([-2, -1, 1, 2])
846 849 xx = numpy.zeros([4, 4])
847 850
848 851 for fil in range(4):
849 852 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
850 853
851 854 xx_inv = numpy.linalg.inv(xx)
852 855 xx_aux = xx_inv[0, :]
853 856
854 857 for ich in range(num_chan):
855 858 yy = jspectra[ich, ind_vel, :]
856 859 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
857 860
858 861 junkid = jspectra[ich, freq_dc, :] <= 0
859 862 cjunkid = sum(junkid)
860 863
861 864 if cjunkid.any():
862 865 jspectra[ich, freq_dc, junkid.nonzero()] = (
863 866 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
864 867
865 868 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
866 869
867 870 return noise
868 871
869 872 @property
870 873 def timeInterval(self):
871 874
872 875 return self.ippSeconds * self.nCohInt * self.nProfiles
873 876
874 877 def splitFunctions(self):
875 878
876 879 pairsList = self.pairsList
877 880 ccf_pairs = []
878 881 acf_pairs = []
879 882 ccf_ind = []
880 883 acf_ind = []
881 884 for l in range(len(pairsList)):
882 885 chan0 = pairsList[l][0]
883 886 chan1 = pairsList[l][1]
884 887
885 888 # Obteniendo pares de Autocorrelacion
886 889 if chan0 == chan1:
887 890 acf_pairs.append(chan0)
888 891 acf_ind.append(l)
889 892 else:
890 893 ccf_pairs.append(pairsList[l])
891 894 ccf_ind.append(l)
892 895
893 896 data_acf = self.data_cf[acf_ind]
894 897 data_ccf = self.data_cf[ccf_ind]
895 898
896 899 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
897 900
898 901 @property
899 902 def normFactor(self):
900 903 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
901 904 acf_pairs = numpy.array(acf_pairs)
902 905 normFactor = numpy.zeros((self.nPairs, self.nHeights))
903 906
904 907 for p in range(self.nPairs):
905 908 pair = self.pairsList[p]
906 909
907 910 ch0 = pair[0]
908 911 ch1 = pair[1]
909 912
910 913 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
911 914 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
912 915 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
913 916
914 917 return normFactor
915 918
916 919
917 920 class Parameters(Spectra):
918 921
919 922 groupList = None # List of Pairs, Groups, etc
920 923 data_param = None # Parameters obtained
921 924 data_pre = None # Data Pre Parametrization
922 925 data_SNR = None # Signal to Noise Ratio
923 926 abscissaList = None # Abscissa, can be velocities, lags or time
924 927 utctimeInit = None # Initial UTC time
925 928 paramInterval = None # Time interval to calculate Parameters in seconds
926 929 useLocalTime = True
927 930 # Fitting
928 931 data_error = None # Error of the estimation
929 932 constants = None
930 933 library = None
931 934 # Output signal
932 935 outputInterval = None # Time interval to calculate output signal in seconds
933 936 data_output = None # Out signal
934 937 nAvg = None
935 938 noise_estimation = None
936 939 GauSPC = None # Fit gaussian SPC
937 940
938 941 data_outlier = None
939 942 data_vdrift = None
940 943 radarControllerHeaderTxt=None #header Controller like text
941 944 txPower = None
942 945 flagProfilesByRange = False
943 946 nProfilesByRange = None
944 947
945 948
946 949 def __init__(self):
947 950 '''
948 951 Constructor
949 952 '''
950 953 self.radarControllerHeaderObj = RadarControllerHeader()
951 954 self.systemHeaderObj = SystemHeader()
952 955 self.processingHeaderObj = ProcessingHeader()
953 956 self.type = "Parameters"
954 957 self.timeZone = 0
955 958
956 959 def getTimeRange1(self, interval):
957 960
958 961 datatime = []
959 962
960 963 if self.useLocalTime:
961 964 time1 = self.utctimeInit - self.timeZone * 60
962 965 else:
963 966 time1 = self.utctimeInit
964 967
965 968 datatime.append(time1)
966 969 datatime.append(time1 + interval)
967 970 datatime = numpy.array(datatime)
968 971
969 972 return datatime
970 973
971 974 @property
972 975 def timeInterval(self):
973 976
974 977 if hasattr(self, 'timeInterval1'):
975 978 return self.timeInterval1
976 979 else:
977 980 return self.paramInterval
978 981
979 982 def setValue(self, value):
980 983
981 984 print("This property should not be initialized")
982 985
983 986 return
984 987
985 988 def getNoise(self):
986 989
987 990 return self.spc_noise
988 991
989 992 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
990 993
991 994
992 995 class PlotterData(object):
993 996 '''
994 997 Object to hold data to be plotted
995 998 '''
996 999
997 1000 MAXNUMX = 200
998 1001 MAXNUMY = 200
999 1002
1000 1003 def __init__(self, code, exp_code, localtime=True):
1001 1004
1002 1005 self.key = code
1003 1006 self.exp_code = exp_code
1004 1007 self.ready = False
1005 1008 self.flagNoData = False
1006 1009 self.localtime = localtime
1007 1010 self.data = {}
1008 1011 self.meta = {}
1009 1012 self.__heights = []
1010 1013
1011 1014 def __str__(self):
1012 1015 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1013 1016 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1014 1017
1015 1018 def __len__(self):
1016 1019 return len(self.data)
1017 1020
1018 1021 def __getitem__(self, key):
1019 1022 if isinstance(key, int):
1020 1023 return self.data[self.times[key]]
1021 1024 elif isinstance(key, str):
1022 1025 ret = numpy.array([self.data[x][key] for x in self.times])
1023 1026 if ret.ndim > 1:
1024 1027 ret = numpy.swapaxes(ret, 0, 1)
1025 1028 return ret
1026 1029
1027 1030 def __contains__(self, key):
1028 1031 return key in self.data[self.min_time]
1029 1032
1030 1033 def setup(self):
1031 1034 '''
1032 1035 Configure object
1033 1036 '''
1034 1037 self.type = ''
1035 1038 self.ready = False
1036 1039 del self.data
1037 1040 self.data = {}
1038 1041 self.__heights = []
1039 1042 self.__all_heights = set()
1040 1043
1041 1044 def shape(self, key):
1042 1045 '''
1043 1046 Get the shape of the one-element data for the given key
1044 1047 '''
1045 1048
1046 1049 if len(self.data[self.min_time][key]):
1047 1050 return self.data[self.min_time][key].shape
1048 1051 return (0,)
1049 1052
1050 1053 def update(self, data, tm, meta={}):
1051 1054 '''
1052 1055 Update data object with new dataOut
1053 1056 '''
1054 1057
1055 1058 self.data[tm] = data
1056 1059
1057 1060 for key, value in meta.items():
1058 1061 setattr(self, key, value)
1059 1062
1060 1063 def normalize_heights(self):
1061 1064 '''
1062 1065 Ensure same-dimension of the data for different heighList
1063 1066 '''
1064 1067
1065 1068 H = numpy.array(list(self.__all_heights))
1066 1069 H.sort()
1067 1070 for key in self.data:
1068 1071 shape = self.shape(key)[:-1] + H.shape
1069 1072 for tm, obj in list(self.data[key].items()):
1070 1073 h = self.__heights[self.times.tolist().index(tm)]
1071 1074 if H.size == h.size:
1072 1075 continue
1073 1076 index = numpy.where(numpy.in1d(H, h))[0]
1074 1077 dummy = numpy.zeros(shape) + numpy.nan
1075 1078 if len(shape) == 2:
1076 1079 dummy[:, index] = obj
1077 1080 else:
1078 1081 dummy[index] = obj
1079 1082 self.data[key][tm] = dummy
1080 1083
1081 1084 self.__heights = [H for tm in self.times]
1082 1085
1083 1086 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1084 1087 '''
1085 1088 Convert data to json
1086 1089 '''
1087 1090
1088 1091 meta = {}
1089 1092 meta['xrange'] = []
1090 1093 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1091 1094 tmp = self.data[tm][self.key]
1092 1095 shape = tmp.shape
1093 1096 if len(shape) == 2:
1094 1097 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1095 1098 elif len(shape) == 3:
1096 1099 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1097 1100 data = self.roundFloats(
1098 1101 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1099 1102 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1100 1103 else:
1101 1104 data = self.roundFloats(self.data[tm][self.key].tolist())
1102 1105
1103 1106 ret = {
1104 1107 'plot': plot_name,
1105 1108 'code': self.exp_code,
1106 1109 'time': float(tm),
1107 1110 'data': data,
1108 1111 }
1109 1112 meta['type'] = plot_type
1110 1113 meta['interval'] = float(self.interval)
1111 1114 meta['localtime'] = self.localtime
1112 1115 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1113 1116 meta.update(self.meta)
1114 1117 ret['metadata'] = meta
1115 1118 return json.dumps(ret)
1116 1119
1117 1120 @property
1118 1121 def times(self):
1119 1122 '''
1120 1123 Return the list of times of the current data
1121 1124 '''
1122 1125
1123 1126 ret = [t for t in self.data]
1124 1127 ret.sort()
1125 1128 return numpy.array(ret)
1126 1129
1127 1130 @property
1128 1131 def min_time(self):
1129 1132 '''
1130 1133 Return the minimun time value
1131 1134 '''
1132 1135
1133 1136 return self.times[0]
1134 1137
1135 1138 @property
1136 1139 def max_time(self):
1137 1140 '''
1138 1141 Return the maximun time value
1139 1142 '''
1140 1143
1141 1144 return self.times[-1]
1142 1145
1143 1146 # @property
1144 1147 # def heights(self):
1145 1148 # '''
1146 1149 # Return the list of heights of the current data
1147 1150 # '''
1148 1151
1149 1152 # return numpy.array(self.__heights[-1])
1150 1153
1151 1154 @staticmethod
1152 1155 def roundFloats(obj):
1153 1156 if isinstance(obj, list):
1154 1157 return list(map(PlotterData.roundFloats, obj))
1155 1158 elif isinstance(obj, float):
1156 1159 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,437 +1,437
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
7 7 from schainpy.utils import log
8 8
9 9 EARTH_RADIUS = 6.3710e3
10 10
11 11
12 12 def ll2xy(lat1, lon1, lat2, lon2):
13 13
14 14 p = 0.017453292519943295
15 15 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
16 16 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
17 17 r = 12742 * numpy.arcsin(numpy.sqrt(a))
18 18 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
19 19 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
20 20 theta = -theta + numpy.pi/2
21 21 return r*numpy.cos(theta), r*numpy.sin(theta)
22 22
23 23
24 24 def km2deg(km):
25 25 '''
26 26 Convert distance in km to degrees
27 27 '''
28 28
29 29 return numpy.rad2deg(km/EARTH_RADIUS)
30 30
31 31
32 32
33 33 class SpectralMomentsPlot(SpectraPlot):
34 34 '''
35 35 Plot for Spectral Moments
36 36 '''
37 37 CODE = 'spc_moments'
38 38 # colormap = 'jet'
39 39 # plot_type = 'pcolor'
40 40
41 41 class DobleGaussianPlot(SpectraPlot):
42 42 '''
43 43 Plot for Double Gaussian Plot
44 44 '''
45 45 CODE = 'gaussian_fit'
46 46 # colormap = 'jet'
47 47 # plot_type = 'pcolor'
48 48
49 49 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
50 50 '''
51 51 Plot SpectraCut with Double Gaussian Fit
52 52 '''
53 53 CODE = 'cut_gaussian_fit'
54 54
55 55
56 56 class SpectralFitObliquePlot(SpectraPlot):
57 57 '''
58 58 Plot for Spectral Oblique
59 59 '''
60 60 CODE = 'spc_moments'
61 61 colormap = 'jet'
62 62 plot_type = 'pcolor'
63 63
64 64
65 65 class SnrPlot(RTIPlot):
66 66 '''
67 67 Plot for SNR Data
68 68 '''
69 69
70 70 CODE = 'snr'
71 71 colormap = 'jet'
72 72
73 73 def update(self, dataOut):
74 74 if len(self.channelList) == 0:
75 75 self.update_list(dataOut)
76 76
77 77 meta = {}
78 78 data = {
79 79 'snr': 10 * numpy.log10(dataOut.data_snr)
80 80 }
81 81 return data, meta
82 82
83 83 class DopplerPlot(RTIPlot):
84 84 '''
85 85 Plot for DOPPLER Data (1st moment)
86 86 '''
87 87
88 88 CODE = 'dop'
89 89 colormap = 'RdBu_r'
90 90
91 91 def update(self, dataOut):
92 92 self.update_list(dataOut)
93 93 data = {
94 94 'dop': dataOut.data_dop
95 95 }
96 96
97 97 return data, {}
98 98
99 99 class PowerPlot(RTIPlot):
100 100 '''
101 101 Plot for Power Data (0 moment)
102 102 '''
103 103
104 104 CODE = 'pow'
105 105 colormap = 'jet'
106 106
107 107 def update(self, dataOut):
108 108 self.update_list(dataOut)
109 109 data = {
110 110 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
111 111 }
112 112 try:
113 113 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
114 114 except:
115 115 pass
116 116 return data, {}
117 117
118 118 class SpectralWidthPlot(RTIPlot):
119 119 '''
120 120 Plot for Spectral Width Data (2nd moment)
121 121 '''
122 122
123 123 CODE = 'width'
124 124 colormap = 'jet'
125 125
126 126 def update(self, dataOut):
127 127 self.update_list(dataOut)
128 128 data = {
129 129 'width': dataOut.data_width
130 130 }
131 131 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
132 132 return data, {}
133 133
134 134 class SkyMapPlot(Plot):
135 135 '''
136 136 Plot for meteors detection data
137 137 '''
138 138
139 139 CODE = 'param'
140 140
141 141 def setup(self):
142 142
143 143 self.ncols = 1
144 144 self.nrows = 1
145 145 self.width = 7.2
146 146 self.height = 7.2
147 147 self.nplots = 1
148 148 self.xlabel = 'Zonal Zenith Angle (deg)'
149 149 self.ylabel = 'Meridional Zenith Angle (deg)'
150 150 self.polar = True
151 151 self.ymin = -180
152 152 self.ymax = 180
153 153 self.colorbar = False
154 154
155 155 def plot(self):
156 156
157 157 arrayParameters = numpy.concatenate(self.data['param'])
158 158 error = arrayParameters[:, -1]
159 159 indValid = numpy.where(error == 0)[0]
160 160 finalMeteor = arrayParameters[indValid, :]
161 161 finalAzimuth = finalMeteor[:, 3]
162 162 finalZenith = finalMeteor[:, 4]
163 163
164 164 x = finalAzimuth * numpy.pi / 180
165 165 y = finalZenith
166 166
167 167 ax = self.axes[0]
168 168
169 169 if ax.firsttime:
170 170 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
171 171 else:
172 172 ax.plot.set_data(x, y)
173 173
174 174 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
175 175 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
176 176 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
177 177 dt2,
178 178 len(x))
179 179 self.titles[0] = title
180 180
181 181 class GenericRTIPlot(Plot):
182 182 '''
183 183 Plot for data_xxxx object
184 184 '''
185 185
186 186 CODE = 'param'
187 187 colormap = 'viridis'
188 188 plot_type = 'pcolorbuffer'
189 189
190 190 def setup(self):
191 191 self.xaxis = 'time'
192 192 self.ncols = 1
193 193 self.nrows = self.data.shape('param')[0]
194 194 self.nplots = self.nrows
195 195 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
196 196
197 197 if not self.xlabel:
198 198 self.xlabel = 'Time'
199 199
200 200 self.ylabel = 'Range [km]'
201 201 if not self.titles:
202 202 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
203 203
204 204 def update(self, dataOut):
205 205
206 206 data = {
207 207 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
208 208 }
209
209
210 210 meta = {}
211 211
212 212 return data, meta
213 213
214 214 def plot(self):
215 215 # self.data.normalize_heights()
216 216 self.x = self.data.times
217 217 self.y = self.data.yrange
218 218 self.z = self.data['param']
219 219
220 220 self.z = numpy.ma.masked_invalid(self.z)
221 221
222 222 if self.decimation is None:
223 223 x, y, z = self.fill_gaps(self.x, self.y, self.z)
224 224 else:
225 225 x, y, z = self.fill_gaps(*self.decimate())
226 226
227 227 for n, ax in enumerate(self.axes):
228 228
229 229 self.zmax = self.zmax if self.zmax is not None else numpy.max(
230 230 self.z[n])
231 231 self.zmin = self.zmin if self.zmin is not None else numpy.min(
232 232 self.z[n])
233 233
234 234 if ax.firsttime:
235 235 if self.zlimits is not None:
236 236 self.zmin, self.zmax = self.zlimits[n]
237 237
238 238 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
239 239 vmin=self.zmin,
240 240 vmax=self.zmax,
241 241 cmap=self.cmaps[n]
242 242 )
243 243 else:
244 244 if self.zlimits is not None:
245 245 self.zmin, self.zmax = self.zlimits[n]
246 246 try:
247 247 ax.collections.remove(ax.collections[0])
248 248 except:
249 249 pass
250 250 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
251 251 vmin=self.zmin,
252 252 vmax=self.zmax,
253 253 cmap=self.cmaps[n]
254 254 )
255 255
256 256
257 257 class PolarMapPlot(Plot):
258 258 '''
259 259 Plot for weather radar
260 260 '''
261 261
262 262 CODE = 'param'
263 263 colormap = 'seismic'
264 264
265 265 def setup(self):
266 266 self.ncols = 1
267 267 self.nrows = 1
268 268 self.width = 9
269 269 self.height = 8
270 270 self.mode = self.data.meta['mode']
271 271 if self.channels is not None:
272 272 self.nplots = len(self.channels)
273 273 self.nrows = len(self.channels)
274 274 else:
275 275 self.nplots = self.data.shape(self.CODE)[0]
276 276 self.nrows = self.nplots
277 277 self.channels = list(range(self.nplots))
278 278 if self.mode == 'E':
279 279 self.xlabel = 'Longitude'
280 280 self.ylabel = 'Latitude'
281 281 else:
282 282 self.xlabel = 'Range (km)'
283 283 self.ylabel = 'Height (km)'
284 284 self.bgcolor = 'white'
285 285 self.cb_labels = self.data.meta['units']
286 286 self.lat = self.data.meta['latitude']
287 287 self.lon = self.data.meta['longitude']
288 288 self.xmin, self.xmax = float(
289 289 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
290 290 self.ymin, self.ymax = float(
291 291 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
292 292 # self.polar = True
293 293
294 294 def plot(self):
295 295
296 296 for n, ax in enumerate(self.axes):
297 297 data = self.data['param'][self.channels[n]]
298 298
299 299 zeniths = numpy.linspace(
300 300 0, self.data.meta['max_range'], data.shape[1])
301 301 if self.mode == 'E':
302 302 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
303 303 r, theta = numpy.meshgrid(zeniths, azimuths)
304 304 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
305 305 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
306 306 x = km2deg(x) + self.lon
307 307 y = km2deg(y) + self.lat
308 308 else:
309 309 azimuths = numpy.radians(self.data.yrange)
310 310 r, theta = numpy.meshgrid(zeniths, azimuths)
311 311 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
312 312 self.y = zeniths
313 313
314 314 if ax.firsttime:
315 315 if self.zlimits is not None:
316 316 self.zmin, self.zmax = self.zlimits[n]
317 317 ax.plt = ax.pcolormesh(# r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
318 318 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
319 319 vmin=self.zmin,
320 320 vmax=self.zmax,
321 321 cmap=self.cmaps[n])
322 322 else:
323 323 if self.zlimits is not None:
324 324 self.zmin, self.zmax = self.zlimits[n]
325 325 ax.collections.remove(ax.collections[0])
326 326 ax.plt = ax.pcolormesh(# r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
327 327 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
328 328 vmin=self.zmin,
329 329 vmax=self.zmax,
330 330 cmap=self.cmaps[n])
331 331
332 332 if self.mode == 'A':
333 333 continue
334 334
335 335 # plot district names
336 336 f = open('/data/workspace/schain_scripts/distrito.csv')
337 337 for line in f:
338 338 label, lon, lat = [s.strip() for s in line.split(',') if s]
339 339 lat = float(lat)
340 340 lon = float(lon)
341 341 # ax.plot(lon, lat, '.b', ms=2)
342 342 ax.text(lon, lat, label.decode('utf8'), ha='center',
343 343 va='bottom', size='8', color='black')
344 344
345 345 # plot limites
346 346 limites = []
347 347 tmp = []
348 348 for line in open('/data/workspace/schain_scripts/lima.csv'):
349 349 if '#' in line:
350 350 if tmp:
351 351 limites.append(tmp)
352 352 tmp = []
353 353 continue
354 354 values = line.strip().split(',')
355 355 tmp.append((float(values[0]), float(values[1])))
356 356 for points in limites:
357 357 ax.add_patch(
358 358 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
359 359
360 360 # plot Cuencas
361 361 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
362 362 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
363 363 values = [line.strip().split(',') for line in f]
364 364 points = [(float(s[0]), float(s[1])) for s in values]
365 365 ax.add_patch(Polygon(points, ec='b', fc='none'))
366 366
367 367 # plot grid
368 368 for r in (15, 30, 45, 60):
369 369 ax.add_artist(plt.Circle((self.lon, self.lat),
370 370 km2deg(r), color='0.6', fill=False, lw=0.2))
371 371 ax.text(
372 372 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
373 373 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
374 374 '{}km'.format(r),
375 375 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
376 376
377 377 if self.mode == 'E':
378 378 title = 'El={}\N{DEGREE SIGN}'.format(self.data.meta['elevation'])
379 379 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
380 380 else:
381 381 title = 'Az={}\N{DEGREE SIGN}'.format(self.data.meta['azimuth'])
382 382 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
383 383
384 384 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
385 385 self.titles = ['{} {}'.format(
386 386 self.data.parameters[x], title) for x in self.channels]
387 387
388 388
389 389
390 390 class TxPowerPlot(Plot):
391 391 '''
392 392 Plot for TX Power from external file
393 393 '''
394 394
395 395 CODE = 'tx_power'
396 396 plot_type = 'scatterbuffer'
397 397
398 398 def setup(self):
399 399 self.xaxis = 'time'
400 400 self.ncols = 1
401 401 self.nrows = 1
402 402 self.nplots = 1
403 403 self.ylabel = 'Power [kW]'
404 404 self.xlabel = 'Time'
405 405 self.titles = ['TX power']
406 406 self.colorbar = False
407 407 self.plots_adjust.update({'right': 0.85 })
408 408 #if not self.titles:
409 409 self.titles = ['TX Power Plot']
410 410
411 411 def update(self, dataOut):
412 412
413 413 data = {}
414 414 meta = {}
415 415
416 416 data['tx_power'] = dataOut.txPower/1000
417 417 meta['yrange'] = numpy.array([])
418 418 #print(dataOut.txPower/1000)
419 419 return data, meta
420 420
421 421 def plot(self):
422 422
423 423 x = self.data.times
424 424 xmin = self.data.min_time
425 425 xmax = xmin + self.xrange * 60 * 60
426 426 Y = self.data['tx_power']
427 427
428 428 if self.axes[0].firsttime:
429 429 if self.ymin is None: self.ymin = 0
430 430 if self.ymax is None: self.ymax = numpy.nanmax(Y) + 5
431 431 if self.ymax == 5:
432 432 self.ymax = 250
433 433 self.ymin = 100
434 434 self.axes[0].plot(x, Y, lw=1, label='Power')
435 435 plt.legend(bbox_to_anchor=(1.18, 1.0))
436 436 else:
437 437 self.axes[0].lines[0].set_data(x, Y) No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1928 +1,1935
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 datetime
12 12
13 13 from schainpy.model.graphics.jroplot_base import Plot, plt, log
14 14 from itertools import combinations
15 15 from matplotlib.ticker import LinearLocator
16 16
17 17 from schainpy.model.utils.BField import BField
18 18 from scipy.interpolate import splrep
19 19 from scipy.interpolate import splev
20 20
21 21 from matplotlib import __version__ as plt_version
22 22
23 23 if plt_version >='3.3.4':
24 24 EXTRA_POINTS = 0
25 25 else:
26 26 EXTRA_POINTS = 1
27 27 class SpectraPlot(Plot):
28 28 '''
29 29 Plot for Spectra data
30 30 '''
31 31
32 32 CODE = 'spc'
33 33 colormap = 'jet'
34 34 plot_type = 'pcolor'
35 35 buffering = False
36 36 channelList = []
37 37 elevationList = []
38 38 azimuthList = []
39 39
40 40 def setup(self):
41 41
42 42 self.nplots = len(self.data.channels)
43 43 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
44 44 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
45 45 self.height = 3.4 * self.nrows
46 46 self.cb_label = 'dB'
47 47 if self.showprofile:
48 48 self.width = 5.2 * self.ncols
49 49 else:
50 50 self.width = 4.2* self.ncols
51 51 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.12})
52 52 self.ylabel = 'Range [km]'
53 53
54 54 def update_list(self,dataOut):
55 55
56 56 if len(self.channelList) == 0:
57 57 self.channelList = dataOut.channelList
58 58 if len(self.elevationList) == 0:
59 59 self.elevationList = dataOut.elevationList
60 60 if len(self.azimuthList) == 0:
61 61 self.azimuthList = dataOut.azimuthList
62 62
63 63 def update(self, dataOut):
64 64
65 65 self.update_list(dataOut)
66 66 data = {}
67 67 meta = {}
68
69 68 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
70 noise = 10*numpy.log10(dataOut.getNoise()/norm)
71 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
72 for ch in range(dataOut.nChannels):
73 if hasattr(dataOut.normFactor,'ndim'):
74 if dataOut.normFactor.ndim > 1:
75 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
69 if dataOut.type == "Parameters":
70 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
71 spc = 10*numpy.log10(dataOut.data_spc/(dataOut.nProfiles))
72 else:
73 noise = 10*numpy.log10(dataOut.getNoise()/norm)
74
75 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
76 for ch in range(dataOut.nChannels):
77 if hasattr(dataOut.normFactor,'ndim'):
78 if dataOut.normFactor.ndim > 1:
79 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
76 80
81 else:
82 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
77 83 else:
78 84 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
79 else:
80 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
81 85
82 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
83 spc = 10*numpy.log10(z)
86 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
87 spc = 10*numpy.log10(z)
84 88
85 89 data['spc'] = spc
86 90 data['rti'] = spc.mean(axis=1)
87 91 data['noise'] = noise
88 92 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
89 93 if self.CODE == 'spc_moments':
90 94 data['moments'] = dataOut.moments
91 95
92 96 return data, meta
93 97
94 98 def plot(self):
95 99
96 100 if self.xaxis == "frequency":
97 101 x = self.data.xrange[0]
98 102 self.xlabel = "Frequency (kHz)"
99 103 elif self.xaxis == "time":
100 104 x = self.data.xrange[1]
101 105 self.xlabel = "Time (ms)"
102 106 else:
103 107 x = self.data.xrange[2]
104 108 self.xlabel = "Velocity (m/s)"
105 109
106 110 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
107 111 x = self.data.xrange[2]
108 112 self.xlabel = "Velocity (m/s)"
109 113
110 114 self.titles = []
111 115
112 116 y = self.data.yrange
113 117 self.y = y
114 118
115 119 data = self.data[-1]
116 120 z = data['spc']
117 121
118 122 for n, ax in enumerate(self.axes):
119 123 noise = self.data['noise'][n][0]
120 124 # noise = data['noise'][n]
121 125
122 126 if self.CODE == 'spc_moments':
123 127 mean = data['moments'][n, 1]
124 128 if self.CODE == 'gaussian_fit':
125 129 gau0 = data['gaussfit'][n][2,:,0]
126 130 gau1 = data['gaussfit'][n][2,:,1]
127 131 if ax.firsttime:
128 132 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
129 133 self.xmin = self.xmin if self.xmin else -self.xmax
130 134 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
131 135 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
132 136 ax.plt = ax.pcolormesh(x, y, z[n].T,
133 137 vmin=self.zmin,
134 138 vmax=self.zmax,
135 139 cmap=plt.get_cmap(self.colormap)
136 140 )
137 141
138 142 if self.showprofile:
139 143 ax.plt_profile = self.pf_axes[n].plot(
140 144 data['rti'][n], y)[0]
141 145 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
142 146 color="k", linestyle="dashed", lw=1)[0]
143 147 if self.CODE == 'spc_moments':
144 148 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
145 149 if self.CODE == 'gaussian_fit':
146 150 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
147 151 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
148 152 else:
149 153 ax.plt.set_array(z[n].T.ravel())
150 154 if self.showprofile:
151 155 ax.plt_profile.set_data(data['rti'][n], y)
152 156 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
153 157 if self.CODE == 'spc_moments':
154 158 ax.plt_mean.set_data(mean, y)
155 159 if self.CODE == 'gaussian_fit':
156 160 ax.plt_gau0.set_data(gau0, y)
157 161 ax.plt_gau1.set_data(gau1, y)
158 162 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
159 163 self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
160 164 else:
161 165 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
162 166
163 167 class SpectraObliquePlot(Plot):
164 168 '''
165 169 Plot for Spectra data
166 170 '''
167 171
168 172 CODE = 'spc_oblique'
169 173 colormap = 'jet'
170 174 plot_type = 'pcolor'
171 175
172 176 def setup(self):
173 177 self.xaxis = "oblique"
174 178 self.nplots = len(self.data.channels)
175 179 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
176 180 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
177 181 self.height = 2.6 * self.nrows
178 182 self.cb_label = 'dB'
179 183 if self.showprofile:
180 184 self.width = 4 * self.ncols
181 185 else:
182 186 self.width = 3.5 * self.ncols
183 187 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
184 188 self.ylabel = 'Range [km]'
185 189
186 190 def update(self, dataOut):
187 191
188 192 data = {}
189 193 meta = {}
190 194
191 195 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
192 196 data['spc'] = spc
193 197 data['rti'] = dataOut.getPower()
194 198 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
195 199 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
196 200
197 201 data['shift1'] = dataOut.Dop_EEJ_T1[0]
198 202 data['shift2'] = dataOut.Dop_EEJ_T2[0]
199 203 data['max_val_2'] = dataOut.Oblique_params[0,-1,:]
200 204 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
201 205 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
202 206
203 207 return data, meta
204 208
205 209 def plot(self):
206 210
207 211 if self.xaxis == "frequency":
208 212 x = self.data.xrange[0]
209 213 self.xlabel = "Frequency (kHz)"
210 214 elif self.xaxis == "time":
211 215 x = self.data.xrange[1]
212 216 self.xlabel = "Time (ms)"
213 217 else:
214 218 x = self.data.xrange[2]
215 219 self.xlabel = "Velocity (m/s)"
216 220
217 221 self.titles = []
218 222
219 223 y = self.data.yrange
220 224 self.y = y
221 225
222 226 data = self.data[-1]
223 227 z = data['spc']
224 228
225 229 for n, ax in enumerate(self.axes):
226 230 noise = self.data['noise'][n][-1]
227 231 shift1 = data['shift1']
228 232 shift2 = data['shift2']
229 233 max_val_2 = data['max_val_2']
230 234 err1 = data['shift1_error']
231 235 err2 = data['shift2_error']
232 236 if ax.firsttime:
233 237
234 238 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
235 239 self.xmin = self.xmin if self.xmin else -self.xmax
236 240 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
237 241 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
238 242 ax.plt = ax.pcolormesh(x, y, z[n].T,
239 243 vmin=self.zmin,
240 244 vmax=self.zmax,
241 245 cmap=plt.get_cmap(self.colormap)
242 246 )
243 247
244 248 if self.showprofile:
245 249 ax.plt_profile = self.pf_axes[n].plot(
246 250 self.data['rti'][n][-1], y)[0]
247 251 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
248 252 color="k", linestyle="dashed", lw=1)[0]
249 253
250 254 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)
251 255 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)
252 256 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)
253 257
254 258 else:
255 259 self.ploterr1.remove()
256 260 self.ploterr2.remove()
257 261 self.ploterr3.remove()
258 262 ax.plt.set_array(z[n].T.ravel())
259 263 if self.showprofile:
260 264 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
261 265 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
262 266 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)
263 267 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)
264 268 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)
265 269
266 270 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
267 271
268 272
269 273 class CrossSpectraPlot(Plot):
270 274
271 275 CODE = 'cspc'
272 276 colormap = 'jet'
273 277 plot_type = 'pcolor'
274 278 zmin_coh = None
275 279 zmax_coh = None
276 280 zmin_phase = None
277 281 zmax_phase = None
278 282 realChannels = None
279 283 crossPairs = None
280 284
281 285 def setup(self):
282 286
283 287 self.ncols = 4
284 288 self.nplots = len(self.data.pairs) * 2
285 289 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
286 290 self.width = 3.1 * self.ncols
287 291 self.height = 2.6 * self.nrows
288 292 self.ylabel = 'Range [km]'
289 293 self.showprofile = False
290 294 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
291 295
292 296 def update(self, dataOut):
293 297
294 298 data = {}
295 299 meta = {}
296 300
297 301 spc = dataOut.data_spc
298 302 cspc = dataOut.data_cspc
299 303 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
300 304 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
301 305 meta['pairs'] = rawPairs
302 306 if self.crossPairs == None:
303 307 self.crossPairs = dataOut.pairsList
304 308 tmp = []
305 309
306 310 for n, pair in enumerate(meta['pairs']):
307 311 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
308 312 coh = numpy.abs(out)
309 313 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
310 314 tmp.append(coh)
311 315 tmp.append(phase)
312 316
313 317 data['cspc'] = numpy.array(tmp)
314 318
315 319 return data, meta
316 320
317 321 def plot(self):
318 322
319 323 if self.xaxis == "frequency":
320 324 x = self.data.xrange[0]
321 325 self.xlabel = "Frequency (kHz)"
322 326 elif self.xaxis == "time":
323 327 x = self.data.xrange[1]
324 328 self.xlabel = "Time (ms)"
325 329 else:
326 330 x = self.data.xrange[2]
327 331 self.xlabel = "Velocity (m/s)"
328 332
329 333 self.titles = []
330 334
331 335 y = self.data.yrange
332 336 self.y = y
333 337
334 338 data = self.data[-1]
335 339 cspc = data['cspc']
336 340
337 341 for n in range(len(self.data.pairs)):
338 342 pair = self.crossPairs[n]
339 343 coh = cspc[n*2]
340 344 phase = cspc[n*2+1]
341 345 ax = self.axes[2 * n]
342 346 if ax.firsttime:
343 347 ax.plt = ax.pcolormesh(x, y, coh.T,
344 348 vmin=self.zmin_coh,
345 349 vmax=self.zmax_coh,
346 350 cmap=plt.get_cmap(self.colormap_coh)
347 351 )
348 352 else:
349 353 ax.plt.set_array(coh.T.ravel())
350 354 self.titles.append(
351 355 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
352 356
353 357 ax = self.axes[2 * n + 1]
354 358 if ax.firsttime:
355 359 ax.plt = ax.pcolormesh(x, y, phase.T,
356 360 vmin=-180,
357 361 vmax=180,
358 362 cmap=plt.get_cmap(self.colormap_phase)
359 363 )
360 364 else:
361 365 ax.plt.set_array(phase.T.ravel())
362 366 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
363 367
364 368
365 369 class CrossSpectra4Plot(Plot):
366 370
367 371 CODE = 'cspc'
368 372 colormap = 'jet'
369 373 plot_type = 'pcolor'
370 374 zmin_coh = None
371 375 zmax_coh = None
372 376 zmin_phase = None
373 377 zmax_phase = None
374 378
375 379 def setup(self):
376 380
377 381 self.ncols = 4
378 382 self.nrows = len(self.data.pairs)
379 383 self.nplots = self.nrows * 4
380 384 self.width = 3.1 * self.ncols
381 385 self.height = 5 * self.nrows
382 386 self.ylabel = 'Range [km]'
383 387 self.showprofile = False
384 388 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
385 389
386 390 def plot(self):
387 391
388 392 if self.xaxis == "frequency":
389 393 x = self.data.xrange[0]
390 394 self.xlabel = "Frequency (kHz)"
391 395 elif self.xaxis == "time":
392 396 x = self.data.xrange[1]
393 397 self.xlabel = "Time (ms)"
394 398 else:
395 399 x = self.data.xrange[2]
396 400 self.xlabel = "Velocity (m/s)"
397 401
398 402 self.titles = []
399 403
400 404
401 405 y = self.data.heights
402 406 self.y = y
403 407 nspc = self.data['spc']
404 408 spc = self.data['cspc'][0]
405 409 cspc = self.data['cspc'][1]
406 410
407 411 for n in range(self.nrows):
408 412 noise = self.data['noise'][:,-1]
409 413 pair = self.data.pairs[n]
410 414
411 415 ax = self.axes[4 * n]
412 416 if ax.firsttime:
413 417 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
414 418 self.xmin = self.xmin if self.xmin else -self.xmax
415 419 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
416 420 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
417 421 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
418 422 vmin=self.zmin,
419 423 vmax=self.zmax,
420 424 cmap=plt.get_cmap(self.colormap)
421 425 )
422 426 else:
423 427
424 428 ax.plt.set_array(nspc[pair[0]].T.ravel())
425 429 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
426 430
427 431 ax = self.axes[4 * n + 1]
428 432
429 433 if ax.firsttime:
430 434 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
431 435 vmin=self.zmin,
432 436 vmax=self.zmax,
433 437 cmap=plt.get_cmap(self.colormap)
434 438 )
435 439 else:
436 440
437 441 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
438 442 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
439 443
440 444 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
441 445 coh = numpy.abs(out)
442 446 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
443 447
444 448 ax = self.axes[4 * n + 2]
445 449 if ax.firsttime:
446 450 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
447 451 vmin=0,
448 452 vmax=1,
449 453 cmap=plt.get_cmap(self.colormap_coh)
450 454 )
451 455 else:
452 456 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
453 457 self.titles.append(
454 458 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
455 459
456 460 ax = self.axes[4 * n + 3]
457 461 if ax.firsttime:
458 462 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
459 463 vmin=-180,
460 464 vmax=180,
461 465 cmap=plt.get_cmap(self.colormap_phase)
462 466 )
463 467 else:
464 468 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
465 469 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
466 470
467 471
468 472 class CrossSpectra2Plot(Plot):
469 473
470 474 CODE = 'cspc'
471 475 colormap = 'jet'
472 476 plot_type = 'pcolor'
473 477 zmin_coh = None
474 478 zmax_coh = None
475 479 zmin_phase = None
476 480 zmax_phase = None
477 481
478 482 def setup(self):
479 483
480 484 self.ncols = 1
481 485 self.nrows = len(self.data.pairs)
482 486 self.nplots = self.nrows * 1
483 487 self.width = 3.1 * self.ncols
484 488 self.height = 5 * self.nrows
485 489 self.ylabel = 'Range [km]'
486 490 self.showprofile = False
487 491 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
488 492
489 493 def plot(self):
490 494
491 495 if self.xaxis == "frequency":
492 496 x = self.data.xrange[0]
493 497 self.xlabel = "Frequency (kHz)"
494 498 elif self.xaxis == "time":
495 499 x = self.data.xrange[1]
496 500 self.xlabel = "Time (ms)"
497 501 else:
498 502 x = self.data.xrange[2]
499 503 self.xlabel = "Velocity (m/s)"
500 504
501 505 self.titles = []
502 506
503 507
504 508 y = self.data.heights
505 509 self.y = y
506 510 cspc = self.data['cspc'][1]
507 511
508 512 for n in range(self.nrows):
509 513 noise = self.data['noise'][:,-1]
510 514 pair = self.data.pairs[n]
511 515 out = cspc[n]
512 516 cross = numpy.abs(out)
513 517 z = cross/self.data.nFactor
514 518 cross = 10*numpy.log10(z)
515 519
516 520 ax = self.axes[1 * n]
517 521 if ax.firsttime:
518 522 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
519 523 self.xmin = self.xmin if self.xmin else -self.xmax
520 524 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
521 525 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
522 526 ax.plt = ax.pcolormesh(x, y, cross.T,
523 527 vmin=self.zmin,
524 528 vmax=self.zmax,
525 529 cmap=plt.get_cmap(self.colormap)
526 530 )
527 531 else:
528 532 ax.plt.set_array(cross.T.ravel())
529 533 self.titles.append(
530 534 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
531 535
532 536
533 537 class CrossSpectra3Plot(Plot):
534 538
535 539 CODE = 'cspc'
536 540 colormap = 'jet'
537 541 plot_type = 'pcolor'
538 542 zmin_coh = None
539 543 zmax_coh = None
540 544 zmin_phase = None
541 545 zmax_phase = None
542 546
543 547 def setup(self):
544 548
545 549 self.ncols = 3
546 550 self.nrows = len(self.data.pairs)
547 551 self.nplots = self.nrows * 3
548 552 self.width = 3.1 * self.ncols
549 553 self.height = 5 * self.nrows
550 554 self.ylabel = 'Range [km]'
551 555 self.showprofile = False
552 556 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
553 557
554 558 def plot(self):
555 559
556 560 if self.xaxis == "frequency":
557 561 x = self.data.xrange[0]
558 562 self.xlabel = "Frequency (kHz)"
559 563 elif self.xaxis == "time":
560 564 x = self.data.xrange[1]
561 565 self.xlabel = "Time (ms)"
562 566 else:
563 567 x = self.data.xrange[2]
564 568 self.xlabel = "Velocity (m/s)"
565 569
566 570 self.titles = []
567 571
568 572
569 573 y = self.data.heights
570 574 self.y = y
571 575
572 576 cspc = self.data['cspc'][1]
573 577
574 578 for n in range(self.nrows):
575 579 noise = self.data['noise'][:,-1]
576 580 pair = self.data.pairs[n]
577 581 out = cspc[n]
578 582
579 583 cross = numpy.abs(out)
580 584 z = cross/self.data.nFactor
581 585 cross = 10*numpy.log10(z)
582 586
583 587 out_r= out.real/self.data.nFactor
584 588
585 589 out_i= out.imag/self.data.nFactor
586 590
587 591 ax = self.axes[3 * n]
588 592 if ax.firsttime:
589 593 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
590 594 self.xmin = self.xmin if self.xmin else -self.xmax
591 595 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
592 596 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
593 597 ax.plt = ax.pcolormesh(x, y, cross.T,
594 598 vmin=self.zmin,
595 599 vmax=self.zmax,
596 600 cmap=plt.get_cmap(self.colormap)
597 601 )
598 602 else:
599 603 ax.plt.set_array(cross.T.ravel())
600 604 self.titles.append(
601 605 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
602 606
603 607 ax = self.axes[3 * n + 1]
604 608 if ax.firsttime:
605 609 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
606 610 self.xmin = self.xmin if self.xmin else -self.xmax
607 611 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
608 612 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
609 613 ax.plt = ax.pcolormesh(x, y, out_r.T,
610 614 vmin=-1.e6,
611 615 vmax=0,
612 616 cmap=plt.get_cmap(self.colormap)
613 617 )
614 618 else:
615 619 ax.plt.set_array(out_r.T.ravel())
616 620 self.titles.append(
617 621 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
618 622
619 623 ax = self.axes[3 * n + 2]
620 624
621 625
622 626 if ax.firsttime:
623 627 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
624 628 self.xmin = self.xmin if self.xmin else -self.xmax
625 629 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
626 630 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
627 631 ax.plt = ax.pcolormesh(x, y, out_i.T,
628 632 vmin=-1.e6,
629 633 vmax=1.e6,
630 634 cmap=plt.get_cmap(self.colormap)
631 635 )
632 636 else:
633 637 ax.plt.set_array(out_i.T.ravel())
634 638 self.titles.append(
635 639 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
636 640
637 641 class RTIPlot(Plot):
638 642 '''
639 643 Plot for RTI data
640 644 '''
641 645
642 646 CODE = 'rti'
643 647 colormap = 'jet'
644 648 plot_type = 'pcolorbuffer'
645 649 titles = None
646 650 channelList = []
647 651 elevationList = []
648 652 azimuthList = []
649 653
650 654 def setup(self):
651 655 self.xaxis = 'time'
652 656 self.ncols = 1
653 657 self.nrows = len(self.data.channels)
654 658 self.nplots = len(self.data.channels)
655 659 self.ylabel = 'Range [km]'
656 660 #self.xlabel = 'Time'
657 661 self.cb_label = 'dB'
658 662 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
659 663 self.titles = ['{} Channel {}'.format(
660 664 self.CODE.upper(), x) for x in range(self.nplots)]
661 665
662 666 def update_list(self,dataOut):
663 667
664 668 if len(self.channelList) == 0:
665 669 self.channelList = dataOut.channelList
666 670 if len(self.elevationList) == 0:
667 671 self.elevationList = dataOut.elevationList
668 672 if len(self.azimuthList) == 0:
669 673 self.azimuthList = dataOut.azimuthList
670 674
671 675
672 676 def update(self, dataOut):
673 677
674 678 if len(self.channelList) == 0:
675 679 self.update_list(dataOut)
676 680 data = {}
677 681 meta = {}
678 682 data['rti'] = dataOut.getPower()
679 683 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
680 684 noise = 10*numpy.log10(dataOut.getNoise()/norm)
681 685 data['noise'] = noise
682 686
683 687 return data, meta
684 688
685 689 def plot(self):
686 690
687 691 self.x = self.data.times
688 692 self.y = self.data.yrange
689 693 self.z = self.data[self.CODE]
690 694 self.z = numpy.array(self.z, dtype=float)
691 695 self.z = numpy.ma.masked_invalid(self.z)
692 696
693 697 try:
694 698 if self.channelList != None:
695 699 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
696 700 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
697 701 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
698 702 else:
699 703 self.titles = ['{} Channel {}'.format(
700 704 self.CODE.upper(), x) for x in self.channelList]
701 705 except:
702 706 if self.channelList.any() != None:
703 707 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
704 708 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
705 709 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
706 710 else:
707 711 self.titles = ['{} Channel {}'.format(
708 712 self.CODE.upper(), x) for x in self.channelList]
709 713
710 714 if self.decimation is None:
711 715 x, y, z = self.fill_gaps(self.x, self.y, self.z)
712 716 else:
713 717 x, y, z = self.fill_gaps(*self.decimate())
714 718
715 719 for n, ax in enumerate(self.axes):
716 720
717 721 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
718 722 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
719 723 data = self.data[-1]
720 724 if ax.firsttime:
721 725 ax.plt = ax.pcolormesh(x, y, z[n].T,
722 726 vmin=self.zmin,
723 727 vmax=self.zmax,
724 728 cmap=plt.get_cmap(self.colormap)
725 729 )
726 730 if self.showprofile:
727 731 ax.plot_profile = self.pf_axes[n].plot(
728 data['rti'][n], self.y)[0]
732 data[self.CODE][n], self.y)[0]
729 733 if "noise" in self.data:
730 734 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
731 735 color="k", linestyle="dashed", lw=1)[0]
732 736 else:
733 737 # ax.collections.remove(ax.collections[0]) # error while running
734 738 ax.plt = ax.pcolormesh(x, y, z[n].T,
735 739 vmin=self.zmin,
736 740 vmax=self.zmax,
737 741 cmap=plt.get_cmap(self.colormap)
738 742 )
739 743 if self.showprofile:
740 ax.plot_profile.set_data(data['rti'][n], self.y)
744 ax.plot_profile.set_data(data[self.CODE][n], self.y)
741 745 if "noise" in self.data:
742 746 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
743 747 color="k", linestyle="dashed", lw=1)[0]
744 748
745 749 class SpectrogramPlot(Plot):
746 750 '''
747 751 Plot for Spectrogram data
748 752 '''
749 753
750 754 CODE = 'Spectrogram_Profile'
751 755 colormap = 'binary'
752 756 plot_type = 'pcolorbuffer'
753 757
754 758 def setup(self):
755 759 self.xaxis = 'time'
756 760 self.ncols = 1
757 761 self.nrows = len(self.data.channels)
758 762 self.nplots = len(self.data.channels)
759 763 self.xlabel = 'Time'
760 764 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
761 765 self.titles = []
762 766
763 767 self.titles = ['{} Channel {}'.format(
764 768 self.CODE.upper(), x) for x in range(self.nrows)]
765 769
766 770
767 771 def update(self, dataOut):
768 772 data = {}
769 773 meta = {}
770 774
771 775 maxHei = 1620#+12000
772 776 indb = numpy.where(dataOut.heightList <= maxHei)
773 777 hei = indb[0][-1]
774 778
775 779 factor = dataOut.nIncohInt
776 780 z = dataOut.data_spc[:,:,hei] / factor
777 781 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
778 782
779 783 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
780 784 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
781 785
782 786 data['hei'] = hei
783 787 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
784 788 data['nProfiles'] = dataOut.nProfiles
785 789
786 790 return data, meta
787 791
788 792 def plot(self):
789 793
790 794 self.x = self.data.times
791 795 self.z = self.data[self.CODE]
792 796 self.y = self.data.xrange[0]
793 797
794 798 hei = self.data['hei'][-1]
795 799 DH = self.data['DH'][-1]
796 800 nProfiles = self.data['nProfiles'][-1]
797 801
798 802 self.ylabel = "Frequency (kHz)"
799 803
800 804 self.z = numpy.ma.masked_invalid(self.z)
801 805
802 806 if self.decimation is None:
803 807 x, y, z = self.fill_gaps(self.x, self.y, self.z)
804 808 else:
805 809 x, y, z = self.fill_gaps(*self.decimate())
806 810
807 811 for n, ax in enumerate(self.axes):
808 812 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
809 813 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
810 814 data = self.data[-1]
811 815 if ax.firsttime:
812 816 ax.plt = ax.pcolormesh(x, y, z[n].T,
813 817 vmin=self.zmin,
814 818 vmax=self.zmax,
815 819 cmap=plt.get_cmap(self.colormap)
816 820 )
817 821 else:
818 822 # ax.collections.remove(ax.collections[0]) # error while running
819 823 ax.plt = ax.pcolormesh(x, y, z[n].T,
820 824 vmin=self.zmin,
821 825 vmax=self.zmax,
822 826 cmap=plt.get_cmap(self.colormap)
823 827 )
824 828
825 829
826 830
827 831 class CoherencePlot(RTIPlot):
828 832 '''
829 833 Plot for Coherence data
830 834 '''
831 835
832 836 CODE = 'coh'
833 837 titles = None
834 838
835 839 def setup(self):
836 840 self.xaxis = 'time'
837 841 self.ncols = 1
838 842 self.nrows = len(self.data.pairs)
839 843 self.nplots = len(self.data.pairs)
840 844 self.ylabel = 'Range [km]'
841 845 self.xlabel = 'Time'
842 846 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
843 847 if self.CODE == 'coh':
844 848 self.cb_label = ''
845 849 self.titles = [
846 850 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
847 851 else:
848 852 self.cb_label = 'Degrees'
849 853 self.titles = [
850 854 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
851 855
852 856 def update(self, dataOut):
853 857
854 858 data = {}
855 859 meta = {}
856 860 data['coh'] = dataOut.getCoherence()
857 861 meta['pairs'] = dataOut.pairsList
858 862
859 863 return data, meta
860 864
861 865 class PhasePlot(CoherencePlot):
862 866 '''
863 867 Plot for Phase map data
864 868 '''
865 869
866 870 CODE = 'phase'
867 871 colormap = 'seismic'
868 872
869 873 def update(self, dataOut):
870 874
871 875 data = {}
872 876 meta = {}
873 877 data['phase'] = dataOut.getCoherence(phase=True)
874 878 meta['pairs'] = dataOut.pairsList
875 879
876 880 return data, meta
877 881
878 882 class NoisePlot(Plot):
879 883 '''
880 884 Plot for noise
881 885 '''
882 886
883 887 CODE = 'noise'
884 888 plot_type = 'scatterbuffer'
885 889
886 890 def setup(self):
887 891 self.xaxis = 'time'
888 892 self.ncols = 1
889 893 self.nrows = 1
890 894 self.nplots = 1
891 895 self.ylabel = 'Intensity [dB]'
892 896 self.xlabel = 'Time'
893 897 self.titles = ['Noise']
894 898 self.colorbar = False
895 899 self.plots_adjust.update({'right': 0.85 })
900 self.titles = ['Noise Plot']
896 901
897 902 def update(self, dataOut):
898 903
899 904 data = {}
900 905 meta = {}
901 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
906 noise = 10*numpy.log10(dataOut.getNoise())
907 noise = noise.reshape(dataOut.nChannels, 1)
908 data['noise'] = noise
902 909 meta['yrange'] = numpy.array([])
903 910
904 911 return data, meta
905 912
906 913 def plot(self):
907 914
908 915 x = self.data.times
909 916 xmin = self.data.min_time
910 917 xmax = xmin + self.xrange * 60 * 60
911 918 Y = self.data['noise']
912 919
913 920 if self.axes[0].firsttime:
914 921 self.ymin = numpy.nanmin(Y) - 5
915 922 self.ymax = numpy.nanmax(Y) + 5
916 923 for ch in self.data.channels:
917 924 y = Y[ch]
918 925 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
919 926 plt.legend(bbox_to_anchor=(1.18, 1.0))
920 927 else:
921 928 for ch in self.data.channels:
922 929 y = Y[ch]
923 930 self.axes[0].lines[ch].set_data(x, y)
924 931
925 932 class PowerProfilePlot(Plot):
926 933
927 934 CODE = 'pow_profile'
928 935 plot_type = 'scatter'
929 936
930 937 def setup(self):
931 938
932 939 self.ncols = 1
933 940 self.nrows = 1
934 941 self.nplots = 1
935 942 self.height = 4
936 943 self.width = 3
937 944 self.ylabel = 'Range [km]'
938 945 self.xlabel = 'Intensity [dB]'
939 946 self.titles = ['Power Profile']
940 947 self.colorbar = False
941 948
942 949 def update(self, dataOut):
943 950
944 951 data = {}
945 952 meta = {}
946 953 data[self.CODE] = dataOut.getPower()
947 954
948 955 return data, meta
949 956
950 957 def plot(self):
951 958
952 959 y = self.data.yrange
953 960 self.y = y
954 961
955 962 x = self.data[-1][self.CODE]
956 963
957 964 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
958 965 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
959 966
960 967 if self.axes[0].firsttime:
961 968 for ch in self.data.channels:
962 969 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
963 970 plt.legend()
964 971 else:
965 972 for ch in self.data.channels:
966 973 self.axes[0].lines[ch].set_data(x[ch], y)
967 974
968 975
969 976 class SpectraCutPlot(Plot):
970 977
971 978 CODE = 'spc_cut'
972 979 plot_type = 'scatter'
973 980 buffering = False
974 981 heights = []
975 982 channelList = []
976 983 maintitle = "Spectra Cuts"
977 984 flag_setIndex = False
978 985
979 986 def setup(self):
980 987
981 988 self.nplots = len(self.data.channels)
982 989 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
983 990 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
984 991 self.width = 4.5 * self.ncols + 2.5
985 992 self.height = 4.8 * self.nrows
986 993 self.ylabel = 'Power [dB]'
987 994 self.colorbar = False
988 995 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
989 996
990 997 if len(self.selectedHeightsList) > 0:
991 998 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
992 999
993 1000
994 1001
995 1002 def update(self, dataOut):
996 1003 if len(self.channelList) == 0:
997 1004 self.channelList = dataOut.channelList
998 1005
999 1006 self.heights = dataOut.heightList
1000 1007 #print("sels: ",self.selectedHeightsList)
1001 1008 if len(self.selectedHeightsList)>0 and not self.flag_setIndex:
1002 1009
1003 1010 for sel_height in self.selectedHeightsList:
1004 1011 index_list = numpy.where(self.heights >= sel_height)
1005 1012 index_list = index_list[0]
1006 1013 self.height_index.append(index_list[0])
1007 1014 #print("sels i:"", self.height_index)
1008 1015 self.flag_setIndex = True
1009 1016 #print(self.height_index)
1010 1017 data = {}
1011 1018 meta = {}
1012 1019
1013 1020 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
1014 1021 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1015 1022 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1016 1023
1017 1024
1018 1025 z = []
1019 1026 for ch in range(dataOut.nChannels):
1020 1027 if hasattr(dataOut.normFactor,'shape'):
1021 1028 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1022 1029 else:
1023 1030 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1024 1031
1025 1032 z = numpy.asarray(z)
1026 1033 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1027 1034 spc = 10*numpy.log10(z)
1028 1035
1029 1036
1030 1037 data['spc'] = spc - noise
1031 1038 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1032 1039
1033 1040 return data, meta
1034 1041
1035 1042 def plot(self):
1036 1043 if self.xaxis == "frequency":
1037 1044 x = self.data.xrange[0][0:]
1038 1045 self.xlabel = "Frequency (kHz)"
1039 1046 elif self.xaxis == "time":
1040 1047 x = self.data.xrange[1]
1041 1048 self.xlabel = "Time (ms)"
1042 1049 else:
1043 1050 x = self.data.xrange[2]
1044 1051 self.xlabel = "Velocity (m/s)"
1045 1052
1046 1053 self.titles = []
1047 1054
1048 1055 y = self.data.yrange
1049 1056 z = self.data[-1]['spc']
1050 1057 #print(z.shape)
1051 1058 if len(self.height_index) > 0:
1052 1059 index = self.height_index
1053 1060 else:
1054 1061 index = numpy.arange(0, len(y), int((len(y))/9))
1055 1062 #print("inde x ", index, self.axes)
1056 1063
1057 1064 for n, ax in enumerate(self.axes):
1058 1065
1059 1066 if ax.firsttime:
1060 1067
1061 1068
1062 1069 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1063 1070 self.xmin = self.xmin if self.xmin else -self.xmax
1064 1071 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
1065 1072 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
1066 1073
1067 1074
1068 1075 ax.plt = ax.plot(x, z[n, :, index].T)
1069 1076 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1070 1077 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
1071 1078 ax.minorticks_on()
1072 1079 ax.grid(which='major', axis='both')
1073 1080 ax.grid(which='minor', axis='x')
1074 1081 else:
1075 1082 for i, line in enumerate(ax.plt):
1076 1083 line.set_data(x, z[n, :, index[i]])
1077 1084
1078 1085
1079 1086 self.titles.append('CH {}'.format(self.channelList[n]))
1080 1087 plt.suptitle(self.maintitle, fontsize=10)
1081 1088
1082 1089
1083 1090 class BeaconPhase(Plot):
1084 1091
1085 1092 __isConfig = None
1086 1093 __nsubplots = None
1087 1094
1088 1095 PREFIX = 'beacon_phase'
1089 1096
1090 1097 def __init__(self):
1091 1098 Plot.__init__(self)
1092 1099 self.timerange = 24*60*60
1093 1100 self.isConfig = False
1094 1101 self.__nsubplots = 1
1095 1102 self.counter_imagwr = 0
1096 1103 self.WIDTH = 800
1097 1104 self.HEIGHT = 400
1098 1105 self.WIDTHPROF = 120
1099 1106 self.HEIGHTPROF = 0
1100 1107 self.xdata = None
1101 1108 self.ydata = None
1102 1109
1103 1110 self.PLOT_CODE = BEACON_CODE
1104 1111
1105 1112 self.FTP_WEI = None
1106 1113 self.EXP_CODE = None
1107 1114 self.SUB_EXP_CODE = None
1108 1115 self.PLOT_POS = None
1109 1116
1110 1117 self.filename_phase = None
1111 1118
1112 1119 self.figfile = None
1113 1120
1114 1121 self.xmin = None
1115 1122 self.xmax = None
1116 1123
1117 1124 def getSubplots(self):
1118 1125
1119 1126 ncol = 1
1120 1127 nrow = 1
1121 1128
1122 1129 return nrow, ncol
1123 1130
1124 1131 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1125 1132
1126 1133 self.__showprofile = showprofile
1127 1134 self.nplots = nplots
1128 1135
1129 1136 ncolspan = 7
1130 1137 colspan = 6
1131 1138 self.__nsubplots = 2
1132 1139
1133 1140 self.createFigure(id = id,
1134 1141 wintitle = wintitle,
1135 1142 widthplot = self.WIDTH+self.WIDTHPROF,
1136 1143 heightplot = self.HEIGHT+self.HEIGHTPROF,
1137 1144 show=show)
1138 1145
1139 1146 nrow, ncol = self.getSubplots()
1140 1147
1141 1148 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1142 1149
1143 1150 def save_phase(self, filename_phase):
1144 1151 f = open(filename_phase,'w+')
1145 1152 f.write('\n\n')
1146 1153 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1147 1154 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1148 1155 f.close()
1149 1156
1150 1157 def save_data(self, filename_phase, data, data_datetime):
1151 1158 f=open(filename_phase,'a')
1152 1159 timetuple_data = data_datetime.timetuple()
1153 1160 day = str(timetuple_data.tm_mday)
1154 1161 month = str(timetuple_data.tm_mon)
1155 1162 year = str(timetuple_data.tm_year)
1156 1163 hour = str(timetuple_data.tm_hour)
1157 1164 minute = str(timetuple_data.tm_min)
1158 1165 second = str(timetuple_data.tm_sec)
1159 1166 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1160 1167 f.close()
1161 1168
1162 1169 def plot(self):
1163 1170 log.warning('TODO: Not yet implemented...')
1164 1171
1165 1172 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1166 1173 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1167 1174 timerange=None,
1168 1175 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1169 1176 server=None, folder=None, username=None, password=None,
1170 1177 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1171 1178
1172 1179 if dataOut.flagNoData:
1173 1180 return dataOut
1174 1181
1175 1182 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1176 1183 return
1177 1184
1178 1185 if pairsList == None:
1179 1186 pairsIndexList = dataOut.pairsIndexList[:10]
1180 1187 else:
1181 1188 pairsIndexList = []
1182 1189 for pair in pairsList:
1183 1190 if pair not in dataOut.pairsList:
1184 1191 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1185 1192 pairsIndexList.append(dataOut.pairsList.index(pair))
1186 1193
1187 1194 if pairsIndexList == []:
1188 1195 return
1189 1196
1190 1197 # if len(pairsIndexList) > 4:
1191 1198 # pairsIndexList = pairsIndexList[0:4]
1192 1199
1193 1200 hmin_index = None
1194 1201 hmax_index = None
1195 1202
1196 1203 if hmin != None and hmax != None:
1197 1204 indexes = numpy.arange(dataOut.nHeights)
1198 1205 hmin_list = indexes[dataOut.heightList >= hmin]
1199 1206 hmax_list = indexes[dataOut.heightList <= hmax]
1200 1207
1201 1208 if hmin_list.any():
1202 1209 hmin_index = hmin_list[0]
1203 1210
1204 1211 if hmax_list.any():
1205 1212 hmax_index = hmax_list[-1]+1
1206 1213
1207 1214 x = dataOut.getTimeRange()
1208 1215
1209 1216 thisDatetime = dataOut.datatime
1210 1217
1211 1218 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1212 1219 xlabel = "Local Time"
1213 1220 ylabel = "Phase (degrees)"
1214 1221
1215 1222 update_figfile = False
1216 1223
1217 1224 nplots = len(pairsIndexList)
1218 1225 phase_beacon = numpy.zeros(len(pairsIndexList))
1219 1226 for i in range(nplots):
1220 1227 pair = dataOut.pairsList[pairsIndexList[i]]
1221 1228 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1222 1229 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1223 1230 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1224 1231 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1225 1232 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1226 1233
1227 1234 if dataOut.beacon_heiIndexList:
1228 1235 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1229 1236 else:
1230 1237 phase_beacon[i] = numpy.average(phase)
1231 1238
1232 1239 if not self.isConfig:
1233 1240
1234 1241 nplots = len(pairsIndexList)
1235 1242
1236 1243 self.setup(id=id,
1237 1244 nplots=nplots,
1238 1245 wintitle=wintitle,
1239 1246 showprofile=showprofile,
1240 1247 show=show)
1241 1248
1242 1249 if timerange != None:
1243 1250 self.timerange = timerange
1244 1251
1245 1252 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1246 1253
1247 1254 if ymin == None: ymin = 0
1248 1255 if ymax == None: ymax = 360
1249 1256
1250 1257 self.FTP_WEI = ftp_wei
1251 1258 self.EXP_CODE = exp_code
1252 1259 self.SUB_EXP_CODE = sub_exp_code
1253 1260 self.PLOT_POS = plot_pos
1254 1261
1255 1262 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1256 1263 self.isConfig = True
1257 1264 self.figfile = figfile
1258 1265 self.xdata = numpy.array([])
1259 1266 self.ydata = numpy.array([])
1260 1267
1261 1268 update_figfile = True
1262 1269
1263 1270 #open file beacon phase
1264 1271 path = '%s%03d' %(self.PREFIX, self.id)
1265 1272 beacon_file = os.path.join(path,'%s.txt'%self.name)
1266 1273 self.filename_phase = os.path.join(figpath,beacon_file)
1267 1274
1268 1275 self.setWinTitle(title)
1269 1276
1270 1277
1271 1278 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1272 1279
1273 1280 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1274 1281
1275 1282 axes = self.axesList[0]
1276 1283
1277 1284 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1278 1285
1279 1286 if len(self.ydata)==0:
1280 1287 self.ydata = phase_beacon.reshape(-1,1)
1281 1288 else:
1282 1289 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1283 1290
1284 1291
1285 1292 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1286 1293 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1287 1294 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1288 1295 XAxisAsTime=True, grid='both'
1289 1296 )
1290 1297
1291 1298 self.draw()
1292 1299
1293 1300 if dataOut.ltctime >= self.xmax:
1294 1301 self.counter_imagwr = wr_period
1295 1302 self.isConfig = False
1296 1303 update_figfile = True
1297 1304
1298 1305 self.save(figpath=figpath,
1299 1306 figfile=figfile,
1300 1307 save=save,
1301 1308 ftp=ftp,
1302 1309 wr_period=wr_period,
1303 1310 thisDatetime=thisDatetime,
1304 1311 update_figfile=update_figfile)
1305 1312
1306 1313 return dataOut
1307 1314
1308 1315 #####################################
1309 1316 class NoiselessSpectraPlot(Plot):
1310 1317 '''
1311 1318 Plot for Spectra data, subtracting
1312 1319 the noise in all channels, using for
1313 1320 amisr-14 data
1314 1321 '''
1315 1322
1316 1323 CODE = 'noiseless_spc'
1317 1324 colormap = 'jet'
1318 1325 plot_type = 'pcolor'
1319 1326 buffering = False
1320 1327 channelList = []
1321 1328 last_noise = None
1322 1329
1323 1330 def setup(self):
1324 1331
1325 1332 self.nplots = len(self.data.channels)
1326 1333 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
1327 1334 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
1328 1335 self.height = 3.5 * self.nrows
1329 1336
1330 1337 self.cb_label = 'dB'
1331 1338 if self.showprofile:
1332 1339 self.width = 5.8 * self.ncols
1333 1340 else:
1334 1341 self.width = 4.8* self.ncols
1335 1342 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12})
1336 1343
1337 1344 self.ylabel = 'Range [km]'
1338 1345
1339 1346
1340 1347 def update_list(self,dataOut):
1341 1348 if len(self.channelList) == 0:
1342 1349 self.channelList = dataOut.channelList
1343 1350
1344 1351 def update(self, dataOut):
1345 1352
1346 1353 self.update_list(dataOut)
1347 1354 data = {}
1348 1355 meta = {}
1349 1356
1350 1357 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1351 1358 n0 = (dataOut.getNoise()/norm)
1352 1359 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1353 1360 noise = 10*numpy.log10(noise)
1354 1361
1355 1362 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
1356 1363 for ch in range(dataOut.nChannels):
1357 1364 if hasattr(dataOut.normFactor,'ndim'):
1358 1365 if dataOut.normFactor.ndim > 1:
1359 1366 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1360 1367 else:
1361 1368 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1362 1369 else:
1363 1370 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1364 1371
1365 1372 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1366 1373 spc = 10*numpy.log10(z)
1367 1374
1368 1375
1369 1376 data['spc'] = spc - noise
1370 1377 #print(spc.shape)
1371 1378 data['rti'] = spc.mean(axis=1)
1372 1379 data['noise'] = noise
1373 1380
1374 1381
1375 1382
1376 1383 # data['noise'] = noise
1377 1384 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1378 1385
1379 1386 return data, meta
1380 1387
1381 1388 def plot(self):
1382 1389 if self.xaxis == "frequency":
1383 1390 x = self.data.xrange[0]
1384 1391 self.xlabel = "Frequency (kHz)"
1385 1392 elif self.xaxis == "time":
1386 1393 x = self.data.xrange[1]
1387 1394 self.xlabel = "Time (ms)"
1388 1395 else:
1389 1396 x = self.data.xrange[2]
1390 1397 self.xlabel = "Velocity (m/s)"
1391 1398
1392 1399 self.titles = []
1393 1400 y = self.data.yrange
1394 1401 self.y = y
1395 1402
1396 1403 data = self.data[-1]
1397 1404 z = data['spc']
1398 1405
1399 1406 for n, ax in enumerate(self.axes):
1400 1407 #noise = data['noise'][n]
1401 1408
1402 1409 if ax.firsttime:
1403 1410 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1404 1411 self.xmin = self.xmin if self.xmin else -self.xmax
1405 1412 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1406 1413 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1407 1414 ax.plt = ax.pcolormesh(x, y, z[n].T,
1408 1415 vmin=self.zmin,
1409 1416 vmax=self.zmax,
1410 1417 cmap=plt.get_cmap(self.colormap)
1411 1418 )
1412 1419
1413 1420 if self.showprofile:
1414 1421 ax.plt_profile = self.pf_axes[n].plot(
1415 1422 data['rti'][n], y)[0]
1416 1423
1417 1424
1418 1425 else:
1419 1426 ax.plt.set_array(z[n].T.ravel())
1420 1427 if self.showprofile:
1421 1428 ax.plt_profile.set_data(data['rti'][n], y)
1422 1429
1423 1430
1424 1431 self.titles.append('CH {}'.format(self.channelList[n]))
1425 1432
1426 1433
1427 1434 class NoiselessRTIPlot(RTIPlot):
1428 1435 '''
1429 1436 Plot for RTI data
1430 1437 '''
1431 1438
1432 1439 CODE = 'noiseless_rti'
1433 1440 colormap = 'jet'
1434 1441 plot_type = 'pcolorbuffer'
1435 1442 titles = None
1436 1443 channelList = []
1437 1444 elevationList = []
1438 1445 azimuthList = []
1439 1446 last_noise = None
1440 1447
1441 1448 def setup(self):
1442 1449 self.xaxis = 'time'
1443 1450 self.ncols = 1
1444 1451 #print("dataChannels ",self.data.channels)
1445 1452 self.nrows = len(self.data.channels)
1446 1453 self.nplots = len(self.data.channels)
1447 1454 self.ylabel = 'Range [km]'
1448 1455 #self.xlabel = 'Time'
1449 1456 self.cb_label = 'dB'
1450 1457 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1451 1458 self.titles = ['{} Channel {}'.format(
1452 1459 self.CODE.upper(), x) for x in range(self.nplots)]
1453 1460
1454 1461 def update_list(self,dataOut):
1455 1462 if len(self.channelList) == 0:
1456 1463 self.channelList = dataOut.channelList
1457 1464 if len(self.elevationList) == 0:
1458 1465 self.elevationList = dataOut.elevationList
1459 1466 if len(self.azimuthList) == 0:
1460 1467 self.azimuthList = dataOut.azimuthList
1461 1468
1462 1469 def update(self, dataOut):
1463 1470 if len(self.channelList) == 0:
1464 1471 self.update_list(dataOut)
1465 1472
1466 1473 data = {}
1467 1474 meta = {}
1468 1475 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1469 1476 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt
1470 1477 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1471 1478 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1472 1479 data['noise'] = n0
1473 1480 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1474 1481 noiseless_data = dataOut.getPower() - noise
1475 1482
1476 1483 #print("power, noise:", dataOut.getPower(), n0)
1477 1484 #print(noise)
1478 1485 #print(noiseless_data)
1479 1486
1480 1487 data['noiseless_rti'] = noiseless_data
1481 1488
1482 1489 return data, meta
1483 1490
1484 1491 def plot(self):
1485 1492 from matplotlib import pyplot as plt
1486 1493 self.x = self.data.times
1487 1494 self.y = self.data.yrange
1488 1495 self.z = self.data['noiseless_rti']
1489 1496 self.z = numpy.array(self.z, dtype=float)
1490 1497 self.z = numpy.ma.masked_invalid(self.z)
1491 1498
1492 1499
1493 1500 try:
1494 1501 if self.channelList != None:
1495 1502 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1496 1503 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1497 1504 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1498 1505 else:
1499 1506 self.titles = ['{} Channel {}'.format(
1500 1507 self.CODE.upper(), x) for x in self.channelList]
1501 1508 except:
1502 1509 if self.channelList.any() != None:
1503 1510 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1504 1511 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1505 1512 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1506 1513 else:
1507 1514 self.titles = ['{} Channel {}'.format(
1508 1515 self.CODE.upper(), x) for x in self.channelList]
1509 1516
1510 1517
1511 1518 if self.decimation is None:
1512 1519 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1513 1520 else:
1514 1521 x, y, z = self.fill_gaps(*self.decimate())
1515 1522
1516 1523 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
1517 1524 #print("plot shapes ", z.shape, x.shape, y.shape)
1518 1525 #print(self.axes)
1519 1526 for n, ax in enumerate(self.axes):
1520 1527
1521 1528
1522 1529 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
1523 1530 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1524 1531 data = self.data[-1]
1525 1532 if ax.firsttime:
1526 1533 if (n+1) == len(self.channelList):
1527 1534 ax.set_xlabel('Time')
1528 1535 ax.plt = ax.pcolormesh(x, y, z[n].T,
1529 1536 vmin=self.zmin,
1530 1537 vmax=self.zmax,
1531 1538 cmap=plt.get_cmap(self.colormap)
1532 1539 )
1533 1540 if self.showprofile:
1534 1541 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1535 1542
1536 1543 else:
1537 1544 # ax.collections.remove(ax.collections[0]) # error while running
1538 1545 ax.plt = ax.pcolormesh(x, y, z[n].T,
1539 1546 vmin=self.zmin,
1540 1547 vmax=self.zmax,
1541 1548 cmap=plt.get_cmap(self.colormap)
1542 1549 )
1543 1550 if self.showprofile:
1544 1551 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1545 1552 # if "noise" in self.data:
1546 1553 # #ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
1547 1554 # ax.plot_noise.set_data(data['noise'][n], self.y)
1548 1555
1549 1556
1550 1557 class OutliersRTIPlot(Plot):
1551 1558 '''
1552 1559 Plot for data_xxxx object
1553 1560 '''
1554 1561
1555 1562 CODE = 'outlier_rtc' # Range Time Counts
1556 1563 colormap = 'cool'
1557 1564 plot_type = 'pcolorbuffer'
1558 1565
1559 1566 def setup(self):
1560 1567 self.xaxis = 'time'
1561 1568 self.ncols = 1
1562 1569 self.nrows = self.data.shape('outlier_rtc')[0]
1563 1570 self.nplots = self.nrows
1564 1571 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1565 1572
1566 1573
1567 1574 if not self.xlabel:
1568 1575 self.xlabel = 'Time'
1569 1576
1570 1577 self.ylabel = 'Height [km]'
1571 1578 if not self.titles:
1572 1579 self.titles = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)]
1573 1580
1574 1581 def update(self, dataOut):
1575 1582
1576 1583 data = {}
1577 1584 data['outlier_rtc'] = dataOut.data_outlier
1578 1585
1579 1586 meta = {}
1580 1587
1581 1588 return data, meta
1582 1589
1583 1590 def plot(self):
1584 1591 # self.data.normalize_heights()
1585 1592 self.x = self.data.times
1586 1593 self.y = self.data.yrange
1587 1594 self.z = self.data['outlier_rtc']
1588 1595
1589 1596 #self.z = numpy.ma.masked_invalid(self.z)
1590 1597
1591 1598 if self.decimation is None:
1592 1599 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1593 1600 else:
1594 1601 x, y, z = self.fill_gaps(*self.decimate())
1595 1602
1596 1603 for n, ax in enumerate(self.axes):
1597 1604
1598 1605 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1599 1606 self.z[n])
1600 1607 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1601 1608 self.z[n])
1602 1609 data = self.data[-1]
1603 1610 if ax.firsttime:
1604 1611 if self.zlimits is not None:
1605 1612 self.zmin, self.zmax = self.zlimits[n]
1606 1613
1607 1614 ax.plt = ax.pcolormesh(x, y, z[n].T,
1608 1615 vmin=self.zmin,
1609 1616 vmax=self.zmax,
1610 1617 cmap=self.cmaps[n]
1611 1618 )
1612 1619 if self.showprofile:
1613 1620 ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0]
1614 1621 self.pf_axes[n].set_xlabel('')
1615 1622 else:
1616 1623 if self.zlimits is not None:
1617 1624 self.zmin, self.zmax = self.zlimits[n]
1618 1625 # ax.collections.remove(ax.collections[0]) # error while running
1619 1626 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1620 1627 vmin=self.zmin,
1621 1628 vmax=self.zmax,
1622 1629 cmap=self.cmaps[n]
1623 1630 )
1624 1631 if self.showprofile:
1625 1632 ax.plot_profile.set_data(data['outlier_rtc'][n], self.y)
1626 1633 self.pf_axes[n].set_xlabel('')
1627 1634
1628 1635 class NIncohIntRTIPlot(Plot):
1629 1636 '''
1630 1637 Plot for data_xxxx object
1631 1638 '''
1632 1639
1633 1640 CODE = 'integrations_rtc' # Range Time Counts
1634 1641 colormap = 'BuGn'
1635 1642 plot_type = 'pcolorbuffer'
1636 1643
1637 1644 def setup(self):
1638 1645 self.xaxis = 'time'
1639 1646 self.ncols = 1
1640 1647 self.nrows = self.data.shape('integrations_rtc')[0]
1641 1648 self.nplots = self.nrows
1642 1649 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1643 1650
1644 1651
1645 1652 if not self.xlabel:
1646 1653 self.xlabel = 'Time'
1647 1654
1648 1655 self.ylabel = 'Height [km]'
1649 1656 if not self.titles:
1650 1657 self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)]
1651 1658
1652 1659 def update(self, dataOut):
1653 1660
1654 1661 data = {}
1655 1662 data['integrations_rtc'] = dataOut.nIncohInt
1656 1663
1657 1664 meta = {}
1658 1665
1659 1666 return data, meta
1660 1667
1661 1668 def plot(self):
1662 1669 # self.data.normalize_heights()
1663 1670 self.x = self.data.times
1664 1671 self.y = self.data.yrange
1665 1672 self.z = self.data['integrations_rtc']
1666 1673
1667 1674 #self.z = numpy.ma.masked_invalid(self.z)
1668 1675
1669 1676 if self.decimation is None:
1670 1677 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1671 1678 else:
1672 1679 x, y, z = self.fill_gaps(*self.decimate())
1673 1680
1674 1681 for n, ax in enumerate(self.axes):
1675 1682
1676 1683 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1677 1684 self.z[n])
1678 1685 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1679 1686 self.z[n])
1680 1687 data = self.data[-1]
1681 1688 if ax.firsttime:
1682 1689 if self.zlimits is not None:
1683 1690 self.zmin, self.zmax = self.zlimits[n]
1684 1691
1685 1692 ax.plt = ax.pcolormesh(x, y, z[n].T,
1686 1693 vmin=self.zmin,
1687 1694 vmax=self.zmax,
1688 1695 cmap=self.cmaps[n]
1689 1696 )
1690 1697 if self.showprofile:
1691 1698 ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0]
1692 1699 self.pf_axes[n].set_xlabel('')
1693 1700 else:
1694 1701 if self.zlimits is not None:
1695 1702 self.zmin, self.zmax = self.zlimits[n]
1696 1703 # ax.collections.remove(ax.collections[0]) # error while running
1697 1704 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1698 1705 vmin=self.zmin,
1699 1706 vmax=self.zmax,
1700 1707 cmap=self.cmaps[n]
1701 1708 )
1702 1709 if self.showprofile:
1703 1710 ax.plot_profile.set_data(data['integrations_rtc'][n], self.y)
1704 1711 self.pf_axes[n].set_xlabel('')
1705 1712
1706 1713
1707 1714
1708 1715 class RTIMapPlot(Plot):
1709 1716 '''
1710 1717 Plot for RTI data
1711 1718
1712 1719 Example:
1713 1720
1714 1721 controllerObj = Project()
1715 1722 controllerObj.setup(id = '11', name='eej_proc', description=desc)
1716 1723 ##.......................................................................................
1717 1724 ##.......................................................................................
1718 1725 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', path=inPath, startDate='2023/05/24',endDate='2023/05/24',
1719 1726 startTime='12:00:00',endTime='12:45:59',walk=1,timezone='lt',margin_days=1,code = code,nCode = nCode,
1720 1727 nBaud = nBaud,nOsamp = nosamp,nChannels=nChannels,nFFT=NFFT,
1721 1728 syncronization=False,shiftChannels=0)
1722 1729
1723 1730 volts_proc = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
1724 1731
1725 1732 opObj01 = volts_proc.addOperation(name='Decoder', optype='other')
1726 1733 opObj01.addParameter(name='code', value=code, format='floatlist')
1727 1734 opObj01.addParameter(name='nCode', value=1, format='int')
1728 1735 opObj01.addParameter(name='nBaud', value=nBaud, format='int')
1729 1736 opObj01.addParameter(name='osamp', value=nosamp, format='int')
1730 1737
1731 1738 opObj12 = volts_proc.addOperation(name='selectHeights', optype='self')
1732 1739 opObj12.addParameter(name='minHei', value='90', format='float')
1733 1740 opObj12.addParameter(name='maxHei', value='150', format='float')
1734 1741
1735 1742 proc_spc = controllerObj.addProcUnit(datatype='SpectraProc', inputId=volts_proc.getId())
1736 1743 proc_spc.addParameter(name='nFFTPoints', value='8', format='int')
1737 1744
1738 1745 opObj11 = proc_spc.addOperation(name='IncohInt', optype='other')
1739 1746 opObj11.addParameter(name='n', value='1', format='int')
1740 1747
1741 1748 beamMapFile = "/home/japaza/Documents/AMISR_sky_mapper/UMET_beamcodes.csv"
1742 1749
1743 1750 opObj12 = proc_spc.addOperation(name='RTIMapPlot', optype='external')
1744 1751 opObj12.addParameter(name='selectedHeightsList', value='95, 100, 105, 110 ', format='int')
1745 1752 opObj12.addParameter(name='bField', value='100', format='int')
1746 1753 opObj12.addParameter(name='filename', value=beamMapFile, format='str')
1747 1754
1748 1755 '''
1749 1756
1750 1757 CODE = 'rti_skymap'
1751 1758
1752 1759 plot_type = 'scatter'
1753 1760 titles = None
1754 1761 colormap = 'jet'
1755 1762 channelList = []
1756 1763 elevationList = []
1757 1764 azimuthList = []
1758 1765 last_noise = None
1759 1766 flag_setIndex = False
1760 1767 heights = []
1761 1768 dcosx = []
1762 1769 dcosy = []
1763 1770 fullDcosy = None
1764 1771 fullDcosy = None
1765 1772 hindex = []
1766 1773 mapFile = False
1767 1774 ##### BField ####
1768 1775 flagBField = False
1769 1776 dcosxB = []
1770 1777 dcosyB = []
1771 1778 Bmarker = ['+','*','D','x','s','>','o','^']
1772 1779
1773 1780
1774 1781 def setup(self):
1775 1782
1776 1783 self.xaxis = 'Range (Km)'
1777 1784 if len(self.selectedHeightsList) > 0:
1778 1785 self.nplots = len(self.selectedHeightsList)
1779 1786 else:
1780 1787 self.nplots = 4
1781 1788 self.ncols = int(numpy.ceil(self.nplots/2))
1782 1789 self.nrows = int(numpy.ceil(self.nplots/self.ncols))
1783 1790 self.ylabel = 'dcosy'
1784 1791 self.xlabel = 'dcosx'
1785 1792 self.colorbar = True
1786 1793 self.width = 6 + 4.1*self.nrows
1787 1794 self.height = 3 + 3.5*self.ncols
1788 1795
1789 1796
1790 1797 if self.extFile!=None:
1791 1798 try:
1792 1799 pointings = numpy.genfromtxt(self.extFile, delimiter=',')
1793 1800 full_azi = pointings[:,1]
1794 1801 full_elev = pointings[:,2]
1795 1802 self.fullDcosx = numpy.cos(numpy.radians(full_elev))*numpy.sin(numpy.radians(full_azi))
1796 1803 self.fullDcosy = numpy.cos(numpy.radians(full_elev))*numpy.cos(numpy.radians(full_azi))
1797 1804 mapFile = True
1798 1805 except Exception as e:
1799 1806 self.extFile = None
1800 1807 print(e)
1801 1808
1802 1809
1803 1810 def update_list(self,dataOut):
1804 1811 if len(self.channelList) == 0:
1805 1812 self.channelList = dataOut.channelList
1806 1813 if len(self.elevationList) == 0:
1807 1814 self.elevationList = dataOut.elevationList
1808 1815 if len(self.azimuthList) == 0:
1809 1816 self.azimuthList = dataOut.azimuthList
1810 1817 a = numpy.radians(numpy.asarray(self.azimuthList))
1811 1818 e = numpy.radians(numpy.asarray(self.elevationList))
1812 1819 self.heights = dataOut.heightList
1813 1820 self.dcosx = numpy.cos(e)*numpy.sin(a)
1814 1821 self.dcosy = numpy.cos(e)*numpy.cos(a)
1815 1822
1816 1823 if len(self.bFieldList)>0:
1817 1824 datetObj = datetime.datetime.fromtimestamp(dataOut.utctime)
1818 1825 doy = datetObj.timetuple().tm_yday
1819 1826 year = datetObj.year
1820 1827 # self.dcosxB, self.dcosyB
1821 1828 ObjB = BField(year=year,doy=doy,site=2,heights=self.bFieldList)
1822 1829 [dcos, alpha, nlon, nlat] = ObjB.getBField()
1823 1830
1824 1831 alpha_location = numpy.zeros((nlon,2,len(self.bFieldList)))
1825 1832 for ih in range(len(self.bFieldList)):
1826 1833 alpha_location[:,0,ih] = dcos[:,0,ih,0]
1827 1834 for ilon in numpy.arange(nlon):
1828 1835 myx = (alpha[ilon,:,ih])[::-1]
1829 1836 myy = (dcos[ilon,:,ih,0])[::-1]
1830 1837 tck = splrep(myx,myy,s=0)
1831 1838 mydcosx = splev(ObjB.alpha_i,tck,der=0)
1832 1839
1833 1840 myx = (alpha[ilon,:,ih])[::-1]
1834 1841 myy = (dcos[ilon,:,ih,1])[::-1]
1835 1842 tck = splrep(myx,myy,s=0)
1836 1843 mydcosy = splev(ObjB.alpha_i,tck,der=0)
1837 1844 alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy])
1838 1845 self.dcosxB.append(alpha_location[:,0,ih])
1839 1846 self.dcosyB.append(alpha_location[:,1,ih])
1840 1847 self.flagBField = True
1841 1848
1842 1849 if len(self.celestialList)>0:
1843 1850 #getBField(self.bFieldList, date)
1844 1851 #pass = kwargs.get('celestial', [])
1845 1852 pass
1846 1853
1847 1854
1848 1855 def update(self, dataOut):
1849 1856
1850 1857 if len(self.channelList) == 0:
1851 1858 self.update_list(dataOut)
1852 1859
1853 1860 if not self.flag_setIndex:
1854 1861 if len(self.selectedHeightsList)>0:
1855 1862 for sel_height in self.selectedHeightsList:
1856 1863 index_list = numpy.where(self.heights >= sel_height)
1857 1864 index_list = index_list[0]
1858 1865 self.hindex.append(index_list[0])
1859 1866 self.flag_setIndex = True
1860 1867
1861 1868 data = {}
1862 1869 meta = {}
1863 1870
1864 1871 data['rti_skymap'] = dataOut.getPower()
1865 1872 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1866 1873 noise = 10*numpy.log10(dataOut.getNoise()/norm)
1867 1874 data['noise'] = noise
1868 1875
1869 1876 return data, meta
1870 1877
1871 1878 def plot(self):
1872 1879
1873 1880 self.x = self.dcosx
1874 1881 self.y = self.dcosy
1875 1882 self.z = self.data[-1]['rti_skymap']
1876 1883 self.z = numpy.array(self.z, dtype=float)
1877 1884
1878 1885 if len(self.hindex) > 0:
1879 1886 index = self.hindex
1880 1887 else:
1881 1888 index = numpy.arange(0, len(self.heights), int((len(self.heights))/4.2))
1882 1889
1883 1890 self.titles = ['Height {:.2f} km '.format(self.heights[i])+" " for i in index]
1884 1891 for n, ax in enumerate(self.axes):
1885 1892
1886 1893 if ax.firsttime:
1887 1894
1888 1895 self.xmax = self.xmax if self.xmax else numpy.nanmax(self.x)
1889 1896 self.xmin = self.xmin if self.xmin else numpy.nanmin(self.x)
1890 1897 self.ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
1891 1898 self.ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
1892 1899 self.zmax = self.zmax if self.zmax else numpy.nanmax(self.z)
1893 1900 self.zmin = self.zmin if self.zmin else numpy.nanmin(self.z)
1894 1901
1895 1902 if self.extFile!=None:
1896 1903 ax.scatter(self.fullDcosx, self.fullDcosy, marker="+", s=20)
1897 1904
1898 1905 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1899 1906 s=60, marker="s", vmax = self.zmax)
1900 1907
1901 1908
1902 1909 ax.minorticks_on()
1903 1910 ax.grid(which='major', axis='both')
1904 1911 ax.grid(which='minor', axis='x')
1905 1912
1906 1913 if self.flagBField :
1907 1914
1908 1915 for ih in range(len(self.bFieldList)):
1909 1916 label = str(self.bFieldList[ih]) + ' km'
1910 1917 ax.plot(self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1911 1918 label=label, linestyle='--', ms=4.0,lw=0.5)
1912 1919 handles, labels = ax.get_legend_handles_labels()
1913 1920 a = -0.05
1914 1921 b = 1.15 - 1.19*(self.nrows)
1915 1922 self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥')
1916 1923
1917 1924 else:
1918 1925
1919 1926 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1920 1927 s=80, marker="s", vmax = self.zmax)
1921 1928
1922 1929 if self.flagBField :
1923 1930 for ih in range(len(self.bFieldList)):
1924 1931 ax.plot (self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1925 1932 linestyle='--', ms=4.0,lw=0.5)
1926 1933
1927 1934
1928 1935
Show file before | Show file after | Hide comments
@@ -1,815 +1,819
1 1 import os
2 2 import time
3 3 import datetime
4 4
5 5 import numpy
6 6 import h5py
7 7
8 8 import schainpy.admin
9 9 from schainpy.model.data.jrodata import *
10 10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 11 from schainpy.model.io.jroIO_base import *
12 12 from schainpy.utils import log
13 13
14 14
15 15 class HDFReader(Reader, ProcessingUnit):
16 16 """Processing unit to read HDF5 format files
17 17
18 18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 20 attributes.
21 21 It is possible to read any HDF5 file by given the structure in the `description`
22 22 parameter, also you can add extra values to metadata with the parameter `extras`.
23 23
24 24 Parameters:
25 25 -----------
26 26 path : str
27 27 Path where files are located.
28 28 startDate : date
29 29 Start date of the files
30 30 endDate : list
31 31 End date of the files
32 32 startTime : time
33 33 Start time of the files
34 34 endTime : time
35 35 End time of the files
36 36 description : dict, optional
37 37 Dictionary with the description of the HDF5 file
38 38 extras : dict, optional
39 39 Dictionary with extra metadata to be be added to `dataOut`
40 40
41 41 Attention: Be carefull, add attribute utcoffset, in the last part of reader in order to work in Local Time without time problems.
42 42
43 43 -----------
44 44 utcoffset='-18000'
45 45
46 46
47 47 Examples
48 48 --------
49 49
50 50 desc = {
51 51 'Data': {
52 52 'data_output': ['u', 'v', 'w'],
53 53 'utctime': 'timestamps',
54 54 } ,
55 55 'Metadata': {
56 56 'heightList': 'heights'
57 57 }
58 58 }
59 59
60 60 desc = {
61 61 'Data': {
62 62 'data_output': 'winds',
63 63 'utctime': 'timestamps'
64 64 },
65 65 'Metadata': {
66 66 'heightList': 'heights'
67 67 }
68 68 }
69 69
70 70 extras = {
71 71 'timeZone': 300
72 72 }
73 73
74 74 reader = project.addReadUnit(
75 75 name='HDFReader',
76 76 path='/path/to/files',
77 77 startDate='2019/01/01',
78 78 endDate='2019/01/31',
79 79 startTime='00:00:00',
80 80 endTime='23:59:59',
81 81 utcoffset='-18000'
82 82 # description=json.dumps(desc),
83 83 # extras=json.dumps(extras),
84 84 )
85 85
86 86 """
87 87
88 88 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
89 89
90 90 def __init__(self):
91 91
92 92 ProcessingUnit.__init__(self)
93 93 self.ext = ".hdf5"
94 94 self.optchar = "D"
95 95 self.meta = {}
96 96 self.data = {}
97 97 self.open_file = h5py.File
98 98 self.open_mode = 'r'
99 99 self.description = {}
100 100 self.extras = {}
101 101 self.filefmt = "*%Y%j***"
102 102 self.folderfmt = "*%Y%j"
103 103 self.utcoffset = 0
104 104 self.flagUpdateDataOut = False
105 105 self.dataOut = Parameters()
106 106 self.dataOut.error=False ## NOTE: Importante definir esto antes inicio
107 107 self.dataOut.flagNoData = True
108 108
109 109 def setup(self, **kwargs):
110 110
111 111 self.set_kwargs(**kwargs)
112 112 if not self.ext.startswith('.'):
113 113 self.ext = '.{}'.format(self.ext)
114 114
115 115 if self.online:
116 116 log.log("Searching files in online mode...", self.name)
117 117
118 118 for nTries in range(self.nTries):
119 119 fullpath = self.searchFilesOnLine(self.path, self.startDate,
120 120 self.endDate, self.expLabel, self.ext, self.walk,
121 121 self.filefmt, self.folderfmt)
122 122 pathname, filename = os.path.split(fullpath)
123 123 try:
124 124 fullpath = next(fullpath)
125 125 except:
126 126 fullpath = None
127 127
128 128 if fullpath:
129 129 break
130 130
131 131 log.warning(
132 132 'Waiting {} sec for a valid file in {}: try {} ...'.format(
133 133 self.delay, self.path, nTries + 1),
134 134 self.name)
135 135 time.sleep(self.delay)
136 136
137 137 if not(fullpath):
138 138 raise schainpy.admin.SchainError(
139 139 'There isn\'t any valid file in {}'.format(self.path))
140 140
141 141 pathname, filename = os.path.split(fullpath)
142 142 self.year = int(filename[1:5])
143 143 self.doy = int(filename[5:8])
144 144 self.set = int(filename[8:11]) - 1
145 145 else:
146 146 log.log("Searching files in {}".format(self.path), self.name)
147 147 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
148 148 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
149 149
150 150 self.setNextFile()
151 151
152 152 return
153 153
154 154 # def readFirstHeader(self):
155 155 # '''Read metadata and data'''
156 156
157 157 # self.__readMetadata()
158 158 # self.__readData()
159 159 # self.__setBlockList()
160 160
161 161 # if 'type' in self.meta:
162 162 # self.dataOut = eval(self.meta['type'])()
163 163
164 164 # for attr in self.meta:
165 165 # setattr(self.dataOut, attr, self.meta[attr])
166 166
167 167 # self.blockIndex = 0
168 168
169 169 # return
170 170
171 171 def readFirstHeader(self):
172 172 '''Read metadata and data'''
173 173
174 174 self.__readMetadata2()
175 175 self.__readData()
176 176 self.__setBlockList()
177 177 if 'type' in self.meta:
178 178 self.dataOut = eval(self.meta['type'])()
179 179
180 180 for attr in self.meta:
181 181 if "processingHeaderObj" in attr:
182 182 self.flagUpdateDataOut=True
183 183 at = attr.split('.')
184 184 if len(at) > 1:
185 185 setattr(eval("self.dataOut."+at[0]),at[1], self.meta[attr])
186 186 else:
187 187 setattr(self.dataOut, attr, self.meta[attr])
188 188 self.blockIndex = 0
189 189
190 190 if self.flagUpdateDataOut:
191 191 self.updateDataOut()
192 192
193 193 return
194 194
195 195 def updateDataOut(self):
196 196
197 197 self.dataOut.azimuthList = self.dataOut.processingHeaderObj.azimuthList
198 198 self.dataOut.elevationList = self.dataOut.processingHeaderObj.elevationList
199 199 self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
200 200 self.dataOut.ippSeconds = self.dataOut.processingHeaderObj.ipp
201 201 self.dataOut.elevationList = self.dataOut.processingHeaderObj.elevationList
202 202 self.dataOut.channelList = self.dataOut.processingHeaderObj.channelList
203 203 self.dataOut.nCohInt = self.dataOut.processingHeaderObj.nCohInt
204 204 self.dataOut.nFFTPoints = self.dataOut.processingHeaderObj.nFFTPoints
205 205 self.flagUpdateDataOut = False
206 206 self.dataOut.frequency = self.dataOut.radarControllerHeaderObj.frequency
207 207 #self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
208 208
209 209 def __setBlockList(self):
210 210 '''
211 211 Selects the data within the times defined
212 212
213 213 self.fp
214 214 self.startTime
215 215 self.endTime
216 216 self.blockList
217 217 self.blocksPerFile
218 218
219 219 '''
220 220
221 221 startTime = self.startTime
222 222 endTime = self.endTime
223 223 thisUtcTime = self.data['utctime'] + self.utcoffset
224 224 # self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
225 225 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
226 226 self.startFileDatetime = thisDatetime
227 227 thisDate = thisDatetime.date()
228 228 thisTime = thisDatetime.time()
229 229 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
230 230 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
231 231 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
232 232
233 233 self.blockList = ind
234 234 self.blocksPerFile = len(ind)
235 235 # self.blocksPerFile = len(thisUtcTime)
236 236 if len(ind)==0:
237 237 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.blockIndex,
238 238 self.blocksPerFile,
239 239 thisDatetime))
240 240 self.setNextFile()
241 241
242 242 return
243 243
244 244 def __readMetadata(self):
245 245 '''
246 246 Reads Metadata
247 247 '''
248 248
249 249 meta = {}
250 250
251 251 if self.description:
252 252 for key, value in self.description['Metadata'].items():
253 253 meta[key] = self.fp[value][()]
254 254 else:
255 255 grp = self.fp['Metadata']
256 256 for name in grp:
257 257 meta[name] = grp[name][()]
258 258
259 259 if self.extras:
260 260 for key, value in self.extras.items():
261 261 meta[key] = value
262 262 self.meta = meta
263 263
264 264 return
265 265
266 266 def __readMetadata2(self):
267 267 '''
268 268 Reads Metadata
269 269 '''
270 270 meta = {}
271 271 if self.description:
272 272 for key, value in self.description['Metadata'].items():
273 273 meta[key] = self.fp[value][()]
274 274 else:
275 275 grp = self.fp['Metadata']
276 276 for item in grp.values():
277 277 name = item.name
278 278 if isinstance(item, h5py.Dataset):
279 279 name = name.split("/")[-1]
280 280 meta[name] = item[()]
281 281 else:
282 282 grp2 = self.fp[name]
283 283 Obj = name.split("/")[-1]
284 284
285 285 for item2 in grp2.values():
286 286 name2 = Obj+"."+item2.name.split("/")[-1]
287 287 meta[name2] = item2[()]
288 288
289 289 if self.extras:
290 290 for key, value in self.extras.items():
291 291 meta[key] = value
292 292 self.meta = meta
293 293
294 294 return
295 295
296 296 def __readData(self):
297 297
298 298 data = {}
299 299
300 300 if self.description:
301 301 for key, value in self.description['Data'].items():
302 302 if isinstance(value, str):
303 303 if isinstance(self.fp[value], h5py.Dataset):
304 304 data[key] = self.fp[value][()]
305 305 elif isinstance(self.fp[value], h5py.Group):
306 306 array = []
307 307 for ch in self.fp[value]:
308 308 array.append(self.fp[value][ch][()])
309 309 data[key] = numpy.array(array)
310 310 elif isinstance(value, list):
311 311 array = []
312 312 for ch in value:
313 313 array.append(self.fp[ch][()])
314 314 data[key] = numpy.array(array)
315 315 else:
316 316 grp = self.fp['Data']
317 317 for name in grp:
318 318 if isinstance(grp[name], h5py.Dataset):
319 319 array = grp[name][()]
320 320 elif isinstance(grp[name], h5py.Group):
321 321 array = []
322 322 for ch in grp[name]:
323 323 array.append(grp[name][ch][()])
324 324 array = numpy.array(array)
325 325 else:
326 326 log.warning('Unknown type: {}'.format(name))
327 327
328 328 if name in self.description:
329 329 key = self.description[name]
330 330 else:
331 331 key = name
332 332 data[key] = array
333 333
334 334 self.data = data
335 335 return
336 336
337 337 def getData(self):
338 338
339 339 if not self.isDateTimeInRange(self.startFileDatetime, self.startDate, self.endDate, self.startTime, self.endTime):
340 340 self.dataOut.flagNoData = True
341 341 self.blockIndex = self.blocksPerFile
342 342 self.dataOut.error = True # TERMINA EL PROGRAMA
343 343 return
344 344 for attr in self.data:
345 345
346 346 if self.data[attr].ndim == 1:
347 347 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
348 348 else:
349 349 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
350 350
351 351
352 352 self.blockIndex += 1
353 353
354 354 if self.blockIndex == 1:
355 355 log.log("Block No. {}/{} -> {}".format(
356 356 self.blockIndex,
357 357 self.blocksPerFile,
358 358 self.dataOut.datatime.ctime()), self.name)
359 359 else:
360 360 log.log("Block No. {}/{} ".format(
361 361 self.blockIndex,
362 362 self.blocksPerFile),self.name)
363 363
364 364 if self.blockIndex == self.blocksPerFile:
365 365 self.setNextFile()
366 366
367 367 self.dataOut.flagNoData = False
368 368
369 369 return
370 370
371 371 def run(self, **kwargs):
372 372
373 373 if not(self.isConfig):
374 374 self.setup(**kwargs)
375 375 self.isConfig = True
376 376
377 377 if self.blockIndex == self.blocksPerFile:
378 378 self.setNextFile()
379 379
380 380 self.getData()
381 381
382 382 return
383 383
384 384 @MPDecorator
385 385 class HDFWriter(Operation):
386 386 """Operation to write HDF5 files.
387 387
388 388 The HDF5 file contains by default two groups Data and Metadata where
389 389 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
390 390 parameters, data attributes are normaly time dependent where the metadata
391 391 are not.
392 392 It is possible to customize the structure of the HDF5 file with the
393 393 optional description parameter see the examples.
394 394
395 395 Parameters:
396 396 -----------
397 397 path : str
398 398 Path where files will be saved.
399 399 blocksPerFile : int
400 400 Number of blocks per file
401 401 metadataList : list
402 402 List of the dataOut attributes that will be saved as metadata
403 403 dataList : int
404 404 List of the dataOut attributes that will be saved as data
405 405 setType : bool
406 406 If True the name of the files corresponds to the timestamp of the data
407 407 description : dict, optional
408 408 Dictionary with the desired description of the HDF5 file
409 409
410 410 Examples
411 411 --------
412 412
413 413 desc = {
414 414 'data_output': {'winds': ['z', 'w', 'v']},
415 415 'utctime': 'timestamps',
416 416 'heightList': 'heights'
417 417 }
418 418 desc = {
419 419 'data_output': ['z', 'w', 'v'],
420 420 'utctime': 'timestamps',
421 421 'heightList': 'heights'
422 422 }
423 423 desc = {
424 424 'Data': {
425 425 'data_output': 'winds',
426 426 'utctime': 'timestamps'
427 427 },
428 428 'Metadata': {
429 429 'heightList': 'heights'
430 430 }
431 431 }
432 432
433 433 writer = proc_unit.addOperation(name='HDFWriter')
434 434 writer.addParameter(name='path', value='/path/to/file')
435 435 writer.addParameter(name='blocksPerFile', value='32')
436 436 writer.addParameter(name='metadataList', value='heightList,timeZone')
437 437 writer.addParameter(name='dataList',value='data_output,utctime')
438 438 # writer.addParameter(name='description',value=json.dumps(desc))
439 439
440 440 """
441 441
442 442 ext = ".hdf5"
443 443 optchar = "D"
444 444 filename = None
445 445 path = None
446 446 setFile = None
447 447 fp = None
448 ds = None
448 449 firsttime = True
449 450 #Configurations
450 451 blocksPerFile = None
451 452 blockIndex = None
452 453 dataOut = None #eval ??????
453 454 #Data Arrays
454 455 dataList = None
455 456 metadataList = None
456 457 currentDay = None
457 458 lastTime = None
458 459 timeZone = "ut"
459 460 hourLimit = 3
460 461 breakDays = True
461 462
462 463 def __init__(self):
463 464
464 465 Operation.__init__(self)
465 466 return
466 467
467 468 def set_kwargs(self, **kwargs):
468 469
469 470 for key, value in kwargs.items():
470 471 setattr(self, key, value)
471 472
472 473 def set_kwargs_obj(self, obj, **kwargs):
473 474
474 475 for key, value in kwargs.items():
475 476 setattr(obj, key, value)
476 477
477 478 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None,
478 479 description={},timeZone = "ut",hourLimit = 3, breakDays=True, **kwargs):
479 480 self.path = path
480 481 self.blocksPerFile = blocksPerFile
481 482 self.metadataList = metadataList
482 483 self.dataList = [s.strip() for s in dataList]
483 484 self.setType = setType
484 485 self.description = description
485 486 self.timeZone = timeZone
486 487 self.hourLimit = hourLimit
487 488 self.breakDays = breakDays
488 489 self.set_kwargs(**kwargs)
489 490
490 491 if self.metadataList is None:
491 492 self.metadataList = self.dataOut.metadata_list
492 493
494 self.metadataList = list(set(self.metadataList))
495
493 496 tableList = []
494 497 dsList = []
495 498
496 499 for i in range(len(self.dataList)):
497 500 dsDict = {}
498 501 if hasattr(self.dataOut, self.dataList[i]):
499 502 dataAux = getattr(self.dataOut, self.dataList[i])
500 503 dsDict['variable'] = self.dataList[i]
501 504 else:
502 505 log.warning('Attribute {} not found in dataOut'.format(self.dataList[i]),self.name)
503 506 continue
504 507
505 508 if dataAux is None:
506 509 continue
507 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
510 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float_)):
508 511 dsDict['nDim'] = 0
509 512 else:
510 513 dsDict['nDim'] = len(dataAux.shape)
511 514 dsDict['shape'] = dataAux.shape
512 515 dsDict['dsNumber'] = dataAux.shape[0]
513 516 dsDict['dtype'] = dataAux.dtype
514 517
515 518 dsList.append(dsDict)
516 519
520 self.blockIndex = 0
517 521 self.dsList = dsList
518 522 self.currentDay = self.dataOut.datatime.date()
519 523
520 524 def timeFlag(self):
521 525 currentTime = self.dataOut.utctime
522 526 timeTuple = None
523 527 if self.timeZone == "lt":
524 528 timeTuple = time.localtime(currentTime)
525 529 else :
526 530 timeTuple = time.gmtime(currentTime)
527 531 dataDay = timeTuple.tm_yday
528 532
529 533 if self.lastTime is None:
530 534 self.lastTime = currentTime
531 535 self.currentDay = dataDay
532 536 return False
533 537
534 538 timeDiff = currentTime - self.lastTime
535 539
536 540 # Si el dia es diferente o si la diferencia entre un
537 541 # dato y otro supera self.hourLimit
538 542 if (dataDay != self.currentDay) and self.breakDays:
539 543 self.currentDay = dataDay
540 544 return True
541 545 elif timeDiff > self.hourLimit*60*60:
542 546 self.lastTime = currentTime
543 547 return True
544 548 else:
545 549 self.lastTime = currentTime
546 550 return False
547 551
548 552 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
549 553 dataList=[], setType=None, description={}, **kwargs):
550 554
551 555 self.dataOut = dataOut
552 556 self.set_kwargs_obj(self.dataOut, **kwargs)
553 557 if not(self.isConfig):
554 558 self.setup(path=path, blocksPerFile=blocksPerFile,
555 559 metadataList=metadataList, dataList=dataList,
556 560 setType=setType, description=description, **kwargs)
557 561
558 562 self.isConfig = True
559 563 self.setNextFile()
560 564
561 565 self.putData()
562 566 return
563 567
564 568 def setNextFile(self):
565 569
566 570 ext = self.ext
567 571 path = self.path
568 572 setFile = self.setFile
569 573 timeTuple = None
570 574 if self.timeZone == "lt":
571 575 timeTuple = time.localtime(self.dataOut.utctime)
572 576 elif self.timeZone == "ut":
573 577 timeTuple = time.gmtime(self.dataOut.utctime)
574 578 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
575 579 fullpath = os.path.join(path, subfolder)
576 580
577 581 if os.path.exists(fullpath):
578 582 filesList = os.listdir(fullpath)
579 583 filesList = [k for k in filesList if k.startswith(self.optchar)]
580 584 if len(filesList) > 0:
581 585 filesList = sorted(filesList, key=str.lower)
582 586 filen = filesList[-1]
583 587 # el filename debera tener el siguiente formato
584 588 # 0 1234 567 89A BCDE (hex)
585 589 # x YYYY DDD SSS .ext
586 590 if isNumber(filen[8:11]):
587 591 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
588 592 else:
589 593 setFile = -1
590 594 else:
591 595 setFile = -1 #inicializo mi contador de seteo
592 596 else:
593 597 os.makedirs(fullpath)
594 598 setFile = -1 #inicializo mi contador de seteo
595 599
596 600 if self.setType is None:
597 601 setFile += 1
598 602 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
599 603 timeTuple.tm_year,
600 604 timeTuple.tm_yday,
601 605 setFile,
602 606 ext)
603 607 else:
604 608 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
605 609 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
606 610 timeTuple.tm_year,
607 611 timeTuple.tm_yday,
608 612 setFile,
609 613 ext)
610 614
611 615 self.filename = os.path.join(path, subfolder, file)
612 616
613 617
614 618
615 619 def getLabel(self, name, x=None):
616 620
617 621 if x is None:
618 622 if 'Data' in self.description:
619 623 data = self.description['Data']
620 624 if 'Metadata' in self.description:
621 625 data.update(self.description['Metadata'])
622 626 else:
623 627 data = self.description
624 628 if name in data:
625 629 if isinstance(data[name], str):
626 630 return data[name]
627 631 elif isinstance(data[name], list):
628 632 return None
629 633 elif isinstance(data[name], dict):
630 634 for key, value in data[name].items():
631 635 return key
632 636 return name
633 637 else:
634 638 if 'Metadata' in self.description:
635 639 meta = self.description['Metadata']
636 640 else:
637 641 meta = self.description
638 642 if name in meta:
639 643 if isinstance(meta[name], list):
640 644 return meta[name][x]
641 645 elif isinstance(meta[name], dict):
642 646 for key, value in meta[name].items():
643 647 return value[x]
644 648 if 'cspc' in name:
645 649 return 'pair{:02d}'.format(x)
646 650 else:
647 651 return 'channel{:02d}'.format(x)
648 652
649 653 def writeMetadata(self, fp):
650 654
651 655 if self.description:
652 656 if 'Metadata' in self.description:
653 657 grp = fp.create_group('Metadata')
654 658 else:
655 659 grp = fp
656 660 else:
657 661 grp = fp.create_group('Metadata')
658 662
659 663 for i in range(len(self.metadataList)):
660 664 if not hasattr(self.dataOut, self.metadataList[i]):
661 665 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
662 666 continue
663 667 value = getattr(self.dataOut, self.metadataList[i])
664 668 if isinstance(value, bool):
665 669 if value is True:
666 670 value = 1
667 671 else:
668 672 value = 0
669 673 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
670 674 return
671 675
672 676 def writeMetadata2(self, fp):
673 677
674 678 if self.description:
675 679 if 'Metadata' in self.description:
676 680 grp = fp.create_group('Metadata')
677 681 else:
678 682 grp = fp
679 683 else:
680 684 grp = fp.create_group('Metadata')
681 685
682 686 for i in range(len(self.metadataList)):
683 687
684 688 attribute = self.metadataList[i]
685 689 attr = attribute.split('.')
686 690 if len(attr) > 1:
687 691 if not hasattr(eval("self.dataOut."+attr[0]),attr[1]):
688 692 log.warning('Metadata: {}.{} not found'.format(attr[0],attr[1]), self.name)
689 693 continue
690 694 value = getattr(eval("self.dataOut."+attr[0]),attr[1])
691 695 if isinstance(value, bool):
692 696 if value is True:
693 697 value = 1
694 698 else:
695 699 value = 0
696 700 if isinstance(value,type(None)):
697 701 log.warning("Invalid value detected, {} is None".format(attribute), self.name)
698 702 value = 0
699 703 grp2 = None
700 704 if not 'Metadata/'+attr[0] in fp:
701 705 grp2 = fp.create_group('Metadata/'+attr[0])
702 706 else:
703 707 grp2 = fp['Metadata/'+attr[0]]
704 708 grp2.create_dataset(attr[1], data=value)
705 709
706 710 else:
707 711 if not hasattr(self.dataOut, attr[0] ):
708 712 log.warning('Metadata: `{}` not found'.format(attribute), self.name)
709 713 continue
710 714 value = getattr(self.dataOut, attr[0])
711 715 if isinstance(value, bool):
712 716 if value is True:
713 717 value = 1
714 718 else:
715 719 value = 0
716 720 if isinstance(value, type(None)):
717 721 log.error("Value {} is None".format(attribute),self.name)
718 722
719 723 grp.create_dataset(self.getLabel(attribute), data=value)
720 724
721 725 return
722 726
723 727 def writeData(self, fp):
724 728
725 729 if self.description:
726 730 if 'Data' in self.description:
727 731 grp = fp.create_group('Data')
728 732 else:
729 733 grp = fp
730 734 else:
731 735 grp = fp.create_group('Data')
732 736
733 737 dtsets = []
734 738 data = []
735 739
736 740 for dsInfo in self.dsList:
737 741 if dsInfo['nDim'] == 0:
738 742 ds = grp.create_dataset(
739 743 self.getLabel(dsInfo['variable']),
740 744 (self.blocksPerFile,),
741 745 chunks=True,
742 746 dtype=numpy.float64)
743 747 dtsets.append(ds)
744 748 data.append((dsInfo['variable'], -1))
745 749 else:
746 750 label = self.getLabel(dsInfo['variable'])
747 751 if label is not None:
748 752 sgrp = grp.create_group(label)
749 753 else:
750 754 sgrp = grp
751 755 for i in range(dsInfo['dsNumber']):
752 756 ds = sgrp.create_dataset(
753 757 self.getLabel(dsInfo['variable'], i),
754 758 (self.blocksPerFile,) + dsInfo['shape'][1:],
755 759 chunks=True,
756 760 dtype=dsInfo['dtype'])
757 761 dtsets.append(ds)
758 762 data.append((dsInfo['variable'], i))
759 763 fp.flush()
760 764
761 765 log.log('Creating file: {}'.format(fp.filename), self.name)
762 766
763 767 self.ds = dtsets
764 768 self.data = data
765 769 self.firsttime = True
766 self.blockIndex = 0
770
767 771 return
768 772
769 773 def putData(self):
770 774
771 775 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
772 776 self.closeFile()
773 777 self.setNextFile()
774 778 self.dataOut.flagNoData = False
775 779 self.blockIndex = 0
776 780
777 781 if self.blockIndex == 0:
778 782 #Setting HDF5 File
779 783 self.fp = h5py.File(self.filename, 'w')
780 784 #write metadata
781 785 self.writeMetadata2(self.fp)
782 786 #Write data
783 787 self.writeData(self.fp)
784 788 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
785 789 elif (self.blockIndex % 10 ==0):
786 790 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
787 791 else:
788 792
789 793 log.log('Block No. {}/{}'.format(self.blockIndex+1, self.blocksPerFile), self.name)
790 794
791 795 for i, ds in enumerate(self.ds):
792 796 attr, ch = self.data[i]
793 797 if ch == -1:
794 798 ds[self.blockIndex] = getattr(self.dataOut, attr)
795 799 else:
796 800 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
797 801
798 802 self.blockIndex += 1
799 803
800 804 self.fp.flush()
801 805 self.dataOut.flagNoData = True
802 806
803 807 def closeFile(self):
804 808
805 809 if self.blockIndex != self.blocksPerFile:
806 810 for ds in self.ds:
807 811 ds.resize(self.blockIndex, axis=0)
808 812
809 813 if self.fp:
810 814 self.fp.flush()
811 815 self.fp.close()
812 816
813 817 def close(self):
814 818
815 819 self.closeFile()
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now