##// END OF EJS Templates
schain-jc
RSS
jrodata: TimeRange added for ACF module
Miguel Valdez -
r771:8c74dbbfdff0
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,1158 +1,1165
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
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12
13 13 def getNumpyDtype(dataTypeCode):
14 14
15 15 if dataTypeCode == 0:
16 16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 17 elif dataTypeCode == 1:
18 18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 19 elif dataTypeCode == 2:
20 20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 21 elif dataTypeCode == 3:
22 22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 23 elif dataTypeCode == 4:
24 24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 25 elif dataTypeCode == 5:
26 26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 27 else:
28 28 raise ValueError, 'dataTypeCode was not defined'
29 29
30 30 return numpyDtype
31 31
32 32 def getDataTypeCode(numpyDtype):
33 33
34 34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 35 datatype = 0
36 36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 37 datatype = 1
38 38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 39 datatype = 2
40 40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 41 datatype = 3
42 42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 43 datatype = 4
44 44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 45 datatype = 5
46 46 else:
47 47 datatype = None
48 48
49 49 return datatype
50 50
51 51 def hildebrand_sekhon(data, navg):
52 52 """
53 53 This method is for the objective determination of the noise level in Doppler spectra. This
54 54 implementation technique is based on the fact that the standard deviation of the spectral
55 55 densities is equal to the mean spectral density for white Gaussian noise
56 56
57 57 Inputs:
58 58 Data : heights
59 59 navg : numbers of averages
60 60
61 61 Return:
62 62 -1 : any error
63 63 anoise : noise's level
64 64 """
65 65
66 66 sortdata = numpy.sort(data,axis=None)
67 67 lenOfData = len(sortdata)
68 68 nums_min = lenOfData*0.2
69 69
70 70 if nums_min <= 5:
71 71 nums_min = 5
72 72
73 73 sump = 0.
74 74
75 75 sumq = 0.
76 76
77 77 j = 0
78 78
79 79 cont = 1
80 80
81 81 while((cont==1)and(j<lenOfData)):
82 82
83 83 sump += sortdata[j]
84 84
85 85 sumq += sortdata[j]**2
86 86
87 87 if j > nums_min:
88 88 rtest = float(j)/(j-1) + 1.0/navg
89 89 if ((sumq*j) > (rtest*sump**2)):
90 90 j = j - 1
91 91 sump = sump - sortdata[j]
92 92 sumq = sumq - sortdata[j]**2
93 93 cont = 0
94 94
95 95 j += 1
96 96
97 97 lnoise = sump /j
98 98 # stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
99 99 return lnoise
100 100
101 101 class Beam:
102 102 def __init__(self):
103 103 self.codeList = []
104 104 self.azimuthList = []
105 105 self.zenithList = []
106 106
107 107 class GenericData(object):
108 108
109 109 flagNoData = True
110 110
111 111 def __init__(self):
112 112
113 113 raise NotImplementedError
114 114
115 115 def copy(self, inputObj=None):
116 116
117 117 if inputObj == None:
118 118 return copy.deepcopy(self)
119 119
120 120 for key in inputObj.__dict__.keys():
121 121
122 122 attribute = inputObj.__dict__[key]
123 123
124 124 #If this attribute is a tuple or list
125 125 if type(inputObj.__dict__[key]) in (tuple, list):
126 126 self.__dict__[key] = attribute[:]
127 127 continue
128 128
129 129 #If this attribute is another object or instance
130 130 if hasattr(attribute, '__dict__'):
131 131 self.__dict__[key] = attribute.copy()
132 132 continue
133 133
134 134 self.__dict__[key] = inputObj.__dict__[key]
135 135
136 136 def deepcopy(self):
137 137
138 138 return copy.deepcopy(self)
139 139
140 140 def isEmpty(self):
141 141
142 142 return self.flagNoData
143 143
144 144 class JROData(GenericData):
145 145
146 146 # m_BasicHeader = BasicHeader()
147 147 # m_ProcessingHeader = ProcessingHeader()
148 148
149 149 systemHeaderObj = SystemHeader()
150 150
151 151 radarControllerHeaderObj = RadarControllerHeader()
152 152
153 153 # data = None
154 154
155 155 type = None
156 156
157 157 datatype = None #dtype but in string
158 158
159 159 # dtype = None
160 160
161 161 # nChannels = None
162 162
163 163 # nHeights = None
164 164
165 165 nProfiles = None
166 166
167 167 heightList = None
168 168
169 169 channelList = None
170 170
171 171 flagDiscontinuousBlock = False
172 172
173 173 useLocalTime = False
174 174
175 175 utctime = None
176 176
177 177 timeZone = None
178 178
179 179 dstFlag = None
180 180
181 181 errorCount = None
182 182
183 183 blocksize = None
184 184
185 185 # nCode = None
186 186 #
187 187 # nBaud = None
188 188 #
189 189 # code = None
190 190
191 191 flagDecodeData = False #asumo q la data no esta decodificada
192 192
193 193 flagDeflipData = False #asumo q la data no esta sin flip
194 194
195 195 flagShiftFFT = False
196 196
197 197 # ippSeconds = None
198 198
199 199 # timeInterval = None
200 200
201 201 nCohInt = None
202 202
203 203 # noise = None
204 204
205 205 windowOfFilter = 1
206 206
207 207 #Speed of ligth
208 208 C = 3e8
209 209
210 210 frequency = 49.92e6
211 211
212 212 realtime = False
213 213
214 214 beacon_heiIndexList = None
215 215
216 216 last_block = None
217 217
218 218 blocknow = None
219 219
220 220 azimuth = None
221 221
222 222 zenith = None
223 223
224 224 beam = Beam()
225 225
226 226 profileIndex = None
227 227
228 228 def __init__(self):
229 229
230 230 raise NotImplementedError
231 231
232 232 def getNoise(self):
233 233
234 234 raise NotImplementedError
235 235
236 236 def getNChannels(self):
237 237
238 238 return len(self.channelList)
239 239
240 240 def getChannelIndexList(self):
241 241
242 242 return range(self.nChannels)
243 243
244 244 def getNHeights(self):
245 245
246 246 return len(self.heightList)
247 247
248 248 def getHeiRange(self, extrapoints=0):
249 249
250 250 heis = self.heightList
251 251 # deltah = self.heightList[1] - self.heightList[0]
252 252 #
253 253 # heis.append(self.heightList[-1])
254 254
255 255 return heis
256 256
257 257 def getDeltaH(self):
258 258
259 259 delta = self.heightList[1] - self.heightList[0]
260 260
261 261 return delta
262 262
263 263 def getltctime(self):
264 264
265 265 if self.useLocalTime:
266 266 return self.utctime - self.timeZone*60
267 267
268 268 return self.utctime
269 269
270 270 def getDatatime(self):
271 271
272 272 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
273 273 return datatimeValue
274 274
275 275 def getTimeRange(self):
276 276
277 277 datatime = []
278 278
279 279 datatime.append(self.ltctime)
280 280 datatime.append(self.ltctime + self.timeInterval+60)
281 281
282 282 datatime = numpy.array(datatime)
283 283
284 284 return datatime
285 285
286 286 def getFmaxTimeResponse(self):
287 287
288 288 period = (10**-6)*self.getDeltaH()/(0.15)
289 289
290 290 PRF = 1./(period * self.nCohInt)
291 291
292 292 fmax = PRF
293 293
294 294 return fmax
295 295
296 296 def getFmax(self):
297 297
298 298 PRF = 1./(self.ippSeconds * self.nCohInt)
299 299
300 300 fmax = PRF
301 301
302 302 return fmax
303 303
304 304 def getVmax(self):
305 305
306 306 _lambda = self.C/self.frequency
307 307
308 308 vmax = self.getFmax() * _lambda
309 309
310 310 return vmax
311 311
312 312 def get_ippSeconds(self):
313 313 '''
314 314 '''
315 315 return self.radarControllerHeaderObj.ippSeconds
316 316
317 317 def set_ippSeconds(self, ippSeconds):
318 318 '''
319 319 '''
320 320
321 321 self.radarControllerHeaderObj.ippSeconds = ippSeconds
322 322
323 323 return
324 324
325 325 def get_dtype(self):
326 326 '''
327 327 '''
328 328 return getNumpyDtype(self.datatype)
329 329
330 330 def set_dtype(self, numpyDtype):
331 331 '''
332 332 '''
333 333
334 334 self.datatype = getDataTypeCode(numpyDtype)
335 335
336 336 def get_code(self):
337 337 '''
338 338 '''
339 339 return self.radarControllerHeaderObj.code
340 340
341 341 def set_code(self, code):
342 342 '''
343 343 '''
344 344 self.radarControllerHeaderObj.code = code
345 345
346 346 return
347 347
348 348 def get_ncode(self):
349 349 '''
350 350 '''
351 351 return self.radarControllerHeaderObj.nCode
352 352
353 353 def set_ncode(self, nCode):
354 354 '''
355 355 '''
356 356 self.radarControllerHeaderObj.nCode = nCode
357 357
358 358 return
359 359
360 360 def get_nbaud(self):
361 361 '''
362 362 '''
363 363 return self.radarControllerHeaderObj.nBaud
364 364
365 365 def set_nbaud(self, nBaud):
366 366 '''
367 367 '''
368 368 self.radarControllerHeaderObj.nBaud = nBaud
369 369
370 370 return
371 371
372 372 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
373 373 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
374 374 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
375 375 #noise = property(getNoise, "I'm the 'nHeights' property.")
376 376 datatime = property(getDatatime, "I'm the 'datatime' property")
377 377 ltctime = property(getltctime, "I'm the 'ltctime' property")
378 378 ippSeconds = property(get_ippSeconds, set_ippSeconds)
379 379 dtype = property(get_dtype, set_dtype)
380 380 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
381 381 code = property(get_code, set_code)
382 382 nCode = property(get_ncode, set_ncode)
383 383 nBaud = property(get_nbaud, set_nbaud)
384 384
385 385 class Voltage(JROData):
386 386
387 387 #data es un numpy array de 2 dmensiones (canales, alturas)
388 388 data = None
389 389
390 390 def __init__(self):
391 391 '''
392 392 Constructor
393 393 '''
394 394
395 395 self.useLocalTime = True
396 396
397 397 self.radarControllerHeaderObj = RadarControllerHeader()
398 398
399 399 self.systemHeaderObj = SystemHeader()
400 400
401 401 self.type = "Voltage"
402 402
403 403 self.data = None
404 404
405 405 # self.dtype = None
406 406
407 407 # self.nChannels = 0
408 408
409 409 # self.nHeights = 0
410 410
411 411 self.nProfiles = None
412 412
413 413 self.heightList = None
414 414
415 415 self.channelList = None
416 416
417 417 # self.channelIndexList = None
418 418
419 419 self.flagNoData = True
420 420
421 421 self.flagDiscontinuousBlock = False
422 422
423 423 self.utctime = None
424 424
425 425 self.timeZone = None
426 426
427 427 self.dstFlag = None
428 428
429 429 self.errorCount = None
430 430
431 431 self.nCohInt = None
432 432
433 433 self.blocksize = None
434 434
435 435 self.flagDecodeData = False #asumo q la data no esta decodificada
436 436
437 437 self.flagDeflipData = False #asumo q la data no esta sin flip
438 438
439 439 self.flagShiftFFT = False
440 440
441 441 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
442 442
443 443 self.profileIndex = 0
444 444
445 445 def getNoisebyHildebrand(self, channel = None):
446 446 """
447 447 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
448 448
449 449 Return:
450 450 noiselevel
451 451 """
452 452
453 453 if channel != None:
454 454 data = self.data[channel]
455 455 nChannels = 1
456 456 else:
457 457 data = self.data
458 458 nChannels = self.nChannels
459 459
460 460 noise = numpy.zeros(nChannels)
461 461 power = data * numpy.conjugate(data)
462 462
463 463 for thisChannel in range(nChannels):
464 464 if nChannels == 1:
465 465 daux = power[:].real
466 466 else:
467 467 daux = power[thisChannel,:].real
468 468 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
469 469
470 470 return noise
471 471
472 472 def getNoise(self, type = 1, channel = None):
473 473
474 474 if type == 1:
475 475 noise = self.getNoisebyHildebrand(channel)
476 476
477 477 return noise
478 478
479 479 def getPower(self, channel = None):
480 480
481 481 if channel != None:
482 482 data = self.data[channel]
483 483 else:
484 484 data = self.data
485 485
486 486 power = data * numpy.conjugate(data)
487 487
488 488 return 10*numpy.log10(power.real)
489 489
490 490 def getTimeInterval(self):
491 491
492 492 timeInterval = self.ippSeconds * self.nCohInt
493 493
494 494 return timeInterval
495 495
496 496 noise = property(getNoise, "I'm the 'nHeights' property.")
497 497 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
498 498
499 499 class Spectra(JROData):
500 500
501 501 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
502 502 data_spc = None
503 503
504 504 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
505 505 data_cspc = None
506 506
507 507 #data es un numpy array de 2 dmensiones (canales, alturas)
508 508 data_dc = None
509 509
510 510 nFFTPoints = None
511 511
512 512 # nPairs = None
513 513
514 514 pairsList = None
515 515
516 516 nIncohInt = None
517 517
518 518 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
519 519
520 520 nCohInt = None #se requiere para determinar el valor de timeInterval
521 521
522 522 ippFactor = None
523 523
524 524 profileIndex = 0
525 525
526 526 def __init__(self):
527 527 '''
528 528 Constructor
529 529 '''
530 530
531 531 self.useLocalTime = True
532 532
533 533 self.radarControllerHeaderObj = RadarControllerHeader()
534 534
535 535 self.systemHeaderObj = SystemHeader()
536 536
537 537 self.type = "Spectra"
538 538
539 539 # self.data = None
540 540
541 541 # self.dtype = None
542 542
543 543 # self.nChannels = 0
544 544
545 545 # self.nHeights = 0
546 546
547 547 self.nProfiles = None
548 548
549 549 self.heightList = None
550 550
551 551 self.channelList = None
552 552
553 553 # self.channelIndexList = None
554 554
555 555 self.pairsList = None
556 556
557 557 self.flagNoData = True
558 558
559 559 self.flagDiscontinuousBlock = False
560 560
561 561 self.utctime = None
562 562
563 563 self.nCohInt = None
564 564
565 565 self.nIncohInt = None
566 566
567 567 self.blocksize = None
568 568
569 569 self.nFFTPoints = None
570 570
571 571 self.wavelength = None
572 572
573 573 self.flagDecodeData = False #asumo q la data no esta decodificada
574 574
575 575 self.flagDeflipData = False #asumo q la data no esta sin flip
576 576
577 577 self.flagShiftFFT = False
578 578
579 579 self.ippFactor = 1
580 580
581 581 #self.noise = None
582 582
583 583 self.beacon_heiIndexList = []
584 584
585 585 self.noise_estimation = None
586 586
587 587
588 588 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
589 589 """
590 590 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
591 591
592 592 Return:
593 593 noiselevel
594 594 """
595 595
596 596 noise = numpy.zeros(self.nChannels)
597 597
598 598 for channel in range(self.nChannels):
599 599 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
600 600 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
601 601
602 602 return noise
603 603
604 604 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
605 605
606 606 if self.noise_estimation is not None:
607 607 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
608 608 else:
609 609 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
610 610 return noise
611 611
612 612 def getFreqRangeTimeResponse(self, extrapoints=0):
613 613
614 614 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
615 615 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
616 616
617 617 return freqrange
618 618
619 def getAcfRange(self, extrapoints=0):
620
621 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
622 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
623
624 return freqrange
625
619 626 def getFreqRange(self, extrapoints=0):
620 627
621 628 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
622 629 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
623 630
624 631 return freqrange
625 632
626 633 def getVelRange(self, extrapoints=0):
627 634
628 635 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
629 636 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
630 637
631 638 return velrange
632 639
633 640 def getNPairs(self):
634 641
635 642 return len(self.pairsList)
636 643
637 644 def getPairsIndexList(self):
638 645
639 646 return range(self.nPairs)
640 647
641 648 def getNormFactor(self):
642 649
643 650 pwcode = 1
644 651
645 652 if self.flagDecodeData:
646 653 pwcode = numpy.sum(self.code[0]**2)
647 654 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
648 655 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
649 656
650 657 return normFactor
651 658
652 659 def getFlagCspc(self):
653 660
654 661 if self.data_cspc is None:
655 662 return True
656 663
657 664 return False
658 665
659 666 def getFlagDc(self):
660 667
661 668 if self.data_dc is None:
662 669 return True
663 670
664 671 return False
665 672
666 673 def getTimeInterval(self):
667 674
668 675 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
669 676
670 677 return timeInterval
671 678
672 679 def setValue(self, value):
673 680
674 681 print "This property should not be initialized"
675 682
676 683 return
677 684
678 685 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
679 686 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
680 687 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
681 688 flag_cspc = property(getFlagCspc, setValue)
682 689 flag_dc = property(getFlagDc, setValue)
683 690 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
684 691 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
685 692
686 693 class SpectraHeis(Spectra):
687 694
688 695 data_spc = None
689 696
690 697 data_cspc = None
691 698
692 699 data_dc = None
693 700
694 701 nFFTPoints = None
695 702
696 703 # nPairs = None
697 704
698 705 pairsList = None
699 706
700 707 nCohInt = None
701 708
702 709 nIncohInt = None
703 710
704 711 def __init__(self):
705 712
706 713 self.radarControllerHeaderObj = RadarControllerHeader()
707 714
708 715 self.systemHeaderObj = SystemHeader()
709 716
710 717 self.type = "SpectraHeis"
711 718
712 719 # self.dtype = None
713 720
714 721 # self.nChannels = 0
715 722
716 723 # self.nHeights = 0
717 724
718 725 self.nProfiles = None
719 726
720 727 self.heightList = None
721 728
722 729 self.channelList = None
723 730
724 731 # self.channelIndexList = None
725 732
726 733 self.flagNoData = True
727 734
728 735 self.flagDiscontinuousBlock = False
729 736
730 737 # self.nPairs = 0
731 738
732 739 self.utctime = None
733 740
734 741 self.blocksize = None
735 742
736 743 self.profileIndex = 0
737 744
738 745 self.nCohInt = 1
739 746
740 747 self.nIncohInt = 1
741 748
742 749 def getNormFactor(self):
743 750 pwcode = 1
744 751 if self.flagDecodeData:
745 752 pwcode = numpy.sum(self.code[0]**2)
746 753
747 754 normFactor = self.nIncohInt*self.nCohInt*pwcode
748 755
749 756 return normFactor
750 757
751 758 def getTimeInterval(self):
752 759
753 760 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
754 761
755 762 return timeInterval
756 763
757 764 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
758 765 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
759 766
760 767 class Fits(JROData):
761 768
762 769 heightList = None
763 770
764 771 channelList = None
765 772
766 773 flagNoData = True
767 774
768 775 flagDiscontinuousBlock = False
769 776
770 777 useLocalTime = False
771 778
772 779 utctime = None
773 780
774 781 timeZone = None
775 782
776 783 # ippSeconds = None
777 784
778 785 # timeInterval = None
779 786
780 787 nCohInt = None
781 788
782 789 nIncohInt = None
783 790
784 791 noise = None
785 792
786 793 windowOfFilter = 1
787 794
788 795 #Speed of ligth
789 796 C = 3e8
790 797
791 798 frequency = 49.92e6
792 799
793 800 realtime = False
794 801
795 802
796 803 def __init__(self):
797 804
798 805 self.type = "Fits"
799 806
800 807 self.nProfiles = None
801 808
802 809 self.heightList = None
803 810
804 811 self.channelList = None
805 812
806 813 # self.channelIndexList = None
807 814
808 815 self.flagNoData = True
809 816
810 817 self.utctime = None
811 818
812 819 self.nCohInt = 1
813 820
814 821 self.nIncohInt = 1
815 822
816 823 self.useLocalTime = True
817 824
818 825 self.profileIndex = 0
819 826
820 827 # self.utctime = None
821 828 # self.timeZone = None
822 829 # self.ltctime = None
823 830 # self.timeInterval = None
824 831 # self.header = None
825 832 # self.data_header = None
826 833 # self.data = None
827 834 # self.datatime = None
828 835 # self.flagNoData = False
829 836 # self.expName = ''
830 837 # self.nChannels = None
831 838 # self.nSamples = None
832 839 # self.dataBlocksPerFile = None
833 840 # self.comments = ''
834 841 #
835 842
836 843
837 844 def getltctime(self):
838 845
839 846 if self.useLocalTime:
840 847 return self.utctime - self.timeZone*60
841 848
842 849 return self.utctime
843 850
844 851 def getDatatime(self):
845 852
846 853 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
847 854 return datatime
848 855
849 856 def getTimeRange(self):
850 857
851 858 datatime = []
852 859
853 860 datatime.append(self.ltctime)
854 861 datatime.append(self.ltctime + self.timeInterval)
855 862
856 863 datatime = numpy.array(datatime)
857 864
858 865 return datatime
859 866
860 867 def getHeiRange(self):
861 868
862 869 heis = self.heightList
863 870
864 871 return heis
865 872
866 873 def getNHeights(self):
867 874
868 875 return len(self.heightList)
869 876
870 877 def getNChannels(self):
871 878
872 879 return len(self.channelList)
873 880
874 881 def getChannelIndexList(self):
875 882
876 883 return range(self.nChannels)
877 884
878 885 def getNoise(self, type = 1):
879 886
880 887 #noise = numpy.zeros(self.nChannels)
881 888
882 889 if type == 1:
883 890 noise = self.getNoisebyHildebrand()
884 891
885 892 if type == 2:
886 893 noise = self.getNoisebySort()
887 894
888 895 if type == 3:
889 896 noise = self.getNoisebyWindow()
890 897
891 898 return noise
892 899
893 900 def getTimeInterval(self):
894 901
895 902 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
896 903
897 904 return timeInterval
898 905
899 906 datatime = property(getDatatime, "I'm the 'datatime' property")
900 907 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
901 908 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
902 909 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
903 910 noise = property(getNoise, "I'm the 'nHeights' property.")
904 911
905 912 ltctime = property(getltctime, "I'm the 'ltctime' property")
906 913 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
907 914
908 915 class Correlation(JROData):
909 916
910 917 noise = None
911 918
912 919 SNR = None
913 920
914 921 pairsAutoCorr = None #Pairs of Autocorrelation
915 922
916 923 #--------------------------------------------------
917 924
918 925 data_corr = None
919 926
920 927 data_volt = None
921 928
922 929 lagT = None # each element value is a profileIndex
923 930
924 931 lagR = None # each element value is in km
925 932
926 933 pairsList = None
927 934
928 935 calculateVelocity = None
929 936
930 937 nPoints = None
931 938
932 939 nAvg = None
933 940
934 941 bufferSize = None
935 942
936 943 def __init__(self):
937 944 '''
938 945 Constructor
939 946 '''
940 947 self.radarControllerHeaderObj = RadarControllerHeader()
941 948
942 949 self.systemHeaderObj = SystemHeader()
943 950
944 951 self.type = "Correlation"
945 952
946 953 self.data = None
947 954
948 955 self.dtype = None
949 956
950 957 self.nProfiles = None
951 958
952 959 self.heightList = None
953 960
954 961 self.channelList = None
955 962
956 963 self.flagNoData = True
957 964
958 965 self.flagDiscontinuousBlock = False
959 966
960 967 self.utctime = None
961 968
962 969 self.timeZone = None
963 970
964 971 self.dstFlag = None
965 972
966 973 self.errorCount = None
967 974
968 975 self.blocksize = None
969 976
970 977 self.flagDecodeData = False #asumo q la data no esta decodificada
971 978
972 979 self.flagDeflipData = False #asumo q la data no esta sin flip
973 980
974 981 self.pairsList = None
975 982
976 983 self.nPoints = None
977 984
978 985 def getLagTRange(self, extrapoints=0):
979 986
980 987 lagTRange = self.lagT
981 988 diff = lagTRange[1] - lagTRange[0]
982 989 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
983 990 lagTRange = numpy.hstack((lagTRange, extra))
984 991
985 992 return lagTRange
986 993
987 994 def getLagRRange(self, extrapoints=0):
988 995
989 996 return self.lagR
990 997
991 998 def getPairsList(self):
992 999
993 1000 return self.pairsList
994 1001
995 1002 def getCalculateVelocity(self):
996 1003
997 1004 return self.calculateVelocity
998 1005
999 1006 def getNPoints(self):
1000 1007
1001 1008 return self.nPoints
1002 1009
1003 1010 def getNAvg(self):
1004 1011
1005 1012 return self.nAvg
1006 1013
1007 1014 def getBufferSize(self):
1008 1015
1009 1016 return self.bufferSize
1010 1017
1011 1018 def getPairsAutoCorr(self):
1012 1019 pairsList = self.pairsList
1013 1020 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
1014 1021
1015 1022 for l in range(len(pairsList)):
1016 1023 firstChannel = pairsList[l][0]
1017 1024 secondChannel = pairsList[l][1]
1018 1025
1019 1026 #Obteniendo pares de Autocorrelacion
1020 1027 if firstChannel == secondChannel:
1021 1028 pairsAutoCorr[firstChannel] = int(l)
1022 1029
1023 1030 pairsAutoCorr = pairsAutoCorr.astype(int)
1024 1031
1025 1032 return pairsAutoCorr
1026 1033
1027 1034 def getNoise(self, mode = 2):
1028 1035
1029 1036 indR = numpy.where(self.lagR == 0)[0][0]
1030 1037 indT = numpy.where(self.lagT == 0)[0][0]
1031 1038
1032 1039 jspectra0 = self.data_corr[:,:,indR,:]
1033 1040 jspectra = copy.copy(jspectra0)
1034 1041
1035 1042 num_chan = jspectra.shape[0]
1036 1043 num_hei = jspectra.shape[2]
1037 1044
1038 1045 freq_dc = jspectra.shape[1]/2
1039 1046 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1040 1047
1041 1048 if ind_vel[0]<0:
1042 1049 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1043 1050
1044 1051 if mode == 1:
1045 1052 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1046 1053
1047 1054 if mode == 2:
1048 1055
1049 1056 vel = numpy.array([-2,-1,1,2])
1050 1057 xx = numpy.zeros([4,4])
1051 1058
1052 1059 for fil in range(4):
1053 1060 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1054 1061
1055 1062 xx_inv = numpy.linalg.inv(xx)
1056 1063 xx_aux = xx_inv[0,:]
1057 1064
1058 1065 for ich in range(num_chan):
1059 1066 yy = jspectra[ich,ind_vel,:]
1060 1067 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1061 1068
1062 1069 junkid = jspectra[ich,freq_dc,:]<=0
1063 1070 cjunkid = sum(junkid)
1064 1071
1065 1072 if cjunkid.any():
1066 1073 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1067 1074
1068 1075 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1069 1076
1070 1077 return noise
1071 1078
1072 1079 def getTimeInterval(self):
1073 1080
1074 1081 timeInterval = self.ippSeconds * self.nCohInt * self.nPoints
1075 1082
1076 1083 return timeInterval
1077 1084
1078 1085 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1079 1086 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
1080 1087 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
1081 1088 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
1082 1089 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
1083 1090 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
1084 1091
1085 1092
1086 1093 class Parameters(JROData):
1087 1094
1088 1095 #Information from previous data
1089 1096
1090 1097 inputUnit = None #Type of data to be processed
1091 1098
1092 1099 operation = None #Type of operation to parametrize
1093 1100
1094 1101 normFactor = None #Normalization Factor
1095 1102
1096 1103 groupList = None #List of Pairs, Groups, etc
1097 1104
1098 1105 #Parameters
1099 1106
1100 1107 data_param = None #Parameters obtained
1101 1108
1102 1109 data_pre = None #Data Pre Parametrization
1103 1110
1104 1111 data_SNR = None #Signal to Noise Ratio
1105 1112
1106 1113 # heightRange = None #Heights
1107 1114
1108 1115 abscissaList = None #Abscissa, can be velocities, lags or time
1109 1116
1110 1117 noise = None #Noise Potency
1111 1118
1112 1119 utctimeInit = None #Initial UTC time
1113 1120
1114 1121 paramInterval = None #Time interval to calculate Parameters in seconds
1115 1122
1116 1123 #Fitting
1117 1124
1118 1125 data_error = None #Error of the estimation
1119 1126
1120 1127 constants = None
1121 1128
1122 1129 library = None
1123 1130
1124 1131 #Output signal
1125 1132
1126 1133 outputInterval = None #Time interval to calculate output signal in seconds
1127 1134
1128 1135 data_output = None #Out signal
1129 1136
1130 1137
1131 1138
1132 1139 def __init__(self):
1133 1140 '''
1134 1141 Constructor
1135 1142 '''
1136 1143 self.radarControllerHeaderObj = RadarControllerHeader()
1137 1144
1138 1145 self.systemHeaderObj = SystemHeader()
1139 1146
1140 1147 self.type = "Parameters"
1141 1148
1142 1149 def getTimeRange1(self):
1143 1150
1144 1151 datatime = []
1145 1152
1146 1153 if self.useLocalTime:
1147 1154 time1 = self.utctimeInit - self.timeZone*60
1148 1155 else:
1149 1156 time1 = utctimeInit
1150 1157
1151 1158 # datatime.append(self.utctimeInit)
1152 1159 # datatime.append(self.utctimeInit + self.outputInterval)
1153 1160 datatime.append(time1)
1154 1161 datatime.append(time1 + self.outputInterval)
1155 1162
1156 1163 datatime = numpy.array(datatime)
1157 1164
1158 1165 return datatime
General Comments 0
You need to be logged in to leave comments. Login now