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