##// END OF EJS Templates
schain
RSS
timeZone verification before ltctime calculation
joabAM -
r1742:9af415573ac7
parent child
Show More
Side by Side Unified
Context file:
r1742:9af415573ac7 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1159 +1,1163
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 return self.utctime - self.timeZone * 60
248
247 if self.timeZone =='lt':
248 return self.utctime - 300 * 60
249 elif self.timeZone =='ut':
250 return self.utctime
251 else:
252 log.error("No valid timeZone detected")
249 253 return self.utctime
250 254
251 255 @property
252 256 def datatime(self):
253 257
254 258 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
255 259 return datatimeValue
256 260
257 261 def getTimeRange(self):
258 262
259 263 datatime = []
260 264
261 265 datatime.append(self.ltctime)
262 266 datatime.append(self.ltctime + self.timeInterval + 1)
263 267
264 268 datatime = numpy.array(datatime)
265 269
266 270 return datatime
267 271
268 272 def getFmaxTimeResponse(self):
269 273
270 274 period = (10**-6) * self.getDeltaH() / (0.15)
271 275
272 276 PRF = 1. / (period * self.nCohInt)
273 277
274 278 fmax = PRF
275 279
276 280 return fmax
277 281
278 282 def getFmax(self):
279 283 PRF = 1. / (self.__ippSeconds * self.nCohInt)
280 284
281 285 fmax = PRF
282 286 return fmax
283 287
284 288 def getVmax(self):
285 289
286 290 _lambda = self.C / self.frequency
287 291
288 292 vmax = self.getFmax() * _lambda / 2
289 293
290 294 return vmax
291 295
292 296 ## Radar Controller Header must be immutable
293 297 @property
294 298 def ippSeconds(self):
295 299 '''
296 300 '''
297 301 #return self.radarControllerHeaderObj.ippSeconds
298 302 return self.__ippSeconds
299 303
300 304 @ippSeconds.setter
301 305 def ippSeconds(self, ippSeconds):
302 306 '''
303 307 '''
304 308 #self.radarControllerHeaderObj.ippSeconds = ippSeconds
305 309 self.__ippSeconds = ippSeconds
306 310 self.__ipp = ippSeconds*SPEED_OF_LIGHT/2000.0
307 311
308 312 @property
309 313 def code(self):
310 314 '''
311 315 '''
312 316 # return self.radarControllerHeaderObj.code
313 317 return self.__code
314 318
315 319 @code.setter
316 320 def code(self, code):
317 321 '''
318 322 '''
319 323 # self.radarControllerHeaderObj.code = code
320 324 self.__code = code
321 325
322 326 @property
323 327 def nCode(self):
324 328 '''
325 329 '''
326 330 # return self.radarControllerHeaderObj.nCode
327 331 return self.__nCode
328 332
329 333 @nCode.setter
330 334 def nCode(self, ncode):
331 335 '''
332 336 '''
333 337 # self.radarControllerHeaderObj.nCode = ncode
334 338 self.__nCode = ncode
335 339
336 340 @property
337 341 def nBaud(self):
338 342 '''
339 343 '''
340 344 # return self.radarControllerHeaderObj.nBaud
341 345 return self.__nBaud
342 346
343 347 @nBaud.setter
344 348 def nBaud(self, nbaud):
345 349 '''
346 350 '''
347 351 # self.radarControllerHeaderObj.nBaud = nbaud
348 352 self.__nBaud = nbaud
349 353
350 354 @property
351 355 def ipp(self):
352 356 '''
353 357 '''
354 358 # return self.radarControllerHeaderObj.ipp
355 359 return self.__ipp
356 360
357 361 @ipp.setter
358 362 def ipp(self, ipp):
359 363 '''
360 364 '''
361 365 # self.radarControllerHeaderObj.ipp = ipp
362 366 self.__ipp = ipp
363 367
364 368 @property
365 369 def metadata(self):
366 370 '''
367 371 '''
368 372
369 373 return {attr: getattr(self, attr) for attr in self.metadata_list}
370 374
371 375
372 376 class Voltage(JROData):
373 377
374 378 dataPP_POW = None
375 379 dataPP_DOP = None
376 380 dataPP_WIDTH = None
377 381 dataPP_SNR = None
378 382
379 383 # To use oper
380 384 flagProfilesByRange = False
381 385 nProfilesByRange = None
382 386 max_nIncohInt = 1
383 387
384 388 def __init__(self):
385 389 '''
386 390 Constructor
387 391 '''
388 392
389 393 self.useLocalTime = True
390 394 self.radarControllerHeaderObj = RadarControllerHeader()
391 395 self.systemHeaderObj = SystemHeader()
392 396 self.processingHeaderObj = ProcessingHeader()
393 397 self.type = "Voltage"
394 398 self.data = None
395 399 self.nProfiles = None
396 400 self.heightList = None
397 401 self.channelList = None
398 402 self.flagNoData = True
399 403 self.flagDiscontinuousBlock = False
400 404 self.utctime = None
401 405 self.timeZone = 0
402 406 self.dstFlag = None
403 407 self.errorCount = None
404 408 self.nCohInt = None
405 409 self.blocksize = None
406 410 self.flagCohInt = False
407 411 self.flagDecodeData = False # asumo q la data no esta decodificada
408 412 self.flagDeflipData = False # asumo q la data no esta sin flip
409 413 self.flagShiftFFT = False
410 414 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
411 415 self.profileIndex = 0
412 416 self.ippFactor=1
413 417 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
414 418 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
415 419
416 420 def getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None):
417 421 """
418 422 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
419 423
420 424 Return:
421 425 noiselevel
422 426 """
423 427
424 428 if channel != None:
425 429 data = self.data[channel,ymin_index:ymax_index]
426 430 nChannels = 1
427 431 else:
428 432 data = self.data[:,ymin_index:ymax_index]
429 433 nChannels = self.nChannels
430 434
431 435 noise = numpy.zeros(nChannels)
432 436 power = data * numpy.conjugate(data)
433 437
434 438 for thisChannel in range(nChannels):
435 439 if nChannels == 1:
436 440 daux = power[:].real
437 441 else:
438 442 daux = power[thisChannel, :].real
439 443 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
440 444
441 445 return noise
442 446
443 447 def getNoise(self, type=1, channel=None,ymin_index=None, ymax_index=None):
444 448
445 449 if type == 1:
446 450 noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index)
447 451
448 452 return noise
449 453
450 454 def getPower(self, channel=None):
451 455
452 456 if channel != None:
453 457 data = self.data[channel]
454 458 else:
455 459 data = self.data
456 460
457 461 power = data * numpy.conjugate(data)
458 462 powerdB = 10 * numpy.log10(power.real)
459 463 powerdB = numpy.squeeze(powerdB)
460 464
461 465 return powerdB
462 466 @property
463 467 def data_pow(self):
464 468 return self.getPower()
465 469
466 470 @property
467 471 def timeInterval(self):
468 472
469 473 return self.ippSeconds * self.nCohInt
470 474
471 475 noise = property(getNoise, "I'm the 'nHeights' property.")
472 476
473 477
474 478 class Spectra(JROData):
475 479
476 480 data_outlier = None
477 481 flagProfilesByRange = False
478 482 nProfilesByRange = None
479 483
480 484 def __init__(self):
481 485 '''
482 486 Constructor
483 487 '''
484 488
485 489 self.data_dc = None
486 490 self.data_spc = None
487 491 self.data_cspc = None
488 492 self.useLocalTime = True
489 493 self.radarControllerHeaderObj = RadarControllerHeader()
490 494 self.systemHeaderObj = SystemHeader()
491 495 self.processingHeaderObj = ProcessingHeader()
492 496 self.type = "Spectra"
493 497 self.timeZone = 0
494 498 self.nProfiles = None
495 499 self.heightList = None
496 500 self.channelList = None
497 501 self.pairsList = None
498 502 self.flagNoData = True
499 503 self.flagDiscontinuousBlock = False
500 504 self.utctime = None
501 505 self.nCohInt = None
502 506 self.nIncohInt = None
503 507 self.blocksize = None
504 508 self.nFFTPoints = None
505 509 self.wavelength = None
506 510 self.flagDecodeData = False # asumo q la data no esta decodificada
507 511 self.flagDeflipData = False # asumo q la data no esta sin flip
508 512 self.flagShiftFFT = False
509 513 self.ippFactor = 1
510 514 self.beacon_heiIndexList = []
511 515 self.noise_estimation = None
512 516 self.codeList = []
513 517 self.azimuthList = []
514 518 self.elevationList = []
515 519 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
516 520 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
517 521
518 522 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
519 523 """
520 524 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
521 525
522 526 Return:
523 527 noiselevel
524 528 """
525 529
526 530 noise = numpy.zeros(self.nChannels)
527 531
528 532 for channel in range(self.nChannels):
529 533 daux = self.data_spc[channel,
530 534 xmin_index:xmax_index, ymin_index:ymax_index]
531 535 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
532 536 noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt[channel])
533 537
534 538 return noise
535 539
536 540 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
537 541
538 542 if self.noise_estimation is not None:
539 543 # this was estimated by getNoise Operation defined in jroproc_spectra.py
540 544 return self.noise_estimation
541 545 else:
542 546 noise = self.getNoisebyHildebrand(
543 547 xmin_index, xmax_index, ymin_index, ymax_index)
544 548 return noise
545 549
546 550 def getFreqRangeTimeResponse(self, extrapoints=0):
547 551
548 552 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
549 553 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
550 554
551 555 return freqrange
552 556
553 557 def getAcfRange(self, extrapoints=0):
554 558
555 559 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
556 560 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
557 561
558 562 return freqrange
559 563
560 564 def getFreqRange(self, extrapoints=0):
561 565
562 566 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
563 567 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
564 568
565 569 return freqrange
566 570
567 571 def getVelRange(self, extrapoints=0):
568 572
569 573 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
570 574 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
571 575
572 576 if self.nmodes:
573 577 return velrange/self.nmodes
574 578 else:
575 579 return velrange
576 580
577 581 @property
578 582 def nPairs(self):
579 583
580 584 return len(self.pairsList)
581 585
582 586 @property
583 587 def pairsIndexList(self):
584 588
585 589 return list(range(self.nPairs))
586 590
587 591 @property
588 592 def normFactor(self):
589 593
590 594 pwcode = 1
591 595 if self.flagDecodeData:
592 596 try:
593 597 pwcode = numpy.sum(self.code[0]**2)
594 598 except Exception as e:
595 599 log.warning("Failed pwcode read, setting to 1")
596 600 pwcode = 1
597 601 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
598 602 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
599 603 if self.flagProfilesByRange:
600 604 normFactor *= (self.nProfilesByRange/self.nProfilesByRange.max())
601 605 return normFactor
602 606
603 607 @property
604 608 def flag_cspc(self):
605 609
606 610 if self.data_cspc is None:
607 611 return True
608 612
609 613 return False
610 614
611 615 @property
612 616 def flag_dc(self):
613 617
614 618 if self.data_dc is None:
615 619 return True
616 620
617 621 return False
618 622
619 623 @property
620 624 def timeInterval(self):
621 625
622 626 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
623 627 if self.nmodes:
624 628 return self.nmodes*timeInterval
625 629 else:
626 630 return timeInterval
627 631
628 632 def getPower(self):
629 633
630 634 factor = self.normFactor
631 635 power = numpy.zeros( (self.nChannels,self.nHeights) )
632 636 for ch in range(self.nChannels):
633 637 z = None
634 638 if hasattr(factor,'shape'):
635 639 if factor.ndim > 1:
636 640 z = self.data_spc[ch]/factor[ch]
637 641 else:
638 642 z = self.data_spc[ch]/factor
639 643 else:
640 644 z = self.data_spc[ch]/factor
641 645 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
642 646 avg = numpy.average(z, axis=0)
643 647 power[ch] = 10 * numpy.log10(avg)
644 648 return power
645 649
646 650 @property
647 651 def max_nIncohInt(self):
648 652
649 653 ints = numpy.zeros(self.nChannels)
650 654 for ch in range(self.nChannels):
651 655 if hasattr(self.nIncohInt,'shape'):
652 656 if self.nIncohInt.ndim > 1:
653 657 ints[ch,] = self.nIncohInt[ch].max()
654 658 else:
655 659 ints[ch,] = self.nIncohInt
656 660 self.nIncohInt = int(self.nIncohInt)
657 661 else:
658 662 ints[ch,] = self.nIncohInt
659 663
660 664 return ints
661 665
662 666 def getCoherence(self, pairsList=None, phase=False):
663 667
664 668 z = []
665 669 if pairsList is None:
666 670 pairsIndexList = self.pairsIndexList
667 671 else:
668 672 pairsIndexList = []
669 673 for pair in pairsList:
670 674 if pair not in self.pairsList:
671 675 raise ValueError("Pair %s is not in dataOut.pairsList" % (
672 676 pair))
673 677 pairsIndexList.append(self.pairsList.index(pair))
674 678 for i in range(len(pairsIndexList)):
675 679 pair = self.pairsList[pairsIndexList[i]]
676 680 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
677 681 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
678 682 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
679 683 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
680 684 if phase:
681 685 data = numpy.arctan2(avgcoherenceComplex.imag,
682 686 avgcoherenceComplex.real) * 180 / numpy.pi
683 687 else:
684 688 data = numpy.abs(avgcoherenceComplex)
685 689
686 690 z.append(data)
687 691
688 692 return numpy.array(z)
689 693
690 694 def setValue(self, value):
691 695
692 696 print("This property should not be initialized", value)
693 697
694 698 return
695 699
696 700 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
697 701
698 702
699 703 class SpectraHeis(Spectra):
700 704
701 705 def __init__(self):
702 706
703 707 self.radarControllerHeaderObj = RadarControllerHeader()
704 708 self.systemHeaderObj = SystemHeader()
705 709 self.type = "SpectraHeis"
706 710 self.nProfiles = None
707 711 self.heightList = None
708 712 self.channelList = None
709 713 self.flagNoData = True
710 714 self.flagDiscontinuousBlock = False
711 715 self.utctime = None
712 716 self.blocksize = None
713 717 self.profileIndex = 0
714 718 self.nCohInt = 1
715 719 self.nIncohInt = 1
716 720
717 721 @property
718 722 def normFactor(self):
719 723 pwcode = 1
720 724 if self.flagDecodeData:
721 725 pwcode = numpy.sum(self.code[0]**2)
722 726
723 727 normFactor = self.nIncohInt * self.nCohInt * pwcode
724 728
725 729 return normFactor
726 730
727 731 @property
728 732 def timeInterval(self):
729 733
730 734 return self.ippSeconds * self.nCohInt * self.nIncohInt
731 735
732 736
733 737 class Fits(JROData):
734 738
735 739 def __init__(self):
736 740
737 741 self.type = "Fits"
738 742 self.nProfiles = None
739 743 self.heightList = None
740 744 self.channelList = None
741 745 self.flagNoData = True
742 746 self.utctime = None
743 747 self.nCohInt = 1
744 748 self.nIncohInt = 1
745 749 self.useLocalTime = True
746 750 self.profileIndex = 0
747 751 self.timeZone = 0
748 752
749 753 def getTimeRange(self):
750 754
751 755 datatime = []
752 756
753 757 datatime.append(self.ltctime)
754 758 datatime.append(self.ltctime + self.timeInterval)
755 759
756 760 datatime = numpy.array(datatime)
757 761
758 762 return datatime
759 763
760 764 def getChannelIndexList(self):
761 765
762 766 return list(range(self.nChannels))
763 767
764 768 def getNoise(self, type=1):
765 769
766 770
767 771 if type == 1:
768 772 noise = self.getNoisebyHildebrand()
769 773
770 774 if type == 2:
771 775 noise = self.getNoisebySort()
772 776
773 777 if type == 3:
774 778 noise = self.getNoisebyWindow()
775 779
776 780 return noise
777 781
778 782 @property
779 783 def timeInterval(self):
780 784
781 785 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
782 786
783 787 return timeInterval
784 788
785 789 @property
786 790 def ippSeconds(self):
787 791 '''
788 792 '''
789 793 return self.ipp_sec
790 794
791 795 noise = property(getNoise, "I'm the 'nHeights' property.")
792 796
793 797
794 798 class Correlation(JROData):
795 799
796 800 def __init__(self):
797 801 '''
798 802 Constructor
799 803 '''
800 804 self.radarControllerHeaderObj = RadarControllerHeader()
801 805 self.systemHeaderObj = SystemHeader()
802 806 self.type = "Correlation"
803 807 self.data = None
804 808 self.dtype = None
805 809 self.nProfiles = None
806 810 self.heightList = None
807 811 self.channelList = None
808 812 self.flagNoData = True
809 813 self.flagDiscontinuousBlock = False
810 814 self.utctime = None
811 815 self.timeZone = 0
812 816 self.dstFlag = None
813 817 self.errorCount = None
814 818 self.blocksize = None
815 819 self.flagDecodeData = False # asumo q la data no esta decodificada
816 820 self.flagDeflipData = False # asumo q la data no esta sin flip
817 821 self.pairsList = None
818 822 self.nPoints = None
819 823
820 824 def getPairsList(self):
821 825
822 826 return self.pairsList
823 827
824 828 def getNoise(self, mode=2):
825 829
826 830 indR = numpy.where(self.lagR == 0)[0][0]
827 831 indT = numpy.where(self.lagT == 0)[0][0]
828 832
829 833 jspectra0 = self.data_corr[:, :, indR, :]
830 834 jspectra = copy.copy(jspectra0)
831 835
832 836 num_chan = jspectra.shape[0]
833 837 num_hei = jspectra.shape[2]
834 838
835 839 freq_dc = jspectra.shape[1] / 2
836 840 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
837 841
838 842 if ind_vel[0] < 0:
839 843 ind_vel[list(range(0, 1))] = ind_vel[list(
840 844 range(0, 1))] + self.num_prof
841 845
842 846 if mode == 1:
843 847 jspectra[:, freq_dc, :] = (
844 848 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
845 849
846 850 if mode == 2:
847 851
848 852 vel = numpy.array([-2, -1, 1, 2])
849 853 xx = numpy.zeros([4, 4])
850 854
851 855 for fil in range(4):
852 856 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
853 857
854 858 xx_inv = numpy.linalg.inv(xx)
855 859 xx_aux = xx_inv[0, :]
856 860
857 861 for ich in range(num_chan):
858 862 yy = jspectra[ich, ind_vel, :]
859 863 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
860 864
861 865 junkid = jspectra[ich, freq_dc, :] <= 0
862 866 cjunkid = sum(junkid)
863 867
864 868 if cjunkid.any():
865 869 jspectra[ich, freq_dc, junkid.nonzero()] = (
866 870 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
867 871
868 872 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
869 873
870 874 return noise
871 875
872 876 @property
873 877 def timeInterval(self):
874 878
875 879 return self.ippSeconds * self.nCohInt * self.nProfiles
876 880
877 881 def splitFunctions(self):
878 882
879 883 pairsList = self.pairsList
880 884 ccf_pairs = []
881 885 acf_pairs = []
882 886 ccf_ind = []
883 887 acf_ind = []
884 888 for l in range(len(pairsList)):
885 889 chan0 = pairsList[l][0]
886 890 chan1 = pairsList[l][1]
887 891
888 892 # Obteniendo pares de Autocorrelacion
889 893 if chan0 == chan1:
890 894 acf_pairs.append(chan0)
891 895 acf_ind.append(l)
892 896 else:
893 897 ccf_pairs.append(pairsList[l])
894 898 ccf_ind.append(l)
895 899
896 900 data_acf = self.data_cf[acf_ind]
897 901 data_ccf = self.data_cf[ccf_ind]
898 902
899 903 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
900 904
901 905 @property
902 906 def normFactor(self):
903 907 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
904 908 acf_pairs = numpy.array(acf_pairs)
905 909 normFactor = numpy.zeros((self.nPairs, self.nHeights))
906 910
907 911 for p in range(self.nPairs):
908 912 pair = self.pairsList[p]
909 913
910 914 ch0 = pair[0]
911 915 ch1 = pair[1]
912 916
913 917 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
914 918 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
915 919 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
916 920
917 921 return normFactor
918 922
919 923
920 924 class Parameters(Spectra):
921 925
922 926 groupList = None # List of Pairs, Groups, etc
923 927 data_param = None # Parameters obtained
924 928 data_pre = None # Data Pre Parametrization
925 929 data_SNR = None # Signal to Noise Ratio
926 930 abscissaList = None # Abscissa, can be velocities, lags or time
927 931 utctimeInit = None # Initial UTC time
928 932 paramInterval = None # Time interval to calculate Parameters in seconds
929 933 useLocalTime = True
930 934 # Fitting
931 935 data_error = None # Error of the estimation
932 936 constants = None
933 937 library = None
934 938 # Output signal
935 939 outputInterval = None # Time interval to calculate output signal in seconds
936 940 data_output = None # Out signal
937 941 nAvg = None
938 942 noise_estimation = None
939 943 GauSPC = None # Fit gaussian SPC
940 944
941 945 data_outlier = None
942 946 data_vdrift = None
943 947 radarControllerHeaderTxt=None #header Controller like text
944 948 txPower = None
945 949 flagProfilesByRange = False
946 950 nProfilesByRange = None
947 951
948 952
949 953 def __init__(self):
950 954 '''
951 955 Constructor
952 956 '''
953 957 self.radarControllerHeaderObj = RadarControllerHeader()
954 958 self.systemHeaderObj = SystemHeader()
955 959 self.processingHeaderObj = ProcessingHeader()
956 960 self.type = "Parameters"
957 961 self.timeZone = 0
958 962
959 963 def getTimeRange1(self, interval):
960 964
961 965 datatime = []
962 966
963 967 if self.useLocalTime:
964 968 time1 = self.utctimeInit - self.timeZone * 60
965 969 else:
966 970 time1 = self.utctimeInit
967 971
968 972 datatime.append(time1)
969 973 datatime.append(time1 + interval)
970 974 datatime = numpy.array(datatime)
971 975
972 976 return datatime
973 977
974 978 @property
975 979 def timeInterval(self):
976 980
977 981 if hasattr(self, 'timeInterval1'):
978 982 return self.timeInterval1
979 983 else:
980 984 return self.paramInterval
981 985
982 986 def setValue(self, value):
983 987
984 988 print("This property should not be initialized")
985 989
986 990 return
987 991
988 992 def getNoise(self):
989 993
990 994 return self.spc_noise
991 995
992 996 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
993 997
994 998
995 999 class PlotterData(object):
996 1000 '''
997 1001 Object to hold data to be plotted
998 1002 '''
999 1003
1000 1004 MAXNUMX = 200
1001 1005 MAXNUMY = 200
1002 1006
1003 1007 def __init__(self, code, exp_code, localtime=True):
1004 1008
1005 1009 self.key = code
1006 1010 self.exp_code = exp_code
1007 1011 self.ready = False
1008 1012 self.flagNoData = False
1009 1013 self.localtime = localtime
1010 1014 self.data = {}
1011 1015 self.meta = {}
1012 1016 self.__heights = []
1013 1017
1014 1018 def __str__(self):
1015 1019 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1016 1020 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1017 1021
1018 1022 def __len__(self):
1019 1023 return len(self.data)
1020 1024
1021 1025 def __getitem__(self, key):
1022 1026 if isinstance(key, int):
1023 1027 return self.data[self.times[key]]
1024 1028 elif isinstance(key, str):
1025 1029 ret = numpy.array([self.data[x][key] for x in self.times])
1026 1030 if ret.ndim > 1:
1027 1031 ret = numpy.swapaxes(ret, 0, 1)
1028 1032 return ret
1029 1033
1030 1034 def __contains__(self, key):
1031 1035 return key in self.data[self.min_time]
1032 1036
1033 1037 def setup(self):
1034 1038 '''
1035 1039 Configure object
1036 1040 '''
1037 1041 self.type = ''
1038 1042 self.ready = False
1039 1043 del self.data
1040 1044 self.data = {}
1041 1045 self.__heights = []
1042 1046 self.__all_heights = set()
1043 1047
1044 1048 def shape(self, key):
1045 1049 '''
1046 1050 Get the shape of the one-element data for the given key
1047 1051 '''
1048 1052
1049 1053 if len(self.data[self.min_time][key]):
1050 1054 return self.data[self.min_time][key].shape
1051 1055 return (0,)
1052 1056
1053 1057 def update(self, data, tm, meta={}):
1054 1058 '''
1055 1059 Update data object with new dataOut
1056 1060 '''
1057 1061
1058 1062 self.data[tm] = data
1059 1063
1060 1064 for key, value in meta.items():
1061 1065 setattr(self, key, value)
1062 1066
1063 1067 def normalize_heights(self):
1064 1068 '''
1065 1069 Ensure same-dimension of the data for different heighList
1066 1070 '''
1067 1071
1068 1072 H = numpy.array(list(self.__all_heights))
1069 1073 H.sort()
1070 1074 for key in self.data:
1071 1075 shape = self.shape(key)[:-1] + H.shape
1072 1076 for tm, obj in list(self.data[key].items()):
1073 1077 h = self.__heights[self.times.tolist().index(tm)]
1074 1078 if H.size == h.size:
1075 1079 continue
1076 1080 index = numpy.where(numpy.in1d(H, h))[0]
1077 1081 dummy = numpy.zeros(shape) + numpy.nan
1078 1082 if len(shape) == 2:
1079 1083 dummy[:, index] = obj
1080 1084 else:
1081 1085 dummy[index] = obj
1082 1086 self.data[key][tm] = dummy
1083 1087
1084 1088 self.__heights = [H for tm in self.times]
1085 1089
1086 1090 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1087 1091 '''
1088 1092 Convert data to json
1089 1093 '''
1090 1094
1091 1095 meta = {}
1092 1096 meta['xrange'] = []
1093 1097 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1094 1098 tmp = self.data[tm][self.key]
1095 1099 shape = tmp.shape
1096 1100 if len(shape) == 2:
1097 1101 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1098 1102 elif len(shape) == 3:
1099 1103 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1100 1104 data = self.roundFloats(
1101 1105 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1102 1106 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1103 1107 else:
1104 1108 data = self.roundFloats(self.data[tm][self.key].tolist())
1105 1109
1106 1110 ret = {
1107 1111 'plot': plot_name,
1108 1112 'code': self.exp_code,
1109 1113 'time': float(tm),
1110 1114 'data': data,
1111 1115 }
1112 1116 meta['type'] = plot_type
1113 1117 meta['interval'] = float(self.interval)
1114 1118 meta['localtime'] = self.localtime
1115 1119 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1116 1120 meta.update(self.meta)
1117 1121 ret['metadata'] = meta
1118 1122 return json.dumps(ret)
1119 1123
1120 1124 @property
1121 1125 def times(self):
1122 1126 '''
1123 1127 Return the list of times of the current data
1124 1128 '''
1125 1129
1126 1130 ret = [t for t in self.data]
1127 1131 ret.sort()
1128 1132 return numpy.array(ret)
1129 1133
1130 1134 @property
1131 1135 def min_time(self):
1132 1136 '''
1133 1137 Return the minimun time value
1134 1138 '''
1135 1139
1136 1140 return self.times[0]
1137 1141
1138 1142 @property
1139 1143 def max_time(self):
1140 1144 '''
1141 1145 Return the maximun time value
1142 1146 '''
1143 1147
1144 1148 return self.times[-1]
1145 1149
1146 1150 # @property
1147 1151 # def heights(self):
1148 1152 # '''
1149 1153 # Return the list of heights of the current data
1150 1154 # '''
1151 1155
1152 1156 # return numpy.array(self.__heights[-1])
1153 1157
1154 1158 @staticmethod
1155 1159 def roundFloats(obj):
1156 1160 if isinstance(obj, list):
1157 1161 return list(map(PlotterData.roundFloats, obj))
1158 1162 elif isinstance(obj, float):
1159 1163 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,779 +1,779
1 1 ''''
2 2 Created on Set 9, 2015
3 3
4 4 @author: roj-idl71 Karim Kuyeng
5 5
6 6 @upgrade: 2021, Joab Apaza
7 7
8 8 '''
9 9
10 10 import os
11 11 import sys
12 12 import glob
13 13 import fnmatch
14 14 import datetime
15 15 import time
16 16 import re
17 17 import h5py
18 18 import numpy
19 19
20 20 try:
21 21 from gevent import sleep
22 22 except:
23 23 from time import sleep
24 24
25 25 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader,ProcessingHeader
26 26 from schainpy.model.data.jrodata import Voltage
27 27 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
28 28 from numpy import imag
29 29 from schainpy.utils import log
30 30
31 31
32 32 class AMISRReader(ProcessingUnit):
33 33 '''
34 34 classdocs
35 35 '''
36 36
37 37 def __init__(self):
38 38 '''
39 39 Constructor
40 40 '''
41 41
42 42 ProcessingUnit.__init__(self)
43 43
44 44 self.set = None
45 45 self.subset = None
46 46 self.extension_file = '.h5'
47 47 self.dtc_str = 'dtc'
48 48 self.dtc_id = 0
49 49 self.status = True
50 50 self.isConfig = False
51 51 self.dirnameList = []
52 52 self.filenameList = []
53 53 self.fileIndex = None
54 54 self.flagNoMoreFiles = False
55 55 self.flagIsNewFile = 0
56 56 self.filename = ''
57 57 self.amisrFilePointer = None
58 58
59 59 self.beamCodeMap = None
60 60 self.azimuthList = []
61 61 self.elevationList = []
62 62 self.dataShape = None
63 63 self.flag_old_beams = False
64 64
65 65 self.flagAsync = False #Use when the experiment has no syncronization
66 66 self.shiftChannels = 0
67 67 self.profileIndex = 0
68 68
69 69
70 70 self.beamCodeByFrame = None
71 71 self.radacTimeByFrame = None
72 72
73 73 self.dataset = None
74 74
75 75 self.__firstFile = True
76 76
77 77 self.buffer = None
78 78
79 79 self.timezone = 'ut'
80 80
81 81 self.__waitForNewFile = 20
82 82 self.__filename_online = None
83 83 #Is really necessary create the output object in the initializer
84 84 self.dataOut = Voltage()
85 85 self.dataOut.error=False
86 86 self.margin_days = 1
87 87 self.flag_ignoreFiles = False #to activate the ignoring Files flag
88 88 self.flag_standby = False # just keep waiting, use when ignoring files
89 89 self.ignStartDateTime=None
90 90 self.ignEndDateTime=None
91 91
92 92 def setup(self,path=None,
93 93 startDate=None,
94 94 endDate=None,
95 95 startTime=None,
96 96 endTime=None,
97 97 walk=True,
98 98 timezone='ut',
99 99 all=0,
100 100 code = 1,
101 101 nCode = 1,
102 102 nBaud = 0,
103 103 nOsamp = 0,
104 104 online=False,
105 105 old_beams=False,
106 106 margin_days=1,
107 107 nFFT = None,
108 108 nChannels = None,
109 109 ignStartDate=None,
110 110 ignEndDate=None,
111 111 ignStartTime=None,
112 112 ignEndTime=None,
113 113 syncronization=True,
114 114 shiftChannels=0
115 115 ):
116 116
117 117
118 118
119 119 self.timezone = timezone
120 120 self.all = all
121 121 self.online = online
122 122 self.flag_old_beams = old_beams
123 123 self.code = code
124 124 self.nCode = int(nCode)
125 125 self.nBaud = int(nBaud)
126 126 self.nOsamp = int(nOsamp)
127 127 self.margin_days = margin_days
128 128 self.__sampleRate = None
129 129 self.flagAsync = not syncronization
130 130 self.shiftChannels = shiftChannels
131 131 self.nFFT = nFFT
132 132 self.nChannels = nChannels
133 133 if ignStartTime!=None and ignEndTime!=None:
134 134 if ignStartDate!=None and ignEndDate!=None:
135 135 self.ignStartDateTime=datetime.datetime.combine(ignStartDate,ignStartTime)
136 136 self.ignEndDateTime=datetime.datetime.combine(ignEndDate,ignEndTime)
137 137 else:
138 138 self.ignStartDateTime=datetime.datetime.combine(startDate,ignStartTime)
139 139 self.ignEndDateTime=datetime.datetime.combine(endDate,ignEndTime)
140 140 self.flag_ignoreFiles = True
141 141
142 142 #self.findFiles()
143 143 if not(online):
144 144 #Busqueda de archivos offline
145 145 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk,)
146 146 else:
147 147 self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
148 148
149 149 if not(self.filenameList):
150 150 raise schainpy.admin.SchainWarning("There is no files into the folder: %s"%(path))
151 151 #sys.exit(0)
152 152 self.dataOut.error = True
153 153
154 154 self.fileIndex = 0
155 155
156 156 self.readNextFile(online)
157 157
158 158 '''
159 159 Add code
160 160 '''
161 161 self.isConfig = True
162 162 # print("Setup Done")
163 163 pass
164 164
165 165
166 166 def readAMISRHeader(self,fp):
167 167
168 168 if self.isConfig and (not self.flagNoMoreFiles):
169 169 newShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
170 170 if self.dataShape != newShape and newShape != None and not self.flag_standby:
171 171 raise schainpy.admin.SchainError("NEW FILE HAS A DIFFERENT SHAPE: ")
172 172 print(self.dataShape,newShape,"\n")
173 173 return 0
174 174 else:
175 175 self.dataShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
176 176
177 177
178 178 header = 'Raw11/Data/RadacHeader'
179 179 if self.nChannels == None:
180 180 expFile = fp['Setup/Experimentfile'][()].decode()
181 181 linesExp = expFile.split("\n")
182 182 a = [line for line in linesExp if "nbeamcodes" in line]
183 183 self.nChannels = int(a[0][11:])
184 184
185 185 if not self.flagAsync: #for experiments with no syncronization
186 186 self.shiftChannels = 0
187 187
188 188
189 189
190 190 self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
191 191
192 192
193 193 if (self.startDate > datetime.date(2021, 7, 15)) or self.flag_old_beams: #Se cambiΓ³ la forma de extracciΓ³n de Apuntes el 17 o forzar con flag de reorganizaciΓ³n
194 194 self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode()
195 195 self.trueBeams = self.beamcodeFile.split("\n")
196 196 self.trueBeams.pop()#remove last
197 197 if self.nFFT == None:
198 198 log.error("FFT or number of repetitions per channels is needed",self.name)
199 199 beams_idx = [k*self.nFFT for k in range(self.nChannels)]
200 200 beams = [self.trueBeams[b] for b in beams_idx]
201 201 self.beamCode = [int(x, 16) for x in beams]
202 202
203 203 if(self.flagAsync and self.shiftChannels == 0):
204 204 initBeam = self.beamCodeByPulse[0, 0]
205 205 self.shiftChannels = numpy.argwhere(self.beamCode ==initBeam)[0,0]
206 206
207 207 else:
208 208 _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes
209 209 self.beamCode = _beamCode[0,:]
210 210
211 211
212 212
213 213
214 214 if self.beamCodeMap == None:
215 215 self.beamCodeMap = fp['Setup/BeamcodeMap']
216 216 for beam in self.beamCode:
217 217 beamAziElev = numpy.where(self.beamCodeMap[:,0]==beam)
218 218 beamAziElev = beamAziElev[0].squeeze()
219 219 self.azimuthList.append(self.beamCodeMap[beamAziElev,1])
220 220 self.elevationList.append(self.beamCodeMap[beamAziElev,2])
221 221 #print("Beamssss: ",self.beamCodeMap[beamAziElev,1],self.beamCodeMap[beamAziElev,2])
222 222 #print(self.beamCode)
223 223 #self.code = fp.get(header+'/Code') # NOT USE FOR THIS
224 224 self.frameCount = fp.get(header+'/FrameCount')# NOT USE FOR THIS
225 225 self.modeGroup = fp.get(header+'/ModeGroup')# NOT USE FOR THIS
226 226 self.nsamplesPulse = fp.get(header+'/NSamplesPulse')# TO GET NSA OR USING DATA FOR THAT
227 227 self.pulseCount = fp.get(header+'/PulseCount')# NOT USE FOR THIS
228 228 self.radacTime = fp.get(header+'/RadacTime')# 1st TIME ON FILE ANDE CALCULATE THE REST WITH IPP*nindexprofile
229 229 self.timeCount = fp.get(header+'/TimeCount')# NOT USE FOR THIS
230 230 self.timeStatus = fp.get(header+'/TimeStatus')# NOT USE FOR THIS
231 231 self.rangeFromFile = fp.get('Raw11/Data/Samples/Range')
232 232 self.frequency = fp.get('Rx/Frequency')
233 233 txAus = fp.get('Raw11/Data/Pulsewidth') #seconds
234 234 self.baud = fp.get('Raw11/Data/TxBaud')
235 235 sampleRate = fp.get('Rx/SampleRate')
236 236 self.__sampleRate = sampleRate[()]
237 237 self.nblocks = self.pulseCount.shape[0] #nblocks
238 238 self.profPerBlockRAW = self.pulseCount.shape[1] #profiles per block in raw data
239 239 self.nprofiles = self.pulseCount.shape[1] #nprofile
240 240 #self.nsa = self.nsamplesPulse[0,0] #ngates
241 241 self.nsa = len(self.rangeFromFile[0])
242 242 self.nchannels = len(self.beamCode)
243 243 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
244 244 #print("IPPS secs: ",self.ippSeconds)
245 245 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
246 246 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
247 247
248 248 #filling radar controller header parameters
249 249 self.__ippKm = self.ippSeconds *.15*1e6 # in km
250 250 #self.__txA = txAus[()]*.15 #(ipp[us]*.15km/1us) in km
251 251 self.__txA = txAus[()] #seconds
252 252 self.__txAKm = self.__txA*1e6*.15
253 253 self.__txB = 0
254 254 nWindows=1
255 255 self.__nSamples = self.nsa
256 256 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
257 257 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
258 258 #print("amisr-ipp:",self.ippSeconds, self.__ippKm)
259 259 #for now until understand why the code saved is different (code included even though code not in tuf file)
260 260 #self.__codeType = 0
261 261 # self.__nCode = None
262 262 # self.__nBaud = None
263 263 self.__code = self.code
264 264 self.__codeType = 0
265 265 if self.code != None:
266 266 self.__codeType = 1
267 267 self.__nCode = self.nCode
268 268 self.__nBaud = self.nBaud
269 269 self.__baudTX = self.__txA/(self.nBaud)
270 270 #self.__code = 0
271 271
272 272 #filling system header parameters
273 273 self.__nSamples = self.nsa
274 274 self.newProfiles = self.nprofiles/self.nchannels
275 275 self.__channelList = [n for n in range(self.nchannels)]
276 276
277 277 self.__frequency = self.frequency[0][0]
278 278
279 279
280 280 return 1
281 281
282 282
283 283 def createBuffers(self):
284 284
285 285 pass
286 286
287 287 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
288 288 self.path = path
289 289 self.startDate = startDate
290 290 self.endDate = endDate
291 291 self.startTime = startTime
292 292 self.endTime = endTime
293 293 self.walk = walk
294 294
295 295
296 296 def __checkPath(self):
297 297 if os.path.exists(self.path):
298 298 self.status = 1
299 299 else:
300 300 self.status = 0
301 301 print('Path:%s does not exists'%self.path)
302 302
303 303 return
304 304
305 305
306 306 def __selDates(self, amisr_dirname_format):
307 307 try:
308 308 year = int(amisr_dirname_format[0:4])
309 309 month = int(amisr_dirname_format[4:6])
310 310 dom = int(amisr_dirname_format[6:8])
311 311 thisDate = datetime.date(year,month,dom)
312 312 #margen de un dΓ­a extra, igual luego se filtra for fecha y hora
313 313 if (thisDate>=(self.startDate - datetime.timedelta(days=self.margin_days)) and thisDate <= (self.endDate)+ datetime.timedelta(days=1)):
314 314 return amisr_dirname_format
315 315 except:
316 316 return None
317 317
318 318
319 319 def __findDataForDates(self,online=False):
320 320
321 321 if not(self.status):
322 322 return None
323 323
324 324 pat = '\d+.\d+'
325 325 dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
326 326 dirnameList = [x for x in dirnameList if x!=None]
327 327 dirnameList = [x.string for x in dirnameList]
328 328 if not(online):
329 329 dirnameList = [self.__selDates(x) for x in dirnameList]
330 330 dirnameList = [x for x in dirnameList if x!=None]
331 331 if len(dirnameList)>0:
332 332 self.status = 1
333 333 self.dirnameList = dirnameList
334 334 self.dirnameList.sort()
335 335 else:
336 336 self.status = 0
337 337 return None
338 338
339 339 def __getTimeFromData(self):
340 340 startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
341 341 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
342 342
343 343 print('Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader))
344 344 print('........................................')
345 345 filter_filenameList = []
346 346 self.filenameList.sort()
347 347 total_files = len(self.filenameList)
348 348 #for i in range(len(self.filenameList)-1):
349 349 for i in range(total_files):
350 350 filename = self.filenameList[i]
351 351 #print("file-> ",filename)
352 352 try:
353 353 fp = h5py.File(filename,'r')
354 354 time_str = fp.get('Time/RadacTimeString')
355 355
356 356 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
357 357 #startDateTimeStr_File = "2019-12-16 09:21:11"
358 358 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
359 359 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
360 360
361 361 #endDateTimeStr_File = "2019-12-16 11:10:11"
362 362 endDateTimeStr_File = time_str[-1][-1].decode('UTF-8').split('.')[0]
363 363 junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
364 364 endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
365 365
366 366 fp.close()
367 367
368 368 #print("check time", startDateTime_File)
369 369 if self.timezone == 'lt':
370 370 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
371 371 endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
372 372 if (startDateTime_File >=startDateTime_Reader and endDateTime_File<=endDateTime_Reader):
373 373 filter_filenameList.append(filename)
374 374
375 375 if (startDateTime_File>endDateTime_Reader):
376 376 break
377 377 except Exception as e:
378 378 log.warning("Error opening file {} -> {}".format(os.path.split(filename)[1],e))
379 379
380 380 filter_filenameList.sort()
381 381 self.filenameList = filter_filenameList
382 382
383 383 return 1
384 384
385 385 def __filterByGlob1(self, dirName):
386 386 filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
387 387 filter_files.sort()
388 388 filterDict = {}
389 389 filterDict.setdefault(dirName)
390 390 filterDict[dirName] = filter_files
391 391 return filterDict
392 392
393 393 def __getFilenameList(self, fileListInKeys, dirList):
394 394 for value in fileListInKeys:
395 395 dirName = list(value.keys())[0]
396 396 for file in value[dirName]:
397 397 filename = os.path.join(dirName, file)
398 398 self.filenameList.append(filename)
399 399
400 400
401 401 def __selectDataForTimes(self, online=False):
402 402 #aun no esta implementado el filtro for tiempo-> implementado en readNextFile
403 403 if not(self.status):
404 404 return None
405 405
406 406 dirList = [os.path.join(self.path,x) for x in self.dirnameList]
407 407 fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
408 408 self.__getFilenameList(fileListInKeys, dirList)
409 409 if not(online):
410 410 #filtro por tiempo
411 411 if not(self.all):
412 412 self.__getTimeFromData()
413 413
414 414 if len(self.filenameList)>0:
415 415 self.status = 1
416 416 self.filenameList.sort()
417 417 else:
418 418 self.status = 0
419 419 return None
420 420
421 421 else:
422 422 #get the last file - 1
423 423 self.filenameList = [self.filenameList[-2]]
424 424 new_dirnameList = []
425 425 for dirname in self.dirnameList:
426 426 junk = numpy.array([dirname in x for x in self.filenameList])
427 427 junk_sum = junk.sum()
428 428 if junk_sum > 0:
429 429 new_dirnameList.append(dirname)
430 430 self.dirnameList = new_dirnameList
431 431 return 1
432 432
433 433 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
434 434 endTime=datetime.time(23,59,59),walk=True):
435 435
436 436 if endDate ==None:
437 437 startDate = datetime.datetime.utcnow().date()
438 438 endDate = datetime.datetime.utcnow().date()
439 439
440 440 self.__setParameters(path=path, startDate=startDate, endDate=endDate,startTime = startTime,endTime=endTime, walk=walk)
441 441
442 442 self.__checkPath()
443 443
444 444 self.__findDataForDates(online=True)
445 445
446 446 self.dirnameList = [self.dirnameList[-1]]
447 447
448 448 self.__selectDataForTimes(online=True)
449 449
450 450 return
451 451
452 452
453 453 def searchFilesOffLine(self,
454 454 path,
455 455 startDate,
456 456 endDate,
457 457 startTime=datetime.time(0,0,0),
458 458 endTime=datetime.time(23,59,59),
459 459 walk=True):
460 460
461 461 self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
462 462
463 463 self.__checkPath()
464 464
465 465 self.__findDataForDates()
466 466
467 467 self.__selectDataForTimes()
468 468
469 469 for i in range(len(self.filenameList)):
470 470 print("%s" %(self.filenameList[i]))
471 471
472 472 return
473 473
474 474 def __setNextFileOffline(self):
475 475
476 476 try:
477 477 self.filename = self.filenameList[self.fileIndex]
478 478 self.amisrFilePointer = h5py.File(self.filename,'r')
479 479 self.fileIndex += 1
480 480 except:
481 481 self.flagNoMoreFiles = 1
482 482 raise schainpy.admin.SchainError('No more files to read')
483 483 return 0
484 484
485 485 self.flagIsNewFile = 1
486 486 print("Setting the file: %s"%self.filename)
487 487
488 488 return 1
489 489
490 490
491 491 def __setNextFileOnline(self):
492 492 filename = self.filenameList[0]
493 493 if self.__filename_online != None:
494 494 self.__selectDataForTimes(online=True)
495 495 filename = self.filenameList[0]
496 496 wait = 0
497 497 self.__waitForNewFile=300 ## DEBUG:
498 498 while self.__filename_online == filename:
499 499 print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
500 500 if wait == 5:
501 501 self.flagNoMoreFiles = 1
502 502 return 0
503 503 sleep(self.__waitForNewFile)
504 504 self.__selectDataForTimes(online=True)
505 505 filename = self.filenameList[0]
506 506 wait += 1
507 507
508 508 self.__filename_online = filename
509 509
510 510 self.amisrFilePointer = h5py.File(filename,'r')
511 511 self.flagIsNewFile = 1
512 512 self.filename = filename
513 513 print("Setting the file: %s"%self.filename)
514 514 return 1
515 515
516 516
517 517 def readData(self):
518 518 buffer = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
519 519 re = buffer[:,:,:,0]
520 520 im = buffer[:,:,:,1]
521 521 dataset = re + im*1j
522 522
523 523 self.radacTime = self.amisrFilePointer.get('Raw11/Data/RadacHeader/RadacTime')
524 524 timeset = self.radacTime[:,0]
525 525
526 526 return dataset,timeset
527 527
528 528 def reshapeData(self):
529 529 #print(self.beamCodeByPulse, self.beamCode, self.nblocks, self.nprofiles, self.nsa)
530 530 channels = self.beamCodeByPulse[0,:]
531 531 nchan = self.nchannels
532 532 #self.newProfiles = self.nprofiles/nchan #must be defined on filljroheader
533 533 nblocks = self.nblocks
534 534 nsamples = self.nsa
535 535 #print("Channels: ",self.nChannels)
536 536 #Dimensions : nChannels, nProfiles, nSamples
537 537 new_block = numpy.empty((nblocks, nchan, numpy.int_(self.newProfiles), nsamples), dtype="complex64")
538 538 ############################################
539 539 profPerCH = int(self.profPerBlockRAW / (self.nFFT* self.nChannels))
540 540 #profPerCH = int(self.profPerBlockRAW / self.nChannels)
541 541 for thisChannel in range(nchan):
542 542
543 543 ich = thisChannel
544 544
545 545 idx_ch = [self.nFFT*(ich + nchan*k) for k in range(profPerCH)]
546 546 #print(idx_ch)
547 547 if self.nFFT > 1:
548 548 aux = [numpy.arange(i, i+self.nFFT) for i in idx_ch]
549 549 idx_ch = None
550 550 idx_ch =aux
551 551 idx_ch = numpy.array(idx_ch, dtype=int).flatten()
552 552 else:
553 553 idx_ch = numpy.array(idx_ch, dtype=int)
554 554
555 555 #print(ich,profPerCH,idx_ch)
556 556 #print(numpy.where(channels==self.beamCode[ich])[0])
557 557 #new_block[:,ich,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[ich])[0],:]
558 558 new_block[:,ich,:,:] = self.dataset[:,idx_ch,:]
559 559
560 560 new_block = numpy.transpose(new_block, (1,0,2,3))
561 561 new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
562 562 if self.flagAsync:
563 563 new_block = numpy.roll(new_block, self.shiftChannels, axis=0)
564 564 return new_block
565 565
566 566 def updateIndexes(self):
567 567
568 568 pass
569 569
570 570 def fillJROHeader(self):
571 571
572 572 #fill radar controller header
573 573
574 574 #fill system header
575 575 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
576 576 nProfiles=self.newProfiles,
577 577 nChannels=len(self.__channelList),
578 578 adcResolution=14,
579 579 pciDioBusWidth=32)
580 580
581 581 self.dataOut.type = "Voltage"
582 582 self.dataOut.data = None
583 583 self.dataOut.dtype = numpy.dtype([('real','<i8'),('imag','<i8')])
584 584 # self.dataOut.nChannels = 0
585 585
586 586 # self.dataOut.nHeights = 0
587 587
588 588 self.dataOut.nProfiles = self.newProfiles*self.nblocks
589 589 #self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
590 590 ranges = numpy.reshape(self.rangeFromFile[()],(-1))
591 591 self.dataOut.heightList = ranges/1000.0 #km
592 592 self.dataOut.channelList = self.__channelList
593 593
594 594 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
595 595
596 596 # self.dataOut.channelIndexList = None
597 597
598 598
599 599 # #self.dataOut.azimuthList = numpy.roll( numpy.array(self.azimuthList) ,self.shiftChannels)
600 600 # #self.dataOut.elevationList = numpy.roll(numpy.array(self.elevationList) ,self.shiftChannels)
601 601 # #self.dataOut.codeList = numpy.roll(numpy.array(self.beamCode), self.shiftChannels)
602 602
603 603 self.dataOut.azimuthList = self.azimuthList
604 604 self.dataOut.elevationList = self.elevationList
605 605 self.dataOut.codeList = self.beamCode
606 606
607 607
608 608
609 609 #print(self.dataOut.elevationList)
610 610 self.dataOut.flagNoData = True
611 611
612 612 #Set to TRUE if the data is discontinuous
613 613 self.dataOut.flagDiscontinuousBlock = False
614 614
615 615 self.dataOut.utctime = None
616 616
617 617 #self.dataOut.timeZone = -5 #self.__timezone/60 #timezone like jroheader, difference in minutes between UTC and localtime
618 618 if self.timezone == 'lt':
619 619 self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
620 620 else:
621 621 self.dataOut.timeZone = 0 #by default time is UTC
622
622
623 623 self.dataOut.dstFlag = 0
624 624 self.dataOut.errorCount = 0
625 625 self.dataOut.nCohInt = 1
626 626 self.dataOut.flagDecodeData = False #asumo que la data esta decodificada
627 627 self.dataOut.flagDeflipData = False #asumo que la data esta sin flip
628 628 self.dataOut.flagShiftFFT = False
629 629 self.dataOut.ippSeconds = self.ippSeconds
630 630 self.dataOut.ipp = self.__ippKm
631 631 self.dataOut.nCode = self.__nCode
632 632 self.dataOut.code = self.__code
633 633 self.dataOut.nBaud = self.__nBaud
634 634
635 635
636 636 self.dataOut.frequency = self.__frequency
637 637 self.dataOut.realtime = self.online
638 638
639 639 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(ipp=self.__ippKm,
640 640 txA=self.__txAKm,
641 641 txB=0,
642 642 nWindows=1,
643 643 nHeights=self.__nSamples,
644 644 firstHeight=self.__firstHeight,
645 645 codeType=self.__codeType,
646 646 nCode=self.__nCode, nBaud=self.__nBaud,
647 647 code = self.__code,
648 648 nOsamp=self.nOsamp,
649 649 frequency = self.__frequency,
650 650 sampleRate= self.__sampleRate,
651 651 fClock=self.__sampleRate)
652 652
653 653
654 654 self.dataOut.radarControllerHeaderObj.heightList = ranges/1000.0 #km
655 655 self.dataOut.radarControllerHeaderObj.heightResolution = self.__deltaHeight
656 656 self.dataOut.radarControllerHeaderObj.rangeIpp = self.__ippKm #km
657 657 self.dataOut.radarControllerHeaderObj.rangeTxA = self.__txA*1e6*.15 #km
658 658 self.dataOut.radarControllerHeaderObj.nChannels = self.nchannels
659 659 self.dataOut.radarControllerHeaderObj.channelList = self.__channelList
660 660 self.dataOut.radarControllerHeaderObj.azimuthList = self.azimuthList
661 661 self.dataOut.radarControllerHeaderObj.elevationList = self.elevationList
662 662 self.dataOut.radarControllerHeaderObj.dtype = "Voltage"
663 663 self.dataOut.ippSeconds = self.ippSeconds
664 664 self.dataOut.ippFactor = self.nFFT
665 665 pass
666 666
667 667 def readNextFile(self,online=False):
668 668
669 669 if not(online):
670 670 newFile = self.__setNextFileOffline()
671 671 else:
672 672 newFile = self.__setNextFileOnline()
673 673
674 674 if not(newFile):
675 675 self.dataOut.error = True
676 676 return 0
677 677
678 678 if not self.readAMISRHeader(self.amisrFilePointer):
679 679 self.dataOut.error = True
680 680 return 0
681 681
682 682 #self.createBuffers()
683 683 self.fillJROHeader()
684 684
685 685 #self.__firstFile = False
686 686
687 687 self.dataset,self.timeset = self.readData()
688 688
689 689 if self.endDate!=None:
690 690 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
691 691 time_str = self.amisrFilePointer.get('Time/RadacTimeString')
692 692 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
693 693 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
694 694 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
695 695 if self.timezone == 'lt':
696 696 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
697 697 if (startDateTime_File>endDateTime_Reader):
698 698 self.flag_standby = False
699 699 return 0
700 700 if self.flag_ignoreFiles and (startDateTime_File >= self.ignStartDateTime and startDateTime_File <= self.ignEndDateTime):
701 701 print("Ignoring...")
702 702 self.flag_standby = True
703 703 return 1
704 704 self.flag_standby = False
705 705
706 706 self.jrodataset = self.reshapeData()
707 707 #----self.updateIndexes()
708 708 self.profileIndex = 0
709 709
710 710 return 1
711 711
712 712
713 713 def __hasNotDataInBuffer(self):
714 714 if self.profileIndex >= (self.newProfiles*self.nblocks):
715 715 return 1
716 716 return 0
717 717
718 718
719 719 def getData(self):
720 720
721 721 if self.flagNoMoreFiles:
722 722 self.dataOut.flagNoData = True
723 723 return 0
724 724
725 725 if self.profileIndex >= (self.newProfiles*self.nblocks): #
726 726 #if self.__hasNotDataInBuffer():
727 727 if not (self.readNextFile(self.online)):
728 728 print("Profile Index break...")
729 729 return 0
730 730
731 731 if self.flag_standby: #Standby mode, if files are being ignoring, just return with no error flag
732 732 return 0
733 733
734 734 if self.dataset is None: # setear esta condicion cuando no hayan datos por leer
735 735 self.dataOut.flagNoData = True
736 736 print("No more data break...")
737 737 return 0
738 738
739 739 #self.dataOut.data = numpy.reshape(self.jrodataset[self.profileIndex,:],(1,-1))
740 740
741 741 self.dataOut.data = self.jrodataset[:,self.profileIndex,:]
742 742
743 743 #print("R_t",self.timeset)
744 744
745 745 #self.dataOut.utctime = self.jrotimeset[self.profileIndex]
746 746 #verificar basic header de jro data y ver si es compatible con este valor
747 747 #self.dataOut.utctime = self.timeset + (self.profileIndex * self.ippSeconds * self.nchannels)
748 748 indexprof = numpy.mod(self.profileIndex, self.newProfiles)
749 749 indexblock = self.profileIndex/self.newProfiles
750 750 #print (indexblock, indexprof)
751 751 diffUTC = 0
752 752 t_comp = (indexprof * self.ippSeconds * self.nchannels) + diffUTC #
753 753
754 754 #print("utc :",indexblock," __ ",t_comp)
755 755 #print(numpy.shape(self.timeset))
756 756 self.dataOut.utctime = self.timeset[numpy.int_(indexblock)] + t_comp
757 757 #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
758 758
759 759 self.dataOut.profileIndex = self.profileIndex
760 760 #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
761 761 self.dataOut.flagNoData = False
762 762 # if indexprof == 0:
763 763 # print("kamisr: ",self.dataOut.utctime)
764 764
765 765 self.profileIndex += 1
766 766
767 767 return self.dataOut.data #retorno necesario??
768 768
769 769
770 770 def run(self, **kwargs):
771 771 '''
772 772 This method will be called many times so here you should put all your code
773 773 '''
774 774 #print("running kamisr")
775 775 if not self.isConfig:
776 776 self.setup(**kwargs)
777 777 self.isConfig = True
778 778
779 779 self.getData()
General Comments 0
You need to be logged in to leave comments. Login now