##// END OF EJS Templates
schain
RSS
jrodata.py: plotting attribute added to SpectraData
Miguel Valdez -
r773:45e0100d584a
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,1165 +1,1167
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 plotting = "spectra"
527
526 528 def __init__(self):
527 529 '''
528 530 Constructor
529 531 '''
530 532
531 533 self.useLocalTime = True
532 534
533 535 self.radarControllerHeaderObj = RadarControllerHeader()
534 536
535 537 self.systemHeaderObj = SystemHeader()
536 538
537 539 self.type = "Spectra"
538 540
539 541 # self.data = None
540 542
541 543 # self.dtype = None
542 544
543 545 # self.nChannels = 0
544 546
545 547 # self.nHeights = 0
546 548
547 549 self.nProfiles = None
548 550
549 551 self.heightList = None
550 552
551 553 self.channelList = None
552 554
553 555 # self.channelIndexList = None
554 556
555 557 self.pairsList = None
556 558
557 559 self.flagNoData = True
558 560
559 561 self.flagDiscontinuousBlock = False
560 562
561 563 self.utctime = None
562 564
563 565 self.nCohInt = None
564 566
565 567 self.nIncohInt = None
566 568
567 569 self.blocksize = None
568 570
569 571 self.nFFTPoints = None
570 572
571 573 self.wavelength = None
572 574
573 575 self.flagDecodeData = False #asumo q la data no esta decodificada
574 576
575 577 self.flagDeflipData = False #asumo q la data no esta sin flip
576 578
577 579 self.flagShiftFFT = False
578 580
579 581 self.ippFactor = 1
580 582
581 583 #self.noise = None
582 584
583 585 self.beacon_heiIndexList = []
584 586
585 587 self.noise_estimation = None
586 588
587 589
588 590 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
589 591 """
590 592 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
591 593
592 594 Return:
593 595 noiselevel
594 596 """
595 597
596 598 noise = numpy.zeros(self.nChannels)
597 599
598 600 for channel in range(self.nChannels):
599 601 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
600 602 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
601 603
602 604 return noise
603 605
604 606 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
605 607
606 608 if self.noise_estimation is not None:
607 609 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
608 610 else:
609 611 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
610 612 return noise
611 613
612 614 def getFreqRangeTimeResponse(self, extrapoints=0):
613 615
614 616 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
615 617 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
616 618
617 619 return freqrange
618 620
619 621 def getAcfRange(self, extrapoints=0):
620 622
621 623 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
622 624 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
623 625
624 626 return freqrange
625 627
626 628 def getFreqRange(self, extrapoints=0):
627 629
628 630 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
629 631 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
630 632
631 633 return freqrange
632 634
633 635 def getVelRange(self, extrapoints=0):
634 636
635 637 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
636 638 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
637 639
638 640 return velrange
639 641
640 642 def getNPairs(self):
641 643
642 644 return len(self.pairsList)
643 645
644 646 def getPairsIndexList(self):
645 647
646 648 return range(self.nPairs)
647 649
648 650 def getNormFactor(self):
649 651
650 652 pwcode = 1
651 653
652 654 if self.flagDecodeData:
653 655 pwcode = numpy.sum(self.code[0]**2)
654 656 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
655 657 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
656 658
657 659 return normFactor
658 660
659 661 def getFlagCspc(self):
660 662
661 663 if self.data_cspc is None:
662 664 return True
663 665
664 666 return False
665 667
666 668 def getFlagDc(self):
667 669
668 670 if self.data_dc is None:
669 671 return True
670 672
671 673 return False
672 674
673 675 def getTimeInterval(self):
674 676
675 677 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
676 678
677 679 return timeInterval
678 680
679 681 def setValue(self, value):
680 682
681 683 print "This property should not be initialized"
682 684
683 685 return
684 686
685 687 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
686 688 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
687 689 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
688 690 flag_cspc = property(getFlagCspc, setValue)
689 691 flag_dc = property(getFlagDc, setValue)
690 692 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
691 693 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
692 694
693 695 class SpectraHeis(Spectra):
694 696
695 697 data_spc = None
696 698
697 699 data_cspc = None
698 700
699 701 data_dc = None
700 702
701 703 nFFTPoints = None
702 704
703 705 # nPairs = None
704 706
705 707 pairsList = None
706 708
707 709 nCohInt = None
708 710
709 711 nIncohInt = None
710 712
711 713 def __init__(self):
712 714
713 715 self.radarControllerHeaderObj = RadarControllerHeader()
714 716
715 717 self.systemHeaderObj = SystemHeader()
716 718
717 719 self.type = "SpectraHeis"
718 720
719 721 # self.dtype = None
720 722
721 723 # self.nChannels = 0
722 724
723 725 # self.nHeights = 0
724 726
725 727 self.nProfiles = None
726 728
727 729 self.heightList = None
728 730
729 731 self.channelList = None
730 732
731 733 # self.channelIndexList = None
732 734
733 735 self.flagNoData = True
734 736
735 737 self.flagDiscontinuousBlock = False
736 738
737 739 # self.nPairs = 0
738 740
739 741 self.utctime = None
740 742
741 743 self.blocksize = None
742 744
743 745 self.profileIndex = 0
744 746
745 747 self.nCohInt = 1
746 748
747 749 self.nIncohInt = 1
748 750
749 751 def getNormFactor(self):
750 752 pwcode = 1
751 753 if self.flagDecodeData:
752 754 pwcode = numpy.sum(self.code[0]**2)
753 755
754 756 normFactor = self.nIncohInt*self.nCohInt*pwcode
755 757
756 758 return normFactor
757 759
758 760 def getTimeInterval(self):
759 761
760 762 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
761 763
762 764 return timeInterval
763 765
764 766 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
765 767 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
766 768
767 769 class Fits(JROData):
768 770
769 771 heightList = None
770 772
771 773 channelList = None
772 774
773 775 flagNoData = True
774 776
775 777 flagDiscontinuousBlock = False
776 778
777 779 useLocalTime = False
778 780
779 781 utctime = None
780 782
781 783 timeZone = None
782 784
783 785 # ippSeconds = None
784 786
785 787 # timeInterval = None
786 788
787 789 nCohInt = None
788 790
789 791 nIncohInt = None
790 792
791 793 noise = None
792 794
793 795 windowOfFilter = 1
794 796
795 797 #Speed of ligth
796 798 C = 3e8
797 799
798 800 frequency = 49.92e6
799 801
800 802 realtime = False
801 803
802 804
803 805 def __init__(self):
804 806
805 807 self.type = "Fits"
806 808
807 809 self.nProfiles = None
808 810
809 811 self.heightList = None
810 812
811 813 self.channelList = None
812 814
813 815 # self.channelIndexList = None
814 816
815 817 self.flagNoData = True
816 818
817 819 self.utctime = None
818 820
819 821 self.nCohInt = 1
820 822
821 823 self.nIncohInt = 1
822 824
823 825 self.useLocalTime = True
824 826
825 827 self.profileIndex = 0
826 828
827 829 # self.utctime = None
828 830 # self.timeZone = None
829 831 # self.ltctime = None
830 832 # self.timeInterval = None
831 833 # self.header = None
832 834 # self.data_header = None
833 835 # self.data = None
834 836 # self.datatime = None
835 837 # self.flagNoData = False
836 838 # self.expName = ''
837 839 # self.nChannels = None
838 840 # self.nSamples = None
839 841 # self.dataBlocksPerFile = None
840 842 # self.comments = ''
841 843 #
842 844
843 845
844 846 def getltctime(self):
845 847
846 848 if self.useLocalTime:
847 849 return self.utctime - self.timeZone*60
848 850
849 851 return self.utctime
850 852
851 853 def getDatatime(self):
852 854
853 855 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
854 856 return datatime
855 857
856 858 def getTimeRange(self):
857 859
858 860 datatime = []
859 861
860 862 datatime.append(self.ltctime)
861 863 datatime.append(self.ltctime + self.timeInterval)
862 864
863 865 datatime = numpy.array(datatime)
864 866
865 867 return datatime
866 868
867 869 def getHeiRange(self):
868 870
869 871 heis = self.heightList
870 872
871 873 return heis
872 874
873 875 def getNHeights(self):
874 876
875 877 return len(self.heightList)
876 878
877 879 def getNChannels(self):
878 880
879 881 return len(self.channelList)
880 882
881 883 def getChannelIndexList(self):
882 884
883 885 return range(self.nChannels)
884 886
885 887 def getNoise(self, type = 1):
886 888
887 889 #noise = numpy.zeros(self.nChannels)
888 890
889 891 if type == 1:
890 892 noise = self.getNoisebyHildebrand()
891 893
892 894 if type == 2:
893 895 noise = self.getNoisebySort()
894 896
895 897 if type == 3:
896 898 noise = self.getNoisebyWindow()
897 899
898 900 return noise
899 901
900 902 def getTimeInterval(self):
901 903
902 904 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
903 905
904 906 return timeInterval
905 907
906 908 datatime = property(getDatatime, "I'm the 'datatime' property")
907 909 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
908 910 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
909 911 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
910 912 noise = property(getNoise, "I'm the 'nHeights' property.")
911 913
912 914 ltctime = property(getltctime, "I'm the 'ltctime' property")
913 915 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
914 916
915 917 class Correlation(JROData):
916 918
917 919 noise = None
918 920
919 921 SNR = None
920 922
921 923 pairsAutoCorr = None #Pairs of Autocorrelation
922 924
923 925 #--------------------------------------------------
924 926
925 927 data_corr = None
926 928
927 929 data_volt = None
928 930
929 931 lagT = None # each element value is a profileIndex
930 932
931 933 lagR = None # each element value is in km
932 934
933 935 pairsList = None
934 936
935 937 calculateVelocity = None
936 938
937 939 nPoints = None
938 940
939 941 nAvg = None
940 942
941 943 bufferSize = None
942 944
943 945 def __init__(self):
944 946 '''
945 947 Constructor
946 948 '''
947 949 self.radarControllerHeaderObj = RadarControllerHeader()
948 950
949 951 self.systemHeaderObj = SystemHeader()
950 952
951 953 self.type = "Correlation"
952 954
953 955 self.data = None
954 956
955 957 self.dtype = None
956 958
957 959 self.nProfiles = None
958 960
959 961 self.heightList = None
960 962
961 963 self.channelList = None
962 964
963 965 self.flagNoData = True
964 966
965 967 self.flagDiscontinuousBlock = False
966 968
967 969 self.utctime = None
968 970
969 971 self.timeZone = None
970 972
971 973 self.dstFlag = None
972 974
973 975 self.errorCount = None
974 976
975 977 self.blocksize = None
976 978
977 979 self.flagDecodeData = False #asumo q la data no esta decodificada
978 980
979 981 self.flagDeflipData = False #asumo q la data no esta sin flip
980 982
981 983 self.pairsList = None
982 984
983 985 self.nPoints = None
984 986
985 987 def getLagTRange(self, extrapoints=0):
986 988
987 989 lagTRange = self.lagT
988 990 diff = lagTRange[1] - lagTRange[0]
989 991 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
990 992 lagTRange = numpy.hstack((lagTRange, extra))
991 993
992 994 return lagTRange
993 995
994 996 def getLagRRange(self, extrapoints=0):
995 997
996 998 return self.lagR
997 999
998 1000 def getPairsList(self):
999 1001
1000 1002 return self.pairsList
1001 1003
1002 1004 def getCalculateVelocity(self):
1003 1005
1004 1006 return self.calculateVelocity
1005 1007
1006 1008 def getNPoints(self):
1007 1009
1008 1010 return self.nPoints
1009 1011
1010 1012 def getNAvg(self):
1011 1013
1012 1014 return self.nAvg
1013 1015
1014 1016 def getBufferSize(self):
1015 1017
1016 1018 return self.bufferSize
1017 1019
1018 1020 def getPairsAutoCorr(self):
1019 1021 pairsList = self.pairsList
1020 1022 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
1021 1023
1022 1024 for l in range(len(pairsList)):
1023 1025 firstChannel = pairsList[l][0]
1024 1026 secondChannel = pairsList[l][1]
1025 1027
1026 1028 #Obteniendo pares de Autocorrelacion
1027 1029 if firstChannel == secondChannel:
1028 1030 pairsAutoCorr[firstChannel] = int(l)
1029 1031
1030 1032 pairsAutoCorr = pairsAutoCorr.astype(int)
1031 1033
1032 1034 return pairsAutoCorr
1033 1035
1034 1036 def getNoise(self, mode = 2):
1035 1037
1036 1038 indR = numpy.where(self.lagR == 0)[0][0]
1037 1039 indT = numpy.where(self.lagT == 0)[0][0]
1038 1040
1039 1041 jspectra0 = self.data_corr[:,:,indR,:]
1040 1042 jspectra = copy.copy(jspectra0)
1041 1043
1042 1044 num_chan = jspectra.shape[0]
1043 1045 num_hei = jspectra.shape[2]
1044 1046
1045 1047 freq_dc = jspectra.shape[1]/2
1046 1048 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1047 1049
1048 1050 if ind_vel[0]<0:
1049 1051 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1050 1052
1051 1053 if mode == 1:
1052 1054 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1053 1055
1054 1056 if mode == 2:
1055 1057
1056 1058 vel = numpy.array([-2,-1,1,2])
1057 1059 xx = numpy.zeros([4,4])
1058 1060
1059 1061 for fil in range(4):
1060 1062 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1061 1063
1062 1064 xx_inv = numpy.linalg.inv(xx)
1063 1065 xx_aux = xx_inv[0,:]
1064 1066
1065 1067 for ich in range(num_chan):
1066 1068 yy = jspectra[ich,ind_vel,:]
1067 1069 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1068 1070
1069 1071 junkid = jspectra[ich,freq_dc,:]<=0
1070 1072 cjunkid = sum(junkid)
1071 1073
1072 1074 if cjunkid.any():
1073 1075 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1074 1076
1075 1077 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1076 1078
1077 1079 return noise
1078 1080
1079 1081 def getTimeInterval(self):
1080 1082
1081 1083 timeInterval = self.ippSeconds * self.nCohInt * self.nPoints
1082 1084
1083 1085 return timeInterval
1084 1086
1085 1087 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1086 1088 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
1087 1089 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
1088 1090 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
1089 1091 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
1090 1092 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
1091 1093
1092 1094
1093 1095 class Parameters(JROData):
1094 1096
1095 1097 #Information from previous data
1096 1098
1097 1099 inputUnit = None #Type of data to be processed
1098 1100
1099 1101 operation = None #Type of operation to parametrize
1100 1102
1101 1103 normFactor = None #Normalization Factor
1102 1104
1103 1105 groupList = None #List of Pairs, Groups, etc
1104 1106
1105 1107 #Parameters
1106 1108
1107 1109 data_param = None #Parameters obtained
1108 1110
1109 1111 data_pre = None #Data Pre Parametrization
1110 1112
1111 1113 data_SNR = None #Signal to Noise Ratio
1112 1114
1113 1115 # heightRange = None #Heights
1114 1116
1115 1117 abscissaList = None #Abscissa, can be velocities, lags or time
1116 1118
1117 1119 noise = None #Noise Potency
1118 1120
1119 1121 utctimeInit = None #Initial UTC time
1120 1122
1121 1123 paramInterval = None #Time interval to calculate Parameters in seconds
1122 1124
1123 1125 #Fitting
1124 1126
1125 1127 data_error = None #Error of the estimation
1126 1128
1127 1129 constants = None
1128 1130
1129 1131 library = None
1130 1132
1131 1133 #Output signal
1132 1134
1133 1135 outputInterval = None #Time interval to calculate output signal in seconds
1134 1136
1135 1137 data_output = None #Out signal
1136 1138
1137 1139
1138 1140
1139 1141 def __init__(self):
1140 1142 '''
1141 1143 Constructor
1142 1144 '''
1143 1145 self.radarControllerHeaderObj = RadarControllerHeader()
1144 1146
1145 1147 self.systemHeaderObj = SystemHeader()
1146 1148
1147 1149 self.type = "Parameters"
1148 1150
1149 1151 def getTimeRange1(self):
1150 1152
1151 1153 datatime = []
1152 1154
1153 1155 if self.useLocalTime:
1154 1156 time1 = self.utctimeInit - self.timeZone*60
1155 1157 else:
1156 1158 time1 = utctimeInit
1157 1159
1158 1160 # datatime.append(self.utctimeInit)
1159 1161 # datatime.append(self.utctimeInit + self.outputInterval)
1160 1162 datatime.append(time1)
1161 1163 datatime.append(time1 + self.outputInterval)
1162 1164
1163 1165 datatime = numpy.array(datatime)
1164 1166
1165 1167 return datatime
General Comments 0
You need to be logged in to leave comments. Login now