##// END OF EJS Templates
schain
RSS
BLTR reader now starts in last block
jespinoza -
r1414:f3c20d0284ef
parent child
Show More
Side by Side Unified
Context file:
r1414:f3c20d0284ef Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1169 +1,1170
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Definition of diferent Data objects for different types of data
6 6
7 7 Here you will find the diferent data objects for the different types
8 8 of data, this data objects must be used as dataIn or dataOut objects in
9 9 processing units and operations. Currently the supported data objects are:
10 10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 11 """
12 12
13 13 import copy
14 14 import numpy
15 15 import datetime
16 16 import json
17 17
18 18 import schainpy.admin
19 19 from schainpy.utils import log
20 20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 21 from schainpy.model.data import _noise
22 22
23 23
24 24 def getNumpyDtype(dataTypeCode):
25 25
26 26 if dataTypeCode == 0:
27 27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 28 elif dataTypeCode == 1:
29 29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 30 elif dataTypeCode == 2:
31 31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 32 elif dataTypeCode == 3:
33 33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 34 elif dataTypeCode == 4:
35 35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 36 elif dataTypeCode == 5:
37 37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 38 else:
39 39 raise ValueError('dataTypeCode was not defined')
40 40
41 41 return numpyDtype
42 42
43 43
44 44 def getDataTypeCode(numpyDtype):
45 45
46 46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 47 datatype = 0
48 48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 49 datatype = 1
50 50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 51 datatype = 2
52 52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 53 datatype = 3
54 54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 55 datatype = 4
56 56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 57 datatype = 5
58 58 else:
59 59 datatype = None
60 60
61 61 return datatype
62 62
63 63
64 64 def hildebrand_sekhon(data, navg):
65 65 """
66 66 This method is for the objective determination of the noise level in Doppler spectra. This
67 67 implementation technique is based on the fact that the standard deviation of the spectral
68 68 densities is equal to the mean spectral density for white Gaussian noise
69 69
70 70 Inputs:
71 71 Data : heights
72 72 navg : numbers of averages
73 73
74 74 Return:
75 75 mean : noise's level
76 76 """
77 77
78 78 sortdata = numpy.sort(data, axis=None)
79 79 #print(numpy.shape(data))
80 80 #exit()
81 81 '''
82 82 lenOfData = len(sortdata)
83 83 nums_min = lenOfData*0.2
84 84
85 85 if nums_min <= 5:
86 86
87 87 nums_min = 5
88 88
89 89 sump = 0.
90 90 sumq = 0.
91 91
92 92 j = 0
93 93 cont = 1
94 94
95 95 while((cont == 1)and(j < lenOfData)):
96 96
97 97 sump += sortdata[j]
98 98 sumq += sortdata[j]**2
99 99
100 100 if j > nums_min:
101 101 rtest = float(j)/(j-1) + 1.0/navg
102 102 if ((sumq*j) > (rtest*sump**2)):
103 103 j = j - 1
104 104 sump = sump - sortdata[j]
105 105 sumq = sumq - sortdata[j]**2
106 106 cont = 0
107 107
108 108 j += 1
109 109
110 110 lnoise = sump / j
111 111 '''
112 112 return _noise.hildebrand_sekhon(sortdata, navg)
113 113
114 114
115 115 class Beam:
116 116
117 117 def __init__(self):
118 118 self.codeList = []
119 119 self.azimuthList = []
120 120 self.zenithList = []
121 121
122 122
123 123 class GenericData(object):
124 124
125 125 flagNoData = True
126 126
127 127 def copy(self, inputObj=None):
128 128
129 129 if inputObj == None:
130 130 return copy.deepcopy(self)
131 131
132 132 for key in list(inputObj.__dict__.keys()):
133 133
134 134 attribute = inputObj.__dict__[key]
135 135
136 136 # If this attribute is a tuple or list
137 137 if type(inputObj.__dict__[key]) in (tuple, list):
138 138 self.__dict__[key] = attribute[:]
139 139 continue
140 140
141 141 # If this attribute is another object or instance
142 142 if hasattr(attribute, '__dict__'):
143 143 self.__dict__[key] = attribute.copy()
144 144 continue
145 145
146 146 self.__dict__[key] = inputObj.__dict__[key]
147 147
148 148 def deepcopy(self):
149 149
150 150 return copy.deepcopy(self)
151 151
152 152 def isEmpty(self):
153 153
154 154 return self.flagNoData
155 155
156 156 def isReady(self):
157 157
158 158 return not self.flagNoData
159 159
160 160
161 161 class JROData(GenericData):
162 162
163 163 systemHeaderObj = SystemHeader()
164 164 radarControllerHeaderObj = RadarControllerHeader()
165 165 type = None
166 166 datatype = None # dtype but in string
167 167 nProfiles = None
168 168 heightList = None
169 169 channelList = None
170 170 flagDiscontinuousBlock = False
171 171 useLocalTime = False
172 172 utctime = None
173 173 timeZone = None
174 174 dstFlag = None
175 175 errorCount = None
176 176 blocksize = None
177 177 flagDecodeData = False # asumo q la data no esta decodificada
178 178 flagDeflipData = False # asumo q la data no esta sin flip
179 179 flagShiftFFT = False
180 180 nCohInt = None
181 181 windowOfFilter = 1
182 182 C = 3e8
183 183 frequency = 49.92e6
184 184 realtime = False
185 185 beacon_heiIndexList = None
186 186 last_block = None
187 187 blocknow = None
188 188 azimuth = None
189 189 zenith = None
190 190 beam = Beam()
191 191 profileIndex = None
192 192 error = None
193 193 data = None
194 194 nmodes = None
195 195 metadata_list = ['heightList', 'timeZone', 'type']
196 196
197 197 def __str__(self):
198 198
199 199 return '{} - {}'.format(self.type, self.datatime())
200 200
201 201 def getNoise(self):
202 202
203 203 raise NotImplementedError
204 204
205 205 @property
206 206 def nChannels(self):
207 207
208 208 return len(self.channelList)
209 209
210 210 @property
211 211 def channelIndexList(self):
212 212
213 213 return list(range(self.nChannels))
214 214
215 215 @property
216 216 def nHeights(self):
217 217
218 218 return len(self.heightList)
219 219
220 220 def getDeltaH(self):
221 221
222 222 return self.heightList[1] - self.heightList[0]
223 223
224 224 @property
225 225 def ltctime(self):
226 226
227 227 if self.useLocalTime:
228 228 return self.utctime - self.timeZone * 60
229 229
230 230 return self.utctime
231 231
232 232 @property
233 233 def datatime(self):
234 234
235 235 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
236 236 return datatimeValue
237 237
238 238 def getTimeRange(self):
239 239
240 240 datatime = []
241 241
242 242 datatime.append(self.ltctime)
243 243 datatime.append(self.ltctime + self.timeInterval + 1)
244 244
245 245 datatime = numpy.array(datatime)
246 246
247 247 return datatime
248 248
249 249 def getFmaxTimeResponse(self):
250 250
251 251 period = (10 ** -6) * self.getDeltaH() / (0.15)
252 252
253 253 PRF = 1. / (period * self.nCohInt)
254 254
255 255 fmax = PRF
256 256
257 257 return fmax
258 258
259 259 def getFmax(self):
260 260 PRF = 1. / (self.ippSeconds * self.nCohInt)
261 261
262 262 fmax = PRF
263 263 return fmax
264 264
265 265 def getVmax(self):
266 266
267 267 _lambda = self.C / self.frequency
268 268
269 269 vmax = self.getFmax() * _lambda / 2
270 270
271 271 return vmax
272 272
273 273 @property
274 274 def ippSeconds(self):
275 275 '''
276 276 '''
277 277 return self.radarControllerHeaderObj.ippSeconds
278 278
279 279 @ippSeconds.setter
280 280 def ippSeconds(self, ippSeconds):
281 281 '''
282 282 '''
283 283 self.radarControllerHeaderObj.ippSeconds = ippSeconds
284 284
285 285 @property
286 286 def code(self):
287 287 '''
288 288 '''
289 289 return self.radarControllerHeaderObj.code
290 290
291 291 @code.setter
292 292 def code(self, code):
293 293 '''
294 294 '''
295 295 self.radarControllerHeaderObj.code = code
296 296
297 297 @property
298 298 def nCode(self):
299 299 '''
300 300 '''
301 301 return self.radarControllerHeaderObj.nCode
302 302
303 303 @nCode.setter
304 304 def nCode(self, ncode):
305 305 '''
306 306 '''
307 307 self.radarControllerHeaderObj.nCode = ncode
308 308
309 309 @property
310 310 def nBaud(self):
311 311 '''
312 312 '''
313 313 return self.radarControllerHeaderObj.nBaud
314 314
315 315 @nBaud.setter
316 316 def nBaud(self, nbaud):
317 317 '''
318 318 '''
319 319 self.radarControllerHeaderObj.nBaud = nbaud
320 320
321 321 @property
322 322 def ipp(self):
323 323 '''
324 324 '''
325 325 return self.radarControllerHeaderObj.ipp
326 326
327 327 @ipp.setter
328 328 def ipp(self, ipp):
329 329 '''
330 330 '''
331 331 self.radarControllerHeaderObj.ipp = ipp
332 332
333 333 @property
334 334 def metadata(self):
335 335 '''
336 336 '''
337 337
338 338 return {attr: getattr(self, attr) for attr in self.metadata_list}
339 339
340 340
341 341 class Voltage(JROData):
342 342
343 343 dataPP_POW = None
344 344 dataPP_DOP = None
345 345 dataPP_WIDTH = None
346 346 dataPP_SNR = None
347 347
348 348 def __init__(self):
349 349 '''
350 350 Constructor
351 351 '''
352 352
353 353 self.useLocalTime = True
354 354 self.radarControllerHeaderObj = RadarControllerHeader()
355 355 self.systemHeaderObj = SystemHeader()
356 356 self.type = "Voltage"
357 357 self.data = None
358 358 self.nProfiles = None
359 359 self.heightList = None
360 360 self.channelList = None
361 361 self.flagNoData = True
362 362 self.flagDiscontinuousBlock = False
363 363 self.utctime = None
364 364 self.timeZone = 0
365 365 self.dstFlag = None
366 366 self.errorCount = None
367 367 self.nCohInt = None
368 368 self.blocksize = None
369 369 self.flagCohInt = False
370 370 self.flagDecodeData = False # asumo q la data no esta decodificada
371 371 self.flagDeflipData = False # asumo q la data no esta sin flip
372 372 self.flagShiftFFT = False
373 373 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
374 374 self.profileIndex = 0
375 375 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
376 376 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
377 377
378 378 def getNoisebyHildebrand(self, channel=None):
379 379 """
380 380 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
381 381
382 382 Return:
383 383 noiselevel
384 384 """
385 385
386 386 if channel != None:
387 387 data = self.data[channel]
388 388 nChannels = 1
389 389 else:
390 390 data = self.data
391 391 nChannels = self.nChannels
392 392
393 393 noise = numpy.zeros(nChannels)
394 394 power = data * numpy.conjugate(data)
395 395
396 396 for thisChannel in range(nChannels):
397 397 if nChannels == 1:
398 398 daux = power[:].real
399 399 else:
400 400 daux = power[thisChannel, :].real
401 401 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
402 402
403 403 return noise
404 404
405 405 def getNoise(self, type=1, channel=None):
406 406
407 407 if type == 1:
408 408 noise = self.getNoisebyHildebrand(channel)
409 409
410 410 return noise
411 411
412 412 def getPower(self, channel=None):
413 413
414 414 if channel != None:
415 415 data = self.data[channel]
416 416 else:
417 417 data = self.data
418 418
419 419 power = data * numpy.conjugate(data)
420 420 powerdB = 10 * numpy.log10(power.real)
421 421 powerdB = numpy.squeeze(powerdB)
422 422
423 423 return powerdB
424 424
425 425 @property
426 426 def timeInterval(self):
427 427
428 428 return self.ippSeconds * self.nCohInt
429 429
430 430 noise = property(getNoise, "I'm the 'nHeights' property.")
431 431
432 432
433 433 class CrossProds(JROData):
434 434
435 435 # data es un numpy array de 2 dmensiones (canales, alturas)
436 436 data = None
437 437
438 438 def __init__(self):
439 439 '''
440 440 Constructor
441 441 '''
442 442
443 443 self.useLocalTime = True
444 444 '''
445 445 self.radarControllerHeaderObj = RadarControllerHeader()
446 446 self.systemHeaderObj = SystemHeader()
447 447 self.type = "Voltage"
448 448 self.data = None
449 449 # self.dtype = None
450 450 # self.nChannels = 0
451 451 # self.nHeights = 0
452 452 self.nProfiles = None
453 453 self.heightList = None
454 454 self.channelList = None
455 455 # self.channelIndexList = None
456 456 self.flagNoData = True
457 457 self.flagDiscontinuousBlock = False
458 458 self.utctime = None
459 459 self.timeZone = None
460 460 self.dstFlag = None
461 461 self.errorCount = None
462 462 self.nCohInt = None
463 463 self.blocksize = None
464 464 self.flagDecodeData = False # asumo q la data no esta decodificada
465 465 self.flagDeflipData = False # asumo q la data no esta sin flip
466 466 self.flagShiftFFT = False
467 467 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
468 468 self.profileIndex = 0
469 469
470 470
471 471 def getNoisebyHildebrand(self, channel=None):
472 472
473 473
474 474 if channel != None:
475 475 data = self.data[channel]
476 476 nChannels = 1
477 477 else:
478 478 data = self.data
479 479 nChannels = self.nChannels
480 480
481 481 noise = numpy.zeros(nChannels)
482 482 power = data * numpy.conjugate(data)
483 483
484 484 for thisChannel in range(nChannels):
485 485 if nChannels == 1:
486 486 daux = power[:].real
487 487 else:
488 488 daux = power[thisChannel, :].real
489 489 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
490 490
491 491 return noise
492 492
493 493 def getNoise(self, type=1, channel=None):
494 494
495 495 if type == 1:
496 496 noise = self.getNoisebyHildebrand(channel)
497 497
498 498 return noise
499 499
500 500 def getPower(self, channel=None):
501 501
502 502 if channel != None:
503 503 data = self.data[channel]
504 504 else:
505 505 data = self.data
506 506
507 507 power = data * numpy.conjugate(data)
508 508 powerdB = 10 * numpy.log10(power.real)
509 509 powerdB = numpy.squeeze(powerdB)
510 510
511 511 return powerdB
512 512
513 513 def getTimeInterval(self):
514 514
515 515 timeInterval = self.ippSeconds * self.nCohInt
516 516
517 517 return timeInterval
518 518
519 519 noise = property(getNoise, "I'm the 'nHeights' property.")
520 520 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
521 521 '''
522 522 def getTimeInterval(self):
523 523
524 524 timeInterval = self.ippSeconds * self.nCohInt
525 525
526 526 return timeInterval
527 527
528 528
529 529
530 530 class Spectra(JROData):
531 531
532 532 def __init__(self):
533 533 '''
534 534 Constructor
535 535 '''
536 536
537 537 self.data_dc = None
538 538 self.data_spc = None
539 539 self.data_cspc = None
540 540 self.useLocalTime = True
541 541 self.radarControllerHeaderObj = RadarControllerHeader()
542 542 self.systemHeaderObj = SystemHeader()
543 543 self.type = "Spectra"
544 544 self.timeZone = 0
545 545 self.nProfiles = None
546 546 self.heightList = None
547 547 self.channelList = None
548 548 self.pairsList = None
549 549 self.flagNoData = True
550 550 self.flagDiscontinuousBlock = False
551 551 self.utctime = None
552 552 self.nCohInt = None
553 553 self.nIncohInt = None
554 554 self.blocksize = None
555 555 self.nFFTPoints = None
556 556 self.wavelength = None
557 557 self.flagDecodeData = False # asumo q la data no esta decodificada
558 558 self.flagDeflipData = False # asumo q la data no esta sin flip
559 559 self.flagShiftFFT = False
560 560 self.ippFactor = 1
561 561 self.beacon_heiIndexList = []
562 562 self.noise_estimation = None
563 563 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
564 564 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp', 'nIncohInt', 'nFFTPoints', 'nProfiles']
565 565
566 566 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
567 567 """
568 568 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
569 569
570 570 Return:
571 571 noiselevel
572 572 """
573 573
574 574 noise = numpy.zeros(self.nChannels)
575 575
576 576 for channel in range(self.nChannels):
577 577 daux = self.data_spc[channel,
578 578 xmin_index:xmax_index, ymin_index:ymax_index]
579 579 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
580 580
581 581 return noise
582 582
583 583 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
584 584
585 585 if self.noise_estimation is not None:
586 586 # this was estimated by getNoise Operation defined in jroproc_spectra.py
587 587 return self.noise_estimation
588 588 else:
589 589 noise = self.getNoisebyHildebrand(
590 590 xmin_index, xmax_index, ymin_index, ymax_index)
591 591 return noise
592 592
593 593 def getFreqRangeTimeResponse(self, extrapoints=0):
594 594
595 595 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
596 596 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
597 597
598 598 return freqrange
599 599
600 600 def getAcfRange(self, extrapoints=0):
601 601
602 602 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
603 603 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
604 604
605 605 return freqrange
606 606
607 607 def getFreqRange(self, extrapoints=0):
608 608
609 609 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
610 610 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
611 611
612 612 return freqrange
613 613
614 614 def getVelRange(self, extrapoints=0):
615 615
616 616 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
617 617 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
618 618
619 619 if self.nmodes:
620 620 return velrange / self.nmodes
621 621 else:
622 622 return velrange
623 623
624 624 @property
625 625 def nPairs(self):
626 626
627 627 return len(self.pairsList)
628 628
629 629 @property
630 630 def pairsIndexList(self):
631 631
632 632 return list(range(self.nPairs))
633 633
634 634 @property
635 635 def normFactor(self):
636 636
637 637 pwcode = 1
638 638
639 639 if self.flagDecodeData:
640 640 pwcode = numpy.sum(self.code[0] ** 2)
641 641 # normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
642 642 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
643 643
644 644 return normFactor
645 645
646 646 @property
647 647 def flag_cspc(self):
648 648
649 649 if self.data_cspc is None:
650 650 return True
651 651
652 652 return False
653 653
654 654 @property
655 655 def flag_dc(self):
656 656
657 657 if self.data_dc is None:
658 658 return True
659 659
660 660 return False
661 661
662 662 @property
663 663 def timeInterval(self):
664 664
665 665 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
666 666 if self.nmodes:
667 667 return self.nmodes * timeInterval
668 668 else:
669 669 return timeInterval
670 670
671 671 def getPower(self):
672 672
673 673 factor = self.normFactor
674 674 z = self.data_spc / factor
675 675 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
676 676 avg = numpy.average(z, axis=1)
677 677
678 678 return 10 * numpy.log10(avg)
679 679
680 680 def getCoherence(self, pairsList=None, phase=False):
681 681
682 682 z = []
683 683 if pairsList is None:
684 684 pairsIndexList = self.pairsIndexList
685 685 else:
686 686 pairsIndexList = []
687 687 for pair in pairsList:
688 688 if pair not in self.pairsList:
689 689 raise ValueError("Pair %s is not in dataOut.pairsList" % (
690 690 pair))
691 691 pairsIndexList.append(self.pairsList.index(pair))
692 692 for i in range(len(pairsIndexList)):
693 693 pair = self.pairsList[pairsIndexList[i]]
694 694 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
695 695 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
696 696 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
697 697 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
698 698 if phase:
699 699 data = numpy.arctan2(avgcoherenceComplex.imag,
700 700 avgcoherenceComplex.real) * 180 / numpy.pi
701 701 else:
702 702 data = numpy.abs(avgcoherenceComplex)
703 703
704 704 z.append(data)
705 705
706 706 return numpy.array(z)
707 707
708 708 def setValue(self, value):
709 709
710 710 print("This property should not be initialized")
711 711
712 712 return
713 713
714 714 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
715 715
716 716
717 717 class SpectraHeis(Spectra):
718 718
719 719 def __init__(self):
720 720
721 721 self.radarControllerHeaderObj = RadarControllerHeader()
722 722 self.systemHeaderObj = SystemHeader()
723 723 self.type = "SpectraHeis"
724 724 self.nProfiles = None
725 725 self.heightList = None
726 726 self.channelList = None
727 727 self.flagNoData = True
728 728 self.flagDiscontinuousBlock = False
729 729 self.utctime = None
730 730 self.blocksize = None
731 731 self.profileIndex = 0
732 732 self.nCohInt = 1
733 733 self.nIncohInt = 1
734 734
735 735 @property
736 736 def normFactor(self):
737 737 pwcode = 1
738 738 if self.flagDecodeData:
739 739 pwcode = numpy.sum(self.code[0] ** 2)
740 740
741 741 normFactor = self.nIncohInt * self.nCohInt * pwcode
742 742
743 743 return normFactor
744 744
745 745 @property
746 746 def timeInterval(self):
747 747
748 748 return self.ippSeconds * self.nCohInt * self.nIncohInt
749 749
750 750
751 751 class Fits(JROData):
752 752
753 753 def __init__(self):
754 754
755 755 self.type = "Fits"
756 756 self.nProfiles = None
757 757 self.heightList = None
758 758 self.channelList = None
759 759 self.flagNoData = True
760 760 self.utctime = None
761 761 self.nCohInt = 1
762 762 self.nIncohInt = 1
763 763 self.useLocalTime = True
764 764 self.profileIndex = 0
765 765 self.timeZone = 0
766 766
767 767 def getTimeRange(self):
768 768
769 769 datatime = []
770 770
771 771 datatime.append(self.ltctime)
772 772 datatime.append(self.ltctime + self.timeInterval)
773 773
774 774 datatime = numpy.array(datatime)
775 775
776 776 return datatime
777 777
778 778 def getChannelIndexList(self):
779 779
780 780 return list(range(self.nChannels))
781 781
782 782 def getNoise(self, type=1):
783 783
784 784
785 785 if type == 1:
786 786 noise = self.getNoisebyHildebrand()
787 787
788 788 if type == 2:
789 789 noise = self.getNoisebySort()
790 790
791 791 if type == 3:
792 792 noise = self.getNoisebyWindow()
793 793
794 794 return noise
795 795
796 796 @property
797 797 def timeInterval(self):
798 798
799 799 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
800 800
801 801 return timeInterval
802 802
803 803 @property
804 804 def ippSeconds(self):
805 805 '''
806 806 '''
807 807 return self.ipp_sec
808 808
809 809 noise = property(getNoise, "I'm the 'nHeights' property.")
810 810
811 811
812 812 class Correlation(JROData):
813 813
814 814 def __init__(self):
815 815 '''
816 816 Constructor
817 817 '''
818 818 self.radarControllerHeaderObj = RadarControllerHeader()
819 819 self.systemHeaderObj = SystemHeader()
820 820 self.type = "Correlation"
821 821 self.data = None
822 822 self.dtype = None
823 823 self.nProfiles = None
824 824 self.heightList = None
825 825 self.channelList = None
826 826 self.flagNoData = True
827 827 self.flagDiscontinuousBlock = False
828 828 self.utctime = None
829 829 self.timeZone = 0
830 830 self.dstFlag = None
831 831 self.errorCount = None
832 832 self.blocksize = None
833 833 self.flagDecodeData = False # asumo q la data no esta decodificada
834 834 self.flagDeflipData = False # asumo q la data no esta sin flip
835 835 self.pairsList = None
836 836 self.nPoints = None
837 837
838 838 def getPairsList(self):
839 839
840 840 return self.pairsList
841 841
842 842 def getNoise(self, mode=2):
843 843
844 844 indR = numpy.where(self.lagR == 0)[0][0]
845 845 indT = numpy.where(self.lagT == 0)[0][0]
846 846
847 847 jspectra0 = self.data_corr[:, :, indR, :]
848 848 jspectra = copy.copy(jspectra0)
849 849
850 850 num_chan = jspectra.shape[0]
851 851 num_hei = jspectra.shape[2]
852 852
853 853 freq_dc = jspectra.shape[1] / 2
854 854 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
855 855
856 856 if ind_vel[0] < 0:
857 857 ind_vel[list(range(0, 1))] = ind_vel[list(
858 858 range(0, 1))] + self.num_prof
859 859
860 860 if mode == 1:
861 861 jspectra[:, freq_dc, :] = (
862 862 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
863 863
864 864 if mode == 2:
865 865
866 866 vel = numpy.array([-2, -1, 1, 2])
867 867 xx = numpy.zeros([4, 4])
868 868
869 869 for fil in range(4):
870 870 xx[fil, :] = vel[fil] ** numpy.asarray(list(range(4)))
871 871
872 872 xx_inv = numpy.linalg.inv(xx)
873 873 xx_aux = xx_inv[0, :]
874 874
875 875 for ich in range(num_chan):
876 876 yy = jspectra[ich, ind_vel, :]
877 877 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
878 878
879 879 junkid = jspectra[ich, freq_dc, :] <= 0
880 880 cjunkid = sum(junkid)
881 881
882 882 if cjunkid.any():
883 883 jspectra[ich, freq_dc, junkid.nonzero()] = (
884 884 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
885 885
886 886 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
887 887
888 888 return noise
889 889
890 890 @property
891 891 def timeInterval(self):
892 892
893 893 return self.ippSeconds * self.nCohInt * self.nProfiles
894 894
895 895 def splitFunctions(self):
896 896
897 897 pairsList = self.pairsList
898 898 ccf_pairs = []
899 899 acf_pairs = []
900 900 ccf_ind = []
901 901 acf_ind = []
902 902 for l in range(len(pairsList)):
903 903 chan0 = pairsList[l][0]
904 904 chan1 = pairsList[l][1]
905 905
906 906 # Obteniendo pares de Autocorrelacion
907 907 if chan0 == chan1:
908 908 acf_pairs.append(chan0)
909 909 acf_ind.append(l)
910 910 else:
911 911 ccf_pairs.append(pairsList[l])
912 912 ccf_ind.append(l)
913 913
914 914 data_acf = self.data_cf[acf_ind]
915 915 data_ccf = self.data_cf[ccf_ind]
916 916
917 917 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
918 918
919 919 @property
920 920 def normFactor(self):
921 921 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
922 922 acf_pairs = numpy.array(acf_pairs)
923 923 normFactor = numpy.zeros((self.nPairs, self.nHeights))
924 924
925 925 for p in range(self.nPairs):
926 926 pair = self.pairsList[p]
927 927
928 928 ch0 = pair[0]
929 929 ch1 = pair[1]
930 930
931 931 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
932 932 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
933 933 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
934 934
935 935 return normFactor
936 936
937 937
938 938 class Parameters(Spectra):
939 939
940 940 groupList = None # List of Pairs, Groups, etc
941 941 data_param = None # Parameters obtained
942 942 data_pre = None # Data Pre Parametrization
943 943 data_SNR = None # Signal to Noise Ratio
944 944 abscissaList = None # Abscissa, can be velocities, lags or time
945 945 utctimeInit = None # Initial UTC time
946 946 paramInterval = None # Time interval to calculate Parameters in seconds
947 947 useLocalTime = True
948 948 # Fitting
949 949 data_error = None # Error of the estimation
950 950 constants = None
951 951 library = None
952 952 # Output signal
953 953 outputInterval = None # Time interval to calculate output signal in seconds
954 954 data_output = None # Out signal
955 955 nAvg = None
956 956 noise_estimation = None
957 957 GauSPC = None # Fit gaussian SPC
958 958
959 959 def __init__(self):
960 960 '''
961 961 Constructor
962 962 '''
963 963 self.radarControllerHeaderObj = RadarControllerHeader()
964 964 self.systemHeaderObj = SystemHeader()
965 965 self.type = "Parameters"
966 966 self.timeZone = 0
967 self.ippFactor = 1
967 968
968 969 def getTimeRange1(self, interval):
969 970
970 971 datatime = []
971 972
972 973 if self.useLocalTime:
973 974 time1 = self.utctimeInit - self.timeZone * 60
974 975 else:
975 976 time1 = self.utctimeInit
976 977
977 978 datatime.append(time1)
978 979 datatime.append(time1 + interval)
979 980 datatime = numpy.array(datatime)
980 981
981 982 return datatime
982 983
983 984 @property
984 985 def timeInterval(self):
985 986
986 987 if hasattr(self, 'timeInterval1'):
987 988 return self.timeInterval1
988 989 else:
989 990 return self.paramInterval
990 991
991 992
992 993 def setValue(self, value):
993 994
994 995 print("This property should not be initialized")
995 996
996 997 return
997 998
998 999 def getNoise(self):
999 1000
1000 1001 return self.spc_noise
1001 1002
1002 1003 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1003 1004
1004 1005
1005 1006 class PlotterData(object):
1006 1007 '''
1007 1008 Object to hold data to be plotted
1008 1009 '''
1009 1010
1010 1011 MAXNUMX = 200
1011 1012 MAXNUMY = 200
1012 1013
1013 1014 def __init__(self, code, exp_code, localtime=True):
1014 1015
1015 1016 self.key = code
1016 1017 self.exp_code = exp_code
1017 1018 self.ready = False
1018 1019 self.flagNoData = False
1019 1020 self.localtime = localtime
1020 1021 self.data = {}
1021 1022 self.meta = {}
1022 1023 self.__heights = []
1023 1024
1024 1025 def __str__(self):
1025 1026 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1026 1027 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1027 1028
1028 1029 def __len__(self):
1029 1030 return len(self.data)
1030 1031
1031 1032 def __getitem__(self, key):
1032 1033 if isinstance(key, int):
1033 1034 return self.data[self.times[key]]
1034 1035 elif isinstance(key, str):
1035 1036 ret = numpy.array([self.data[x][key] for x in self.times])
1036 1037 if ret.ndim > 1:
1037 1038 ret = numpy.swapaxes(ret, 0, 1)
1038 1039 return ret
1039 1040
1040 1041 def __contains__(self, key):
1041 1042 return key in self.data[self.min_time]
1042 1043
1043 1044 def setup(self):
1044 1045 '''
1045 1046 Configure object
1046 1047 '''
1047 1048 self.type = ''
1048 1049 self.ready = False
1049 1050 del self.data
1050 1051 self.data = {}
1051 1052 self.__heights = []
1052 1053 self.__all_heights = set()
1053 1054
1054 1055 def shape(self, key):
1055 1056 '''
1056 1057 Get the shape of the one-element data for the given key
1057 1058 '''
1058 1059
1059 1060 if len(self.data[self.min_time][key]):
1060 1061 return self.data[self.min_time][key].shape
1061 1062 return (0,)
1062 1063
1063 1064 def update(self, data, tm, meta={}):
1064 1065 '''
1065 1066 Update data object with new dataOut
1066 1067 '''
1067 1068
1068 1069 self.data[tm] = data
1069 1070
1070 1071 for key, value in meta.items():
1071 1072 setattr(self, key, value)
1072 1073
1073 1074 def normalize_heights(self):
1074 1075 '''
1075 1076 Ensure same-dimension of the data for different heighList
1076 1077 '''
1077 1078
1078 1079 H = numpy.array(list(self.__all_heights))
1079 1080 H.sort()
1080 1081 for key in self.data:
1081 1082 shape = self.shape(key)[:-1] + H.shape
1082 1083 for tm, obj in list(self.data[key].items()):
1083 1084 h = self.__heights[self.times.tolist().index(tm)]
1084 1085 if H.size == h.size:
1085 1086 continue
1086 1087 index = numpy.where(numpy.in1d(H, h))[0]
1087 1088 dummy = numpy.zeros(shape) + numpy.nan
1088 1089 if len(shape) == 2:
1089 1090 dummy[:, index] = obj
1090 1091 else:
1091 1092 dummy[index] = obj
1092 1093 self.data[key][tm] = dummy
1093 1094
1094 1095 self.__heights = [H for tm in self.times]
1095 1096
1096 1097 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1097 1098 '''
1098 1099 Convert data to json
1099 1100 '''
1100 1101
1101 1102 meta = {}
1102 1103 meta['xrange'] = []
1103 1104 dy = int(len(self.yrange) / self.MAXNUMY) + 1
1104 1105 tmp = self.data[tm][self.key]
1105 1106 shape = tmp.shape
1106 1107 if len(shape) == 2:
1107 1108 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1108 1109 elif len(shape) == 3:
1109 1110 dx = int(self.data[tm][self.key].shape[1] / self.MAXNUMX) + 1
1110 1111 data = self.roundFloats(
1111 1112 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1112 1113 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1113 1114 else:
1114 1115 data = self.roundFloats(self.data[tm][self.key].tolist())
1115 1116
1116 1117 ret = {
1117 1118 'plot': plot_name,
1118 1119 'code': self.exp_code,
1119 1120 'time': float(tm),
1120 1121 'data': data,
1121 1122 }
1122 1123 meta['type'] = plot_type
1123 1124 meta['interval'] = float(self.interval)
1124 1125 meta['localtime'] = self.localtime
1125 1126 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1126 1127 meta.update(self.meta)
1127 1128 ret['metadata'] = meta
1128 1129 return json.dumps(ret)
1129 1130
1130 1131 @property
1131 1132 def times(self):
1132 1133 '''
1133 1134 Return the list of times of the current data
1134 1135 '''
1135 1136
1136 1137 ret = [t for t in self.data]
1137 1138 ret.sort()
1138 1139 return numpy.array(ret)
1139 1140
1140 1141 @property
1141 1142 def min_time(self):
1142 1143 '''
1143 1144 Return the minimun time value
1144 1145 '''
1145 1146
1146 1147 return self.times[0]
1147 1148
1148 1149 @property
1149 1150 def max_time(self):
1150 1151 '''
1151 1152 Return the maximun time value
1152 1153 '''
1153 1154
1154 1155 return self.times[-1]
1155 1156
1156 1157 # @property
1157 1158 # def heights(self):
1158 1159 # '''
1159 1160 # Return the list of heights of the current data
1160 1161 # '''
1161 1162
1162 1163 # return numpy.array(self.__heights[-1])
1163 1164
1164 1165 @staticmethod
1165 1166 def roundFloats(obj):
1166 1167 if isinstance(obj, list):
1167 1168 return list(map(PlotterData.roundFloats, obj))
1168 1169 elif isinstance(obj, float):
1169 1170 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,909 +1,908
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROHeaderIO.py 151 2012-10-31 19:00:51Z murco $
5 5 '''
6 6 import sys
7 7 import numpy
8 8 import copy
9 9 import datetime
10 10 import inspect
11 11 from schainpy.utils import log
12 12
13 13 SPEED_OF_LIGHT = 299792458
14 14 SPEED_OF_LIGHT = 3e8
15 15
16 16 BASIC_STRUCTURE = numpy.dtype([
17 17 ('nSize', '<u4'),
18 18 ('nVersion', '<u2'),
19 19 ('nDataBlockId', '<u4'),
20 20 ('nUtime', '<u4'),
21 21 ('nMilsec', '<u2'),
22 22 ('nTimezone', '<i2'),
23 23 ('nDstflag', '<i2'),
24 24 ('nErrorCount', '<u4')
25 25 ])
26 26
27 27 SYSTEM_STRUCTURE = numpy.dtype([
28 28 ('nSize', '<u4'),
29 29 ('nNumSamples', '<u4'),
30 30 ('nNumProfiles', '<u4'),
31 31 ('nNumChannels', '<u4'),
32 32 ('nADCResolution', '<u4'),
33 33 ('nPCDIOBusWidth', '<u4'),
34 34 ])
35 35
36 36 RADAR_STRUCTURE = numpy.dtype([
37 37 ('nSize', '<u4'),
38 38 ('nExpType', '<u4'),
39 39 ('nNTx', '<u4'),
40 40 ('fIpp', '<f4'),
41 41 ('fTxA', '<f4'),
42 42 ('fTxB', '<f4'),
43 43 ('nNumWindows', '<u4'),
44 44 ('nNumTaus', '<u4'),
45 45 ('nCodeType', '<u4'),
46 46 ('nLine6Function', '<u4'),
47 47 ('nLine5Function', '<u4'),
48 48 ('fClock', '<f4'),
49 49 ('nPrePulseBefore', '<u4'),
50 50 ('nPrePulseAfter', '<u4'),
51 51 ('sRangeIPP', '<a20'),
52 52 ('sRangeTxA', '<a20'),
53 53 ('sRangeTxB', '<a20'),
54 54 ])
55 55
56 56 SAMPLING_STRUCTURE = numpy.dtype(
57 57 [('h0', '<f4'), ('dh', '<f4'), ('nsa', '<u4')])
58 58
59 59
60 60 PROCESSING_STRUCTURE = numpy.dtype([
61 61 ('nSize', '<u4'),
62 62 ('nDataType', '<u4'),
63 63 ('nSizeOfDataBlock', '<u4'),
64 64 ('nProfilesperBlock', '<u4'),
65 65 ('nDataBlocksperFile', '<u4'),
66 66 ('nNumWindows', '<u4'),
67 67 ('nProcessFlags', '<u4'),
68 68 ('nCoherentIntegrations', '<u4'),
69 69 ('nIncoherentIntegrations', '<u4'),
70 70 ('nTotalSpectra', '<u4')
71 71 ])
72 72
73 73
74 74 class Header(object):
75 75
76 76 def __init__(self):
77 77 raise NotImplementedError
78 78
79 79 def copy(self):
80 80 return copy.deepcopy(self)
81 81
82 82 def read(self):
83 83
84 84 raise NotImplementedError
85 85
86 86 def write(self):
87 87
88 88 raise NotImplementedError
89 89
90 90 def getAllowedArgs(self):
91 91 args = inspect.getargspec(self.__init__).args
92 92 try:
93 93 args.remove('self')
94 94 except:
95 95 pass
96 96 return args
97 97
98 98 def getAsDict(self):
99 99 args = self.getAllowedArgs()
100 100 asDict = {}
101 101 for x in args:
102 102 asDict[x] = self[x]
103 103 return asDict
104 104
105 105 def __getitem__(self, name):
106 106 return getattr(self, name)
107 107
108 108 def printInfo(self):
109 109
110 110 message = "#" * 50 + "\n"
111 111 message += self.__class__.__name__.upper() + "\n"
112 112 message += "#" * 50 + "\n"
113 113
114 114 keyList = list(self.__dict__.keys())
115 115 keyList.sort()
116 116
117 117 for key in keyList:
118 118 message += "%s = %s" % (key, self.__dict__[key]) + "\n"
119 119
120 120 if "size" not in keyList:
121 121 attr = getattr(self, "size")
122 122
123 123 if attr:
124 124 message += "%s = %s" % ("size", attr) + "\n"
125 125
126 126 print(message)
127 127
128 128
129 129 class BasicHeader(Header):
130 130
131 131 size = None
132 132 version = None
133 133 dataBlock = None
134 134 utc = None
135 135 ltc = None
136 136 miliSecond = None
137 137 timeZone = None
138 138 dstFlag = None
139 139 errorCount = None
140 140 F = None
141 141 structure = BASIC_STRUCTURE
142 142 __LOCALTIME = None
143 143
144 144 def __init__(self, useLocalTime=True):
145 145
146 146 self.size = 24
147 147 self.version = 0
148 148 self.dataBlock = 0
149 149 self.utc = 0
150 150 self.miliSecond = 0
151 151 self.timeZone = 0
152 152 self.dstFlag = 0
153 153 self.errorCount = 0
154
155 154 self.useLocalTime = useLocalTime
156 155
157 156 def read(self, fp):
158 157
159 158 self.length = 0
160 159 try:
161 160 if hasattr(fp, 'read'):
162 161 header = numpy.fromfile(fp, BASIC_STRUCTURE, 1)
163 162 else:
164 163 header = numpy.fromstring(fp, BASIC_STRUCTURE, 1)
165 164 except Exception as e:
166 165 print("BasicHeader: ")
167 166 print(e)
168 167 return 0
169 168
170 169 self.size = int(header['nSize'][0])
171 170 self.version = int(header['nVersion'][0])
172 171 self.dataBlock = int(header['nDataBlockId'][0])
173 172 self.utc = int(header['nUtime'][0])
174 173 self.miliSecond = int(header['nMilsec'][0])
175 174 self.timeZone = int(header['nTimezone'][0])
176 175 self.dstFlag = int(header['nDstflag'][0])
177 176 self.errorCount = int(header['nErrorCount'][0])
178 177
179 178 if self.size < 24:
180 179 return 0
181 180
182 181 self.length = header.nbytes
183 182 return 1
184 183
185 184 def write(self, fp):
186 185
187 186 headerTuple = (self.size, self.version, self.dataBlock, self.utc,
188 187 self.miliSecond, self.timeZone, self.dstFlag, self.errorCount)
189 188 header = numpy.array(headerTuple, BASIC_STRUCTURE)
190 189 header.tofile(fp)
191 190
192 191 return 1
193 192
194 193 def get_ltc(self):
195 194
196 195 return self.utc - self.timeZone * 60
197 196
198 197 def set_ltc(self, value):
199 198
200 199 self.utc = value + self.timeZone * 60
201 200
202 201 def get_datatime(self):
203 202
204 203 return datetime.datetime.utcfromtimestamp(self.ltc)
205 204
206 205 ltc = property(get_ltc, set_ltc)
207 206 datatime = property(get_datatime)
208 207
209 208
210 209 class SystemHeader(Header):
211 210
212 211 size = None
213 212 nSamples = None
214 213 nProfiles = None
215 214 nChannels = None
216 215 adcResolution = None
217 216 pciDioBusWidth = None
218 217 structure = SYSTEM_STRUCTURE
219 218
220 219 def __init__(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWidth=0):
221 220
222 221 self.size = 24
223 222 self.nSamples = nSamples
224 223 self.nProfiles = nProfiles
225 224 self.nChannels = nChannels
226 225 self.adcResolution = adcResolution
227 226 self.pciDioBusWidth = pciDioBusWidth
228 227
229 228 def read(self, fp):
230 229 self.length = 0
231 230 try:
232 231 startFp = fp.tell()
233 232 except Exception as e:
234 233 startFp = None
235 234 pass
236 235
237 236 try:
238 237 if hasattr(fp, 'read'):
239 238 header = numpy.fromfile(fp, SYSTEM_STRUCTURE, 1)
240 239 else:
241 240 header = numpy.fromstring(fp, SYSTEM_STRUCTURE, 1)
242 241 except Exception as e:
243 242 print("System Header: " + str(e))
244 243 return 0
245 244
246 245 self.size = header['nSize'][0]
247 246 self.nSamples = header['nNumSamples'][0]
248 247 self.nProfiles = header['nNumProfiles'][0]
249 248 self.nChannels = header['nNumChannels'][0]
250 249 self.adcResolution = header['nADCResolution'][0]
251 250 self.pciDioBusWidth = header['nPCDIOBusWidth'][0]
252 251
253 252 if startFp is not None:
254 253 endFp = self.size + startFp
255 254
256 255 if fp.tell() > endFp:
257 256 sys.stderr.write(
258 257 "Warning %s: Size value read from System Header is lower than it has to be\n" % fp.name)
259 258 return 0
260 259
261 260 if fp.tell() < endFp:
262 261 sys.stderr.write(
263 262 "Warning %s: Size value read from System Header size is greater than it has to be\n" % fp.name)
264 263 return 0
265 264
266 265 self.length = header.nbytes
267 266 return 1
268 267
269 268 def write(self, fp):
270 269
271 270 headerTuple = (self.size, self.nSamples, self.nProfiles,
272 271 self.nChannels, self.adcResolution, self.pciDioBusWidth)
273 272 header = numpy.array(headerTuple, SYSTEM_STRUCTURE)
274 273 header.tofile(fp)
275 274
276 275 return 1
277 276
278 277
279 278 class RadarControllerHeader(Header):
280 279
281 280 expType = None
282 281 nTx = None
283 282 ipp = None
284 283 txA = None
285 284 txB = None
286 285 nWindows = None
287 286 numTaus = None
288 287 codeType = None
289 288 line6Function = None
290 289 line5Function = None
291 290 fClock = None
292 291 prePulseBefore = None
293 292 prePulseAfter = None
294 293 rangeIpp = None
295 294 rangeTxA = None
296 295 rangeTxB = None
297 296 structure = RADAR_STRUCTURE
298 297 __size = None
299 298
300 299 def __init__(self, expType=2, nTx=1,
301 300 ipp=None, txA=0, txB=0,
302 301 nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None,
303 302 numTaus=0, line6Function=0, line5Function=0, fClock=None,
304 303 prePulseBefore=0, prePulseAfter=0,
305 304 codeType=0, nCode=0, nBaud=0, code=[],
306 305 flip1=0, flip2=0):
307 306
308 307 # self.size = 116
309 308 self.expType = expType
310 309 self.nTx = nTx
311 310 self.ipp = ipp
312 311 self.txA = txA
313 312 self.txB = txB
314 313 self.rangeIpp = ipp
315 314 self.rangeTxA = txA
316 315 self.rangeTxB = txB
317 316
318 317 self.nWindows = nWindows
319 318 self.numTaus = numTaus
320 319 self.codeType = codeType
321 320 self.line6Function = line6Function
322 321 self.line5Function = line5Function
323 322 self.fClock = fClock
324 323 self.prePulseBefore = prePulseBefore
325 324 self.prePulseAfter = prePulseAfter
326 325
327 326 self.nHeights = nHeights
328 327 self.firstHeight = firstHeight
329 328 self.deltaHeight = deltaHeight
330 329 self.samplesWin = nHeights
331 330
332 331 self.nCode = nCode
333 332 self.nBaud = nBaud
334 333 self.code = code
335 334 self.flip1 = flip1
336 335 self.flip2 = flip2
337 336
338 337 self.code_size = int(numpy.ceil(self.nBaud / 32.)) * self.nCode * 4
339 338 # self.dynamic = numpy.array([],numpy.dtype('byte'))
340 339
341 340 if self.fClock is None and self.deltaHeight is not None:
342 341 self.fClock = 0.15 / (deltaHeight * 1e-6) # 0.15Km / (height * 1u)
343 342
344 343 def read(self, fp):
345 344 self.length = 0
346 345 try:
347 346 startFp = fp.tell()
348 347 except Exception as e:
349 348 startFp = None
350 349 pass
351 350
352 351 try:
353 352 if hasattr(fp, 'read'):
354 353 header = numpy.fromfile(fp, RADAR_STRUCTURE, 1)
355 354 else:
356 355 header = numpy.fromstring(fp, RADAR_STRUCTURE, 1)
357 356 self.length += header.nbytes
358 357 except Exception as e:
359 358 print("RadarControllerHeader: " + str(e))
360 359 return 0
361 360
362 361 size = int(header['nSize'][0])
363 362 self.expType = int(header['nExpType'][0])
364 363 self.nTx = int(header['nNTx'][0])
365 364 self.ipp = float(header['fIpp'][0])
366 365 #print(self.ipp)
367 366 self.txA = float(header['fTxA'][0])
368 367 self.txB = float(header['fTxB'][0])
369 368 self.nWindows = int(header['nNumWindows'][0])
370 369 self.numTaus = int(header['nNumTaus'][0])
371 370 self.codeType = int(header['nCodeType'][0])
372 371 self.line6Function = int(header['nLine6Function'][0])
373 372 self.line5Function = int(header['nLine5Function'][0])
374 373 self.fClock = float(header['fClock'][0])
375 374 self.prePulseBefore = int(header['nPrePulseBefore'][0])
376 375 self.prePulseAfter = int(header['nPrePulseAfter'][0])
377 376 self.rangeIpp = header['sRangeIPP'][0]
378 377 self.rangeTxA = header['sRangeTxA'][0]
379 378 self.rangeTxB = header['sRangeTxB'][0]
380 379
381 380 try:
382 381 if hasattr(fp, 'read'):
383 382 samplingWindow = numpy.fromfile(
384 383 fp, SAMPLING_STRUCTURE, self.nWindows)
385 384 else:
386 385 samplingWindow = numpy.fromstring(
387 386 fp[self.length:], SAMPLING_STRUCTURE, self.nWindows)
388 387 self.length += samplingWindow.nbytes
389 388 except Exception as e:
390 389 print("RadarControllerHeader: " + str(e))
391 390 return 0
392 391 self.nHeights = int(numpy.sum(samplingWindow['nsa']))
393 392 self.firstHeight = samplingWindow['h0']
394 393 self.deltaHeight = samplingWindow['dh']
395 394 self.samplesWin = samplingWindow['nsa']
396 395
397 396 try:
398 397 if hasattr(fp, 'read'):
399 398 self.Taus = numpy.fromfile(fp, '<f4', self.numTaus)
400 399 else:
401 400 self.Taus = numpy.fromstring(
402 401 fp[self.length:], '<f4', self.numTaus)
403 402 self.length += self.Taus.nbytes
404 403 except Exception as e:
405 404 print("RadarControllerHeader: " + str(e))
406 405 return 0
407 406
408 407 self.code_size = 0
409 408 if self.codeType != 0:
410 409
411 410 try:
412 411 if hasattr(fp, 'read'):
413 412 self.nCode = numpy.fromfile(fp, '<u4', 1)[0]
414 413 self.length += self.nCode.nbytes
415 414 self.nBaud = numpy.fromfile(fp, '<u4', 1)[0]
416 415 self.length += self.nBaud.nbytes
417 416 else:
418 417 self.nCode = numpy.fromstring(
419 418 fp[self.length:], '<u4', 1)[0]
420 419 self.length += self.nCode.nbytes
421 420 self.nBaud = numpy.fromstring(
422 421 fp[self.length:], '<u4', 1)[0]
423 422 self.length += self.nBaud.nbytes
424 423 except Exception as e:
425 424 print("RadarControllerHeader: " + str(e))
426 425 return 0
427 426 code = numpy.empty([self.nCode, self.nBaud], dtype='i1')
428 427
429 428 for ic in range(self.nCode):
430 429 try:
431 430 if hasattr(fp, 'read'):
432 431 temp = numpy.fromfile(fp, 'u4', int(
433 432 numpy.ceil(self.nBaud / 32.)))
434 433 else:
435 434 temp = numpy.fromstring(
436 435 fp, 'u4', int(numpy.ceil(self.nBaud / 32.)))
437 436 self.length += temp.nbytes
438 437 except Exception as e:
439 438 print("RadarControllerHeader: " + str(e))
440 439 return 0
441 440
442 441 for ib in range(self.nBaud - 1, -1, -1):
443 442 code[ic, ib] = temp[int(ib / 32)] % 2
444 443 temp[int(ib / 32)] = temp[int(ib / 32)] / 2
445 444
446 445 self.code = 2.0 * code - 1.0
447 446 self.code_size = int(numpy.ceil(self.nBaud / 32.)) * self.nCode * 4
448 447
449 448 # if self.line5Function == RCfunction.FLIP:
450 449 # self.flip1 = numpy.fromfile(fp,'<u4',1)
451 450 #
452 451 # if self.line6Function == RCfunction.FLIP:
453 452 # self.flip2 = numpy.fromfile(fp,'<u4',1)
454 453 if startFp is not None:
455 454 endFp = size + startFp
456 455
457 456 if fp.tell() != endFp:
458 457 # fp.seek(endFp)
459 458 print("%s: Radar Controller Header size is not consistent: from data [%d] != from header field [%d]" % (fp.name, fp.tell() - startFp, size))
460 459 # return 0
461 460
462 461 if fp.tell() > endFp:
463 462 sys.stderr.write(
464 463 "Warning %s: Size value read from Radar Controller header is lower than it has to be\n" % fp.name)
465 464 # return 0
466 465
467 466 if fp.tell() < endFp:
468 467 sys.stderr.write(
469 468 "Warning %s: Size value read from Radar Controller header is greater than it has to be\n" % fp.name)
470 469
471 470 return 1
472 471
473 472 def write(self, fp):
474 473
475 474 headerTuple = (self.size,
476 475 self.expType,
477 476 self.nTx,
478 477 self.ipp,
479 478 self.txA,
480 479 self.txB,
481 480 self.nWindows,
482 481 self.numTaus,
483 482 self.codeType,
484 483 self.line6Function,
485 484 self.line5Function,
486 485 self.fClock,
487 486 self.prePulseBefore,
488 487 self.prePulseAfter,
489 488 self.rangeIpp,
490 489 self.rangeTxA,
491 490 self.rangeTxB)
492 491
493 492 header = numpy.array(headerTuple, RADAR_STRUCTURE)
494 493 header.tofile(fp)
495 494
496 495 sampleWindowTuple = (
497 496 self.firstHeight, self.deltaHeight, self.samplesWin)
498 497 samplingWindow = numpy.array(sampleWindowTuple, SAMPLING_STRUCTURE)
499 498 samplingWindow.tofile(fp)
500 499
501 500 if self.numTaus > 0:
502 501 self.Taus.tofile(fp)
503 502
504 503 if self.codeType != 0:
505 504 nCode = numpy.array(self.nCode, '<u4')
506 505 nCode.tofile(fp)
507 506 nBaud = numpy.array(self.nBaud, '<u4')
508 507 nBaud.tofile(fp)
509 508 code1 = (self.code + 1.0) / 2.
510 509
511 510 for ic in range(self.nCode):
512 511 tempx = numpy.zeros(int(numpy.ceil(self.nBaud / 32.)))
513 512 start = 0
514 513 end = 32
515 514 for i in range(len(tempx)):
516 515 code_selected = code1[ic, start:end]
517 516 for j in range(len(code_selected) - 1, -1, -1):
518 517 if code_selected[j] == 1:
519 518 tempx[i] = tempx[i] + \
520 519 2 ** (len(code_selected) - 1 - j)
521 520 start = start + 32
522 521 end = end + 32
523 522
524 523 tempx = tempx.astype('u4')
525 524 tempx.tofile(fp)
526 525
527 526 # if self.line5Function == RCfunction.FLIP:
528 527 # self.flip1.tofile(fp)
529 528 #
530 529 # if self.line6Function == RCfunction.FLIP:
531 530 # self.flip2.tofile(fp)
532 531
533 532 return 1
534 533
535 534 def get_ippSeconds(self):
536 535 '''
537 536 '''
538 537
539 538 ippSeconds = 2.0 * 1000 * self.ipp / SPEED_OF_LIGHT
540 539
541 540 return ippSeconds
542 541
543 542 def set_ippSeconds(self, ippSeconds):
544 543 '''
545 544 '''
546 545
547 546 self.ipp = ippSeconds * SPEED_OF_LIGHT / (2.0 * 1000)
548 547
549 548 return
550 549
551 550 def get_size(self):
552 551
553 552 self.__size = 116 + 12 * self.nWindows + 4 * self.numTaus
554 553
555 554 if self.codeType != 0:
556 555 self.__size += 4 + 4 + 4 * self.nCode * \
557 556 numpy.ceil(self.nBaud / 32.)
558 557
559 558 return self.__size
560 559
561 560 def set_size(self, value):
562 561
563 562 raise IOError("size is a property and it cannot be set, just read")
564 563
565 564 return
566 565
567 566 ippSeconds = property(get_ippSeconds, set_ippSeconds)
568 567 size = property(get_size, set_size)
569 568
570 569
571 570 class ProcessingHeader(Header):
572 571
573 572 # size = None
574 573 dtype = None
575 574 blockSize = None
576 575 profilesPerBlock = None
577 576 dataBlocksPerFile = None
578 577 nWindows = None
579 578 processFlags = None
580 579 nCohInt = None
581 580 nIncohInt = None
582 581 totalSpectra = None
583 582 structure = PROCESSING_STRUCTURE
584 583 flag_dc = None
585 584 flag_cspc = None
586 585
587 586 def __init__(self, dtype=0, blockSize=0, profilesPerBlock=0, dataBlocksPerFile=0, nWindows=0, processFlags=0, nCohInt=0,
588 587 nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0, deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0,
589 588 code=0, nBaud=None, shif_fft=False, flag_dc=False, flag_cspc=False, flag_decode=False, flag_deflip=False
590 589 ):
591 590
592 591 # self.size = 0
593 592 self.dtype = dtype
594 593 self.blockSize = blockSize
595 594 self.profilesPerBlock = 0
596 595 self.dataBlocksPerFile = 0
597 596 self.nWindows = 0
598 597 self.processFlags = 0
599 598 self.nCohInt = 0
600 599 self.nIncohInt = 0
601 600 self.totalSpectra = 0
602 601
603 602 self.nHeights = 0
604 603 self.firstHeight = 0
605 604 self.deltaHeight = 0
606 605 self.samplesWin = 0
607 606 self.spectraComb = 0
608 607 self.nCode = None
609 608 self.code = None
610 609 self.nBaud = None
611 610
612 611 self.shif_fft = False
613 612 self.flag_dc = False
614 613 self.flag_cspc = False
615 614 self.flag_decode = False
616 615 self.flag_deflip = False
617 616 self.length = 0
618 617
619 618 def read(self, fp):
620 619 self.length = 0
621 620 try:
622 621 startFp = fp.tell()
623 622 except Exception as e:
624 623 startFp = None
625 624 pass
626 625
627 626 try:
628 627 if hasattr(fp, 'read'):
629 628 header = numpy.fromfile(fp, PROCESSING_STRUCTURE, 1)
630 629 else:
631 630 header = numpy.fromstring(fp, PROCESSING_STRUCTURE, 1)
632 631 self.length += header.nbytes
633 632 except Exception as e:
634 633 print("ProcessingHeader: " + str(e))
635 634 return 0
636 635
637 636 size = int(header['nSize'][0])
638 637 self.dtype = int(header['nDataType'][0])
639 638 self.blockSize = int(header['nSizeOfDataBlock'][0])
640 639 self.profilesPerBlock = int(header['nProfilesperBlock'][0])
641 640 self.dataBlocksPerFile = int(header['nDataBlocksperFile'][0])
642 641 self.nWindows = int(header['nNumWindows'][0])
643 642 self.processFlags = header['nProcessFlags']
644 643 self.nCohInt = int(header['nCoherentIntegrations'][0])
645 644
646 645 self.nIncohInt = int(header['nIncoherentIntegrations'][0])
647 646 self.totalSpectra = int(header['nTotalSpectra'][0])
648 647
649 648 try:
650 649 if hasattr(fp, 'read'):
651 650 samplingWindow = numpy.fromfile(
652 651 fp, SAMPLING_STRUCTURE, self.nWindows)
653 652 else:
654 653 samplingWindow = numpy.fromstring(
655 654 fp[self.length:], SAMPLING_STRUCTURE, self.nWindows)
656 655 self.length += samplingWindow.nbytes
657 656 except Exception as e:
658 657 print("ProcessingHeader: " + str(e))
659 658 return 0
660 659
661 660 self.nHeights = int(numpy.sum(samplingWindow['nsa']))
662 661 self.firstHeight = float(samplingWindow['h0'][0])
663 662 self.deltaHeight = float(samplingWindow['dh'][0])
664 663 self.samplesWin = samplingWindow['nsa'][0]
665 664
666 665 try:
667 666 if hasattr(fp, 'read'):
668 667 self.spectraComb = numpy.fromfile(
669 668 fp, 'u1', 2 * self.totalSpectra)
670 669 else:
671 670 self.spectraComb = numpy.fromstring(
672 671 fp[self.length:], 'u1', 2 * self.totalSpectra)
673 672 self.length += self.spectraComb.nbytes
674 673 except Exception as e:
675 674 print("ProcessingHeader: " + str(e))
676 675 return 0
677 676
678 677 if ((self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE) == PROCFLAG.DEFINE_PROCESS_CODE):
679 678 self.nCode = int(numpy.fromfile(fp, '<u4', 1))
680 679 self.nBaud = int(numpy.fromfile(fp, '<u4', 1))
681 680 self.code = numpy.fromfile(
682 681 fp, '<f4', self.nCode * self.nBaud).reshape(self.nCode, self.nBaud)
683 682
684 683 if ((self.processFlags & PROCFLAG.EXP_NAME_ESP) == PROCFLAG.EXP_NAME_ESP):
685 684 exp_name_len = int(numpy.fromfile(fp, '<u4', 1))
686 685 exp_name = numpy.fromfile(fp, 'u1', exp_name_len + 1)
687 686
688 687 if ((self.processFlags & PROCFLAG.SHIFT_FFT_DATA) == PROCFLAG.SHIFT_FFT_DATA):
689 688 self.shif_fft = True
690 689 else:
691 690 self.shif_fft = False
692 691
693 692 if ((self.processFlags & PROCFLAG.SAVE_CHANNELS_DC) == PROCFLAG.SAVE_CHANNELS_DC):
694 693 self.flag_dc = True
695 694 else:
696 695 self.flag_dc = False
697 696
698 697 if ((self.processFlags & PROCFLAG.DECODE_DATA) == PROCFLAG.DECODE_DATA):
699 698 self.flag_decode = True
700 699 else:
701 700 self.flag_decode = False
702 701
703 702 if ((self.processFlags & PROCFLAG.DEFLIP_DATA) == PROCFLAG.DEFLIP_DATA):
704 703 self.flag_deflip = True
705 704 else:
706 705 self.flag_deflip = False
707 706
708 707 nChannels = 0
709 708 nPairs = 0
710 709 pairList = []
711 710
712 711 for i in range(0, self.totalSpectra * 2, 2):
713 712 if self.spectraComb[i] == self.spectraComb[i + 1]:
714 713 nChannels = nChannels + 1 # par de canales iguales
715 714 else:
716 715 nPairs = nPairs + 1 # par de canales diferentes
717 716 pairList.append((self.spectraComb[i], self.spectraComb[i + 1]))
718 717
719 718 self.flag_cspc = False
720 719 if nPairs > 0:
721 720 self.flag_cspc = True
722 721
723 722 if startFp is not None:
724 723 endFp = size + startFp
725 724 if fp.tell() > endFp:
726 725 sys.stderr.write(
727 726 "Warning: Processing header size is lower than it has to be")
728 727 return 0
729 728
730 729 if fp.tell() < endFp:
731 730 sys.stderr.write(
732 731 "Warning: Processing header size is greater than it is considered")
733 732
734 733 return 1
735 734
736 735 def write(self, fp):
737 736 # Clear DEFINE_PROCESS_CODE
738 737 self.processFlags = self.processFlags & (~PROCFLAG.DEFINE_PROCESS_CODE)
739 738
740 739 headerTuple = (self.size,
741 740 self.dtype,
742 741 self.blockSize,
743 742 self.profilesPerBlock,
744 743 self.dataBlocksPerFile,
745 744 self.nWindows,
746 745 self.processFlags,
747 746 self.nCohInt,
748 747 self.nIncohInt,
749 748 self.totalSpectra)
750 749
751 750 header = numpy.array(headerTuple, PROCESSING_STRUCTURE)
752 751 header.tofile(fp)
753 752
754 753 if self.nWindows != 0:
755 754 sampleWindowTuple = (
756 755 self.firstHeight, self.deltaHeight, self.samplesWin)
757 756 samplingWindow = numpy.array(sampleWindowTuple, SAMPLING_STRUCTURE)
758 757 samplingWindow.tofile(fp)
759 758
760 759 if self.totalSpectra != 0:
761 760 # spectraComb = numpy.array([],numpy.dtype('u1'))
762 761 spectraComb = self.spectraComb
763 762 spectraComb.tofile(fp)
764 763
765 764 # if self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE == PROCFLAG.DEFINE_PROCESS_CODE:
766 765 # nCode = numpy.array([self.nCode], numpy.dtype('u4')) #Probar con un dato que almacene codigo, hasta el momento no se hizo la prueba
767 766 # nCode.tofile(fp)
768 767 #
769 768 # nBaud = numpy.array([self.nBaud], numpy.dtype('u4'))
770 769 # nBaud.tofile(fp)
771 770 #
772 771 # code = self.code.reshape(self.nCode*self.nBaud)
773 772 # code = code.astype(numpy.dtype('<f4'))
774 773 # code.tofile(fp)
775 774
776 775 return 1
777 776
778 777 def get_size(self):
779 778
780 779 self.__size = 40 + 12 * self.nWindows + 2 * self.totalSpectra
781 780
782 781 # if self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE == PROCFLAG.DEFINE_PROCESS_CODE:
783 782 # self.__size += 4 + 4 + 4*self.nCode*numpy.ceil(self.nBaud/32.)
784 783 # self.__size += 4 + 4 + 4 * self.nCode * self.nBaud
785 784
786 785 return self.__size
787 786
788 787 def set_size(self, value):
789 788
790 789 raise IOError("size is a property and it cannot be set, just read")
791 790
792 791 return
793 792
794 793 size = property(get_size, set_size)
795 794
796 795
797 796 class RCfunction:
798 797 NONE = 0
799 798 FLIP = 1
800 799 CODE = 2
801 800 SAMPLING = 3
802 801 LIN6DIV256 = 4
803 802 SYNCHRO = 5
804 803
805 804
806 805 class nCodeType:
807 806 NONE = 0
808 807 USERDEFINE = 1
809 808 BARKER2 = 2
810 809 BARKER3 = 3
811 810 BARKER4 = 4
812 811 BARKER5 = 5
813 812 BARKER7 = 6
814 813 BARKER11 = 7
815 814 BARKER13 = 8
816 815 AC128 = 9
817 816 COMPLEMENTARYCODE2 = 10
818 817 COMPLEMENTARYCODE4 = 11
819 818 COMPLEMENTARYCODE8 = 12
820 819 COMPLEMENTARYCODE16 = 13
821 820 COMPLEMENTARYCODE32 = 14
822 821 COMPLEMENTARYCODE64 = 15
823 822 COMPLEMENTARYCODE128 = 16
824 823 CODE_BINARY28 = 17
825 824
826 825
827 826 class PROCFLAG:
828 827
829 828 COHERENT_INTEGRATION = numpy.uint32(0x00000001)
830 829 DECODE_DATA = numpy.uint32(0x00000002)
831 830 SPECTRA_CALC = numpy.uint32(0x00000004)
832 831 INCOHERENT_INTEGRATION = numpy.uint32(0x00000008)
833 832 POST_COHERENT_INTEGRATION = numpy.uint32(0x00000010)
834 833 SHIFT_FFT_DATA = numpy.uint32(0x00000020)
835 834
836 835 DATATYPE_CHAR = numpy.uint32(0x00000040)
837 836 DATATYPE_SHORT = numpy.uint32(0x00000080)
838 837 DATATYPE_LONG = numpy.uint32(0x00000100)
839 838 DATATYPE_INT64 = numpy.uint32(0x00000200)
840 839 DATATYPE_FLOAT = numpy.uint32(0x00000400)
841 840 DATATYPE_DOUBLE = numpy.uint32(0x00000800)
842 841
843 842 DATAARRANGE_CONTIGUOUS_CH = numpy.uint32(0x00001000)
844 843 DATAARRANGE_CONTIGUOUS_H = numpy.uint32(0x00002000)
845 844 DATAARRANGE_CONTIGUOUS_P = numpy.uint32(0x00004000)
846 845
847 846 SAVE_CHANNELS_DC = numpy.uint32(0x00008000)
848 847 DEFLIP_DATA = numpy.uint32(0x00010000)
849 848 DEFINE_PROCESS_CODE = numpy.uint32(0x00020000)
850 849
851 850 ACQ_SYS_NATALIA = numpy.uint32(0x00040000)
852 851 ACQ_SYS_ECHOTEK = numpy.uint32(0x00080000)
853 852 ACQ_SYS_ADRXD = numpy.uint32(0x000C0000)
854 853 ACQ_SYS_JULIA = numpy.uint32(0x00100000)
855 854 ACQ_SYS_XXXXXX = numpy.uint32(0x00140000)
856 855
857 856 EXP_NAME_ESP = numpy.uint32(0x00200000)
858 857 CHANNEL_NAMES_ESP = numpy.uint32(0x00400000)
859 858
860 859 OPERATION_MASK = numpy.uint32(0x0000003F)
861 860 DATATYPE_MASK = numpy.uint32(0x00000FC0)
862 861 DATAARRANGE_MASK = numpy.uint32(0x00007000)
863 862 ACQ_SYS_MASK = numpy.uint32(0x001C0000)
864 863
865 864
866 865 dtype0 = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
867 866 dtype1 = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
868 867 dtype2 = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
869 868 dtype3 = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
870 869 dtype4 = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
871 870 dtype5 = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
872 871
873 872 NUMPY_DTYPE_LIST = [dtype0, dtype1, dtype2, dtype3, dtype4, dtype5]
874 873
875 874 PROCFLAG_DTYPE_LIST = [PROCFLAG.DATATYPE_CHAR,
876 875 PROCFLAG.DATATYPE_SHORT,
877 876 PROCFLAG.DATATYPE_LONG,
878 877 PROCFLAG.DATATYPE_INT64,
879 878 PROCFLAG.DATATYPE_FLOAT,
880 879 PROCFLAG.DATATYPE_DOUBLE]
881 880
882 881 DTYPE_WIDTH = [1, 2, 4, 8, 4, 8]
883 882
884 883
885 884 def get_dtype_index(numpy_dtype):
886 885
887 886 index = None
888 887
889 888 for i in range(len(NUMPY_DTYPE_LIST)):
890 889 if numpy_dtype == NUMPY_DTYPE_LIST[i]:
891 890 index = i
892 891 break
893 892
894 893 return index
895 894
896 895
897 896 def get_numpy_dtype(index):
898 897
899 898 return NUMPY_DTYPE_LIST[index]
900 899
901 900
902 901 def get_procflag_dtype(index):
903 902
904 903 return PROCFLAG_DTYPE_LIST[index]
905 904
906 905
907 906 def get_dtype_width(index):
908 907
909 908 return DTYPE_WIDTH[index]
Show file before | Show file after | Hide comments
@@ -1,355 +1,361
1 1 '''
2 2 Created on Nov 9, 2016
3 3
4 4 @author: roj- LouVD
5 5 '''
6 6
7 7
8 8 import os
9 9 import sys
10 10 import time
11 11 import glob
12 12 import datetime
13 13
14 14 import numpy
15 15
16 16 import schainpy.admin
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator
18 18 from schainpy.model.data.jrodata import Parameters
19 19 from schainpy.model.io.jroIO_base import Reader
20 20 from schainpy.utils import log
21 21
22 22 FILE_HEADER_STRUCTURE = numpy.dtype([
23 23 ('FMN', '<u4'),
24 24 ('nrec', '<u4'),
25 25 ('fr_offset', '<u4'),
26 26 ('id', '<u4'),
27 27 ('site', 'u1', (32,))
28 28 ])
29 29
30 30 REC_HEADER_STRUCTURE = numpy.dtype([
31 31 ('rmn', '<u4'),
32 32 ('rcounter', '<u4'),
33 33 ('nr_offset', '<u4'),
34 34 ('tr_offset', '<u4'),
35 35 ('time', '<u4'),
36 36 ('time_msec', '<u4'),
37 37 ('tag', 'u1', (32,)),
38 38 ('comments', 'u1', (32,)),
39 39 ('lat', '<f4'),
40 40 ('lon', '<f4'),
41 41 ('gps_status', '<u4'),
42 42 ('freq', '<u4'),
43 43 ('freq0', '<u4'),
44 44 ('nchan', '<u4'),
45 45 ('delta_r', '<u4'),
46 46 ('nranges', '<u4'),
47 47 ('r0', '<u4'),
48 48 ('prf', '<u4'),
49 49 ('ncoh', '<u4'),
50 50 ('npoints', '<u4'),
51 51 ('polarization', '<i4'),
52 52 ('rx_filter', '<u4'),
53 53 ('nmodes', '<u4'),
54 54 ('dmode_index', '<u4'),
55 55 ('dmode_rngcorr', '<u4'),
56 56 ('nrxs', '<u4'),
57 57 ('acf_length', '<u4'),
58 58 ('acf_lags', '<u4'),
59 59 ('sea_to_atmos', '<f4'),
60 60 ('sea_notch', '<u4'),
61 61 ('lh_sea', '<u4'),
62 62 ('hh_sea', '<u4'),
63 63 ('nbins_sea', '<u4'),
64 64 ('min_snr', '<f4'),
65 65 ('min_cc', '<f4'),
66 66 ('max_time_diff', '<f4')
67 67 ])
68 68
69 69 DATA_STRUCTURE = numpy.dtype([
70 70 ('range', '<u4'),
71 71 ('status', '<u4'),
72 72 ('zonal', '<f4'),
73 73 ('meridional', '<f4'),
74 74 ('vertical', '<f4'),
75 75 ('zonal_a', '<f4'),
76 76 ('meridional_a', '<f4'),
77 77 ('corrected_fading', '<f4'), # seconds
78 78 ('uncorrected_fading', '<f4'), # seconds
79 79 ('time_diff', '<f4'),
80 80 ('major_axis', '<f4'),
81 81 ('axial_ratio', '<f4'),
82 82 ('orientation', '<f4'),
83 83 ('sea_power', '<u4'),
84 84 ('sea_algorithm', '<u4')
85 85 ])
86 86
87 87
88 88 class BLTRParamReader(Reader, ProcessingUnit):
89 89 '''
90 90 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR
91 91 from *.sswma files
92 92 '''
93 93
94 94 ext = '.sswma'
95 95
96 96 def __init__(self):
97 97
98 98 ProcessingUnit.__init__(self)
99 99
100 100 self.dataOut = Parameters()
101 101 self.dataOut.timezone = 300
102 102 self.counter_records = 0
103 103 self.flagNoMoreFiles = 0
104 104 self.isConfig = False
105 105 self.filename = None
106 106 self.status_value = 0
107 107 self.datatime = datetime.datetime(1900, 1, 1)
108 108 self.filefmt = "*********%Y%m%d******"
109 109
110 110 def setup(self, **kwargs):
111 111
112 112 self.set_kwargs(**kwargs)
113 113
114 114 if self.path is None:
115 115 raise ValueError("The path is not valid")
116 116
117 117 if self.online:
118 118 log.log("Searching files in online mode...", self.name)
119 119
120 120 for nTries in range(self.nTries):
121 121 fullpath = self.searchFilesOnLine(self.path, self.startDate,
122 122 self.endDate, self.expLabel, self.ext, self.walk,
123 123 self.filefmt, self.folderfmt)
124 124 try:
125 125 fullpath = next(fullpath)
126 126 except:
127 127 fullpath = None
128 128
129 129 if fullpath:
130 130 self.fileSize = os.path.getsize(fullpath)
131 131 self.filename = fullpath
132 132 self.flagIsNewFile = 1
133 133 if self.fp != None:
134 134 self.fp.close()
135 135 self.fp = self.open_file(fullpath, self.open_mode)
136 136 self.flagNoMoreFiles = 0
137 137 break
138 138
139 139 log.warning(
140 140 'Waiting {} sec for a valid file in {}: try {} ...'.format(
141 141 self.delay, self.path, nTries + 1),
142 142 self.name)
143 143 time.sleep(self.delay)
144 144
145 145 if not(fullpath):
146 146 raise schainpy.admin.SchainError(
147 147 'There isn\'t any valid file in {}'.format(self.path))
148 148 self.readFirstHeader()
149 149 else:
150 150 log.log("Searching files in {}".format(self.path), self.name)
151 151 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
152 152 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
153 153 self.setNextFile()
154 154
155 155 def checkForRealPath(self, nextFile, nextDay):
156 156 '''
157 157 '''
158 158
159 159 dt = self.datatime + datetime.timedelta(1)
160 160 filename = '{}.{}{}'.format(self.siteFile, dt.strftime('%Y%m%d'), self.ext)
161 161 fullfilename = os.path.join(self.path, filename)
162 162 if os.path.exists(fullfilename):
163 163 return fullfilename, filename
164 164 return None, filename
165 165
166 166
167 167 def readFirstHeader(self):
168 168 '''
169 169 '''
170 170
171 171 # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya
172 172 self.siteFile = self.filename.split('/')[-1].split('.')[0]
173 173 self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
174 174 self.nrecords = self.header_file['nrec'][0]
175 175 self.counter_records = 0
176 176 self.flagIsNewFile = 0
177 self.fileIndex += 1
177 self.fileIndex += 1
178 178
179 179 def readNextBlock(self):
180 180
181 181 while True:
182 182 if not self.online and self.counter_records == self.nrecords:
183 183 self.flagIsNewFile = 1
184 184 if not self.setNextFile():
185 185 return 0
186 186 try:
187 pointer = self.fp.tell()
187 if self.online and self.counter_records == 0:
188 pos = int(self.fileSize / (38512))
189 self.counter_records = pos*2 - 2
190 pointer = 38512 * (pos-1) + 48
191 self.fp.seek(pointer)
192 else:
193 pointer = self.fp.tell()
188 194 self.readBlock()
189 195 except:
190 196 if self.online and self.waitDataBlock(pointer, 38512) == 1:
191 197 continue
192 198 else:
193 199 if not self.setNextFile():
194 200 return 0
195 201
196 202 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
197 203 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
198 204 log.warning(
199 205 'Reading Record No. {}/{} -> {} [Skipping]'.format(
200 206 self.counter_records,
201 207 self.nrecords,
202 208 self.datatime.ctime()),
203 209 'BLTRParamReader')
204 210 continue
205 211 break
206 212
207 213 log.log('Reading Record No. {} -> {}'.format(
208 214 self.counter_records,
209 215 self.datatime.ctime()), 'BLTRParamReader')
210 216
211 217 return 1
212 218
213 219 def readBlock(self):
214 220
215 221 pointer = self.fp.tell()
216 222 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
217 223 self.nchannels = int(header_rec['nchan'][0] / 2)
218 224 self.kchan = header_rec['nrxs'][0]
219 225 self.nmodes = header_rec['nmodes'][0]
220 226 self.nranges = header_rec['nranges'][0]
221 227 self.fp.seek(pointer)
222 228 self.height = numpy.empty((self.nmodes, self.nranges))
223 229 self.snr = numpy.empty((self.nmodes, int(self.nchannels), self.nranges))
224 230 self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
225 231 self.flagDiscontinuousBlock = 0
226 232
227 233 for mode in range(self.nmodes):
228 234 self.readHeader()
229 235 data = self.readData()
230 236 self.height[mode] = (data[0] - self.correction) / 1000.
231 237 self.buffer[mode] = data[1]
232 238 self.snr[mode] = data[2]
233 239
234 240 self.counter_records = self.counter_records + self.nmodes
235 241
236 242 return
237 243
238 244 def readHeader(self):
239 245 '''
240 246 RecordHeader of BLTR rawdata file
241 247 '''
242 248
243 249 header_structure = numpy.dtype(
244 250 REC_HEADER_STRUCTURE.descr + [
245 251 ('antenna_coord', 'f4', (2, int(self.nchannels))),
246 252 ('rx_gains', 'u4', (int(self.nchannels),)),
247 253 ('rx_analysis', 'u4', (int(self.nchannels),))
248 254 ]
249 255 )
250 256
251 257 self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
252 258 self.lat = self.header_rec['lat'][0]
253 259 self.lon = self.header_rec['lon'][0]
254 260 self.delta = self.header_rec['delta_r'][0]
255 261 self.correction = self.header_rec['dmode_rngcorr'][0]
256 262 self.imode = self.header_rec['dmode_index'][0]
257 263 self.antenna = self.header_rec['antenna_coord']
258 264 self.rx_gains = self.header_rec['rx_gains']
259 265 self.time = self.header_rec['time'][0]
260 266 dt = datetime.datetime.utcfromtimestamp(self.time)
261 267 if dt.date() > self.datatime.date():
262 268 self.flagDiscontinuousBlock = 1
263 269 self.datatime = dt
264 270
265 271 def readData(self):
266 272 '''
267 273 Reading and filtering data block record of BLTR rawdata file,
268 274 filtering is according to status_value.
269 275
270 276 Input:
271 277 status_value - Array data is set to NAN for values that are not
272 278 equal to status_value
273 279
274 280 '''
275 281 self.nchannels = int(self.nchannels)
276 282
277 283 data_structure = numpy.dtype(
278 284 DATA_STRUCTURE.descr + [
279 285 ('rx_saturation', 'u4', (self.nchannels,)),
280 286 ('chan_offset', 'u4', (2 * self.nchannels,)),
281 287 ('rx_amp', 'u4', (self.nchannels,)),
282 288 ('rx_snr', 'f4', (self.nchannels,)),
283 289 ('cross_snr', 'f4', (self.kchan,)),
284 290 ('sea_power_relative', 'f4', (self.kchan,))]
285 291 )
286 292
287 293 data = numpy.fromfile(self.fp, data_structure, self.nranges)
288 294
289 295 height = data['range']
290 296 winds = numpy.array(
291 297 (data['zonal'], data['meridional'], data['vertical']))
292 298 snr = data['rx_snr'].T
293 299
294 300 winds[numpy.where(winds == -9999.)] = numpy.nan
295 301 winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
296 302 snr[numpy.where(snr == -9999.)] = numpy.nan
297 303 snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
298 304 snr = numpy.power(10, snr / 10)
299 305
300 306 return height, winds, snr
301 307
302 308 def set_output(self):
303 309 '''
304 310 Storing data from databuffer to dataOut object
305 311 '''
306 312
307 313 self.dataOut.data_snr = self.snr
308 314 self.dataOut.height = self.height
309 315 self.dataOut.data = self.buffer
310 316 self.dataOut.utctimeInit = self.time
311 317 self.dataOut.utctime = self.dataOut.utctimeInit
312 318 self.dataOut.useLocalTime = False
313 319 self.dataOut.paramInterval = 157
314 320 self.dataOut.site = self.siteFile
315 321 self.dataOut.nrecords = self.nrecords / self.nmodes
316 322 self.dataOut.lat = self.lat
317 323 self.dataOut.lon = self.lon
318 324 self.dataOut.channelList = list(range(self.nchannels))
319 325 self.dataOut.kchan = self.kchan
320 326 self.dataOut.delta = self.delta
321 327 self.dataOut.correction = self.correction
322 328 self.dataOut.nmodes = self.nmodes
323 329 self.dataOut.imode = self.imode
324 330 self.dataOut.antenna = self.antenna
325 331 self.dataOut.rx_gains = self.rx_gains
326 332 self.dataOut.flagNoData = False
327 333 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
328 334
329 335 def getData(self):
330 336 '''
331 337 Storing data from databuffer to dataOut object
332 338 '''
333 339 if self.flagNoMoreFiles:
334 340 self.dataOut.flagNoData = True
335 341 return 0
336 342
337 343 if not self.readNextBlock():
338 344 self.dataOut.flagNoData = True
339 345 return 0
340 346
341 347 self.set_output()
342 348
343 349 return 1
344 350
345 351 def run(self, **kwargs):
346 352 '''
347 353 '''
348 354
349 355 if not(self.isConfig):
350 356 self.setup(**kwargs)
351 357 self.isConfig = True
352 358
353 359 self.getData()
354 360
355 361 return
General Comments 0
You need to be logged in to leave comments. Login now