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