##// END OF EJS Templates
schain-jc
RSS
-Corrections to Vmax calculation...
Julio Valdez -
r850:5a02664c0d53
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,1183 +1,1183
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12
13 13 def getNumpyDtype(dataTypeCode):
14 14
15 15 if dataTypeCode == 0:
16 16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 17 elif dataTypeCode == 1:
18 18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 19 elif dataTypeCode == 2:
20 20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 21 elif dataTypeCode == 3:
22 22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 23 elif dataTypeCode == 4:
24 24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 25 elif dataTypeCode == 5:
26 26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 27 else:
28 28 raise ValueError, 'dataTypeCode was not defined'
29 29
30 30 return numpyDtype
31 31
32 32 def getDataTypeCode(numpyDtype):
33 33
34 34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 35 datatype = 0
36 36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 37 datatype = 1
38 38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 39 datatype = 2
40 40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 41 datatype = 3
42 42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 43 datatype = 4
44 44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 45 datatype = 5
46 46 else:
47 47 datatype = None
48 48
49 49 return datatype
50 50
51 51 def hildebrand_sekhon(data, navg):
52 52 """
53 53 This method is for the objective determination of the noise level in Doppler spectra. This
54 54 implementation technique is based on the fact that the standard deviation of the spectral
55 55 densities is equal to the mean spectral density for white Gaussian noise
56 56
57 57 Inputs:
58 58 Data : heights
59 59 navg : numbers of averages
60 60
61 61 Return:
62 62 -1 : any error
63 63 anoise : noise's level
64 64 """
65 65
66 66 sortdata = numpy.sort(data,axis=None)
67 67 lenOfData = len(sortdata)
68 68 nums_min = lenOfData*0.2
69 69
70 70 if nums_min <= 5:
71 71 nums_min = 5
72 72
73 73 sump = 0.
74 74
75 75 sumq = 0.
76 76
77 77 j = 0
78 78
79 79 cont = 1
80 80
81 81 while((cont==1)and(j<lenOfData)):
82 82
83 83 sump += sortdata[j]
84 84
85 85 sumq += sortdata[j]**2
86 86
87 87 if j > nums_min:
88 88 rtest = float(j)/(j-1) + 1.0/navg
89 89 if ((sumq*j) > (rtest*sump**2)):
90 90 j = j - 1
91 91 sump = sump - sortdata[j]
92 92 sumq = sumq - sortdata[j]**2
93 93 cont = 0
94 94
95 95 j += 1
96 96
97 97 lnoise = sump /j
98 98 # stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
99 99 return lnoise
100 100
101 101 class Beam:
102 102 def __init__(self):
103 103 self.codeList = []
104 104 self.azimuthList = []
105 105 self.zenithList = []
106 106
107 107 class GenericData(object):
108 108
109 109 flagNoData = True
110 110
111 111 def __init__(self):
112 112
113 113 raise NotImplementedError
114 114
115 115 def copy(self, inputObj=None):
116 116
117 117 if inputObj == None:
118 118 return copy.deepcopy(self)
119 119
120 120 for key in inputObj.__dict__.keys():
121 121
122 122 attribute = inputObj.__dict__[key]
123 123
124 124 #If this attribute is a tuple or list
125 125 if type(inputObj.__dict__[key]) in (tuple, list):
126 126 self.__dict__[key] = attribute[:]
127 127 continue
128 128
129 129 #If this attribute is another object or instance
130 130 if hasattr(attribute, '__dict__'):
131 131 self.__dict__[key] = attribute.copy()
132 132 continue
133 133
134 134 self.__dict__[key] = inputObj.__dict__[key]
135 135
136 136 def deepcopy(self):
137 137
138 138 return copy.deepcopy(self)
139 139
140 140 def isEmpty(self):
141 141
142 142 return self.flagNoData
143 143
144 144 class JROData(GenericData):
145 145
146 146 # m_BasicHeader = BasicHeader()
147 147 # m_ProcessingHeader = ProcessingHeader()
148 148
149 149 systemHeaderObj = SystemHeader()
150 150
151 151 radarControllerHeaderObj = RadarControllerHeader()
152 152
153 153 # data = None
154 154
155 155 type = None
156 156
157 157 datatype = None #dtype but in string
158 158
159 159 # dtype = None
160 160
161 161 # nChannels = None
162 162
163 163 # nHeights = None
164 164
165 165 nProfiles = None
166 166
167 167 heightList = None
168 168
169 169 channelList = None
170 170
171 171 flagDiscontinuousBlock = False
172 172
173 173 useLocalTime = False
174 174
175 175 utctime = None
176 176
177 177 timeZone = None
178 178
179 179 dstFlag = None
180 180
181 181 errorCount = None
182 182
183 183 blocksize = None
184 184
185 185 # nCode = None
186 186 #
187 187 # nBaud = None
188 188 #
189 189 # code = None
190 190
191 191 flagDecodeData = False #asumo q la data no esta decodificada
192 192
193 193 flagDeflipData = False #asumo q la data no esta sin flip
194 194
195 195 flagShiftFFT = False
196 196
197 197 # ippSeconds = None
198 198
199 199 # timeInterval = None
200 200
201 201 nCohInt = None
202 202
203 203 # noise = None
204 204
205 205 windowOfFilter = 1
206 206
207 207 #Speed of ligth
208 208 C = 3e8
209 209
210 210 frequency = 49.92e6
211 211
212 212 realtime = False
213 213
214 214 beacon_heiIndexList = None
215 215
216 216 last_block = None
217 217
218 218 blocknow = None
219 219
220 220 azimuth = None
221 221
222 222 zenith = None
223 223
224 224 beam = Beam()
225 225
226 226 profileIndex = None
227 227
228 228 def __init__(self):
229 229
230 230 raise NotImplementedError
231 231
232 232 def getNoise(self):
233 233
234 234 raise NotImplementedError
235 235
236 236 def getNChannels(self):
237 237
238 238 return len(self.channelList)
239 239
240 240 def getChannelIndexList(self):
241 241
242 242 return range(self.nChannels)
243 243
244 244 def getNHeights(self):
245 245
246 246 return len(self.heightList)
247 247
248 248 def getHeiRange(self, extrapoints=0):
249 249
250 250 heis = self.heightList
251 251 # deltah = self.heightList[1] - self.heightList[0]
252 252 #
253 253 # heis.append(self.heightList[-1])
254 254
255 255 return heis
256 256
257 257 def getDeltaH(self):
258 258
259 259 delta = self.heightList[1] - self.heightList[0]
260 260
261 261 return delta
262 262
263 263 def getltctime(self):
264 264
265 265 if self.useLocalTime:
266 266 return self.utctime - self.timeZone*60
267 267
268 268 return self.utctime
269 269
270 270 def getDatatime(self):
271 271
272 272 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
273 273 return datatimeValue
274 274
275 275 def getTimeRange(self):
276 276
277 277 datatime = []
278 278
279 279 datatime.append(self.ltctime)
280 280 datatime.append(self.ltctime + self.timeInterval+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 vmax = self.getFmax() * _lambda
308 vmax = self.getFmax() * _lambda/2
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 powerdB = 10*numpy.log10(power.real)
488 powerdB = numpy.squeeze(powerdB)
487 489
488 return 10*numpy.log10(power.real)
490 return powerdB
489 491
490 492 def getTimeInterval(self):
491 493
492 494 timeInterval = self.ippSeconds * self.nCohInt
493 495
494 496 return timeInterval
495 497
496 498 noise = property(getNoise, "I'm the 'nHeights' property.")
497 499 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
498 500
499 501 class Spectra(JROData):
500 502
501 503 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
502 504 data_spc = None
503 505
504 506 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
505 507 data_cspc = None
506 508
507 509 #data dc es un numpy array de 2 dmensiones (canales, alturas)
508 510 data_dc = None
509 511
510 512 #data power
511 513 data_pwr = None
512 514
513 515 nFFTPoints = None
514 516
515 517 # nPairs = None
516 518
517 519 pairsList = None
518 520
519 521 nIncohInt = None
520 522
521 523 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
522 524
523 525 nCohInt = None #se requiere para determinar el valor de timeInterval
524 526
525 527 ippFactor = None
526 528
527 529 profileIndex = 0
528 530
529 531 plotting = "spectra"
530 532
531 533 def __init__(self):
532 534 '''
533 535 Constructor
534 536 '''
535 537
536 538 self.useLocalTime = True
537 539
538 540 self.radarControllerHeaderObj = RadarControllerHeader()
539 541
540 542 self.systemHeaderObj = SystemHeader()
541 543
542 544 self.type = "Spectra"
543 545
544 546 # self.data = None
545 547
546 548 # self.dtype = None
547 549
548 550 # self.nChannels = 0
549 551
550 552 # self.nHeights = 0
551 553
552 554 self.nProfiles = None
553 555
554 556 self.heightList = None
555 557
556 558 self.channelList = None
557 559
558 560 # self.channelIndexList = None
559 561
560 562 self.pairsList = None
561 563
562 564 self.flagNoData = True
563 565
564 566 self.flagDiscontinuousBlock = False
565 567
566 568 self.utctime = None
567 569
568 570 self.nCohInt = None
569 571
570 572 self.nIncohInt = None
571 573
572 574 self.blocksize = None
573 575
574 576 self.nFFTPoints = None
575 577
576 578 self.wavelength = None
577 579
578 580 self.flagDecodeData = False #asumo q la data no esta decodificada
579 581
580 582 self.flagDeflipData = False #asumo q la data no esta sin flip
581 583
582 584 self.flagShiftFFT = False
583 585
584 586 self.ippFactor = 1
585 587
586 588 #self.noise = None
587 589
588 590 self.beacon_heiIndexList = []
589 591
590 592 self.noise_estimation = None
591 593
592 594
593 595 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
594 596 """
595 597 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
596 598
597 599 Return:
598 600 noiselevel
599 601 """
600 602
601 603 noise = numpy.zeros(self.nChannels)
602 604
603 605 for channel in range(self.nChannels):
604 606 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
605 607 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
606 608
607 609 return noise
608 610
609 611 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
610 612
611 613 if self.noise_estimation is not None:
612 614 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
613 615 else:
614 616 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
615 617 return noise
616 618
617 619 def getFreqRangeTimeResponse(self, extrapoints=0):
618 620
619 621 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
620 622 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
621 623
622 624 return freqrange
623 625
624 626 def getAcfRange(self, extrapoints=0):
625 627
626 628 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
627 629 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
628 630
629 631 return freqrange
630 632
631 633 def getFreqRange(self, extrapoints=0):
632 634
633 635 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
634 636 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
635 637
636 638 return freqrange
637 639
638 640 def getVelRange(self, extrapoints=0):
639 641
640 642 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
641 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
643 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
642 644
643 645 return velrange
644 646
645 647 def getNPairs(self):
646 648
647 649 return len(self.pairsList)
648 650
649 651 def getPairsIndexList(self):
650 652
651 653 return range(self.nPairs)
652 654
653 655 def getNormFactor(self):
654 656
655 657 pwcode = 1
656 658
657 659 if self.flagDecodeData:
658 660 pwcode = numpy.sum(self.code[0]**2)
659 661 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
660 662 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
661 663
662 664 return normFactor
663 665
664 666 def getFlagCspc(self):
665 667
666 668 if self.data_cspc is None:
667 669 return True
668 670
669 671 return False
670 672
671 673 def getFlagDc(self):
672 674
673 675 if self.data_dc is None:
674 676 return True
675 677
676 678 return False
677 679
678 680 def getTimeInterval(self):
679 681
680 682 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
681 683
682 684 return timeInterval
683 685
684 686 def getPower(self):
685 687
686 688 factor = self.normFactor
687 689 z = self.data_spc/factor
688 690 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
689 691 avg = numpy.average(z, axis=1)
690 692
691 693 return 10*numpy.log10(avg)
692 694
693 695 def setValue(self, value):
694 696
695 697 print "This property should not be initialized"
696 698
697 699 return
698 700
699 701 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
700 702 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
701 703 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
702 704 flag_cspc = property(getFlagCspc, setValue)
703 705 flag_dc = property(getFlagDc, setValue)
704 706 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
705 707 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
706 708
707 709 class SpectraHeis(Spectra):
708 710
709 711 data_spc = None
710 712
711 713 data_cspc = None
712 714
713 715 data_dc = None
714 716
715 717 nFFTPoints = None
716 718
717 719 # nPairs = None
718 720
719 721 pairsList = None
720 722
721 723 nCohInt = None
722 724
723 725 nIncohInt = None
724 726
725 727 def __init__(self):
726 728
727 729 self.radarControllerHeaderObj = RadarControllerHeader()
728 730
729 731 self.systemHeaderObj = SystemHeader()
730 732
731 733 self.type = "SpectraHeis"
732 734
733 735 # self.dtype = None
734 736
735 737 # self.nChannels = 0
736 738
737 739 # self.nHeights = 0
738 740
739 741 self.nProfiles = None
740 742
741 743 self.heightList = None
742 744
743 745 self.channelList = None
744 746
745 747 # self.channelIndexList = None
746 748
747 749 self.flagNoData = True
748 750
749 751 self.flagDiscontinuousBlock = False
750 752
751 753 # self.nPairs = 0
752 754
753 755 self.utctime = None
754 756
755 757 self.blocksize = None
756 758
757 759 self.profileIndex = 0
758 760
759 761 self.nCohInt = 1
760 762
761 763 self.nIncohInt = 1
762 764
763 765 def getNormFactor(self):
764 766 pwcode = 1
765 767 if self.flagDecodeData:
766 768 pwcode = numpy.sum(self.code[0]**2)
767 769
768 770 normFactor = self.nIncohInt*self.nCohInt*pwcode
769 771
770 772 return normFactor
771 773
772 774 def getTimeInterval(self):
773 775
774 776 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
775 777
776 778 return timeInterval
777 779
778 780 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
779 781 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
780 782
781 783 class Fits(JROData):
782 784
783 785 heightList = None
784 786
785 787 channelList = None
786 788
787 789 flagNoData = True
788 790
789 791 flagDiscontinuousBlock = False
790 792
791 793 useLocalTime = False
792 794
793 795 utctime = None
794 796
795 797 timeZone = None
796 798
797 799 # ippSeconds = None
798 800
799 801 # timeInterval = None
800 802
801 803 nCohInt = None
802 804
803 805 nIncohInt = None
804 806
805 807 noise = None
806 808
807 809 windowOfFilter = 1
808 810
809 811 #Speed of ligth
810 812 C = 3e8
811 813
812 814 frequency = 49.92e6
813 815
814 816 realtime = False
815 817
816 818
817 819 def __init__(self):
818 820
819 821 self.type = "Fits"
820 822
821 823 self.nProfiles = None
822 824
823 825 self.heightList = None
824 826
825 827 self.channelList = None
826 828
827 829 # self.channelIndexList = None
828 830
829 831 self.flagNoData = True
830 832
831 833 self.utctime = None
832 834
833 835 self.nCohInt = 1
834 836
835 837 self.nIncohInt = 1
836 838
837 839 self.useLocalTime = True
838 840
839 841 self.profileIndex = 0
840 842
841 843 # self.utctime = None
842 844 # self.timeZone = None
843 845 # self.ltctime = None
844 846 # self.timeInterval = None
845 847 # self.header = None
846 848 # self.data_header = None
847 849 # self.data = None
848 850 # self.datatime = None
849 851 # self.flagNoData = False
850 852 # self.expName = ''
851 853 # self.nChannels = None
852 854 # self.nSamples = None
853 855 # self.dataBlocksPerFile = None
854 856 # self.comments = ''
855 857 #
856 858
857 859
858 860 def getltctime(self):
859 861
860 862 if self.useLocalTime:
861 863 return self.utctime - self.timeZone*60
862 864
863 865 return self.utctime
864 866
865 867 def getDatatime(self):
866 868
867 869 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
868 870 return datatime
869 871
870 872 def getTimeRange(self):
871 873
872 874 datatime = []
873 875
874 876 datatime.append(self.ltctime)
875 877 datatime.append(self.ltctime + self.timeInterval)
876 878
877 879 datatime = numpy.array(datatime)
878 880
879 881 return datatime
880 882
881 883 def getHeiRange(self):
882 884
883 885 heis = self.heightList
884 886
885 887 return heis
886 888
887 889 def getNHeights(self):
888 890
889 891 return len(self.heightList)
890 892
891 893 def getNChannels(self):
892 894
893 895 return len(self.channelList)
894 896
895 897 def getChannelIndexList(self):
896 898
897 899 return range(self.nChannels)
898 900
899 901 def getNoise(self, type = 1):
900 902
901 903 #noise = numpy.zeros(self.nChannels)
902 904
903 905 if type == 1:
904 906 noise = self.getNoisebyHildebrand()
905 907
906 908 if type == 2:
907 909 noise = self.getNoisebySort()
908 910
909 911 if type == 3:
910 912 noise = self.getNoisebyWindow()
911 913
912 914 return noise
913 915
914 916 def getTimeInterval(self):
915 917
916 918 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
917 919
918 920 return timeInterval
919 921
920 922 datatime = property(getDatatime, "I'm the 'datatime' property")
921 923 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
922 924 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
923 925 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
924 926 noise = property(getNoise, "I'm the 'nHeights' property.")
925 927
926 928 ltctime = property(getltctime, "I'm the 'ltctime' property")
927 929 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
928 930
931
929 932 class Correlation(JROData):
930 933
931 934 noise = None
932 935
933 936 SNR = None
934 937
935 pairsAutoCorr = None #Pairs of Autocorrelation
936
937 938 #--------------------------------------------------
938 939
939 data_corr = None
940 mode = None
940 941
941 data_volt = None
942 split = False
943
944 data_cf = None
942 945
943 lagT = None # each element value is a profileIndex
946 lags = None
944 947
945 lagR = None # each element value is in km
948 lagRange = None
946 949
947 950 pairsList = None
948 951
949 calculateVelocity = None
952 normFactor = None
950 953
951 nPoints = None
954 #--------------------------------------------------
952 955
953 nAvg = None
956 # calculateVelocity = None
957
958 nLags = None
954 959
955 bufferSize = None
960 nPairs = None
961
962 nAvg = None
963
956 964
957 965 def __init__(self):
958 966 '''
959 967 Constructor
960 968 '''
961 969 self.radarControllerHeaderObj = RadarControllerHeader()
962 970
963 971 self.systemHeaderObj = SystemHeader()
964 972
965 973 self.type = "Correlation"
966 974
967 975 self.data = None
968 976
969 977 self.dtype = None
970 978
971 979 self.nProfiles = None
972 980
973 981 self.heightList = None
974 982
975 983 self.channelList = None
976 984
977 985 self.flagNoData = True
978 986
979 987 self.flagDiscontinuousBlock = False
980 988
981 989 self.utctime = None
982 990
983 991 self.timeZone = None
984 992
985 993 self.dstFlag = None
986 994
987 995 self.errorCount = None
988 996
989 997 self.blocksize = None
990 998
991 999 self.flagDecodeData = False #asumo q la data no esta decodificada
992 1000
993 1001 self.flagDeflipData = False #asumo q la data no esta sin flip
994 1002
995 1003 self.pairsList = None
996 1004
997 1005 self.nPoints = None
998
999 def getLagTRange(self, extrapoints=0):
1000
1001 lagTRange = self.lagT
1002 diff = lagTRange[1] - lagTRange[0]
1003 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
1004 lagTRange = numpy.hstack((lagTRange, extra))
1005
1006 return lagTRange
1007
1008 def getLagRRange(self, extrapoints=0):
1009
1010 return self.lagR
1011
1006
1012 1007 def getPairsList(self):
1013 1008
1014 1009 return self.pairsList
1015
1016 def getCalculateVelocity(self):
1017
1018 return self.calculateVelocity
1019
1020 def getNPoints(self):
1021
1022 return self.nPoints
1023
1024 def getNAvg(self):
1025
1026 return self.nAvg
1027
1028 def getBufferSize(self):
1029
1030 return self.bufferSize
1031
1032 def getPairsAutoCorr(self):
1033 pairsList = self.pairsList
1034 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
1035
1036 for l in range(len(pairsList)):
1037 firstChannel = pairsList[l][0]
1038 secondChannel = pairsList[l][1]
1039
1040 #Obteniendo pares de Autocorrelacion
1041 if firstChannel == secondChannel:
1042 pairsAutoCorr[firstChannel] = int(l)
1043
1044 pairsAutoCorr = pairsAutoCorr.astype(int)
1045
1046 return pairsAutoCorr
1047
1010
1048 1011 def getNoise(self, mode = 2):
1049 1012
1050 1013 indR = numpy.where(self.lagR == 0)[0][0]
1051 1014 indT = numpy.where(self.lagT == 0)[0][0]
1052 1015
1053 1016 jspectra0 = self.data_corr[:,:,indR,:]
1054 1017 jspectra = copy.copy(jspectra0)
1055 1018
1056 1019 num_chan = jspectra.shape[0]
1057 1020 num_hei = jspectra.shape[2]
1058 1021
1059 1022 freq_dc = jspectra.shape[1]/2
1060 1023 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1061 1024
1062 1025 if ind_vel[0]<0:
1063 1026 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1064 1027
1065 1028 if mode == 1:
1066 1029 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1067 1030
1068 1031 if mode == 2:
1069 1032
1070 1033 vel = numpy.array([-2,-1,1,2])
1071 1034 xx = numpy.zeros([4,4])
1072 1035
1073 1036 for fil in range(4):
1074 1037 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1075 1038
1076 1039 xx_inv = numpy.linalg.inv(xx)
1077 1040 xx_aux = xx_inv[0,:]
1078 1041
1079 1042 for ich in range(num_chan):
1080 1043 yy = jspectra[ich,ind_vel,:]
1081 1044 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1082 1045
1083 1046 junkid = jspectra[ich,freq_dc,:]<=0
1084 1047 cjunkid = sum(junkid)
1085 1048
1086 1049 if cjunkid.any():
1087 1050 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1088 1051
1089 1052 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1090 1053
1091 1054 return noise
1092 1055
1093 1056 def getTimeInterval(self):
1094 1057
1095 timeInterval = self.ippSeconds * self.nCohInt * self.nPoints
1058 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1096 1059
1097 1060 return timeInterval
1098 1061
1062 def splitFunctions(self):
1063
1064 pairsList = self.pairsList
1065 ccf_pairs = []
1066 acf_pairs = []
1067 ccf_ind = []
1068 acf_ind = []
1069 for l in range(len(pairsList)):
1070 chan0 = pairsList[l][0]
1071 chan1 = pairsList[l][1]
1072
1073 #Obteniendo pares de Autocorrelacion
1074 if chan0 == chan1:
1075 acf_pairs.append(chan0)
1076 acf_ind.append(l)
1077 else:
1078 ccf_pairs.append(pairsList[l])
1079 ccf_ind.append(l)
1080
1081 data_acf = self.data_cf[acf_ind]
1082 data_ccf = self.data_cf[ccf_ind]
1083
1084 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1085
1086 def getNormFactor(self):
1087 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1088 acf_pairs = numpy.array(acf_pairs)
1089 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1090
1091 for p in range(self.nPairs):
1092 pair = self.pairsList[p]
1093
1094 ch0 = pair[0]
1095 ch1 = pair[1]
1096
1097 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1098 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1099 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1100
1101 return normFactor
1102
1099 1103 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1100 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
1101 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
1102 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
1103 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
1104 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
1105
1104 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1106 1105
1107 1106 class Parameters(JROData):
1108 1107
1109 1108 experimentInfo = None #Information about the experiment
1110 1109
1111 1110 #Information from previous data
1112 1111
1113 1112 inputUnit = None #Type of data to be processed
1114 1113
1115 1114 operation = None #Type of operation to parametrize
1116 1115
1117 1116 normFactor = None #Normalization Factor
1118 1117
1119 1118 groupList = None #List of Pairs, Groups, etc
1120 1119
1121 1120 #Parameters
1122 1121
1123 1122 data_param = None #Parameters obtained
1124 1123
1125 1124 data_pre = None #Data Pre Parametrization
1126 1125
1127 1126 data_SNR = None #Signal to Noise Ratio
1128 1127
1129 1128 # heightRange = None #Heights
1130 1129
1131 1130 abscissaList = None #Abscissa, can be velocities, lags or time
1132 1131
1133 1132 noise = None #Noise Potency
1134 1133
1135 1134 utctimeInit = None #Initial UTC time
1136 1135
1137 1136 paramInterval = None #Time interval to calculate Parameters in seconds
1138 1137
1139 1138 useLocalTime = True
1140 1139
1141 1140 #Fitting
1142 1141
1143 1142 data_error = None #Error of the estimation
1144 1143
1145 1144 constants = None
1146 1145
1147 1146 library = None
1148 1147
1149 1148 #Output signal
1150 1149
1151 1150 outputInterval = None #Time interval to calculate output signal in seconds
1152 1151
1153 1152 data_output = None #Out signal
1154 1153
1154 nAvg = None
1155 1155
1156 1156
1157 1157 def __init__(self):
1158 1158 '''
1159 1159 Constructor
1160 1160 '''
1161 1161 self.radarControllerHeaderObj = RadarControllerHeader()
1162 1162
1163 1163 self.systemHeaderObj = SystemHeader()
1164 1164
1165 1165 self.type = "Parameters"
1166 1166
1167 1167 def getTimeRange1(self, interval):
1168 1168
1169 1169 datatime = []
1170 1170
1171 1171 if self.useLocalTime:
1172 1172 time1 = self.utctimeInit - self.timeZone*60
1173 1173 else:
1174 1174 time1 = self.utctimeInit
1175 1175
1176 1176 # datatime.append(self.utctimeInit)
1177 1177 # datatime.append(self.utctimeInit + self.outputInterval)
1178 1178 datatime.append(time1)
1179 1179 datatime.append(time1 + interval)
1180 1180
1181 1181 datatime = numpy.array(datatime)
1182 1182
1183 1183 return datatime
General Comments 0
You need to be logged in to leave comments. Login now