##// END OF EJS Templates
schain
RSS
Fix BLTR spectra reader
Juan C. Espinoza -
r1211:2e26717528e3
parent child
Show More
Side by Side Unified
Context file:
r1211:2e26717528e3 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1358 +1,1363
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10 import json
11 11
12 12 from schainpy.utils import log
13 13 from .jroheaderIO import SystemHeader, RadarControllerHeader
14 14
15 15
16 16 def getNumpyDtype(dataTypeCode):
17 17
18 18 if dataTypeCode == 0:
19 19 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
20 20 elif dataTypeCode == 1:
21 21 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
22 22 elif dataTypeCode == 2:
23 23 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
24 24 elif dataTypeCode == 3:
25 25 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
26 26 elif dataTypeCode == 4:
27 27 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
28 28 elif dataTypeCode == 5:
29 29 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
30 30 else:
31 31 raise ValueError('dataTypeCode was not defined')
32 32
33 33 return numpyDtype
34 34
35 35
36 36 def getDataTypeCode(numpyDtype):
37 37
38 38 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
39 39 datatype = 0
40 40 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
41 41 datatype = 1
42 42 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
43 43 datatype = 2
44 44 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
45 45 datatype = 3
46 46 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
47 47 datatype = 4
48 48 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
49 49 datatype = 5
50 50 else:
51 51 datatype = None
52 52
53 53 return datatype
54 54
55 55
56 56 def hildebrand_sekhon(data, navg):
57 57 """
58 58 This method is for the objective determination of the noise level in Doppler spectra. This
59 59 implementation technique is based on the fact that the standard deviation of the spectral
60 60 densities is equal to the mean spectral density for white Gaussian noise
61 61
62 62 Inputs:
63 63 Data : heights
64 64 navg : numbers of averages
65 65
66 66 Return:
67 67 mean : noise's level
68 68 """
69 69
70 70 sortdata = numpy.sort(data, axis=None)
71 71 lenOfData = len(sortdata)
72 72 nums_min = lenOfData*0.2
73 73
74 74 if nums_min <= 5:
75 75
76 76 nums_min = 5
77 77
78 78 sump = 0.
79 79 sumq = 0.
80 80
81 81 j = 0
82 82 cont = 1
83 83
84 84 while((cont == 1)and(j < lenOfData)):
85 85
86 86 sump += sortdata[j]
87 87 sumq += sortdata[j]**2
88 88
89 89 if j > nums_min:
90 90 rtest = float(j)/(j-1) + 1.0/navg
91 91 if ((sumq*j) > (rtest*sump**2)):
92 92 j = j - 1
93 93 sump = sump - sortdata[j]
94 94 sumq = sumq - sortdata[j]**2
95 95 cont = 0
96 96
97 97 j += 1
98 98
99 99 lnoise = sump / j
100 100
101 101 return lnoise
102 102
103 103
104 104 class Beam:
105 105
106 106 def __init__(self):
107 107 self.codeList = []
108 108 self.azimuthList = []
109 109 self.zenithList = []
110 110
111 111
112 112 class GenericData(object):
113 113
114 114 flagNoData = True
115 115
116 116 def copy(self, inputObj=None):
117 117
118 118 if inputObj == None:
119 119 return copy.deepcopy(self)
120 120
121 121 for key in list(inputObj.__dict__.keys()):
122 122
123 123 attribute = inputObj.__dict__[key]
124 124
125 125 # If this attribute is a tuple or list
126 126 if type(inputObj.__dict__[key]) in (tuple, list):
127 127 self.__dict__[key] = attribute[:]
128 128 continue
129 129
130 130 # If this attribute is another object or instance
131 131 if hasattr(attribute, '__dict__'):
132 132 self.__dict__[key] = attribute.copy()
133 133 continue
134 134
135 135 self.__dict__[key] = inputObj.__dict__[key]
136 136
137 137 def deepcopy(self):
138 138
139 139 return copy.deepcopy(self)
140 140
141 141 def isEmpty(self):
142 142
143 143 return self.flagNoData
144 144
145 145
146 146 class JROData(GenericData):
147 147
148 148 # m_BasicHeader = BasicHeader()
149 149 # m_ProcessingHeader = ProcessingHeader()
150 150
151 151 systemHeaderObj = SystemHeader()
152 152 radarControllerHeaderObj = RadarControllerHeader()
153 153 # data = None
154 154 type = None
155 155 datatype = None # dtype but in string
156 156 # dtype = None
157 157 # nChannels = None
158 158 # nHeights = None
159 159 nProfiles = None
160 160 heightList = None
161 161 channelList = None
162 162 flagDiscontinuousBlock = False
163 163 useLocalTime = False
164 164 utctime = None
165 165 timeZone = None
166 166 dstFlag = None
167 167 errorCount = None
168 168 blocksize = None
169 169 # nCode = None
170 170 # nBaud = None
171 171 # code = None
172 172 flagDecodeData = False # asumo q la data no esta decodificada
173 173 flagDeflipData = False # asumo q la data no esta sin flip
174 174 flagShiftFFT = False
175 175 # ippSeconds = None
176 176 # timeInterval = None
177 177 nCohInt = None
178 178 # noise = None
179 179 windowOfFilter = 1
180 180 # Speed of ligth
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
195 195 def __str__(self):
196 196
197 197 return '{} - {}'.format(self.type, self.getDatatime())
198 198
199 199 def getNoise(self):
200 200
201 201 raise NotImplementedError
202 202
203 203 def getNChannels(self):
204 204
205 205 return len(self.channelList)
206 206
207 207 def getChannelIndexList(self):
208 208
209 209 return list(range(self.nChannels))
210 210
211 211 def getNHeights(self):
212 212
213 213 return len(self.heightList)
214 214
215 215 def getHeiRange(self, extrapoints=0):
216 216
217 217 heis = self.heightList
218 218 # deltah = self.heightList[1] - self.heightList[0]
219 219 #
220 220 # heis.append(self.heightList[-1])
221 221
222 222 return heis
223 223
224 224 def getDeltaH(self):
225 225
226 226 delta = self.heightList[1] - self.heightList[0]
227 227
228 228 return delta
229 229
230 230 def getltctime(self):
231 231
232 232 if self.useLocalTime:
233 233 return self.utctime - self.timeZone * 60
234 234
235 235 return self.utctime
236 236
237 237 def getDatatime(self):
238 238
239 239 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
240 240 return datatimeValue
241 241
242 242 def getTimeRange(self):
243 243
244 244 datatime = []
245 245
246 246 datatime.append(self.ltctime)
247 247 datatime.append(self.ltctime + self.timeInterval + 1)
248 248
249 249 datatime = numpy.array(datatime)
250 250
251 251 return datatime
252 252
253 253 def getFmaxTimeResponse(self):
254 254
255 255 period = (10**-6) * self.getDeltaH() / (0.15)
256 256
257 257 PRF = 1. / (period * self.nCohInt)
258 258
259 259 fmax = PRF
260 260
261 261 return fmax
262 262
263 263 def getFmax(self):
264 264 PRF = 1. / (self.ippSeconds * self.nCohInt)
265 265
266 266 fmax = PRF
267 267 return fmax
268 268
269 269 def getVmax(self):
270 270
271 271 _lambda = self.C / self.frequency
272 272
273 273 vmax = self.getFmax() * _lambda / 2
274 274
275 275 return vmax
276 276
277 277 def get_ippSeconds(self):
278 278 '''
279 279 '''
280 280 return self.radarControllerHeaderObj.ippSeconds
281 281
282 282 def set_ippSeconds(self, ippSeconds):
283 283 '''
284 284 '''
285 285
286 286 self.radarControllerHeaderObj.ippSeconds = ippSeconds
287 287
288 288 return
289 289
290 290 def get_dtype(self):
291 291 '''
292 292 '''
293 293 return getNumpyDtype(self.datatype)
294 294
295 295 def set_dtype(self, numpyDtype):
296 296 '''
297 297 '''
298 298
299 299 self.datatype = getDataTypeCode(numpyDtype)
300 300
301 301 def get_code(self):
302 302 '''
303 303 '''
304 304 return self.radarControllerHeaderObj.code
305 305
306 306 def set_code(self, code):
307 307 '''
308 308 '''
309 309 self.radarControllerHeaderObj.code = code
310 310
311 311 return
312 312
313 313 def get_ncode(self):
314 314 '''
315 315 '''
316 316 return self.radarControllerHeaderObj.nCode
317 317
318 318 def set_ncode(self, nCode):
319 319 '''
320 320 '''
321 321 self.radarControllerHeaderObj.nCode = nCode
322 322
323 323 return
324 324
325 325 def get_nbaud(self):
326 326 '''
327 327 '''
328 328 return self.radarControllerHeaderObj.nBaud
329 329
330 330 def set_nbaud(self, nBaud):
331 331 '''
332 332 '''
333 333 self.radarControllerHeaderObj.nBaud = nBaud
334 334
335 335 return
336 336
337 337 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
338 338 channelIndexList = property(
339 339 getChannelIndexList, "I'm the 'channelIndexList' property.")
340 340 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
341 341 #noise = property(getNoise, "I'm the 'nHeights' property.")
342 342 datatime = property(getDatatime, "I'm the 'datatime' property")
343 343 ltctime = property(getltctime, "I'm the 'ltctime' property")
344 344 ippSeconds = property(get_ippSeconds, set_ippSeconds)
345 345 dtype = property(get_dtype, set_dtype)
346 346 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
347 347 code = property(get_code, set_code)
348 348 nCode = property(get_ncode, set_ncode)
349 349 nBaud = property(get_nbaud, set_nbaud)
350 350
351 351
352 352 class Voltage(JROData):
353 353
354 354 # data es un numpy array de 2 dmensiones (canales, alturas)
355 355 data = None
356 356
357 357 def __init__(self):
358 358 '''
359 359 Constructor
360 360 '''
361 361
362 362 self.useLocalTime = True
363 363 self.radarControllerHeaderObj = RadarControllerHeader()
364 364 self.systemHeaderObj = SystemHeader()
365 365 self.type = "Voltage"
366 366 self.data = None
367 367 # self.dtype = None
368 368 # self.nChannels = 0
369 369 # self.nHeights = 0
370 370 self.nProfiles = None
371 371 self.heightList = None
372 372 self.channelList = None
373 373 # self.channelIndexList = None
374 374 self.flagNoData = True
375 375 self.flagDiscontinuousBlock = False
376 376 self.utctime = None
377 377 self.timeZone = None
378 378 self.dstFlag = None
379 379 self.errorCount = None
380 380 self.nCohInt = None
381 381 self.blocksize = None
382 382 self.flagDecodeData = False # asumo q la data no esta decodificada
383 383 self.flagDeflipData = False # asumo q la data no esta sin flip
384 384 self.flagShiftFFT = False
385 385 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
386 386 self.profileIndex = 0
387 387
388 388 def getNoisebyHildebrand(self, channel=None):
389 389 """
390 390 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
391 391
392 392 Return:
393 393 noiselevel
394 394 """
395 395
396 396 if channel != None:
397 397 data = self.data[channel]
398 398 nChannels = 1
399 399 else:
400 400 data = self.data
401 401 nChannels = self.nChannels
402 402
403 403 noise = numpy.zeros(nChannels)
404 404 power = data * numpy.conjugate(data)
405 405
406 406 for thisChannel in range(nChannels):
407 407 if nChannels == 1:
408 408 daux = power[:].real
409 409 else:
410 410 daux = power[thisChannel, :].real
411 411 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
412 412
413 413 return noise
414 414
415 415 def getNoise(self, type=1, channel=None):
416 416
417 417 if type == 1:
418 418 noise = self.getNoisebyHildebrand(channel)
419 419
420 420 return noise
421 421
422 422 def getPower(self, channel=None):
423 423
424 424 if channel != None:
425 425 data = self.data[channel]
426 426 else:
427 427 data = self.data
428 428
429 429 power = data * numpy.conjugate(data)
430 430 powerdB = 10 * numpy.log10(power.real)
431 431 powerdB = numpy.squeeze(powerdB)
432 432
433 433 return powerdB
434 434
435 435 def getTimeInterval(self):
436 436
437 437 timeInterval = self.ippSeconds * self.nCohInt
438 438
439 439 return timeInterval
440 440
441 441 noise = property(getNoise, "I'm the 'nHeights' property.")
442 442 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
443 443
444 444
445 445 class Spectra(JROData):
446 446
447 447 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
448 448 data_spc = None
449 449 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
450 450 data_cspc = None
451 451 # data dc es un numpy array de 2 dmensiones (canales, alturas)
452 452 data_dc = None
453 453 # data power
454 454 data_pwr = None
455 455 nFFTPoints = None
456 456 # nPairs = None
457 457 pairsList = None
458 458 nIncohInt = None
459 459 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
460 460 nCohInt = None # se requiere para determinar el valor de timeInterval
461 461 ippFactor = None
462 462 profileIndex = 0
463 463 plotting = "spectra"
464 464
465 465 def __init__(self):
466 466 '''
467 467 Constructor
468 468 '''
469 469
470 470 self.useLocalTime = True
471 471 self.radarControllerHeaderObj = RadarControllerHeader()
472 472 self.systemHeaderObj = SystemHeader()
473 473 self.type = "Spectra"
474 474 # self.data = None
475 475 # self.dtype = None
476 476 # self.nChannels = 0
477 477 # self.nHeights = 0
478 478 self.nProfiles = None
479 479 self.heightList = None
480 480 self.channelList = None
481 481 # self.channelIndexList = None
482 482 self.pairsList = None
483 483 self.flagNoData = True
484 484 self.flagDiscontinuousBlock = False
485 485 self.utctime = None
486 486 self.nCohInt = None
487 487 self.nIncohInt = None
488 488 self.blocksize = None
489 489 self.nFFTPoints = None
490 490 self.wavelength = None
491 491 self.flagDecodeData = False # asumo q la data no esta decodificada
492 492 self.flagDeflipData = False # asumo q la data no esta sin flip
493 493 self.flagShiftFFT = False
494 494 self.ippFactor = 1
495 495 #self.noise = None
496 496 self.beacon_heiIndexList = []
497 497 self.noise_estimation = None
498 498
499 499 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
500 500 """
501 501 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
502 502
503 503 Return:
504 504 noiselevel
505 505 """
506 506
507 507 noise = numpy.zeros(self.nChannels)
508 508
509 509 for channel in range(self.nChannels):
510 510 daux = self.data_spc[channel,
511 511 xmin_index:xmax_index, ymin_index:ymax_index]
512 512 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
513 513
514 514 return noise
515 515
516 516 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
517 517
518 518 if self.noise_estimation is not None:
519 519 # this was estimated by getNoise Operation defined in jroproc_spectra.py
520 520 return self.noise_estimation
521 521 else:
522 522 noise = self.getNoisebyHildebrand(
523 523 xmin_index, xmax_index, ymin_index, ymax_index)
524 524 return noise
525 525
526 526 def getFreqRangeTimeResponse(self, extrapoints=0):
527 527
528 528 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
529 529 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
530 530
531 531 return freqrange
532 532
533 533 def getAcfRange(self, extrapoints=0):
534 534
535 535 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
536 536 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
537 537
538 538 return freqrange
539 539
540 540 def getFreqRange(self, extrapoints=0):
541 541
542 542 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
543 543 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
544 544
545 545 return freqrange
546 546
547 547 def getVelRange(self, extrapoints=0):
548 548
549 549 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
550 550 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
551 551
552 552 if self.nmodes:
553 553 return velrange/self.nmodes
554 554 else:
555 555 return velrange
556 556
557 557 def getNPairs(self):
558 558
559 559 return len(self.pairsList)
560 560
561 561 def getPairsIndexList(self):
562 562
563 563 return list(range(self.nPairs))
564 564
565 565 def getNormFactor(self):
566 566
567 567 pwcode = 1
568 568
569 569 if self.flagDecodeData:
570 570 pwcode = numpy.sum(self.code[0]**2)
571 571 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
572 572 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
573 573
574 574 return normFactor
575 575
576 576 def getFlagCspc(self):
577 577
578 578 if self.data_cspc is None:
579 579 return True
580 580
581 581 return False
582 582
583 583 def getFlagDc(self):
584 584
585 585 if self.data_dc is None:
586 586 return True
587 587
588 588 return False
589 589
590 590 def getTimeInterval(self):
591 591
592 592 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
593
593 if self.nmodes:
594 return self.nmodes*timeInterval
595 else:
594 596 return timeInterval
595 597
596 598 def getPower(self):
597 599
598 600 factor = self.normFactor
599 601 z = self.data_spc / factor
600 602 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
601 603 avg = numpy.average(z, axis=1)
602 604
603 605 return 10 * numpy.log10(avg)
604 606
605 607 def getCoherence(self, pairsList=None, phase=False):
606 608
607 609 z = []
608 610 if pairsList is None:
609 611 pairsIndexList = self.pairsIndexList
610 612 else:
611 613 pairsIndexList = []
612 614 for pair in pairsList:
613 615 if pair not in self.pairsList:
614 616 raise ValueError("Pair %s is not in dataOut.pairsList" % (
615 617 pair))
616 618 pairsIndexList.append(self.pairsList.index(pair))
617 619 for i in range(len(pairsIndexList)):
618 620 pair = self.pairsList[pairsIndexList[i]]
619 621 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
620 622 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
621 623 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
622 624 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
623 625 if phase:
624 626 data = numpy.arctan2(avgcoherenceComplex.imag,
625 627 avgcoherenceComplex.real) * 180 / numpy.pi
626 628 else:
627 629 data = numpy.abs(avgcoherenceComplex)
628 630
629 631 z.append(data)
630 632
631 633 return numpy.array(z)
632 634
633 635 def setValue(self, value):
634 636
635 637 print("This property should not be initialized")
636 638
637 639 return
638 640
639 641 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
640 642 pairsIndexList = property(
641 643 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
642 644 normFactor = property(getNormFactor, setValue,
643 645 "I'm the 'getNormFactor' property.")
644 646 flag_cspc = property(getFlagCspc, setValue)
645 647 flag_dc = property(getFlagDc, setValue)
646 648 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
647 649 timeInterval = property(getTimeInterval, setValue,
648 650 "I'm the 'timeInterval' property")
649 651
650 652
651 653 class SpectraHeis(Spectra):
652 654
653 655 data_spc = None
654 656 data_cspc = None
655 657 data_dc = None
656 658 nFFTPoints = None
657 659 # nPairs = None
658 660 pairsList = None
659 661 nCohInt = None
660 662 nIncohInt = None
661 663
662 664 def __init__(self):
663 665
664 666 self.radarControllerHeaderObj = RadarControllerHeader()
665 667
666 668 self.systemHeaderObj = SystemHeader()
667 669
668 670 self.type = "SpectraHeis"
669 671
670 672 # self.dtype = None
671 673
672 674 # self.nChannels = 0
673 675
674 676 # self.nHeights = 0
675 677
676 678 self.nProfiles = None
677 679
678 680 self.heightList = None
679 681
680 682 self.channelList = None
681 683
682 684 # self.channelIndexList = None
683 685
684 686 self.flagNoData = True
685 687
686 688 self.flagDiscontinuousBlock = False
687 689
688 690 # self.nPairs = 0
689 691
690 692 self.utctime = None
691 693
692 694 self.blocksize = None
693 695
694 696 self.profileIndex = 0
695 697
696 698 self.nCohInt = 1
697 699
698 700 self.nIncohInt = 1
699 701
700 702 def getNormFactor(self):
701 703 pwcode = 1
702 704 if self.flagDecodeData:
703 705 pwcode = numpy.sum(self.code[0]**2)
704 706
705 707 normFactor = self.nIncohInt * self.nCohInt * pwcode
706 708
707 709 return normFactor
708 710
709 711 def getTimeInterval(self):
710 712
711 713 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
712 714
713 715 return timeInterval
714 716
715 717 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
716 718 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
717 719
718 720
719 721 class Fits(JROData):
720 722
721 723 heightList = None
722 724 channelList = None
723 725 flagNoData = True
724 726 flagDiscontinuousBlock = False
725 727 useLocalTime = False
726 728 utctime = None
727 729 timeZone = None
728 730 # ippSeconds = None
729 731 # timeInterval = None
730 732 nCohInt = None
731 733 nIncohInt = None
732 734 noise = None
733 735 windowOfFilter = 1
734 736 # Speed of ligth
735 737 C = 3e8
736 738 frequency = 49.92e6
737 739 realtime = False
738 740
739 741 def __init__(self):
740 742
741 743 self.type = "Fits"
742 744
743 745 self.nProfiles = None
744 746
745 747 self.heightList = None
746 748
747 749 self.channelList = None
748 750
749 751 # self.channelIndexList = None
750 752
751 753 self.flagNoData = True
752 754
753 755 self.utctime = None
754 756
755 757 self.nCohInt = 1
756 758
757 759 self.nIncohInt = 1
758 760
759 761 self.useLocalTime = True
760 762
761 763 self.profileIndex = 0
762 764
763 765 # self.utctime = None
764 766 # self.timeZone = None
765 767 # self.ltctime = None
766 768 # self.timeInterval = None
767 769 # self.header = None
768 770 # self.data_header = None
769 771 # self.data = None
770 772 # self.datatime = None
771 773 # self.flagNoData = False
772 774 # self.expName = ''
773 775 # self.nChannels = None
774 776 # self.nSamples = None
775 777 # self.dataBlocksPerFile = None
776 778 # self.comments = ''
777 779 #
778 780
779 781 def getltctime(self):
780 782
781 783 if self.useLocalTime:
782 784 return self.utctime - self.timeZone * 60
783 785
784 786 return self.utctime
785 787
786 788 def getDatatime(self):
787 789
788 790 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
789 791 return datatime
790 792
791 793 def getTimeRange(self):
792 794
793 795 datatime = []
794 796
795 797 datatime.append(self.ltctime)
796 798 datatime.append(self.ltctime + self.timeInterval)
797 799
798 800 datatime = numpy.array(datatime)
799 801
800 802 return datatime
801 803
802 804 def getHeiRange(self):
803 805
804 806 heis = self.heightList
805 807
806 808 return heis
807 809
808 810 def getNHeights(self):
809 811
810 812 return len(self.heightList)
811 813
812 814 def getNChannels(self):
813 815
814 816 return len(self.channelList)
815 817
816 818 def getChannelIndexList(self):
817 819
818 820 return list(range(self.nChannels))
819 821
820 822 def getNoise(self, type=1):
821 823
822 824 #noise = numpy.zeros(self.nChannels)
823 825
824 826 if type == 1:
825 827 noise = self.getNoisebyHildebrand()
826 828
827 829 if type == 2:
828 830 noise = self.getNoisebySort()
829 831
830 832 if type == 3:
831 833 noise = self.getNoisebyWindow()
832 834
833 835 return noise
834 836
835 837 def getTimeInterval(self):
836 838
837 839 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
838 840
839 841 return timeInterval
840 842
841 843 def get_ippSeconds(self):
842 844 '''
843 845 '''
844 846 return self.ipp_sec
845 847
846 848
847 849 datatime = property(getDatatime, "I'm the 'datatime' property")
848 850 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
849 851 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
850 852 channelIndexList = property(
851 853 getChannelIndexList, "I'm the 'channelIndexList' property.")
852 854 noise = property(getNoise, "I'm the 'nHeights' property.")
853 855
854 856 ltctime = property(getltctime, "I'm the 'ltctime' property")
855 857 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
856 858 ippSeconds = property(get_ippSeconds, '')
857 859
858 860 class Correlation(JROData):
859 861
860 862 noise = None
861 863 SNR = None
862 864 #--------------------------------------------------
863 865 mode = None
864 866 split = False
865 867 data_cf = None
866 868 lags = None
867 869 lagRange = None
868 870 pairsList = None
869 871 normFactor = None
870 872 #--------------------------------------------------
871 873 # calculateVelocity = None
872 874 nLags = None
873 875 nPairs = None
874 876 nAvg = None
875 877
876 878 def __init__(self):
877 879 '''
878 880 Constructor
879 881 '''
880 882 self.radarControllerHeaderObj = RadarControllerHeader()
881 883
882 884 self.systemHeaderObj = SystemHeader()
883 885
884 886 self.type = "Correlation"
885 887
886 888 self.data = None
887 889
888 890 self.dtype = None
889 891
890 892 self.nProfiles = None
891 893
892 894 self.heightList = None
893 895
894 896 self.channelList = None
895 897
896 898 self.flagNoData = True
897 899
898 900 self.flagDiscontinuousBlock = False
899 901
900 902 self.utctime = None
901 903
902 904 self.timeZone = None
903 905
904 906 self.dstFlag = None
905 907
906 908 self.errorCount = None
907 909
908 910 self.blocksize = None
909 911
910 912 self.flagDecodeData = False # asumo q la data no esta decodificada
911 913
912 914 self.flagDeflipData = False # asumo q la data no esta sin flip
913 915
914 916 self.pairsList = None
915 917
916 918 self.nPoints = None
917 919
918 920 def getPairsList(self):
919 921
920 922 return self.pairsList
921 923
922 924 def getNoise(self, mode=2):
923 925
924 926 indR = numpy.where(self.lagR == 0)[0][0]
925 927 indT = numpy.where(self.lagT == 0)[0][0]
926 928
927 929 jspectra0 = self.data_corr[:, :, indR, :]
928 930 jspectra = copy.copy(jspectra0)
929 931
930 932 num_chan = jspectra.shape[0]
931 933 num_hei = jspectra.shape[2]
932 934
933 935 freq_dc = jspectra.shape[1] / 2
934 936 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
935 937
936 938 if ind_vel[0] < 0:
937 939 ind_vel[list(range(0, 1))] = ind_vel[list(
938 940 range(0, 1))] + self.num_prof
939 941
940 942 if mode == 1:
941 943 jspectra[:, freq_dc, :] = (
942 944 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
943 945
944 946 if mode == 2:
945 947
946 948 vel = numpy.array([-2, -1, 1, 2])
947 949 xx = numpy.zeros([4, 4])
948 950
949 951 for fil in range(4):
950 952 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
951 953
952 954 xx_inv = numpy.linalg.inv(xx)
953 955 xx_aux = xx_inv[0, :]
954 956
955 957 for ich in range(num_chan):
956 958 yy = jspectra[ich, ind_vel, :]
957 959 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
958 960
959 961 junkid = jspectra[ich, freq_dc, :] <= 0
960 962 cjunkid = sum(junkid)
961 963
962 964 if cjunkid.any():
963 965 jspectra[ich, freq_dc, junkid.nonzero()] = (
964 966 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
965 967
966 968 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
967 969
968 970 return noise
969 971
970 972 def getTimeInterval(self):
971 973
972 974 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
973 975
974 976 return timeInterval
975 977
976 978 def splitFunctions(self):
977 979
978 980 pairsList = self.pairsList
979 981 ccf_pairs = []
980 982 acf_pairs = []
981 983 ccf_ind = []
982 984 acf_ind = []
983 985 for l in range(len(pairsList)):
984 986 chan0 = pairsList[l][0]
985 987 chan1 = pairsList[l][1]
986 988
987 989 # Obteniendo pares de Autocorrelacion
988 990 if chan0 == chan1:
989 991 acf_pairs.append(chan0)
990 992 acf_ind.append(l)
991 993 else:
992 994 ccf_pairs.append(pairsList[l])
993 995 ccf_ind.append(l)
994 996
995 997 data_acf = self.data_cf[acf_ind]
996 998 data_ccf = self.data_cf[ccf_ind]
997 999
998 1000 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
999 1001
1000 1002 def getNormFactor(self):
1001 1003 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1002 1004 acf_pairs = numpy.array(acf_pairs)
1003 1005 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1004 1006
1005 1007 for p in range(self.nPairs):
1006 1008 pair = self.pairsList[p]
1007 1009
1008 1010 ch0 = pair[0]
1009 1011 ch1 = pair[1]
1010 1012
1011 1013 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1012 1014 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1013 1015 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1014 1016
1015 1017 return normFactor
1016 1018
1017 1019 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1018 1020 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1019 1021
1020 1022
1021 1023 class Parameters(Spectra):
1022 1024
1023 1025 experimentInfo = None # Information about the experiment
1024 1026 # Information from previous data
1025 1027 inputUnit = None # Type of data to be processed
1026 1028 operation = None # Type of operation to parametrize
1027 1029 # normFactor = None #Normalization Factor
1028 1030 groupList = None # List of Pairs, Groups, etc
1029 1031 # Parameters
1030 1032 data_param = None # Parameters obtained
1031 1033 data_pre = None # Data Pre Parametrization
1032 1034 data_SNR = None # Signal to Noise Ratio
1033 1035 # heightRange = None #Heights
1034 1036 abscissaList = None # Abscissa, can be velocities, lags or time
1035 1037 # noise = None #Noise Potency
1036 1038 utctimeInit = None # Initial UTC time
1037 1039 paramInterval = None # Time interval to calculate Parameters in seconds
1038 1040 useLocalTime = True
1039 1041 # Fitting
1040 1042 data_error = None # Error of the estimation
1041 1043 constants = None
1042 1044 library = None
1043 1045 # Output signal
1044 1046 outputInterval = None # Time interval to calculate output signal in seconds
1045 1047 data_output = None # Out signal
1046 1048 nAvg = None
1047 1049 noise_estimation = None
1048 1050 GauSPC = None # Fit gaussian SPC
1049 1051
1050 1052 def __init__(self):
1051 1053 '''
1052 1054 Constructor
1053 1055 '''
1054 1056 self.radarControllerHeaderObj = RadarControllerHeader()
1055 1057
1056 1058 self.systemHeaderObj = SystemHeader()
1057 1059
1058 1060 self.type = "Parameters"
1059 1061
1060 1062 def getTimeRange1(self, interval):
1061 1063
1062 1064 datatime = []
1063 1065
1064 1066 if self.useLocalTime:
1065 1067 time1 = self.utctimeInit - self.timeZone * 60
1066 1068 else:
1067 1069 time1 = self.utctimeInit
1068 1070
1069 1071 datatime.append(time1)
1070 1072 datatime.append(time1 + interval)
1071 1073 datatime = numpy.array(datatime)
1072 1074
1073 1075 return datatime
1074 1076
1075 1077 def getTimeInterval(self):
1076 1078
1077 1079 if hasattr(self, 'timeInterval1'):
1078 1080 return self.timeInterval1
1079 1081 else:
1080 1082 return self.paramInterval
1081 1083
1082 1084 def setValue(self, value):
1083 1085
1084 1086 print("This property should not be initialized")
1085 1087
1086 1088 return
1087 1089
1088 1090 def getNoise(self):
1089 1091
1090 1092 return self.spc_noise
1091 1093
1092 1094 timeInterval = property(getTimeInterval)
1093 1095 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1094 1096
1095 1097
1096 1098 class PlotterData(object):
1097 1099 '''
1098 1100 Object to hold data to be plotted
1099 1101 '''
1100 1102
1101 1103 MAXNUMX = 100
1102 1104 MAXNUMY = 100
1103 1105
1104 1106 def __init__(self, code, throttle_value, exp_code, buffering=True):
1105 1107
1106 1108 self.throttle = throttle_value
1107 1109 self.exp_code = exp_code
1108 1110 self.buffering = buffering
1109 1111 self.ready = False
1110 1112 self.localtime = False
1111 1113 self.data = {}
1112 1114 self.meta = {}
1113 1115 self.__times = []
1114 1116 self.__heights = []
1115 1117
1116 1118 if 'snr' in code:
1117 1119 self.plottypes = ['snr']
1118 1120 elif code == 'spc':
1119 1121 self.plottypes = ['spc', 'noise', 'rti']
1120 1122 elif code == 'rti':
1121 1123 self.plottypes = ['noise', 'rti']
1122 1124 else:
1123 1125 self.plottypes = [code]
1124 1126
1125 1127 for plot in self.plottypes:
1126 1128 self.data[plot] = {}
1127 1129
1128 1130 def __str__(self):
1129 1131 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1130 1132 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1131 1133
1132 1134 def __len__(self):
1133 1135 return len(self.__times)
1134 1136
1135 1137 def __getitem__(self, key):
1136 1138
1137 1139 if key not in self.data:
1138 1140 raise KeyError(log.error('Missing key: {}'.format(key)))
1139 1141 if 'spc' in key or not self.buffering:
1140 1142 ret = self.data[key]
1141 1143 elif 'scope' in key:
1142 1144 ret = numpy.array(self.data[key][float(self.tm)])
1143 1145 else:
1144 1146 ret = numpy.array([self.data[key][x] for x in self.times])
1145 1147 if ret.ndim > 1:
1146 1148 ret = numpy.swapaxes(ret, 0, 1)
1147 1149 return ret
1148 1150
1149 1151 def __contains__(self, key):
1150 1152 return key in self.data
1151 1153
1152 1154 def setup(self):
1153 1155 '''
1154 1156 Configure object
1155 1157 '''
1156 1158
1157 1159 self.type = ''
1158 1160 self.ready = False
1159 1161 self.data = {}
1160 1162 self.__times = []
1161 1163 self.__heights = []
1162 1164 self.__all_heights = set()
1163 1165 for plot in self.plottypes:
1164 1166 if 'snr' in plot:
1165 1167 plot = 'snr'
1166 1168 elif 'spc_moments' == plot:
1167 1169 plot = 'moments'
1168 1170 self.data[plot] = {}
1169 1171
1170 1172 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1171 1173 self.data['noise'] = {}
1174 self.data['rti'] = {}
1172 1175 if 'noise' not in self.plottypes:
1173 1176 self.plottypes.append('noise')
1177 if 'rti' not in self.plottypes:
1178 self.plottypes.append('rti')
1174 1179
1175 1180 def shape(self, key):
1176 1181 '''
1177 1182 Get the shape of the one-element data for the given key
1178 1183 '''
1179 1184
1180 1185 if len(self.data[key]):
1181 1186 if 'spc' in key or not self.buffering:
1182 1187 return self.data[key].shape
1183 1188 return self.data[key][self.__times[0]].shape
1184 1189 return (0,)
1185 1190
1186 1191 def update(self, dataOut, tm):
1187 1192 '''
1188 1193 Update data object with new dataOut
1189 1194 '''
1190 1195
1191 1196 if tm in self.__times:
1192 1197 return
1193 1198 self.profileIndex = dataOut.profileIndex
1194 1199 self.tm = tm
1195 1200 self.type = dataOut.type
1196 1201 self.parameters = getattr(dataOut, 'parameters', [])
1197 1202 if hasattr(dataOut, 'pairsList'):
1198 1203 self.pairs = dataOut.pairsList
1199 1204 if hasattr(dataOut, 'meta'):
1200 1205 self.meta = dataOut.meta
1201 1206 self.channels = dataOut.channelList
1202 1207 self.interval = dataOut.getTimeInterval()
1203 1208 self.localtime = dataOut.useLocalTime
1204 1209 if 'spc' in self.plottypes or 'cspc' in self.plottypes or 'spc_moments' in self.plottypes:
1205 1210 self.xrange = (dataOut.getFreqRange(1)/1000.,
1206 1211 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1207 1212 self.factor = dataOut.normFactor
1208 1213 self.__heights.append(dataOut.heightList)
1209 1214 self.__all_heights.update(dataOut.heightList)
1210 1215 self.__times.append(tm)
1211 1216
1212 1217 for plot in self.plottypes:
1213 1218 if plot in ('spc', 'spc_moments'):
1214 1219 z = dataOut.data_spc/dataOut.normFactor
1215 1220 buffer = 10*numpy.log10(z)
1216 1221 if plot == 'cspc':
1217 1222 z = dataOut.data_spc/dataOut.normFactor
1218 1223 buffer = (dataOut.data_spc, dataOut.data_cspc)
1219 1224 if plot == 'noise':
1220 1225 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1221 1226 if plot == 'rti':
1222 1227 buffer = dataOut.getPower()
1223 1228 if plot == 'snr_db':
1224 1229 buffer = dataOut.data_SNR
1225 1230 if plot == 'snr':
1226 1231 buffer = 10*numpy.log10(dataOut.data_SNR)
1227 1232 if plot == 'dop':
1228 1233 buffer = 10*numpy.log10(dataOut.data_DOP)
1229 1234 if plot == 'mean':
1230 1235 buffer = dataOut.data_MEAN
1231 1236 if plot == 'std':
1232 1237 buffer = dataOut.data_STD
1233 1238 if plot == 'coh':
1234 1239 buffer = dataOut.getCoherence()
1235 1240 if plot == 'phase':
1236 1241 buffer = dataOut.getCoherence(phase=True)
1237 1242 if plot == 'output':
1238 1243 buffer = dataOut.data_output
1239 1244 if plot == 'param':
1240 1245 buffer = dataOut.data_param
1241 1246 if plot == 'scope':
1242 1247 buffer = dataOut.data
1243 1248 self.flagDataAsBlock = dataOut.flagDataAsBlock
1244 1249 self.nProfiles = dataOut.nProfiles
1245 1250
1246 1251 if plot == 'spc':
1247 self.data[plot] = buffer
1252 self.data['spc'] = buffer
1248 1253 elif plot == 'cspc':
1249 1254 self.data['spc'] = buffer[0]
1250 1255 self.data['cspc'] = buffer[1]
1251 1256 elif plot == 'spc_moments':
1252 1257 self.data['spc'] = buffer
1253 1258 self.data['moments'][tm] = dataOut.moments
1254 1259 else:
1255 1260 if self.buffering:
1256 1261 self.data[plot][tm] = buffer
1257 1262 else:
1258 1263 self.data[plot] = buffer
1259 1264
1260 1265 def normalize_heights(self):
1261 1266 '''
1262 1267 Ensure same-dimension of the data for different heighList
1263 1268 '''
1264 1269
1265 1270 H = numpy.array(list(self.__all_heights))
1266 1271 H.sort()
1267 1272 for key in self.data:
1268 1273 shape = self.shape(key)[:-1] + H.shape
1269 1274 for tm, obj in list(self.data[key].items()):
1270 1275 h = self.__heights[self.__times.index(tm)]
1271 1276 if H.size == h.size:
1272 1277 continue
1273 1278 index = numpy.where(numpy.in1d(H, h))[0]
1274 1279 dummy = numpy.zeros(shape) + numpy.nan
1275 1280 if len(shape) == 2:
1276 1281 dummy[:, index] = obj
1277 1282 else:
1278 1283 dummy[index] = obj
1279 1284 self.data[key][tm] = dummy
1280 1285
1281 1286 self.__heights = [H for tm in self.__times]
1282 1287
1283 1288 def jsonify(self, decimate=False):
1284 1289 '''
1285 1290 Convert data to json
1286 1291 '''
1287 1292
1288 1293 data = {}
1289 1294 tm = self.times[-1]
1290 1295 dy = int(self.heights.size/self.MAXNUMY) + 1
1291 1296 for key in self.data:
1292 1297 if key in ('spc', 'cspc') or not self.buffering:
1293 1298 dx = int(self.data[key].shape[1]/self.MAXNUMX) + 1
1294 1299 data[key] = self.roundFloats(
1295 1300 self.data[key][::, ::dx, ::dy].tolist())
1296 1301 else:
1297 1302 data[key] = self.roundFloats(self.data[key][tm].tolist())
1298 1303
1299 1304 ret = {'data': data}
1300 1305 ret['exp_code'] = self.exp_code
1301 1306 ret['time'] = float(tm)
1302 1307 ret['interval'] = float(self.interval)
1303 1308 ret['localtime'] = self.localtime
1304 1309 ret['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1305 1310 if 'spc' in self.data or 'cspc' in self.data:
1306 1311 ret['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1307 1312 else:
1308 1313 ret['xrange'] = []
1309 1314 if hasattr(self, 'pairs'):
1310 1315 ret['pairs'] = [(int(p[0]), int(p[1])) for p in self.pairs]
1311 1316 else:
1312 1317 ret['pairs'] = []
1313 1318
1314 1319 for key, value in list(self.meta.items()):
1315 1320 ret[key] = value
1316 1321
1317 1322 return json.dumps(ret)
1318 1323
1319 1324 @property
1320 1325 def times(self):
1321 1326 '''
1322 1327 Return the list of times of the current data
1323 1328 '''
1324 1329
1325 1330 ret = numpy.array(self.__times)
1326 1331 ret.sort()
1327 1332 return ret
1328 1333
1329 1334 @property
1330 1335 def min_time(self):
1331 1336 '''
1332 1337 Return the minimun time value
1333 1338 '''
1334 1339
1335 1340 return self.times[0]
1336 1341
1337 1342 @property
1338 1343 def max_time(self):
1339 1344 '''
1340 1345 Return the maximun time value
1341 1346 '''
1342 1347
1343 1348 return self.times[-1]
1344 1349
1345 1350 @property
1346 1351 def heights(self):
1347 1352 '''
1348 1353 Return the list of heights of the current data
1349 1354 '''
1350 1355
1351 1356 return numpy.array(self.__heights[-1])
1352 1357
1353 1358 @staticmethod
1354 1359 def roundFloats(obj):
1355 1360 if isinstance(obj, list):
1356 1361 return list(map(PlotterData.roundFloats, obj))
1357 1362 elif isinstance(obj, float):
1358 1363 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,466 +1,466
1 1 import os
2 2 import sys
3 3 import glob
4 4 import fnmatch
5 5 import datetime
6 6 import time
7 7 import re
8 8 import h5py
9 9 import numpy
10 10
11 11 import pylab as plb
12 12 from scipy.optimize import curve_fit
13 13 from scipy import asarray as ar, exp
14 14
15 15 SPEED_OF_LIGHT = 299792458
16 16 SPEED_OF_LIGHT = 3e8
17 17
18 18 try:
19 19 from gevent import sleep
20 20 except:
21 21 from time import sleep
22 22
23 23 from .utils import folder_in_range
24 24
25 25 from schainpy.model.data.jrodata import Spectra
26 26 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
27 27 from schainpy.utils import log
28 28 from schainpy.model.io.jroIO_base import JRODataReader
29 29
30 30 def pol2cart(rho, phi):
31 31 x = rho * numpy.cos(phi)
32 32 y = rho * numpy.sin(phi)
33 33 return(x, y)
34 34
35 35 FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes
36 36 ('FileMgcNumber', '<u4'), # 0x23020100
37 37 ('nFDTdataRecors', '<u4'),
38 38 ('OffsetStartHeader', '<u4'),
39 39 ('RadarUnitId', '<u4'),
40 40 ('SiteName', 'S32'), # Null terminated
41 41 ])
42 42
43 43
44 44 class FileHeaderBLTR():
45 45
46 46 def __init__(self, fo):
47 47
48 48 self.fo = fo
49 49 self.size = 48
50 50 self.read()
51 51
52 52 def read(self):
53 53
54 54 header = numpy.fromfile(self.fo, FILE_STRUCTURE, 1)
55 55 self.FileMgcNumber = hex(header['FileMgcNumber'][0])
56 56 self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
57 57 self.RadarUnitId = int(header['RadarUnitId'][0])
58 58 self.OffsetStartHeader = int(header['OffsetStartHeader'][0])
59 59 self.SiteName = header['SiteName'][0]
60 60
61 61 def write(self, fp):
62 62
63 63 headerTuple = (self.FileMgcNumber,
64 64 self.nFDTdataRecors,
65 65 self.RadarUnitId,
66 66 self.SiteName,
67 67 self.size)
68 68
69 69 header = numpy.array(headerTuple, FILE_STRUCTURE)
70 70 header.tofile(fp)
71 71 ''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
72 72
73 73 fid : file or str
74 74 An open file object, or a string containing a filename.
75 75
76 76 sep : str
77 77 Separator between array items for text output. If "" (empty), a binary file is written,
78 78 equivalent to file.write(a.tobytes()).
79 79
80 80 format : str
81 81 Format string for text file output. Each entry in the array is formatted to text by
82 82 first converting it to the closest Python type, and then using "format" % item.
83 83
84 84 '''
85 85
86 86 return 1
87 87
88 88
89 89 RECORD_STRUCTURE = numpy.dtype([ # RECORD HEADER 180+20N bytes
90 90 ('RecMgcNumber', '<u4'), # 0x23030001
91 91 ('RecCounter', '<u4'), # Record counter(0,1, ...)
92 92 # Offset to start of next record form start of this record
93 93 ('Off2StartNxtRec', '<u4'),
94 94 # Offset to start of data from start of this record
95 95 ('Off2StartData', '<u4'),
96 96 # Epoch time stamp of start of acquisition (seconds)
97 97 ('nUtime', '<i4'),
98 98 # Millisecond component of time stamp (0,...,999)
99 99 ('nMilisec', '<u4'),
100 100 # Experiment tag name (null terminated)
101 101 ('ExpTagName', 'S32'),
102 102 # Experiment comment (null terminated)
103 103 ('ExpComment', 'S32'),
104 104 # Site latitude (from GPS) in degrees (positive implies North)
105 105 ('SiteLatDegrees', '<f4'),
106 106 # Site longitude (from GPS) in degrees (positive implies East)
107 107 ('SiteLongDegrees', '<f4'),
108 108 # RTC GPS engine status (0=SEEK, 1=LOCK, 2=NOT FITTED, 3=UNAVAILABLE)
109 109 ('RTCgpsStatus', '<u4'),
110 110 ('TransmitFrec', '<u4'), # Transmit frequency (Hz)
111 111 ('ReceiveFrec', '<u4'), # Receive frequency
112 112 # First local oscillator frequency (Hz)
113 113 ('FirstOsciFrec', '<u4'),
114 114 # (0="O", 1="E", 2="linear 1", 3="linear2")
115 115 ('Polarisation', '<u4'),
116 116 # Receiver filter settings (0,1,2,3)
117 117 ('ReceiverFiltSett', '<u4'),
118 118 # Number of modes in use (1 or 2)
119 119 ('nModesInUse', '<u4'),
120 120 # Dual Mode index number for these data (0 or 1)
121 121 ('DualModeIndex', '<u4'),
122 122 # Dual Mode range correction for these data (m)
123 123 ('DualModeRange', '<u4'),
124 124 # Number of digital channels acquired (2*N)
125 125 ('nDigChannels', '<u4'),
126 126 # Sampling resolution (meters)
127 127 ('SampResolution', '<u4'),
128 128 # Number of range gates sampled
129 129 ('nHeights', '<u4'),
130 130 # Start range of sampling (meters)
131 131 ('StartRangeSamp', '<u4'),
132 132 ('PRFhz', '<u4'), # PRF (Hz)
133 133 ('nCohInt', '<u4'), # Integrations
134 134 # Number of data points transformed
135 135 ('nProfiles', '<u4'),
136 136 # Number of receive beams stored in file (1 or N)
137 137 ('nChannels', '<u4'),
138 138 ('nIncohInt', '<u4'), # Number of spectral averages
139 139 # FFT windowing index (0 = no window)
140 140 ('FFTwindowingInd', '<u4'),
141 141 # Beam steer angle (azimuth) in degrees (clockwise from true North)
142 142 ('BeamAngleAzim', '<f4'),
143 143 # Beam steer angle (zenith) in degrees (0=> vertical)
144 144 ('BeamAngleZen', '<f4'),
145 145 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
146 146 ('AntennaCoord0', '<f4'),
147 147 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
148 148 ('AntennaAngl0', '<f4'),
149 149 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
150 150 ('AntennaCoord1', '<f4'),
151 151 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
152 152 ('AntennaAngl1', '<f4'),
153 153 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
154 154 ('AntennaCoord2', '<f4'),
155 155 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
156 156 ('AntennaAngl2', '<f4'),
157 157 # Receiver phase calibration (degrees) - N values
158 158 ('RecPhaseCalibr0', '<f4'),
159 159 # Receiver phase calibration (degrees) - N values
160 160 ('RecPhaseCalibr1', '<f4'),
161 161 # Receiver phase calibration (degrees) - N values
162 162 ('RecPhaseCalibr2', '<f4'),
163 163 # Receiver amplitude calibration (ratio relative to receiver one) - N values
164 164 ('RecAmpCalibr0', '<f4'),
165 165 # Receiver amplitude calibration (ratio relative to receiver one) - N values
166 166 ('RecAmpCalibr1', '<f4'),
167 167 # Receiver amplitude calibration (ratio relative to receiver one) - N values
168 168 ('RecAmpCalibr2', '<f4'),
169 169 # Receiver gains in dB - N values
170 170 ('ReceiverGaindB0', '<i4'),
171 171 # Receiver gains in dB - N values
172 172 ('ReceiverGaindB1', '<i4'),
173 173 # Receiver gains in dB - N values
174 174 ('ReceiverGaindB2', '<i4'),
175 175 ])
176 176
177 177
178 178 class RecordHeaderBLTR():
179 179
180 180 def __init__(self, fo):
181 181
182 182 self.fo = fo
183 183 self.OffsetStartHeader = 48
184 184 self.Off2StartNxtRec = 811248
185 185
186 186 def read(self, block):
187 187 OffRHeader = self.OffsetStartHeader + block * self.Off2StartNxtRec
188 188 self.fo.seek(OffRHeader, os.SEEK_SET)
189 189 header = numpy.fromfile(self.fo, RECORD_STRUCTURE, 1)
190 190 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) # 0x23030001
191 191 self.RecCounter = int(header['RecCounter'][0])
192 192 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
193 193 self.Off2StartData = int(header['Off2StartData'][0])
194 194 self.nUtime = header['nUtime'][0]
195 195 self.nMilisec = header['nMilisec'][0]
196 196 self.ExpTagName = '' # str(header['ExpTagName'][0])
197 197 self.ExpComment = '' # str(header['ExpComment'][0])
198 198 self.SiteLatDegrees = header['SiteLatDegrees'][0]
199 199 self.SiteLongDegrees = header['SiteLongDegrees'][0]
200 200 self.RTCgpsStatus = header['RTCgpsStatus'][0]
201 201 self.TransmitFrec = header['TransmitFrec'][0]
202 202 self.ReceiveFrec = header['ReceiveFrec'][0]
203 203 self.FirstOsciFrec = header['FirstOsciFrec'][0]
204 204 self.Polarisation = header['Polarisation'][0]
205 205 self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
206 206 self.nModesInUse = header['nModesInUse'][0]
207 207 self.DualModeIndex = header['DualModeIndex'][0]
208 208 self.DualModeRange = header['DualModeRange'][0]
209 209 self.nDigChannels = header['nDigChannels'][0]
210 210 self.SampResolution = header['SampResolution'][0]
211 211 self.nHeights = header['nHeights'][0]
212 212 self.StartRangeSamp = header['StartRangeSamp'][0]
213 213 self.PRFhz = header['PRFhz'][0]
214 214 self.nCohInt = header['nCohInt'][0]
215 215 self.nProfiles = header['nProfiles'][0]
216 216 self.nChannels = header['nChannels'][0]
217 217 self.nIncohInt = header['nIncohInt'][0]
218 218 self.FFTwindowingInd = header['FFTwindowingInd'][0]
219 219 self.BeamAngleAzim = header['BeamAngleAzim'][0]
220 220 self.BeamAngleZen = header['BeamAngleZen'][0]
221 221 self.AntennaCoord0 = header['AntennaCoord0'][0]
222 222 self.AntennaAngl0 = header['AntennaAngl0'][0]
223 223 self.AntennaCoord1 = header['AntennaCoord1'][0]
224 224 self.AntennaAngl1 = header['AntennaAngl1'][0]
225 225 self.AntennaCoord2 = header['AntennaCoord2'][0]
226 226 self.AntennaAngl2 = header['AntennaAngl2'][0]
227 227 self.RecPhaseCalibr0 = header['RecPhaseCalibr0'][0]
228 228 self.RecPhaseCalibr1 = header['RecPhaseCalibr1'][0]
229 229 self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0]
230 230 self.RecAmpCalibr0 = header['RecAmpCalibr0'][0]
231 231 self.RecAmpCalibr1 = header['RecAmpCalibr1'][0]
232 232 self.RecAmpCalibr2 = header['RecAmpCalibr2'][0]
233 233 self.ReceiverGaindB0 = header['ReceiverGaindB0'][0]
234 234 self.ReceiverGaindB1 = header['ReceiverGaindB1'][0]
235 235 self.ReceiverGaindB2 = header['ReceiverGaindB2'][0]
236 236 self.ipp = 0.5 * (SPEED_OF_LIGHT / self.PRFhz)
237 237 self.RHsize = 180 + 20 * self.nChannels
238 238 self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
239 239 endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
240 240
241 241
242 242 if OffRHeader > endFp:
243 243 sys.stderr.write(
244 244 "Warning %s: Size value read from System Header is lower than it has to be\n" % fp)
245 245 return 0
246 246
247 247 if OffRHeader < endFp:
248 248 sys.stderr.write(
249 249 "Warning %s: Size value read from System Header size is greater than it has to be\n" % fp)
250 250 return 0
251 251
252 252 return 1
253 253
254 254 @MPDecorator
255 255 class BLTRSpectraReader (ProcessingUnit):
256 256
257 257 def __init__(self):
258 258
259 259 ProcessingUnit.__init__(self)
260 260
261 261 self.ext = ".fdt"
262 262 self.optchar = "P"
263 263 self.fpFile = None
264 264 self.fp = None
265 265 self.BlockCounter = 0
266 266 self.fileSizeByHeader = None
267 267 self.filenameList = []
268 268 self.fileSelector = 0
269 269 self.Off2StartNxtRec = 0
270 270 self.RecCounter = 0
271 271 self.flagNoMoreFiles = 0
272 272 self.data_spc = None
273 273 self.data_cspc = None
274 274 self.path = None
275 275 self.OffsetStartHeader = 0
276 276 self.Off2StartData = 0
277 277 self.ipp = 0
278 278 self.nFDTdataRecors = 0
279 279 self.blocksize = 0
280 280 self.dataOut = Spectra()
281 281 self.dataOut.flagNoData = False
282 282
283 283 def search_files(self):
284 284 '''
285 285 Function that indicates the number of .fdt files that exist in the folder to be read.
286 286 It also creates an organized list with the names of the files to read.
287 287 '''
288 288
289 289 files = glob.glob(os.path.join(self.path, '*{}'.format(self.ext)))
290 290 files = sorted(files)
291 291 for f in files:
292 292 filename = f.split('/')[-1]
293 293 if folder_in_range(filename.split('.')[1], self.startDate, self.endDate, '%Y%m%d'):
294 294 self.filenameList.append(f)
295 295
296 296 def run(self, **kwargs):
297 297 '''
298 298 This method will be the one that will initiate the data entry, will be called constantly.
299 299 You should first verify that your Setup () is set up and then continue to acquire
300 300 the data to be processed with getData ().
301 301 '''
302 302 if not self.isConfig:
303 303 self.setup(**kwargs)
304 304 self.isConfig = True
305 305
306 306 self.getData()
307 307
308 308 def setup(self,
309 309 path=None,
310 310 startDate=None,
311 311 endDate=None,
312 312 startTime=None,
313 313 endTime=None,
314 314 walk=True,
315 315 code=None,
316 316 online=False,
317 317 mode=None,
318 318 **kwargs):
319 319
320 320 self.isConfig = True
321 321
322 322 self.path = path
323 323 self.startDate = startDate
324 324 self.endDate = endDate
325 325 self.startTime = startTime
326 326 self.endTime = endTime
327 327 self.walk = walk
328 328 self.mode = int(mode)
329 329 self.search_files()
330 330 if self.filenameList:
331 331 self.readFile()
332 332
333 333 def getData(self):
334 334 '''
335 335 Before starting this function, you should check that there is still an unread file,
336 336 If there are still blocks to read or if the data block is empty.
337 337
338 338 You should call the file "read".
339 339
340 340 '''
341 341
342 342 if self.flagNoMoreFiles:
343 343 self.dataOut.flagNoData = True
344 344 self.dataOut.error = 'No more files'
345 345
346 346 self.readBlock()
347 347
348 348 def readFile(self):
349 349 '''
350 350 You must indicate if you are reading in Online or Offline mode and load the
351 351 The parameters for this file reading mode.
352 352
353 353 Then you must do 2 actions:
354 354
355 355 1. Get the BLTR FileHeader.
356 356 2. Start reading the first block.
357 357 '''
358 358
359 359 if self.fileSelector < len(self.filenameList):
360 360 log.success('Opening file: {}'.format(self.filenameList[self.fileSelector]), self.name)
361 361 self.fp = open(self.filenameList[self.fileSelector])
362 362 self.fheader = FileHeaderBLTR(self.fp)
363 363 self.rheader = RecordHeaderBLTR(self.fp)
364 364 self.nFDTdataRecors = self.fheader.nFDTdataRecors
365 365 self.fileSelector += 1
366 366 self.BlockCounter = 0
367 367 return 1
368 368 else:
369 369 self.flagNoMoreFiles = True
370 370 self.dataOut.flagNoData = True
371 371 return 0
372 372
373 373 def readBlock(self):
374 374 '''
375 375 It should be checked if the block has data, if it is not passed to the next file.
376 376
377 377 Then the following is done:
378 378
379 379 1. Read the RecordHeader
380 380 2. Fill the buffer with the current block number.
381 381
382 382 '''
383 383
384 384 if self.BlockCounter == self.nFDTdataRecors:
385 385 if not self.readFile():
386 386 return
387 387
388 388 if self.mode == 1:
389 389 self.rheader.read(self.BlockCounter+1)
390 390 elif self.mode == 0:
391 391 self.rheader.read(self.BlockCounter)
392 392
393 393 self.RecCounter = self.rheader.RecCounter
394 394 self.OffsetStartHeader = self.rheader.OffsetStartHeader
395 395 self.Off2StartNxtRec = self.rheader.Off2StartNxtRec
396 396 self.Off2StartData = self.rheader.Off2StartData
397 397 self.nProfiles = self.rheader.nProfiles
398 398 self.nChannels = self.rheader.nChannels
399 399 self.nHeights = self.rheader.nHeights
400 400 self.frequency = self.rheader.TransmitFrec
401 401 self.DualModeIndex = self.rheader.DualModeIndex
402 402 self.pairsList = [(0, 1), (0, 2), (1, 2)]
403 403 self.dataOut.pairsList = self.pairsList
404 404 self.nRdPairs = len(self.dataOut.pairsList)
405 405 self.dataOut.nRdPairs = self.nRdPairs
406 406 self.dataOut.heightList = (self.rheader.StartRangeSamp + numpy.arange(self.nHeights) * self.rheader.SampResolution) / 1000.
407 407 self.dataOut.channelList = range(self.nChannels)
408 408 self.dataOut.nProfiles=self.rheader.nProfiles
409 409 self.dataOut.nIncohInt=self.rheader.nIncohInt
410 410 self.dataOut.nCohInt=self.rheader.nCohInt
411 411 self.dataOut.ippSeconds= 1/float(self.rheader.PRFhz)
412 412 self.dataOut.PRF=self.rheader.PRFhz
413 413 self.dataOut.nFFTPoints=self.rheader.nProfiles
414 self.dataOut.utctime=self.rheader.nUtime
414 self.dataOut.utctime = self.rheader.nUtime + self.rheader.nMilisec/1000.
415 415 self.dataOut.timeZone = 0
416 416 self.dataOut.useLocalTime = False
417 417 self.dataOut.nmodes = 2
418 418 log.log('Reading block {} - {}'.format(self.BlockCounter, self.dataOut.datatime), self.name)
419 419 OffDATA = self.OffsetStartHeader + self.RecCounter * \
420 420 self.Off2StartNxtRec + self.Off2StartData
421 421 self.fp.seek(OffDATA, os.SEEK_SET)
422 422
423 423 self.data_fft = numpy.fromfile(self.fp, [('complex','<c8')], self.nProfiles*self.nChannels*self.nHeights )
424 424 self.data_fft = self.data_fft.astype(numpy.dtype('complex'))
425 425 self.data_block = numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles))
426 426 self.data_block = numpy.transpose(self.data_block, (1,2,0))
427 427 copy = self.data_block.copy()
428 428 spc = copy * numpy.conjugate(copy)
429 429 self.data_spc = numpy.absolute(spc) # valor absoluto o magnitud
430 430 self.dataOut.data_spc = self.data_spc
431 431
432 432 cspc = self.data_block.copy()
433 433 self.data_cspc = self.data_block.copy()
434 434
435 435 for i in range(self.nRdPairs):
436 436
437 437 chan_index0 = self.dataOut.pairsList[i][0]
438 438 chan_index1 = self.dataOut.pairsList[i][1]
439 439
440 440 self.data_cspc[i, :, :] = cspc[chan_index0, :, :] * numpy.conjugate(cspc[chan_index1, :, :])
441 441
442 442 '''Getting Eij and Nij'''
443 443 (AntennaX0, AntennaY0) = pol2cart(
444 444 self.rheader.AntennaCoord0, self.rheader.AntennaAngl0 * numpy.pi / 180)
445 445 (AntennaX1, AntennaY1) = pol2cart(
446 446 self.rheader.AntennaCoord1, self.rheader.AntennaAngl1 * numpy.pi / 180)
447 447 (AntennaX2, AntennaY2) = pol2cart(
448 448 self.rheader.AntennaCoord2, self.rheader.AntennaAngl2 * numpy.pi / 180)
449 449
450 450 E01 = AntennaX0 - AntennaX1
451 451 N01 = AntennaY0 - AntennaY1
452 452
453 453 E02 = AntennaX0 - AntennaX2
454 454 N02 = AntennaY0 - AntennaY2
455 455
456 456 E12 = AntennaX1 - AntennaX2
457 457 N12 = AntennaY1 - AntennaY2
458 458
459 459 self.ChanDist = numpy.array(
460 460 [[E01, N01], [E02, N02], [E12, N12]])
461 461
462 462 self.dataOut.ChanDist = self.ChanDist
463 463
464 464 self.BlockCounter += 2
465 465 self.dataOut.data_spc = self.data_spc
466 466 self.dataOut.data_cspc =self.data_cspc
General Comments 0
You need to be logged in to leave comments. Login now