##// END OF EJS Templates
schain
RSS
Fix HDFReader, update metadata list for spectra
Juan C. Espinoza -
r1552:634b09aabe58
parent child
Show More
Side by Side Unified
Context file:
r1552:634b09aabe58 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1170 +1,1169
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Definition of diferent Data objects for different types of data
6 6
7 7 Here you will find the diferent data objects for the different types
8 8 of data, this data objects must be used as dataIn or dataOut objects in
9 9 processing units and operations. Currently the supported data objects are:
10 10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 11 """
12 12
13 13 import copy
14 14 import numpy
15 15 import datetime
16 16 import json
17 17
18 18 import schainpy.admin
19 19 from schainpy.utils import log
20 20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 21 from schainpy.model.data import _noise
22 22
23 23
24 24 def getNumpyDtype(dataTypeCode):
25 25
26 26 if dataTypeCode == 0:
27 27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 28 elif dataTypeCode == 1:
29 29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 30 elif dataTypeCode == 2:
31 31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 32 elif dataTypeCode == 3:
33 33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 34 elif dataTypeCode == 4:
35 35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 36 elif dataTypeCode == 5:
37 37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 38 else:
39 39 raise ValueError('dataTypeCode was not defined')
40 40
41 41 return numpyDtype
42 42
43 43
44 44 def getDataTypeCode(numpyDtype):
45 45
46 46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 47 datatype = 0
48 48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 49 datatype = 1
50 50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 51 datatype = 2
52 52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 53 datatype = 3
54 54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 55 datatype = 4
56 56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 57 datatype = 5
58 58 else:
59 59 datatype = None
60 60
61 61 return datatype
62 62
63 63
64 64 def hildebrand_sekhon(data, navg):
65 65 """
66 66 This method is for the objective determination of the noise level in Doppler spectra. This
67 67 implementation technique is based on the fact that the standard deviation of the spectral
68 68 densities is equal to the mean spectral density for white Gaussian noise
69 69
70 70 Inputs:
71 71 Data : heights
72 72 navg : numbers of averages
73 73
74 74 Return:
75 75 mean : noise's level
76 76 """
77 77
78 78 sortdata = numpy.sort(data, axis=None)
79 79 #print(numpy.shape(data))
80 80 #exit()
81 81 '''
82 82 lenOfData = len(sortdata)
83 83 nums_min = lenOfData*0.2
84 84
85 85 if nums_min <= 5:
86 86
87 87 nums_min = 5
88 88
89 89 sump = 0.
90 90 sumq = 0.
91 91
92 92 j = 0
93 93 cont = 1
94 94
95 95 while((cont == 1)and(j < lenOfData)):
96 96
97 97 sump += sortdata[j]
98 98 sumq += sortdata[j]**2
99 99
100 100 if j > nums_min:
101 101 rtest = float(j)/(j-1) + 1.0/navg
102 102 if ((sumq*j) > (rtest*sump**2)):
103 103 j = j - 1
104 104 sump = sump - sortdata[j]
105 105 sumq = sumq - sortdata[j]**2
106 106 cont = 0
107 107
108 108 j += 1
109 109
110 110 lnoise = sump / j
111 111 '''
112 112 return _noise.hildebrand_sekhon(sortdata, navg)
113 113
114 114
115 115 class Beam:
116 116
117 117 def __init__(self):
118 118 self.codeList = []
119 119 self.azimuthList = []
120 120 self.zenithList = []
121 121
122 122
123 123 class GenericData(object):
124 124
125 125 flagNoData = True
126 126
127 127 def copy(self, inputObj=None):
128 128
129 129 if inputObj == None:
130 130 return copy.deepcopy(self)
131 131
132 132 for key in list(inputObj.__dict__.keys()):
133 133
134 134 attribute = inputObj.__dict__[key]
135 135
136 136 # If this attribute is a tuple or list
137 137 if type(inputObj.__dict__[key]) in (tuple, list):
138 138 self.__dict__[key] = attribute[:]
139 139 continue
140 140
141 141 # If this attribute is another object or instance
142 142 if hasattr(attribute, '__dict__'):
143 143 self.__dict__[key] = attribute.copy()
144 144 continue
145 145
146 146 self.__dict__[key] = inputObj.__dict__[key]
147 147
148 148 def deepcopy(self):
149 149
150 150 return copy.deepcopy(self)
151 151
152 152 def isEmpty(self):
153 153
154 154 return self.flagNoData
155 155
156 156 def isReady(self):
157 157
158 158 return not self.flagNoData
159 159
160 160
161 161 class JROData(GenericData):
162 162
163 163 systemHeaderObj = SystemHeader()
164 164 radarControllerHeaderObj = RadarControllerHeader()
165 165 type = None
166 166 datatype = None # dtype but in string
167 167 nProfiles = None
168 168 heightList = None
169 169 channelList = None
170 170 flagDiscontinuousBlock = False
171 171 useLocalTime = False
172 172 utctime = None
173 173 timeZone = None
174 174 dstFlag = None
175 175 errorCount = None
176 176 blocksize = None
177 177 flagDecodeData = False # asumo q la data no esta decodificada
178 178 flagDeflipData = False # asumo q la data no esta sin flip
179 179 flagShiftFFT = False
180 180 nCohInt = None
181 181 windowOfFilter = 1
182 182 C = 3e8
183 183 frequency = 49.92e6
184 184 realtime = False
185 185 beacon_heiIndexList = None
186 186 last_block = None
187 187 blocknow = None
188 188 azimuth = None
189 189 zenith = None
190 190 beam = Beam()
191 191 profileIndex = None
192 192 error = None
193 193 data = None
194 194 nmodes = None
195 195 metadata_list = ['heightList', 'timeZone', 'type']
196 196
197 197 def __str__(self):
198 198
199 return '{} - {}'.format(self.type, self.datatime())
199 return '{} - {}'.format(self.type, self.datatime)
200 200
201 201 def getNoise(self):
202 202
203 203 raise NotImplementedError
204 204
205 205 @property
206 206 def nChannels(self):
207 207
208 208 return len(self.channelList)
209 209
210 210 @property
211 211 def channelIndexList(self):
212 212
213 213 return list(range(self.nChannels))
214 214
215 215 @property
216 216 def nHeights(self):
217 217
218 218 return len(self.heightList)
219 219
220 220 def getDeltaH(self):
221 221
222 222 return self.heightList[1] - self.heightList[0]
223 223
224 224 @property
225 225 def ltctime(self):
226 226
227 227 if self.useLocalTime:
228 228 return self.utctime - self.timeZone * 60
229 229
230 230 return self.utctime
231 231
232 232 @property
233 233 def datatime(self):
234 234
235 235 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
236 236 return datatimeValue
237 237
238 238 def getTimeRange(self):
239 239
240 240 datatime = []
241 241
242 242 datatime.append(self.ltctime)
243 243 datatime.append(self.ltctime + self.timeInterval + 1)
244 244
245 245 datatime = numpy.array(datatime)
246 246
247 247 return datatime
248 248
249 249 def getFmaxTimeResponse(self):
250 250
251 251 period = (10 ** -6) * self.getDeltaH() / (0.15)
252 252
253 253 PRF = 1. / (period * self.nCohInt)
254 254
255 255 fmax = PRF
256 256
257 257 return fmax
258 258
259 259 def getFmax(self):
260 260 PRF = 1. / (self.ippSeconds * self.nCohInt)
261 261
262 262 fmax = PRF
263 263 return fmax
264 264
265 265 def getVmax(self):
266 266
267 267 _lambda = self.C / self.frequency
268 268
269 269 vmax = self.getFmax() * _lambda / 2
270 270
271 271 return vmax
272 272
273 273 @property
274 274 def ippSeconds(self):
275 275 '''
276 276 '''
277 277 return self.radarControllerHeaderObj.ippSeconds
278 278
279 279 @ippSeconds.setter
280 280 def ippSeconds(self, ippSeconds):
281 281 '''
282 282 '''
283 283 self.radarControllerHeaderObj.ippSeconds = ippSeconds
284 284
285 285 @property
286 286 def code(self):
287 287 '''
288 288 '''
289 289 return self.radarControllerHeaderObj.code
290 290
291 291 @code.setter
292 292 def code(self, code):
293 293 '''
294 294 '''
295 295 self.radarControllerHeaderObj.code = code
296 296
297 297 @property
298 298 def nCode(self):
299 299 '''
300 300 '''
301 301 return self.radarControllerHeaderObj.nCode
302 302
303 303 @nCode.setter
304 304 def nCode(self, ncode):
305 305 '''
306 306 '''
307 307 self.radarControllerHeaderObj.nCode = ncode
308 308
309 309 @property
310 310 def nBaud(self):
311 311 '''
312 312 '''
313 313 return self.radarControllerHeaderObj.nBaud
314 314
315 315 @nBaud.setter
316 316 def nBaud(self, nbaud):
317 317 '''
318 318 '''
319 319 self.radarControllerHeaderObj.nBaud = nbaud
320 320
321 321 @property
322 322 def ipp(self):
323 323 '''
324 324 '''
325 325 return self.radarControllerHeaderObj.ipp
326 326
327 327 @ipp.setter
328 328 def ipp(self, ipp):
329 329 '''
330 330 '''
331 331 self.radarControllerHeaderObj.ipp = ipp
332 332
333 333 @property
334 334 def metadata(self):
335 335 '''
336 336 '''
337 337
338 338 return {attr: getattr(self, attr) for attr in self.metadata_list}
339 339
340 340
341 341 class Voltage(JROData):
342 342
343 343 dataPP_POW = None
344 344 dataPP_DOP = None
345 345 dataPP_WIDTH = None
346 346 dataPP_SNR = None
347 347
348 348 def __init__(self):
349 349 '''
350 350 Constructor
351 351 '''
352 352
353 353 self.useLocalTime = True
354 354 self.radarControllerHeaderObj = RadarControllerHeader()
355 355 self.systemHeaderObj = SystemHeader()
356 356 self.type = "Voltage"
357 357 self.data = None
358 358 self.nProfiles = None
359 359 self.heightList = None
360 360 self.channelList = None
361 361 self.flagNoData = True
362 362 self.flagDiscontinuousBlock = False
363 363 self.utctime = None
364 364 self.timeZone = 0
365 365 self.dstFlag = None
366 366 self.errorCount = None
367 367 self.nCohInt = None
368 368 self.blocksize = None
369 369 self.flagCohInt = False
370 370 self.flagDecodeData = False # asumo q la data no esta decodificada
371 371 self.flagDeflipData = False # asumo q la data no esta sin flip
372 372 self.flagShiftFFT = False
373 373 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
374 374 self.profileIndex = 0
375 375 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
376 376 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
377 377
378 378 def getNoisebyHildebrand(self, channel=None):
379 379 """
380 380 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
381 381
382 382 Return:
383 383 noiselevel
384 384 """
385 385
386 386 if channel != None:
387 387 data = self.data[channel]
388 388 nChannels = 1
389 389 else:
390 390 data = self.data
391 391 nChannels = self.nChannels
392 392
393 393 noise = numpy.zeros(nChannels)
394 394 power = data * numpy.conjugate(data)
395 395
396 396 for thisChannel in range(nChannels):
397 397 if nChannels == 1:
398 398 daux = power[:].real
399 399 else:
400 400 daux = power[thisChannel, :].real
401 401 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
402 402
403 403 return noise
404 404
405 405 def getNoise(self, type=1, channel=None):
406 406
407 407 if type == 1:
408 408 noise = self.getNoisebyHildebrand(channel)
409 409
410 410 return noise
411 411
412 412 def getPower(self, channel=None):
413 413
414 414 if channel != None:
415 415 data = self.data[channel]
416 416 else:
417 417 data = self.data
418 418
419 419 power = data * numpy.conjugate(data)
420 420 powerdB = 10 * numpy.log10(power.real)
421 421 powerdB = numpy.squeeze(powerdB)
422 422
423 423 return powerdB
424 424
425 425 @property
426 426 def timeInterval(self):
427 427
428 428 return self.ippSeconds * self.nCohInt
429 429
430 430 noise = property(getNoise, "I'm the 'nHeights' property.")
431 431
432 432
433 433 class CrossProds(JROData):
434 434
435 435 # data es un numpy array de 2 dmensiones (canales, alturas)
436 436 data = None
437 437
438 438 def __init__(self):
439 439 '''
440 440 Constructor
441 441 '''
442 442
443 443 self.useLocalTime = True
444 444 '''
445 445 self.radarControllerHeaderObj = RadarControllerHeader()
446 446 self.systemHeaderObj = SystemHeader()
447 447 self.type = "Voltage"
448 448 self.data = None
449 449 # self.dtype = None
450 450 # self.nChannels = 0
451 451 # self.nHeights = 0
452 452 self.nProfiles = None
453 453 self.heightList = None
454 454 self.channelList = None
455 455 # self.channelIndexList = None
456 456 self.flagNoData = True
457 457 self.flagDiscontinuousBlock = False
458 458 self.utctime = None
459 459 self.timeZone = None
460 460 self.dstFlag = None
461 461 self.errorCount = None
462 462 self.nCohInt = None
463 463 self.blocksize = None
464 464 self.flagDecodeData = False # asumo q la data no esta decodificada
465 465 self.flagDeflipData = False # asumo q la data no esta sin flip
466 466 self.flagShiftFFT = False
467 467 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
468 468 self.profileIndex = 0
469 469
470 470
471 471 def getNoisebyHildebrand(self, channel=None):
472 472
473 473
474 474 if channel != None:
475 475 data = self.data[channel]
476 476 nChannels = 1
477 477 else:
478 478 data = self.data
479 479 nChannels = self.nChannels
480 480
481 481 noise = numpy.zeros(nChannels)
482 482 power = data * numpy.conjugate(data)
483 483
484 484 for thisChannel in range(nChannels):
485 485 if nChannels == 1:
486 486 daux = power[:].real
487 487 else:
488 488 daux = power[thisChannel, :].real
489 489 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
490 490
491 491 return noise
492 492
493 493 def getNoise(self, type=1, channel=None):
494 494
495 495 if type == 1:
496 496 noise = self.getNoisebyHildebrand(channel)
497 497
498 498 return noise
499 499
500 500 def getPower(self, channel=None):
501 501
502 502 if channel != None:
503 503 data = self.data[channel]
504 504 else:
505 505 data = self.data
506 506
507 507 power = data * numpy.conjugate(data)
508 508 powerdB = 10 * numpy.log10(power.real)
509 509 powerdB = numpy.squeeze(powerdB)
510 510
511 511 return powerdB
512 512
513 513 def getTimeInterval(self):
514 514
515 515 timeInterval = self.ippSeconds * self.nCohInt
516 516
517 517 return timeInterval
518 518
519 519 noise = property(getNoise, "I'm the 'nHeights' property.")
520 520 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
521 521 '''
522 522 def getTimeInterval(self):
523 523
524 524 timeInterval = self.ippSeconds * self.nCohInt
525 525
526 526 return timeInterval
527 527
528 528
529 529
530 530 class Spectra(JROData):
531 531
532 532 def __init__(self):
533 533 '''
534 534 Constructor
535 535 '''
536 536
537 537 self.data_dc = None
538 538 self.data_spc = None
539 539 self.data_cspc = None
540 540 self.useLocalTime = True
541 541 self.radarControllerHeaderObj = RadarControllerHeader()
542 542 self.systemHeaderObj = SystemHeader()
543 543 self.type = "Spectra"
544 544 self.timeZone = 0
545 545 self.nProfiles = None
546 546 self.heightList = None
547 547 self.channelList = None
548 548 self.pairsList = None
549 549 self.flagNoData = True
550 550 self.flagDiscontinuousBlock = False
551 551 self.utctime = None
552 552 self.nCohInt = None
553 553 self.nIncohInt = None
554 554 self.blocksize = None
555 555 self.nFFTPoints = None
556 556 self.wavelength = None
557 557 self.flagDecodeData = False # asumo q la data no esta decodificada
558 558 self.flagDeflipData = False # asumo q la data no esta sin flip
559 559 self.flagShiftFFT = False
560 560 self.ippFactor = 1
561 561 self.beacon_heiIndexList = []
562 562 self.noise_estimation = None
563 self.spc_noise = None
563 564 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
564 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp', 'nIncohInt', 'nFFTPoints', 'nProfiles']
565 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp', 'nIncohInt', 'nFFTPoints', 'nProfiles', 'flagDecodeData']
565 566
566 567 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
567 568 """
568 569 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
569 570
570 571 Return:
571 572 noiselevel
572 573 """
573 574
574 575 noise = numpy.zeros(self.nChannels)
575 576
576 577 for channel in range(self.nChannels):
577 578 daux = self.data_spc[channel,
578 579 xmin_index:xmax_index, ymin_index:ymax_index]
579 580 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
580 581
581 582 return noise
582 583
583 584 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
584 585
585 if self.noise_estimation is not None:
586 if self.spc_noise is not None:
587 # this was estimated by getNoise Operation defined in jroproc_parameters.py
588 return self.spc_noise
589 elif self.noise_estimation is not None:
586 590 # this was estimated by getNoise Operation defined in jroproc_spectra.py
587 591 return self.noise_estimation
588 592 else:
589 593 noise = self.getNoisebyHildebrand(
590 594 xmin_index, xmax_index, ymin_index, ymax_index)
591 595 return noise
592 596
593 597 def getFreqRangeTimeResponse(self, extrapoints=0):
594 598
595 599 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
596 600 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
597 601
598 602 return freqrange
599 603
600 604 def getAcfRange(self, extrapoints=0):
601 605
602 606 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
603 607 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
604 608
605 609 return freqrange
606 610
607 611 def getFreqRange(self, extrapoints=0):
608 612
609 613 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
610 614 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
611 615
612 616 return freqrange
613 617
614 618 def getVelRange(self, extrapoints=0):
615 619
616 620 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
617 621 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
618 622
619 623 if self.nmodes:
620 624 return velrange / self.nmodes
621 625 else:
622 626 return velrange
623 627
624 628 @property
625 629 def nPairs(self):
626 630
627 631 return len(self.pairsList)
628 632
629 633 @property
630 634 def pairsIndexList(self):
631 635
632 636 return list(range(self.nPairs))
633 637
634 638 @property
635 639 def normFactor(self):
636 640
637 641 pwcode = 1
638 642
639 643 if self.flagDecodeData:
640 644 pwcode = numpy.sum(self.code[0] ** 2)
641 645 # normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
642 646 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
643 647
644 648 return normFactor
645 649
646 650 @property
647 651 def flag_cspc(self):
648 652
649 653 if self.data_cspc is None:
650 654 return True
651 655
652 656 return False
653 657
654 658 @property
655 659 def flag_dc(self):
656 660
657 661 if self.data_dc is None:
658 662 return True
659 663
660 664 return False
661 665
662 666 @property
663 667 def timeInterval(self):
664 668
665 669 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
666 670 if self.nmodes:
667 671 return self.nmodes * timeInterval
668 672 else:
669 673 return timeInterval
670 674
671 675 def getPower(self):
672 676
673 677 factor = self.normFactor
674 678 z = self.data_spc / factor
675 679 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
676 680 avg = numpy.average(z, axis=1)
677 681
678 682 return 10 * numpy.log10(avg)
679 683
680 684 def getCoherence(self, pairsList=None, phase=False):
681 685
682 686 z = []
683 687 if pairsList is None:
684 688 pairsIndexList = self.pairsIndexList
685 689 else:
686 690 pairsIndexList = []
687 691 for pair in pairsList:
688 692 if pair not in self.pairsList:
689 693 raise ValueError("Pair %s is not in dataOut.pairsList" % (
690 694 pair))
691 695 pairsIndexList.append(self.pairsList.index(pair))
692 696 for i in range(len(pairsIndexList)):
693 697 pair = self.pairsList[pairsIndexList[i]]
694 698 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
695 699 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
696 700 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
697 701 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
698 702 if phase:
699 703 data = numpy.arctan2(avgcoherenceComplex.imag,
700 704 avgcoherenceComplex.real) * 180 / numpy.pi
701 705 else:
702 706 data = numpy.abs(avgcoherenceComplex)
703 707
704 708 z.append(data)
705 709
706 710 return numpy.array(z)
707 711
708 712 def setValue(self, value):
709 713
710 714 print("This property should not be initialized")
711 715
712 716 return
713 717
714 718 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
715 719
716 720
717 721 class SpectraHeis(Spectra):
718 722
719 723 def __init__(self):
720 724
721 725 self.radarControllerHeaderObj = RadarControllerHeader()
722 726 self.systemHeaderObj = SystemHeader()
723 727 self.type = "SpectraHeis"
724 728 self.nProfiles = None
725 729 self.heightList = None
726 730 self.channelList = None
727 731 self.flagNoData = True
728 732 self.flagDiscontinuousBlock = False
729 733 self.utctime = None
730 734 self.blocksize = None
731 735 self.profileIndex = 0
732 736 self.nCohInt = 1
733 737 self.nIncohInt = 1
734 738
735 739 @property
736 740 def normFactor(self):
737 741 pwcode = 1
738 742 if self.flagDecodeData:
739 743 pwcode = numpy.sum(self.code[0] ** 2)
740 744
741 745 normFactor = self.nIncohInt * self.nCohInt * pwcode
742 746
743 747 return normFactor
744 748
745 749 @property
746 750 def timeInterval(self):
747 751
748 752 return self.ippSeconds * self.nCohInt * self.nIncohInt
749 753
750 754
751 755 class Fits(JROData):
752 756
753 757 def __init__(self):
754 758
755 759 self.type = "Fits"
756 760 self.nProfiles = None
757 761 self.heightList = None
758 762 self.channelList = None
759 763 self.flagNoData = True
760 764 self.utctime = None
761 765 self.nCohInt = 1
762 766 self.nIncohInt = 1
763 767 self.useLocalTime = True
764 768 self.profileIndex = 0
765 769 self.timeZone = 0
766 770
767 771 def getTimeRange(self):
768 772
769 773 datatime = []
770 774
771 775 datatime.append(self.ltctime)
772 776 datatime.append(self.ltctime + self.timeInterval)
773 777
774 778 datatime = numpy.array(datatime)
775 779
776 780 return datatime
777 781
778 782 def getChannelIndexList(self):
779 783
780 784 return list(range(self.nChannels))
781 785
782 786 def getNoise(self, type=1):
783 787
784 788
785 789 if type == 1:
786 790 noise = self.getNoisebyHildebrand()
787 791
788 792 if type == 2:
789 793 noise = self.getNoisebySort()
790 794
791 795 if type == 3:
792 796 noise = self.getNoisebyWindow()
793 797
794 798 return noise
795 799
796 800 @property
797 801 def timeInterval(self):
798 802
799 803 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
800 804
801 805 return timeInterval
802 806
803 807 @property
804 808 def ippSeconds(self):
805 809 '''
806 810 '''
807 811 return self.ipp_sec
808 812
809 813 noise = property(getNoise, "I'm the 'nHeights' property.")
810 814
811 815
812 816 class Correlation(JROData):
813 817
814 818 def __init__(self):
815 819 '''
816 820 Constructor
817 821 '''
818 822 self.radarControllerHeaderObj = RadarControllerHeader()
819 823 self.systemHeaderObj = SystemHeader()
820 824 self.type = "Correlation"
821 825 self.data = None
822 826 self.dtype = None
823 827 self.nProfiles = None
824 828 self.heightList = None
825 829 self.channelList = None
826 830 self.flagNoData = True
827 831 self.flagDiscontinuousBlock = False
828 832 self.utctime = None
829 833 self.timeZone = 0
830 834 self.dstFlag = None
831 835 self.errorCount = None
832 836 self.blocksize = None
833 837 self.flagDecodeData = False # asumo q la data no esta decodificada
834 838 self.flagDeflipData = False # asumo q la data no esta sin flip
835 839 self.pairsList = None
836 840 self.nPoints = None
837 841
838 842 def getPairsList(self):
839 843
840 844 return self.pairsList
841 845
842 846 def getNoise(self, mode=2):
843 847
844 848 indR = numpy.where(self.lagR == 0)[0][0]
845 849 indT = numpy.where(self.lagT == 0)[0][0]
846 850
847 851 jspectra0 = self.data_corr[:, :, indR, :]
848 852 jspectra = copy.copy(jspectra0)
849 853
850 854 num_chan = jspectra.shape[0]
851 855 num_hei = jspectra.shape[2]
852 856
853 857 freq_dc = jspectra.shape[1] / 2
854 858 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
855 859
856 860 if ind_vel[0] < 0:
857 861 ind_vel[list(range(0, 1))] = ind_vel[list(
858 862 range(0, 1))] + self.num_prof
859 863
860 864 if mode == 1:
861 865 jspectra[:, freq_dc, :] = (
862 866 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
863 867
864 868 if mode == 2:
865 869
866 870 vel = numpy.array([-2, -1, 1, 2])
867 871 xx = numpy.zeros([4, 4])
868 872
869 873 for fil in range(4):
870 874 xx[fil, :] = vel[fil] ** numpy.asarray(list(range(4)))
871 875
872 876 xx_inv = numpy.linalg.inv(xx)
873 877 xx_aux = xx_inv[0, :]
874 878
875 879 for ich in range(num_chan):
876 880 yy = jspectra[ich, ind_vel, :]
877 881 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
878 882
879 883 junkid = jspectra[ich, freq_dc, :] <= 0
880 884 cjunkid = sum(junkid)
881 885
882 886 if cjunkid.any():
883 887 jspectra[ich, freq_dc, junkid.nonzero()] = (
884 888 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
885 889
886 890 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
887 891
888 892 return noise
889 893
890 894 @property
891 895 def timeInterval(self):
892 896
893 897 return self.ippSeconds * self.nCohInt * self.nProfiles
894 898
895 899 def splitFunctions(self):
896 900
897 901 pairsList = self.pairsList
898 902 ccf_pairs = []
899 903 acf_pairs = []
900 904 ccf_ind = []
901 905 acf_ind = []
902 906 for l in range(len(pairsList)):
903 907 chan0 = pairsList[l][0]
904 908 chan1 = pairsList[l][1]
905 909
906 910 # Obteniendo pares de Autocorrelacion
907 911 if chan0 == chan1:
908 912 acf_pairs.append(chan0)
909 913 acf_ind.append(l)
910 914 else:
911 915 ccf_pairs.append(pairsList[l])
912 916 ccf_ind.append(l)
913 917
914 918 data_acf = self.data_cf[acf_ind]
915 919 data_ccf = self.data_cf[ccf_ind]
916 920
917 921 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
918 922
919 923 @property
920 924 def normFactor(self):
921 925 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
922 926 acf_pairs = numpy.array(acf_pairs)
923 927 normFactor = numpy.zeros((self.nPairs, self.nHeights))
924 928
925 929 for p in range(self.nPairs):
926 930 pair = self.pairsList[p]
927 931
928 932 ch0 = pair[0]
929 933 ch1 = pair[1]
930 934
931 935 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
932 936 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
933 937 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
934 938
935 939 return normFactor
936 940
937 941
938 942 class Parameters(Spectra):
939 943
940 944 groupList = None # List of Pairs, Groups, etc
941 945 data_param = None # Parameters obtained
942 946 data_pre = None # Data Pre Parametrization
943 947 data_SNR = None # Signal to Noise Ratio
944 948 abscissaList = None # Abscissa, can be velocities, lags or time
945 949 utctimeInit = None # Initial UTC time
946 950 paramInterval = None # Time interval to calculate Parameters in seconds
947 951 useLocalTime = True
948 952 # Fitting
949 953 data_error = None # Error of the estimation
950 954 constants = None
951 955 library = None
952 956 # Output signal
953 957 outputInterval = None # Time interval to calculate output signal in seconds
954 958 data_output = None # Out signal
955 959 nAvg = None
956 960 noise_estimation = None
957 961 GauSPC = None # Fit gaussian SPC
962 spc_noise = None
958 963
959 964 def __init__(self):
960 965 '''
961 966 Constructor
962 967 '''
963 968 self.radarControllerHeaderObj = RadarControllerHeader()
964 969 self.systemHeaderObj = SystemHeader()
965 970 self.type = "Parameters"
966 971 self.timeZone = 0
967 972 self.ippFactor = 1
968 973
969 974 def getTimeRange1(self, interval):
970 975
971 976 datatime = []
972 977
973 978 if self.useLocalTime:
974 979 time1 = self.utctimeInit - self.timeZone * 60
975 980 else:
976 981 time1 = self.utctimeInit
977 982
978 983 datatime.append(time1)
979 984 datatime.append(time1 + interval)
980 985 datatime = numpy.array(datatime)
981 986
982 987 return datatime
983 988
984 989 @property
985 990 def timeInterval(self):
986 991
987 992 if hasattr(self, 'timeInterval1'):
988 993 return self.timeInterval1
989 994 else:
990 995 return self.paramInterval
991 996
992 997
993 998 def setValue(self, value):
994 999
995 1000 print("This property should not be initialized")
996 1001
997 1002 return
998 1003
999 def getNoise(self):
1000
1001 return self.spc_noise
1002
1003 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1004
1005 1004
1006 1005 class PlotterData(object):
1007 1006 '''
1008 1007 Object to hold data to be plotted
1009 1008 '''
1010 1009
1011 1010 MAXNUMX = 200
1012 1011 MAXNUMY = 200
1013 1012
1014 1013 def __init__(self, code, exp_code, localtime=True):
1015 1014
1016 1015 self.key = code
1017 1016 self.exp_code = exp_code
1018 1017 self.ready = False
1019 1018 self.flagNoData = False
1020 1019 self.localtime = localtime
1021 1020 self.data = {}
1022 1021 self.meta = {}
1023 1022 self.__heights = []
1024 1023
1025 1024 def __str__(self):
1026 1025 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1027 1026 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1028 1027
1029 1028 def __len__(self):
1030 1029 return len(self.data)
1031 1030
1032 1031 def __getitem__(self, key):
1033 1032 if isinstance(key, int):
1034 1033 return self.data[self.times[key]]
1035 1034 elif isinstance(key, str):
1036 1035 ret = numpy.array([self.data[x][key] for x in self.times])
1037 1036 if ret.ndim > 1:
1038 1037 ret = numpy.swapaxes(ret, 0, 1)
1039 1038 return ret
1040 1039
1041 1040 def __contains__(self, key):
1042 1041 return key in self.data[self.min_time]
1043 1042
1044 1043 def setup(self):
1045 1044 '''
1046 1045 Configure object
1047 1046 '''
1048 1047 self.type = ''
1049 1048 self.ready = False
1050 1049 del self.data
1051 1050 self.data = {}
1052 1051 self.__heights = []
1053 1052 self.__all_heights = set()
1054 1053
1055 1054 def shape(self, key):
1056 1055 '''
1057 1056 Get the shape of the one-element data for the given key
1058 1057 '''
1059 1058
1060 1059 if len(self.data[self.min_time][key]):
1061 1060 return self.data[self.min_time][key].shape
1062 1061 return (0,)
1063 1062
1064 1063 def update(self, data, tm, meta={}):
1065 1064 '''
1066 1065 Update data object with new dataOut
1067 1066 '''
1068 1067
1069 1068 self.data[tm] = data
1070 1069
1071 1070 for key, value in meta.items():
1072 1071 setattr(self, key, value)
1073 1072
1074 1073 def normalize_heights(self):
1075 1074 '''
1076 1075 Ensure same-dimension of the data for different heighList
1077 1076 '''
1078 1077
1079 1078 H = numpy.array(list(self.__all_heights))
1080 1079 H.sort()
1081 1080 for key in self.data:
1082 1081 shape = self.shape(key)[:-1] + H.shape
1083 1082 for tm, obj in list(self.data[key].items()):
1084 1083 h = self.__heights[self.times.tolist().index(tm)]
1085 1084 if H.size == h.size:
1086 1085 continue
1087 1086 index = numpy.where(numpy.in1d(H, h))[0]
1088 1087 dummy = numpy.zeros(shape) + numpy.nan
1089 1088 if len(shape) == 2:
1090 1089 dummy[:, index] = obj
1091 1090 else:
1092 1091 dummy[index] = obj
1093 1092 self.data[key][tm] = dummy
1094 1093
1095 1094 self.__heights = [H for tm in self.times]
1096 1095
1097 1096 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1098 1097 '''
1099 1098 Convert data to json
1100 1099 '''
1101 1100
1102 1101 meta = {}
1103 1102 meta['xrange'] = []
1104 1103 dy = int(len(self.yrange) / self.MAXNUMY) + 1
1105 1104 tmp = self.data[tm][self.key]
1106 1105 shape = tmp.shape
1107 1106 if len(shape) == 2:
1108 1107 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1109 1108 elif len(shape) == 3:
1110 1109 dx = int(self.data[tm][self.key].shape[1] / self.MAXNUMX) + 1
1111 1110 data = self.roundFloats(
1112 1111 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1113 1112 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1114 1113 else:
1115 1114 data = self.roundFloats(self.data[tm][self.key].tolist())
1116 1115
1117 1116 ret = {
1118 1117 'plot': plot_name,
1119 1118 'code': self.exp_code,
1120 1119 'time': float(tm),
1121 1120 'data': data,
1122 1121 }
1123 1122 meta['type'] = plot_type
1124 1123 meta['interval'] = float(self.interval)
1125 1124 meta['localtime'] = self.localtime
1126 1125 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1127 1126 meta.update(self.meta)
1128 1127 ret['metadata'] = meta
1129 1128 return json.dumps(ret)
1130 1129
1131 1130 @property
1132 1131 def times(self):
1133 1132 '''
1134 1133 Return the list of times of the current data
1135 1134 '''
1136 1135
1137 1136 ret = [t for t in self.data]
1138 1137 ret.sort()
1139 1138 return numpy.array(ret)
1140 1139
1141 1140 @property
1142 1141 def min_time(self):
1143 1142 '''
1144 1143 Return the minimun time value
1145 1144 '''
1146 1145
1147 1146 return self.times[0]
1148 1147
1149 1148 @property
1150 1149 def max_time(self):
1151 1150 '''
1152 1151 Return the maximun time value
1153 1152 '''
1154 1153
1155 1154 return self.times[-1]
1156 1155
1157 1156 # @property
1158 1157 # def heights(self):
1159 1158 # '''
1160 1159 # Return the list of heights of the current data
1161 1160 # '''
1162 1161
1163 1162 # return numpy.array(self.__heights[-1])
1164 1163
1165 1164 @staticmethod
1166 1165 def roundFloats(obj):
1167 1166 if isinstance(obj, list):
1168 1167 return list(map(PlotterData.roundFloats, obj))
1169 1168 elif isinstance(obj, float):
1170 1169 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,873 +1,873
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 Examples
42 42 --------
43 43
44 44 desc = {
45 45 'Data': {
46 46 'data_output': ['u', 'v', 'w'],
47 47 'utctime': 'timestamps',
48 48 } ,
49 49 'Metadata': {
50 50 'heightList': 'heights'
51 51 }
52 52 }
53 53
54 54 desc = {
55 55 'Data': {
56 56 'data_output': 'winds',
57 57 'utctime': 'timestamps'
58 58 },
59 59 'Metadata': {
60 60 'heightList': 'heights'
61 61 }
62 62 }
63 63
64 64 extras = {
65 65 'timeZone': 300
66 66 }
67 67
68 68 reader = project.addReadUnit(
69 69 name='HDFReader',
70 70 path='/path/to/files',
71 71 startDate='2019/01/01',
72 72 endDate='2019/01/31',
73 73 startTime='00:00:00',
74 74 endTime='23:59:59',
75 75 # description=json.dumps(desc),
76 76 # extras=json.dumps(extras),
77 77 )
78 78
79 79 """
80 80
81 81 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
82 82
83 83 def __init__(self):
84 84 ProcessingUnit.__init__(self)
85 85 self.dataOut = Parameters()
86 86 self.ext = ".hdf5"
87 87 self.optchar = "D"
88 88 self.meta = {}
89 89 self.data = {}
90 90 self.open_file = h5py.File
91 91 self.open_mode = 'r'
92 92 self.description = {}
93 93 self.extras = {}
94 94 self.filefmt = "*%Y%j***"
95 95 self.folderfmt = "*%Y%j"
96 96 self.utcoffset = 0
97 97
98 98 def setup(self, **kwargs):
99 99
100 100 self.set_kwargs(**kwargs)
101 101 if not self.ext.startswith('.'):
102 102 self.ext = '.{}'.format(self.ext)
103 103
104 104 if self.online:
105 105 log.log("Searching files in online mode...", self.name)
106 106
107 107 for nTries in range(self.nTries):
108 108 fullpath = self.searchFilesOnLine(self.path, self.startDate,
109 109 self.endDate, self.expLabel, self.ext, self.walk,
110 110 self.filefmt, self.folderfmt)
111 111 try:
112 112 fullpath = next(fullpath)
113 113 except:
114 114 fullpath = None
115 115
116 116 if fullpath:
117 117 break
118 118
119 119 log.warning(
120 120 'Waiting {} sec for a valid file in {}: try {} ...'.format(
121 121 self.delay, self.path, nTries + 1),
122 122 self.name)
123 123 time.sleep(self.delay)
124 124
125 125 if not(fullpath):
126 126 raise schainpy.admin.SchainError(
127 127 'There isn\'t any valid file in {}'.format(self.path))
128 128
129 129 pathname, filename = os.path.split(fullpath)
130 130 self.year = int(filename[1:5])
131 131 self.doy = int(filename[5:8])
132 132 self.set = int(filename[8:11]) - 1
133 133 else:
134 134 log.log("Searching files in {}".format(self.path), self.name)
135 135 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
136 136 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
137 137
138 138 self.setNextFile()
139 139
140 140 return
141 141
142 142 def readFirstHeader(self):
143 143 '''Read metadata and data'''
144 144
145 145 self.__readMetadata()
146 146 self.__readData()
147 self.__setBlockList()
148
149 if 'type' in self.meta:
150 self.dataOut = eval(self.meta['type'])()
147 self.__setBlockList()
151 148
152 149 for attr in self.meta:
153 150 setattr(self.dataOut, attr, self.meta[attr])
154 151
155 152 self.blockIndex = 0
156 153
157 154 return
158 155
159 156 def __setBlockList(self):
160 157 '''
161 158 Selects the data within the times defined
162 159
163 160 self.fp
164 161 self.startTime
165 162 self.endTime
166 163 self.blockList
167 164 self.blocksPerFile
168 165
169 166 '''
170 167
171 168 startTime = self.startTime
172 169 endTime = self.endTime
173 170 thisUtcTime = self.data['utctime'] + self.utcoffset
174 171 self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
175 172 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
176 173
177 174 thisDate = thisDatetime.date()
178 175 thisTime = thisDatetime.time()
179 176
180 177 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
181 178 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
182 179
183 180 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
184 181
185 182 self.blockList = ind
186 183 self.blocksPerFile = len(ind)
187 184 return
188 185
189 186 def __readMetadata(self):
190 187 '''
191 188 Reads Metadata
192 189 '''
193 190
194 191 meta = {}
195 192
196 193 if self.description:
197 194 for key, value in self.description['Metadata'].items():
198 195 meta[key] = self.fp[value][()]
199 196 else:
200 197 grp = self.fp['Metadata']
201 198 for name in grp:
202 199 meta[name] = grp[name][()]
203 200
204 201 if self.extras:
205 202 for key, value in self.extras.items():
206 203 meta[key] = value
207 204 self.meta = meta
208 205
209 206 return
210 207
211 208 def __readData(self):
212 209
213 210 data = {}
214 211
215 212 if self.description:
216 213 for key, value in self.description['Data'].items():
217 214 if isinstance(value, str):
218 215 if isinstance(self.fp[value], h5py.Dataset):
219 216 data[key] = self.fp[value][()]
220 217 elif isinstance(self.fp[value], h5py.Group):
221 218 array = []
222 219 for ch in self.fp[value]:
223 220 array.append(self.fp[value][ch][()])
224 221 data[key] = numpy.array(array)
225 222 elif isinstance(value, list):
226 223 array = []
227 224 for ch in value:
228 225 array.append(self.fp[ch][()])
229 226 data[key] = numpy.array(array)
230 227 else:
231 228 grp = self.fp['Data']
232 229 for name in grp:
233 230 if isinstance(grp[name], h5py.Dataset):
234 231 array = grp[name][()]
235 232 elif isinstance(grp[name], h5py.Group):
236 233 array = []
237 234 for ch in grp[name]:
238 235 array.append(grp[name][ch][()])
239 236 array = numpy.array(array)
240 237 else:
241 238 log.warning('Unknown type: {}'.format(name))
242 239
243 240 if name in self.description:
244 241 key = self.description[name]
245 242 else:
246 243 key = name
247 244 data[key] = array
248 245
249 246 self.data = data
250 247 return
251 248
252 249 def getData(self):
253 250
254 251 for attr in self.data:
255 252 if self.data[attr].ndim == 1:
256 253 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
257 254 else:
258 255 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
259 256
260 257 self.dataOut.flagNoData = False
261 258 self.blockIndex += 1
262 259
263 260 log.log("Block No. {}/{} -> {}".format(
264 261 self.blockIndex,
265 262 self.blocksPerFile,
266 263 self.dataOut.datatime.ctime()), self.name)
267 264
268 265 return
269 266
270 267 def run(self, **kwargs):
271 268
272 269 if not(self.isConfig):
273 270 self.setup(**kwargs)
274 271 self.isConfig = True
275 272
276 273 if self.blockIndex == self.blocksPerFile:
277 274 self.setNextFile()
278 275
279 276 self.getData()
280 277
278 if 'type' in self.meta:
279 self.dataOut.type = self.meta['type'].decode('utf-8')
280
281 281 return
282 282
283 283 @MPDecorator
284 284 class HDFWriter(Operation):
285 285 """Operation to write HDF5 files.
286 286
287 287 The HDF5 file contains by default two groups Data and Metadata where
288 288 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
289 289 parameters, data attributes are normaly time dependent where the metadata
290 290 are not.
291 291 It is possible to customize the structure of the HDF5 file with the
292 292 optional description parameter see the examples.
293 293
294 294 Parameters:
295 295 -----------
296 296 path : str
297 297 Path where files will be saved.
298 298 blocksPerFile : int
299 299 Number of blocks per file
300 300 metadataList : list
301 301 List of the dataOut attributes that will be saved as metadata
302 302 dataList : int
303 303 List of the dataOut attributes that will be saved as data
304 304 setType : bool
305 305 If True the name of the files corresponds to the timestamp of the data
306 306 description : dict, optional
307 307 Dictionary with the desired description of the HDF5 file
308 308
309 309 Examples
310 310 --------
311 311
312 312 desc = {
313 313 'data_output': {'winds': ['z', 'w', 'v']},
314 314 'utctime': 'timestamps',
315 315 'heightList': 'heights'
316 316 }
317 317 desc = {
318 318 'data_output': ['z', 'w', 'v'],
319 319 'utctime': 'timestamps',
320 320 'heightList': 'heights'
321 321 }
322 322 desc = {
323 323 'Data': {
324 324 'data_output': 'winds',
325 325 'utctime': 'timestamps'
326 326 },
327 327 'Metadata': {
328 328 'heightList': 'heights'
329 329 }
330 330 }
331 331
332 332 writer = proc_unit.addOperation(name='HDFWriter')
333 333 writer.addParameter(name='path', value='/path/to/file')
334 334 writer.addParameter(name='blocksPerFile', value='32')
335 335 writer.addParameter(name='metadataList', value='heightList,timeZone')
336 336 writer.addParameter(name='dataList',value='data_output,utctime')
337 337 # writer.addParameter(name='description',value=json.dumps(desc))
338 338
339 339 """
340 340
341 341 ext = ".hdf5"
342 342 optchar = "D"
343 343 filename = None
344 344 path = None
345 345 setFile = None
346 346 fp = None
347 347 firsttime = True
348 348 # Configurations
349 349 blocksPerFile = None
350 350 blockIndex = None
351 351 dataOut = None
352 352 # Data Arrays
353 353 dataList = None
354 354 metadataList = None
355 355 currentDay = None
356 356 lastTime = None
357 357
358 358 def __init__(self):
359 359
360 360 Operation.__init__(self)
361 361 return
362 362
363 363 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None):
364 364 self.path = path
365 365 self.blocksPerFile = blocksPerFile
366 366 self.metadataList = metadataList
367 367 self.dataList = [s.strip() for s in dataList]
368 368 self.setType = setType
369 369 self.description = description
370 370
371 371 if self.metadataList is None:
372 372 self.metadataList = self.dataOut.metadata_list
373 373
374 374 tableList = []
375 375 dsList = []
376 376
377 377 for i in range(len(self.dataList)):
378 378 dsDict = {}
379 379 if hasattr(self.dataOut, self.dataList[i]):
380 380 dataAux = getattr(self.dataOut, self.dataList[i])
381 381 dsDict['variable'] = self.dataList[i]
382 382 else:
383 383 log.warning('Attribute {} not found in dataOut', self.name)
384 384 continue
385 385
386 386 if dataAux is None:
387 387 continue
388 388 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
389 389 dsDict['nDim'] = 0
390 390 else:
391 391 dsDict['nDim'] = len(dataAux.shape)
392 392 dsDict['shape'] = dataAux.shape
393 393 dsDict['dsNumber'] = dataAux.shape[0]
394 394 dsDict['dtype'] = dataAux.dtype
395 395
396 396 dsList.append(dsDict)
397 397
398 398 self.dsList = dsList
399 399 self.currentDay = self.dataOut.datatime.date()
400 400
401 401 def timeFlag(self):
402 402 currentTime = self.dataOut.utctime
403 403 timeTuple = time.localtime(currentTime)
404 404 dataDay = timeTuple.tm_yday
405 405
406 406 if self.lastTime is None:
407 407 self.lastTime = currentTime
408 408 self.currentDay = dataDay
409 409 return False
410 410
411 411 timeDiff = currentTime - self.lastTime
412 412
413 413 # Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
414 414 if dataDay != self.currentDay:
415 415 self.currentDay = dataDay
416 416 return True
417 417 elif timeDiff > 3 * 60 * 60:
418 418 self.lastTime = currentTime
419 419 return True
420 420 else:
421 421 self.lastTime = currentTime
422 422 return False
423 423
424 424 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
425 425 dataList=[], setType=None, description={}):
426 426
427 427 self.dataOut = dataOut
428 428 if not(self.isConfig):
429 429 self.setup(path=path, blocksPerFile=blocksPerFile,
430 430 metadataList=metadataList, dataList=dataList,
431 431 setType=setType, description=description)
432 432
433 433 self.isConfig = True
434 434 self.setNextFile()
435 435
436 436 self.putData()
437 437 return
438 438
439 439 def setNextFile(self):
440 440
441 441 ext = self.ext
442 442 path = self.path
443 443 setFile = self.setFile
444 444
445 445 timeTuple = time.localtime(self.dataOut.utctime)
446 446 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year, timeTuple.tm_yday)
447 447 fullpath = os.path.join(path, subfolder)
448 448
449 449 if os.path.exists(fullpath):
450 450 filesList = os.listdir(fullpath)
451 451 filesList = [k for k in filesList if k.startswith(self.optchar)]
452 452 if len(filesList) > 0:
453 453 filesList = sorted(filesList, key=str.lower)
454 454 filen = filesList[-1]
455 455 # el filename debera tener el siguiente formato
456 456 # 0 1234 567 89A BCDE (hex)
457 457 # x YYYY DDD SSS .ext
458 458 if isNumber(filen[8:11]):
459 459 setFile = int(filen[8:11]) # inicializo mi contador de seteo al seteo del ultimo file
460 460 else:
461 461 setFile = -1
462 462 else:
463 463 setFile = -1 # inicializo mi contador de seteo
464 464 else:
465 465 os.makedirs(fullpath)
466 466 setFile = -1 # inicializo mi contador de seteo
467 467
468 468 if self.setType is None:
469 469 setFile += 1
470 470 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
471 471 timeTuple.tm_year,
472 472 timeTuple.tm_yday,
473 473 setFile,
474 474 ext)
475 475 else:
476 476 setFile = timeTuple.tm_hour * 60 + timeTuple.tm_min
477 477 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
478 478 timeTuple.tm_year,
479 479 timeTuple.tm_yday,
480 480 setFile,
481 481 ext)
482 482
483 483 self.filename = os.path.join(path, subfolder, file)
484 484
485 485 # Setting HDF5 File
486 486 self.fp = h5py.File(self.filename, 'w')
487 487 # write metadata
488 488 self.writeMetadata(self.fp)
489 489 # Write data
490 490 self.writeData(self.fp)
491 491
492 492 def getLabel(self, name, x=None):
493 493
494 494 if x is None:
495 495 if 'Data' in self.description:
496 496 data = self.description['Data']
497 497 if 'Metadata' in self.description:
498 498 data.update(self.description['Metadata'])
499 499 else:
500 500 data = self.description
501 501 if name in data:
502 502 if isinstance(data[name], str):
503 503 return data[name]
504 504 elif isinstance(data[name], list):
505 505 return None
506 506 elif isinstance(data[name], dict):
507 507 for key, value in data[name].items():
508 508 return key
509 509 return name
510 510 else:
511 511 if 'Metadata' in self.description:
512 512 meta = self.description['Metadata']
513 513 else:
514 514 meta = self.description
515 515 if name in meta:
516 516 if isinstance(meta[name], list):
517 517 return meta[name][x]
518 518 elif isinstance(meta[name], dict):
519 519 for key, value in meta[name].items():
520 520 return value[x]
521 521 if 'cspc' in name:
522 522 return 'pair{:02d}'.format(x)
523 523 else:
524 524 return 'channel{:02d}'.format(x)
525 525
526 526 def writeMetadata(self, fp):
527 527
528 528 if self.description:
529 529 if 'Metadata' in self.description:
530 530 grp = fp.create_group('Metadata')
531 531 else:
532 532 grp = fp
533 533 else:
534 534 grp = fp.create_group('Metadata')
535 535
536 536 for i in range(len(self.metadataList)):
537 537 if not hasattr(self.dataOut, self.metadataList[i]):
538 538 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
539 539 continue
540 540 value = getattr(self.dataOut, self.metadataList[i])
541 541 if isinstance(value, bool):
542 542 if value is True:
543 543 value = 1
544 544 else:
545 545 value = 0
546 546 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
547 547 return
548 548
549 549 def writeData(self, fp):
550 550
551 551 if self.description:
552 552 if 'Data' in self.description:
553 553 grp = fp.create_group('Data')
554 554 else:
555 555 grp = fp
556 556 else:
557 557 grp = fp.create_group('Data')
558 558
559 559 dtsets = []
560 560 data = []
561 561
562 562 for dsInfo in self.dsList:
563 563 if dsInfo['nDim'] == 0:
564 564 ds = grp.create_dataset(
565 565 self.getLabel(dsInfo['variable']),
566 566 (self.blocksPerFile,),
567 567 chunks=True,
568 568 dtype=numpy.float64)
569 569 dtsets.append(ds)
570 570 data.append((dsInfo['variable'], -1))
571 571 else:
572 572 label = self.getLabel(dsInfo['variable'])
573 573 if label is not None:
574 574 sgrp = grp.create_group(label)
575 575 else:
576 576 sgrp = grp
577 577 for i in range(dsInfo['dsNumber']):
578 578 ds = sgrp.create_dataset(
579 579 self.getLabel(dsInfo['variable'], i),
580 580 (self.blocksPerFile,) + dsInfo['shape'][1:],
581 581 chunks=True,
582 582 dtype=dsInfo['dtype'])
583 583 dtsets.append(ds)
584 584 data.append((dsInfo['variable'], i))
585 585 fp.flush()
586 586
587 587 log.log('Creating file: {}'.format(fp.filename), self.name)
588 588
589 589 self.ds = dtsets
590 590 self.data = data
591 591 self.firsttime = True
592 592 self.blockIndex = 0
593 593 return
594 594
595 595 def putData(self):
596 596
597 597 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
598 598 self.closeFile()
599 599 self.setNextFile()
600 600
601 601 for i, ds in enumerate(self.ds):
602 602 attr, ch = self.data[i]
603 603 if ch == -1:
604 604 ds[self.blockIndex] = getattr(self.dataOut, attr)
605 605 else:
606 606 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
607 607
608 608 self.fp.flush()
609 609 self.blockIndex += 1
610 610 log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
611 611
612 612 return
613 613
614 614 def closeFile(self):
615 615
616 616 if self.blockIndex != self.blocksPerFile:
617 617 for ds in self.ds:
618 618 ds.resize(self.blockIndex, axis=0)
619 619
620 620 if self.fp:
621 621 self.fp.flush()
622 622 self.fp.close()
623 623
624 624 def close(self):
625 625
626 626 self.closeFile()
627 627
628 628
629 629 @MPDecorator
630 630 class ASCIIWriter(Operation):
631 631 """Operation to write data in ascii files.
632 632
633 633 Parameters:
634 634 -----------
635 635 path : str
636 636 Path where files will be saved.
637 637 blocksPerFile : int
638 638 Number of blocks per file
639 639 metadataList : list
640 640 List of the dataOut attributes that will be saved as metadata
641 641 dataDict : dict
642 642 Dictionary with the varaibles to be saved
643 643 setType : bool
644 644 If True the name of the files corresponds to the timestamp of the data
645 645
646 646 Examples
647 647 --------
648 648
649 649 data = {
650 650 'data_output': ['z', 'w', 'v'],
651 651 'utctime': 'time',
652 652 'heightList': 'height'
653 653 }
654 654
655 655 writer = proc_unit.addOperation(name='ASCIIWriter')
656 656 writer.addParameter(name='path', value='/path/to/file')
657 657 writer.addParameter(name='blocksPerFile', value='32')
658 658 writer.addParameter(name='dataDict',value=json.dumps(data))
659 659
660 660 """
661 661
662 662 ext = ".txt"
663 663 optchar = "D"
664 664 filename = None
665 665 path = None
666 666 setFile = None
667 667 fp = None
668 668 firsttime = True
669 669 # Configurations
670 670 blocksPerFile = None
671 671 blockIndex = None
672 672 dataOut = None
673 673 # Data Arrays
674 674 dataDict = None
675 675 metadataList = None
676 676 currentDay = None
677 677 lastTime = None
678 678 localtime = True
679 679
680 680 def __init__(self):
681 681
682 682 Operation.__init__(self)
683 683 return
684 684
685 685 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataDict=None, setType=None, localtime=True):
686 686 self.path = path
687 687 self.blocksPerFile = blocksPerFile
688 688 self.metadataList = metadataList
689 689 self.dataDict = dataDict
690 690 self.setType = setType
691 691 self.localtime = localtime
692 692
693 693 if self.metadataList is None:
694 694 self.metadataList = self.dataOut.metadata_list
695 695
696 696 dsList = []
697 697
698 698 for key, value in self.dataDict.items():
699 699 dsDict = {}
700 700 if hasattr(self.dataOut, key):
701 701 dataAux = getattr(self.dataOut, key)
702 702 dsDict['variable'] = key
703 703 else:
704 704 log.warning('Attribute {} not found in dataOut', self.name)
705 705 continue
706 706
707 707 if dataAux is None:
708 708 continue
709 709 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
710 710 dsDict['nDim'] = 0
711 711 else:
712 712 dsDict['nDim'] = len(dataAux.shape)
713 713 dsDict['shape'] = dataAux.shape
714 714 dsDict['dsNumber'] = dataAux.shape[0]
715 715 dsDict['dtype'] = dataAux.dtype
716 716
717 717 dsList.append(dsDict)
718 718 self.dsList = dsList
719 719 self.currentDay = self.dataOut.datatime.date()
720 720
721 721 def timeFlag(self):
722 722 currentTime = self.dataOut.utctime
723 723 if self.localtime:
724 724 timeTuple = time.localtime(currentTime)
725 725 else:
726 726 timeTuple = time.gmtime(currentTime)
727 727
728 728 dataDay = timeTuple.tm_yday
729 729
730 730 if self.lastTime is None:
731 731 self.lastTime = currentTime
732 732 self.currentDay = dataDay
733 733 return False
734 734
735 735 timeDiff = currentTime - self.lastTime
736 736
737 737 # Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
738 738 if dataDay != self.currentDay:
739 739 self.currentDay = dataDay
740 740 return True
741 741 elif timeDiff > 3 * 60 * 60:
742 742 self.lastTime = currentTime
743 743 return True
744 744 else:
745 745 self.lastTime = currentTime
746 746 return False
747 747
748 748 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
749 749 dataDict={}, setType=None, localtime=True):
750 750
751 751 self.dataOut = dataOut
752 752 if not(self.isConfig):
753 753 self.setup(path=path, blocksPerFile=blocksPerFile,
754 754 metadataList=metadataList, dataDict=dataDict,
755 755 setType=setType, localtime=localtime)
756 756
757 757 self.isConfig = True
758 758 self.setNextFile()
759 759
760 760 self.putData()
761 761 return
762 762
763 763 def setNextFile(self):
764 764
765 765 ext = self.ext
766 766 path = self.path
767 767 setFile = self.setFile
768 768 if self.localtime:
769 769 timeTuple = time.localtime(self.dataOut.utctime)
770 770 else:
771 771 timeTuple = time.gmtime(self.dataOut.utctime)
772 772 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year, timeTuple.tm_yday)
773 773 fullpath = os.path.join(path, subfolder)
774 774
775 775 if os.path.exists(fullpath):
776 776 filesList = os.listdir(fullpath)
777 777 filesList = [k for k in filesList if k.startswith(self.optchar)]
778 778 if len(filesList) > 0:
779 779 filesList = sorted(filesList, key=str.lower)
780 780 filen = filesList[-1]
781 781 # el filename debera tener el siguiente formato
782 782 # 0 1234 567 89A BCDE (hex)
783 783 # x YYYY DDD SSS .ext
784 784 if isNumber(filen[8:11]):
785 785 setFile = int(filen[8:11]) # inicializo mi contador de seteo al seteo del ultimo file
786 786 else:
787 787 setFile = -1
788 788 else:
789 789 setFile = -1 # inicializo mi contador de seteo
790 790 else:
791 791 os.makedirs(fullpath)
792 792 setFile = -1 # inicializo mi contador de seteo
793 793
794 794 if self.setType is None:
795 795 setFile += 1
796 796 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
797 797 timeTuple.tm_year,
798 798 timeTuple.tm_yday,
799 799 setFile,
800 800 ext)
801 801 else:
802 802 setFile = timeTuple.tm_hour * 60 + timeTuple.tm_min
803 803 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
804 804 timeTuple.tm_year,
805 805 timeTuple.tm_yday,
806 806 setFile,
807 807 ext)
808 808
809 809 self.filename = os.path.join(path, subfolder, file)
810 810
811 811 # Setting HDF5 File
812 812 self.fp = open(self.filename, 'w')
813 813 # write metadata
814 814 self.writeMetadata(self.fp)
815 815 # Write data
816 816 self.writeData(self.fp)
817 817
818 818 def writeMetadata(self, fp):
819 819
820 820 line = ''
821 821 for d in self.dsList:
822 822 par = self.dataDict[d['variable']]
823 823 if isinstance(par, (list,tuple)):
824 824 for p in par:
825 825 line += '{:>16}'.format(p)
826 826 else:
827 827 line += '{:>16}'.format(par)
828 828
829 829 line += '\n'
830 830 fp.write(line)
831 831
832 832 def writeData(self, fp):
833 833
834 834 log.log('Creating file: {}'.format(self.filename), self.name)
835 835
836 836 self.firsttime = True
837 837 self.blockIndex = 0
838 838 return
839 839
840 840 def putData(self):
841 841
842 842 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
843 843 self.closeFile()
844 844 self.setNextFile()
845 845
846 846 line = ''
847 847 for j in range(len(self.dataOut.heightList)):
848 848 for ds in self.dsList:
849 849 par = self.dataDict[ds['variable']]
850 850 if ds['nDim'] == 2:
851 851 for i in range(len(par)):
852 852 line += '{:>16}'.format('%8.2f' % getattr(self.dataOut, ds['variable'])[i][j])
853 853 elif ds['nDim'] == 1:
854 854 line += '{:>16}'.format('%8.2f' % getattr(self.dataOut, ds['variable'])[j])
855 855 else:
856 856 line += '{:>16}'.format('%8.2f' % getattr(self.dataOut, ds['variable']))
857 857
858 858 line += '\n'
859 859 self.fp.write(line)
860 860
861 861 self.blockIndex += 1
862 862 log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
863 863
864 864 return
865 865
866 866 def closeFile(self):
867 867
868 868 if self.fp:
869 869 self.fp.close()
870 870
871 871 def close(self):
872 872
873 873 self.closeFile()
General Comments 0
You need to be logged in to leave comments. Login now