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