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