##// END OF EJS Templates
schain-jc
RSS
Coherent Integration by Blocks: Bug fixed, nProfiles is divided by number of integrations
Miguel Valdez -
r720:7c23a7f6f4e8
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,34
1 #!/usr/bin/env python
2 '''
3 Created on Jul 7, 2014
4
5 @author: roj-idl71
6 '''
7 import os, sys
8
9 schainpy_path = os.path.dirname(os.getcwd())
10 source_path = os.path.dirname(schainpy_path)
11 sys.path.insert(0, source_path)
12
13 from schainpy.controller_api import ControllerThread
14
15 def main():
16 desc = "Segundo Test"
17 filename = "/Users/miguel/Downloads/mst_blocks.xml"
18
19 controllerObj = ControllerThread()
20 controllerObj.readXml(filename)
21
22 #Configure use of external plotter before start
23 plotterObj = controllerObj.useExternalPlotter()
24 ########################################
25
26 controllerObj.start()
27 plotterObj.start()
28
29
30 if __name__ == '__main__':
31 import time
32 start_time = time.time()
33 main()
34 print("--- %s seconds ---" % (time.time() - start_time)) No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1123 +1,1123
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12
13 13 def getNumpyDtype(dataTypeCode):
14 14
15 15 if dataTypeCode == 0:
16 16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 17 elif dataTypeCode == 1:
18 18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 19 elif dataTypeCode == 2:
20 20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 21 elif dataTypeCode == 3:
22 22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 23 elif dataTypeCode == 4:
24 24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 25 elif dataTypeCode == 5:
26 26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 27 else:
28 28 raise ValueError, 'dataTypeCode was not defined'
29 29
30 30 return numpyDtype
31 31
32 32 def getDataTypeCode(numpyDtype):
33 33
34 34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 35 datatype = 0
36 36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 37 datatype = 1
38 38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 39 datatype = 2
40 40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 41 datatype = 3
42 42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 43 datatype = 4
44 44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 45 datatype = 5
46 46 else:
47 47 datatype = None
48 48
49 49 return datatype
50 50
51 51 def hildebrand_sekhon(data, navg):
52 52 """
53 53 This method is for the objective determination of the noise level in Doppler spectra. This
54 54 implementation technique is based on the fact that the standard deviation of the spectral
55 55 densities is equal to the mean spectral density for white Gaussian noise
56 56
57 57 Inputs:
58 58 Data : heights
59 59 navg : numbers of averages
60 60
61 61 Return:
62 62 -1 : any error
63 63 anoise : noise's level
64 64 """
65 65
66 66 sortdata = numpy.sort(data,axis=None)
67 67 lenOfData = len(sortdata)
68 68 nums_min = lenOfData*0.2
69 69
70 70 if nums_min <= 5:
71 71 nums_min = 5
72 72
73 73 sump = 0.
74 74
75 75 sumq = 0.
76 76
77 77 j = 0
78 78
79 79 cont = 1
80 80
81 81 while((cont==1)and(j<lenOfData)):
82 82
83 83 sump += sortdata[j]
84 84
85 85 sumq += sortdata[j]**2
86 86
87 87 if j > nums_min:
88 88 rtest = float(j)/(j-1) + 1.0/navg
89 89 if ((sumq*j) > (rtest*sump**2)):
90 90 j = j - 1
91 91 sump = sump - sortdata[j]
92 92 sumq = sumq - sortdata[j]**2
93 93 cont = 0
94 94
95 95 j += 1
96 96
97 97 lnoise = sump /j
98 98 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
99 99 return lnoise
100 100
101 101 class Beam:
102 102 def __init__(self):
103 103 self.codeList = []
104 104 self.azimuthList = []
105 105 self.zenithList = []
106 106
107 107 class GenericData(object):
108 108
109 109 flagNoData = True
110 110
111 111 def __init__(self):
112 112
113 113 raise NotImplementedError
114 114
115 115 def copy(self, inputObj=None):
116 116
117 117 if inputObj == None:
118 118 return copy.deepcopy(self)
119 119
120 120 for key in inputObj.__dict__.keys():
121 121 self.__dict__[key] = inputObj.__dict__[key]
122 122
123 123 def deepcopy(self):
124 124
125 125 return copy.deepcopy(self)
126 126
127 127 def isEmpty(self):
128 128
129 129 return self.flagNoData
130 130
131 131 class JROData(GenericData):
132 132
133 133 # m_BasicHeader = BasicHeader()
134 134 # m_ProcessingHeader = ProcessingHeader()
135 135
136 136 systemHeaderObj = SystemHeader()
137 137
138 138 radarControllerHeaderObj = RadarControllerHeader()
139 139
140 140 # data = None
141 141
142 142 type = None
143 143
144 144 datatype = None #dtype but in string
145 145
146 146 # dtype = None
147 147
148 148 # nChannels = None
149 149
150 150 # nHeights = None
151 151
152 152 nProfiles = None
153 153
154 154 heightList = None
155 155
156 156 channelList = None
157 157
158 158 flagDiscontinuousBlock = False
159 159
160 160 useLocalTime = False
161 161
162 162 utctime = None
163 163
164 164 timeZone = None
165 165
166 166 dstFlag = None
167 167
168 168 errorCount = None
169 169
170 170 blocksize = None
171 171
172 172 # nCode = None
173 173 #
174 174 # nBaud = None
175 175 #
176 176 # code = None
177 177
178 178 flagDecodeData = False #asumo q la data no esta decodificada
179 179
180 180 flagDeflipData = False #asumo q la data no esta sin flip
181 181
182 182 flagShiftFFT = False
183 183
184 184 # ippSeconds = None
185 185
186 186 # timeInterval = None
187 187
188 188 nCohInt = None
189 189
190 190 # noise = None
191 191
192 192 windowOfFilter = 1
193 193
194 194 #Speed of ligth
195 195 C = 3e8
196 196
197 197 frequency = 49.92e6
198 198
199 199 realtime = False
200 200
201 201 beacon_heiIndexList = None
202 202
203 203 last_block = None
204 204
205 205 blocknow = None
206 206
207 207 azimuth = None
208 208
209 209 zenith = None
210 210
211 211 beam = Beam()
212 212
213 213 profileIndex = None
214 214
215 215 def __init__(self):
216 216
217 217 raise NotImplementedError
218 218
219 219 def getNoise(self):
220 220
221 221 raise NotImplementedError
222 222
223 223 def getNChannels(self):
224 224
225 225 return len(self.channelList)
226 226
227 227 def getChannelIndexList(self):
228 228
229 229 return range(self.nChannels)
230 230
231 231 def getNHeights(self):
232 232
233 233 return len(self.heightList)
234 234
235 235 def getHeiRange(self, extrapoints=0):
236 236
237 237 heis = self.heightList
238 238 # deltah = self.heightList[1] - self.heightList[0]
239 239 #
240 240 # heis.append(self.heightList[-1])
241 241
242 242 return heis
243 243
244 244 def getltctime(self):
245 245
246 246 if self.useLocalTime:
247 247 return self.utctime - self.timeZone*60
248 248
249 249 return self.utctime
250 250
251 251 def getDatatime(self):
252 252
253 253 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
254 254 return datatimeValue
255 255
256 256 def getTimeRange(self):
257 257
258 258 datatime = []
259 259
260 260 datatime.append(self.ltctime)
261 261 datatime.append(self.ltctime + self.timeInterval+60)
262 262
263 263 datatime = numpy.array(datatime)
264 264
265 265 return datatime
266 266
267 267 def getFmax(self):
268 268
269 269 PRF = 1./(self.ippSeconds * self.nCohInt)
270 270
271 271 fmax = PRF/2.
272 272
273 273 return fmax
274 274
275 275 def getVmax(self):
276 276
277 277 _lambda = self.C/self.frequency
278 278
279 279 vmax = self.getFmax() * _lambda
280 280
281 281 return vmax
282 282
283 283 def get_ippSeconds(self):
284 284 '''
285 285 '''
286 286 return self.radarControllerHeaderObj.ippSeconds
287 287
288 288 def set_ippSeconds(self, ippSeconds):
289 289 '''
290 290 '''
291 291
292 292 self.radarControllerHeaderObj.ippSeconds = ippSeconds
293 293
294 294 return
295 295
296 296 def get_dtype(self):
297 297 '''
298 298 '''
299 299 return getNumpyDtype(self.datatype)
300 300
301 301 def set_dtype(self, numpyDtype):
302 302 '''
303 303 '''
304 304
305 305 self.datatype = getDataTypeCode(numpyDtype)
306 306
307 307 def get_code(self):
308 308 '''
309 309 '''
310 310 return self.radarControllerHeaderObj.code
311 311
312 312 def set_code(self, code):
313 313 '''
314 314 '''
315 315 self.radarControllerHeaderObj.code = code
316 316
317 317 return
318 318
319 319 def get_ncode(self):
320 320 '''
321 321 '''
322 322 return self.radarControllerHeaderObj.nCode
323 323
324 324 def set_ncode(self, nCode):
325 325 '''
326 326 '''
327 327 self.radarControllerHeaderObj.nCode = nCode
328 328
329 329 return
330 330
331 331 def get_nbaud(self):
332 332 '''
333 333 '''
334 334 return self.radarControllerHeaderObj.nBaud
335 335
336 336 def set_nbaud(self, nBaud):
337 337 '''
338 338 '''
339 339 self.radarControllerHeaderObj.nBaud = nBaud
340 340
341 341 return
342 342
343 343 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
344 344 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
345 345 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
346 346 #noise = property(getNoise, "I'm the 'nHeights' property.")
347 347 datatime = property(getDatatime, "I'm the 'datatime' property")
348 348 ltctime = property(getltctime, "I'm the 'ltctime' property")
349 349 ippSeconds = property(get_ippSeconds, set_ippSeconds)
350 350 dtype = property(get_dtype, set_dtype)
351 351 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
352 352 code = property(get_code, set_code)
353 353 nCode = property(get_ncode, set_ncode)
354 354 nBaud = property(get_nbaud, set_nbaud)
355 355
356 356 class Voltage(JROData):
357 357
358 358 #data es un numpy array de 2 dmensiones (canales, alturas)
359 359 data = None
360 360
361 361 def __init__(self):
362 362 '''
363 363 Constructor
364 364 '''
365 365
366 366 self.useLocalTime = True
367 367
368 368 self.radarControllerHeaderObj = RadarControllerHeader()
369 369
370 370 self.systemHeaderObj = SystemHeader()
371 371
372 372 self.type = "Voltage"
373 373
374 374 self.data = None
375 375
376 376 # self.dtype = None
377 377
378 378 # self.nChannels = 0
379 379
380 380 # self.nHeights = 0
381 381
382 382 self.nProfiles = None
383 383
384 384 self.heightList = None
385 385
386 386 self.channelList = None
387 387
388 388 # self.channelIndexList = None
389 389
390 390 self.flagNoData = True
391 391
392 392 self.flagDiscontinuousBlock = False
393 393
394 394 self.utctime = None
395 395
396 396 self.timeZone = None
397 397
398 398 self.dstFlag = None
399 399
400 400 self.errorCount = None
401 401
402 402 self.nCohInt = None
403 403
404 404 self.blocksize = None
405 405
406 406 self.flagDecodeData = False #asumo q la data no esta decodificada
407 407
408 408 self.flagDeflipData = False #asumo q la data no esta sin flip
409 409
410 410 self.flagShiftFFT = False
411 411
412 412 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
413 413
414 414 self.profileIndex = 0
415 415
416 416 def getNoisebyHildebrand(self, channel = None):
417 417 """
418 418 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
419 419
420 420 Return:
421 421 noiselevel
422 422 """
423 423
424 424 if channel != None:
425 425 data = self.data[channel]
426 426 nChannels = 1
427 427 else:
428 428 data = self.data
429 429 nChannels = self.nChannels
430 430
431 431 noise = numpy.zeros(nChannels)
432 432 power = data * numpy.conjugate(data)
433 433
434 434 for thisChannel in range(nChannels):
435 435 if nChannels == 1:
436 436 daux = power[:].real
437 437 else:
438 438 daux = power[thisChannel,:].real
439 439 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
440 440
441 441 return noise
442 442
443 443 def getNoise(self, type = 1, channel = None):
444 444
445 445 if type == 1:
446 446 noise = self.getNoisebyHildebrand(channel)
447 447
448 return 10*numpy.log10(noise)
448 return noise
449 449
450 450 def getPower(self, channel = None):
451 451
452 452 if channel != None:
453 453 data = self.data[channel]
454 454 else:
455 455 data = self.data
456 456
457 457 power = data * numpy.conjugate(data)
458 458
459 459 return 10*numpy.log10(power.real)
460 460
461 461 def getTimeInterval(self):
462 462
463 463 timeInterval = self.ippSeconds * self.nCohInt
464 464
465 465 return timeInterval
466 466
467 467 noise = property(getNoise, "I'm the 'nHeights' property.")
468 468 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
469 469
470 470 class Spectra(JROData):
471 471
472 472 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
473 473 data_spc = None
474 474
475 475 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
476 476 data_cspc = None
477 477
478 478 #data es un numpy array de 2 dmensiones (canales, alturas)
479 479 data_dc = None
480 480
481 481 nFFTPoints = None
482 482
483 483 # nPairs = None
484 484
485 485 pairsList = None
486 486
487 487 nIncohInt = None
488 488
489 489 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
490 490
491 491 nCohInt = None #se requiere para determinar el valor de timeInterval
492 492
493 493 ippFactor = None
494 494
495 495 profileIndex = 0
496 496
497 497 def __init__(self):
498 498 '''
499 499 Constructor
500 500 '''
501 501
502 502 self.useLocalTime = True
503 503
504 504 self.radarControllerHeaderObj = RadarControllerHeader()
505 505
506 506 self.systemHeaderObj = SystemHeader()
507 507
508 508 self.type = "Spectra"
509 509
510 510 # self.data = None
511 511
512 512 # self.dtype = None
513 513
514 514 # self.nChannels = 0
515 515
516 516 # self.nHeights = 0
517 517
518 518 self.nProfiles = None
519 519
520 520 self.heightList = None
521 521
522 522 self.channelList = None
523 523
524 524 # self.channelIndexList = None
525 525
526 526 self.pairsList = None
527 527
528 528 self.flagNoData = True
529 529
530 530 self.flagDiscontinuousBlock = False
531 531
532 532 self.utctime = None
533 533
534 534 self.nCohInt = None
535 535
536 536 self.nIncohInt = None
537 537
538 538 self.blocksize = None
539 539
540 540 self.nFFTPoints = None
541 541
542 542 self.wavelength = None
543 543
544 544 self.flagDecodeData = False #asumo q la data no esta decodificada
545 545
546 546 self.flagDeflipData = False #asumo q la data no esta sin flip
547 547
548 548 self.flagShiftFFT = False
549 549
550 550 self.ippFactor = 1
551 551
552 552 #self.noise = None
553 553
554 554 self.beacon_heiIndexList = []
555 555
556 556 self.noise_estimation = None
557 557
558 558
559 559 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
560 560 """
561 561 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
562 562
563 563 Return:
564 564 noiselevel
565 565 """
566 566
567 567 noise = numpy.zeros(self.nChannels)
568 568
569 569 for channel in range(self.nChannels):
570 570 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
571 571 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
572 572
573 573 return noise
574 574
575 575 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
576 576
577 577 if self.noise_estimation != None:
578 578 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
579 579 else:
580 580 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
581 581 return noise
582 582
583 583
584 584 def getFreqRange(self, extrapoints=0):
585 585
586 586 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
587 587 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
588 588
589 589 return freqrange
590 590
591 591 def getVelRange(self, extrapoints=0):
592 592
593 593 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
594 594 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
595 595
596 596 return velrange
597 597
598 598 def getNPairs(self):
599 599
600 600 return len(self.pairsList)
601 601
602 602 def getPairsIndexList(self):
603 603
604 604 return range(self.nPairs)
605 605
606 606 def getNormFactor(self):
607 607
608 608 pwcode = 1
609 609
610 610 if self.flagDecodeData:
611 611 pwcode = numpy.sum(self.code[0]**2)
612 612 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
613 613 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
614 614
615 615 return normFactor
616 616
617 617 def getFlagCspc(self):
618 618
619 619 if self.data_cspc is None:
620 620 return True
621 621
622 622 return False
623 623
624 624 def getFlagDc(self):
625 625
626 626 if self.data_dc is None:
627 627 return True
628 628
629 629 return False
630 630
631 631 def getTimeInterval(self):
632 632
633 633 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
634 634
635 635 return timeInterval
636 636
637 637 def setValue(self, value):
638 638
639 639 print "This property should not be initialized"
640 640
641 641 return
642 642
643 643 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
644 644 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
645 645 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
646 646 flag_cspc = property(getFlagCspc, setValue)
647 647 flag_dc = property(getFlagDc, setValue)
648 648 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
649 649 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
650 650
651 651 class SpectraHeis(Spectra):
652 652
653 653 data_spc = None
654 654
655 655 data_cspc = None
656 656
657 657 data_dc = None
658 658
659 659 nFFTPoints = None
660 660
661 661 # nPairs = None
662 662
663 663 pairsList = None
664 664
665 665 nCohInt = None
666 666
667 667 nIncohInt = None
668 668
669 669 def __init__(self):
670 670
671 671 self.radarControllerHeaderObj = RadarControllerHeader()
672 672
673 673 self.systemHeaderObj = SystemHeader()
674 674
675 675 self.type = "SpectraHeis"
676 676
677 677 # self.dtype = None
678 678
679 679 # self.nChannels = 0
680 680
681 681 # self.nHeights = 0
682 682
683 683 self.nProfiles = None
684 684
685 685 self.heightList = None
686 686
687 687 self.channelList = None
688 688
689 689 # self.channelIndexList = None
690 690
691 691 self.flagNoData = True
692 692
693 693 self.flagDiscontinuousBlock = False
694 694
695 695 # self.nPairs = 0
696 696
697 697 self.utctime = None
698 698
699 699 self.blocksize = None
700 700
701 701 self.profileIndex = 0
702 702
703 703 self.nCohInt = 1
704 704
705 705 self.nIncohInt = 1
706 706
707 707 def getNormFactor(self):
708 708 pwcode = 1
709 709 if self.flagDecodeData:
710 710 pwcode = numpy.sum(self.code[0]**2)
711 711
712 712 normFactor = self.nIncohInt*self.nCohInt*pwcode
713 713
714 714 return normFactor
715 715
716 716 def getTimeInterval(self):
717 717
718 718 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
719 719
720 720 return timeInterval
721 721
722 722 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
723 723 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
724 724
725 725 class Fits(JROData):
726 726
727 727 heightList = None
728 728
729 729 channelList = None
730 730
731 731 flagNoData = True
732 732
733 733 flagDiscontinuousBlock = False
734 734
735 735 useLocalTime = False
736 736
737 737 utctime = None
738 738
739 739 timeZone = None
740 740
741 741 # ippSeconds = None
742 742
743 743 # timeInterval = None
744 744
745 745 nCohInt = None
746 746
747 747 nIncohInt = None
748 748
749 749 noise = None
750 750
751 751 windowOfFilter = 1
752 752
753 753 #Speed of ligth
754 754 C = 3e8
755 755
756 756 frequency = 49.92e6
757 757
758 758 realtime = False
759 759
760 760
761 761 def __init__(self):
762 762
763 763 self.type = "Fits"
764 764
765 765 self.nProfiles = None
766 766
767 767 self.heightList = None
768 768
769 769 self.channelList = None
770 770
771 771 # self.channelIndexList = None
772 772
773 773 self.flagNoData = True
774 774
775 775 self.utctime = None
776 776
777 777 self.nCohInt = 1
778 778
779 779 self.nIncohInt = 1
780 780
781 781 self.useLocalTime = True
782 782
783 783 self.profileIndex = 0
784 784
785 785 # self.utctime = None
786 786 # self.timeZone = None
787 787 # self.ltctime = None
788 788 # self.timeInterval = None
789 789 # self.header = None
790 790 # self.data_header = None
791 791 # self.data = None
792 792 # self.datatime = None
793 793 # self.flagNoData = False
794 794 # self.expName = ''
795 795 # self.nChannels = None
796 796 # self.nSamples = None
797 797 # self.dataBlocksPerFile = None
798 798 # self.comments = ''
799 799 #
800 800
801 801
802 802 def getltctime(self):
803 803
804 804 if self.useLocalTime:
805 805 return self.utctime - self.timeZone*60
806 806
807 807 return self.utctime
808 808
809 809 def getDatatime(self):
810 810
811 811 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
812 812 return datatime
813 813
814 814 def getTimeRange(self):
815 815
816 816 datatime = []
817 817
818 818 datatime.append(self.ltctime)
819 819 datatime.append(self.ltctime + self.timeInterval)
820 820
821 821 datatime = numpy.array(datatime)
822 822
823 823 return datatime
824 824
825 825 def getHeiRange(self):
826 826
827 827 heis = self.heightList
828 828
829 829 return heis
830 830
831 831 def getNHeights(self):
832 832
833 833 return len(self.heightList)
834 834
835 835 def getNChannels(self):
836 836
837 837 return len(self.channelList)
838 838
839 839 def getChannelIndexList(self):
840 840
841 841 return range(self.nChannels)
842 842
843 843 def getNoise(self, type = 1):
844 844
845 845 #noise = numpy.zeros(self.nChannels)
846 846
847 847 if type == 1:
848 848 noise = self.getNoisebyHildebrand()
849 849
850 850 if type == 2:
851 851 noise = self.getNoisebySort()
852 852
853 853 if type == 3:
854 854 noise = self.getNoisebyWindow()
855 855
856 856 return noise
857 857
858 858 def getTimeInterval(self):
859 859
860 860 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
861 861
862 862 return timeInterval
863 863
864 864 datatime = property(getDatatime, "I'm the 'datatime' property")
865 865 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
866 866 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
867 867 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
868 868 noise = property(getNoise, "I'm the 'nHeights' property.")
869 869
870 870 ltctime = property(getltctime, "I'm the 'ltctime' property")
871 871 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
872 872
873 873 class Correlation(JROData):
874 874
875 875 noise = None
876 876
877 877 SNR = None
878 878
879 879 pairsAutoCorr = None #Pairs of Autocorrelation
880 880
881 881 #--------------------------------------------------
882 882
883 883 data_corr = None
884 884
885 885 data_volt = None
886 886
887 887 lagT = None # each element value is a profileIndex
888 888
889 889 lagR = None # each element value is in km
890 890
891 891 pairsList = None
892 892
893 893 calculateVelocity = None
894 894
895 895 nPoints = None
896 896
897 897 nAvg = None
898 898
899 899 bufferSize = None
900 900
901 901 def __init__(self):
902 902 '''
903 903 Constructor
904 904 '''
905 905 self.radarControllerHeaderObj = RadarControllerHeader()
906 906
907 907 self.systemHeaderObj = SystemHeader()
908 908
909 909 self.type = "Correlation"
910 910
911 911 self.data = None
912 912
913 913 self.dtype = None
914 914
915 915 self.nProfiles = None
916 916
917 917 self.heightList = None
918 918
919 919 self.channelList = None
920 920
921 921 self.flagNoData = True
922 922
923 923 self.flagDiscontinuousBlock = False
924 924
925 925 self.utctime = None
926 926
927 927 self.timeZone = None
928 928
929 929 self.dstFlag = None
930 930
931 931 self.errorCount = None
932 932
933 933 self.blocksize = None
934 934
935 935 self.flagDecodeData = False #asumo q la data no esta decodificada
936 936
937 937 self.flagDeflipData = False #asumo q la data no esta sin flip
938 938
939 939 self.pairsList = None
940 940
941 941 self.nPoints = None
942 942
943 943 def getLagTRange(self, extrapoints=0):
944 944
945 945 lagTRange = self.lagT
946 946 diff = lagTRange[1] - lagTRange[0]
947 947 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
948 948 lagTRange = numpy.hstack((lagTRange, extra))
949 949
950 950 return lagTRange
951 951
952 952 def getLagRRange(self, extrapoints=0):
953 953
954 954 return self.lagR
955 955
956 956 def getPairsList(self):
957 957
958 958 return self.pairsList
959 959
960 960 def getCalculateVelocity(self):
961 961
962 962 return self.calculateVelocity
963 963
964 964 def getNPoints(self):
965 965
966 966 return self.nPoints
967 967
968 968 def getNAvg(self):
969 969
970 970 return self.nAvg
971 971
972 972 def getBufferSize(self):
973 973
974 974 return self.bufferSize
975 975
976 976 def getPairsAutoCorr(self):
977 977 pairsList = self.pairsList
978 978 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
979 979
980 980 for l in range(len(pairsList)):
981 981 firstChannel = pairsList[l][0]
982 982 secondChannel = pairsList[l][1]
983 983
984 984 #Obteniendo pares de Autocorrelacion
985 985 if firstChannel == secondChannel:
986 986 pairsAutoCorr[firstChannel] = int(l)
987 987
988 988 pairsAutoCorr = pairsAutoCorr.astype(int)
989 989
990 990 return pairsAutoCorr
991 991
992 992 def getNoise(self, mode = 2):
993 993
994 994 indR = numpy.where(self.lagR == 0)[0][0]
995 995 indT = numpy.where(self.lagT == 0)[0][0]
996 996
997 997 jspectra0 = self.data_corr[:,:,indR,:]
998 998 jspectra = copy.copy(jspectra0)
999 999
1000 1000 num_chan = jspectra.shape[0]
1001 1001 num_hei = jspectra.shape[2]
1002 1002
1003 1003 freq_dc = jspectra.shape[1]/2
1004 1004 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1005 1005
1006 1006 if ind_vel[0]<0:
1007 1007 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1008 1008
1009 1009 if mode == 1:
1010 1010 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1011 1011
1012 1012 if mode == 2:
1013 1013
1014 1014 vel = numpy.array([-2,-1,1,2])
1015 1015 xx = numpy.zeros([4,4])
1016 1016
1017 1017 for fil in range(4):
1018 1018 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1019 1019
1020 1020 xx_inv = numpy.linalg.inv(xx)
1021 1021 xx_aux = xx_inv[0,:]
1022 1022
1023 1023 for ich in range(num_chan):
1024 1024 yy = jspectra[ich,ind_vel,:]
1025 1025 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1026 1026
1027 1027 junkid = jspectra[ich,freq_dc,:]<=0
1028 1028 cjunkid = sum(junkid)
1029 1029
1030 1030 if cjunkid.any():
1031 1031 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1032 1032
1033 1033 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1034 1034
1035 1035 return noise
1036 1036
1037 1037 def getTimeInterval(self):
1038 1038
1039 1039 timeInterval = self.ippSeconds * self.nCohInt * self.nPoints
1040 1040
1041 1041 return timeInterval
1042 1042
1043 1043 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1044 1044 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
1045 1045 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
1046 1046 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
1047 1047 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
1048 1048 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
1049 1049
1050 1050
1051 1051 class Parameters(JROData):
1052 1052
1053 1053 #Information from previous data
1054 1054
1055 1055 inputUnit = None #Type of data to be processed
1056 1056
1057 1057 operation = None #Type of operation to parametrize
1058 1058
1059 1059 normFactor = None #Normalization Factor
1060 1060
1061 1061 groupList = None #List of Pairs, Groups, etc
1062 1062
1063 1063 #Parameters
1064 1064
1065 1065 data_param = None #Parameters obtained
1066 1066
1067 1067 data_pre = None #Data Pre Parametrization
1068 1068
1069 1069 data_SNR = None #Signal to Noise Ratio
1070 1070
1071 1071 # heightRange = None #Heights
1072 1072
1073 1073 abscissaList = None #Abscissa, can be velocities, lags or time
1074 1074
1075 1075 noise = None #Noise Potency
1076 1076
1077 1077 utctimeInit = None #Initial UTC time
1078 1078
1079 1079 paramInterval = None #Time interval to calculate Parameters in seconds
1080 1080
1081 1081 #Fitting
1082 1082
1083 1083 data_error = None #Error of the estimation
1084 1084
1085 1085 constants = None
1086 1086
1087 1087 library = None
1088 1088
1089 1089 #Output signal
1090 1090
1091 1091 outputInterval = None #Time interval to calculate output signal in seconds
1092 1092
1093 1093 data_output = None #Out signal
1094 1094
1095 1095
1096 1096
1097 1097 def __init__(self):
1098 1098 '''
1099 1099 Constructor
1100 1100 '''
1101 1101 self.radarControllerHeaderObj = RadarControllerHeader()
1102 1102
1103 1103 self.systemHeaderObj = SystemHeader()
1104 1104
1105 1105 self.type = "Parameters"
1106 1106
1107 1107 def getTimeRange1(self):
1108 1108
1109 1109 datatime = []
1110 1110
1111 1111 if self.useLocalTime:
1112 1112 time1 = self.utctimeInit - self.timeZone*60
1113 1113 else:
1114 1114 time1 = utctimeInit
1115 1115
1116 1116 # datatime.append(self.utctimeInit)
1117 1117 # datatime.append(self.utctimeInit + self.outputInterval)
1118 1118 datatime.append(time1)
1119 1119 datatime.append(time1 + self.outputInterval)
1120 1120
1121 1121 datatime = numpy.array(datatime)
1122 1122
1123 1123 return datatime
Show file before | Show file after | Hide comments
@@ -1,875 +1,878
1 1 import numpy
2 2 import math
3 3
4 4 from jroproc_base import ProcessingUnit, Operation
5 5 from schainpy.model.data.jrodata import Spectra
6 6 from schainpy.model.data.jrodata import hildebrand_sekhon
7 7
8 8 class SpectraProc(ProcessingUnit):
9 9
10 10 def __init__(self):
11 11
12 12 ProcessingUnit.__init__(self)
13 13
14 14 self.buffer = None
15 15 self.firstdatatime = None
16 16 self.profIndex = 0
17 17 self.dataOut = Spectra()
18 18 self.id_min = None
19 19 self.id_max = None
20 20
21 21 def __updateSpecFromVoltage(self):
22 22
23 23 self.dataOut.timeZone = self.dataIn.timeZone
24 24 self.dataOut.dstFlag = self.dataIn.dstFlag
25 25 self.dataOut.errorCount = self.dataIn.errorCount
26 26 self.dataOut.useLocalTime = self.dataIn.useLocalTime
27 27
28 28 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
29 29 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
30 30 self.dataOut.channelList = self.dataIn.channelList
31 31 self.dataOut.heightList = self.dataIn.heightList
32 32 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
33 33
34 34 self.dataOut.nBaud = self.dataIn.nBaud
35 35 self.dataOut.nCode = self.dataIn.nCode
36 36 self.dataOut.code = self.dataIn.code
37 37 self.dataOut.nProfiles = self.dataOut.nFFTPoints
38 38
39 39 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
40 40 self.dataOut.utctime = self.firstdatatime
41 41 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
42 42 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
43 43 self.dataOut.flagShiftFFT = False
44 44
45 45 self.dataOut.nCohInt = self.dataIn.nCohInt
46 46 self.dataOut.nIncohInt = 1
47 47
48 48 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
49 49
50 50 self.dataOut.frequency = self.dataIn.frequency
51 51 self.dataOut.realtime = self.dataIn.realtime
52 52
53 53 self.dataOut.azimuth = self.dataIn.azimuth
54 54 self.dataOut.zenith = self.dataIn.zenith
55 55
56 56 self.dataOut.beam.codeList = self.dataIn.beam.codeList
57 57 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
58 58 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
59 59
60 60 def __getFft(self):
61 61 """
62 62 Convierte valores de Voltaje a Spectra
63 63
64 64 Affected:
65 65 self.dataOut.data_spc
66 66 self.dataOut.data_cspc
67 67 self.dataOut.data_dc
68 68 self.dataOut.heightList
69 69 self.profIndex
70 70 self.buffer
71 71 self.dataOut.flagNoData
72 72 """
73 73 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
74 74 fft_volt = fft_volt.astype(numpy.dtype('complex'))
75 75 dc = fft_volt[:,0,:]
76 76
77 77 #calculo de self-spectra
78 78 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
79 79 spc = fft_volt * numpy.conjugate(fft_volt)
80 80 spc = spc.real
81 81
82 82 blocksize = 0
83 83 blocksize += dc.size
84 84 blocksize += spc.size
85 85
86 86 cspc = None
87 87 pairIndex = 0
88 88 if self.dataOut.pairsList != None:
89 89 #calculo de cross-spectra
90 90 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
91 91 for pair in self.dataOut.pairsList:
92 92 if pair[0] not in self.dataOut.channelList:
93 93 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
94 94 if pair[1] not in self.dataOut.channelList:
95 95 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
96 96
97 97 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
98 98 pairIndex += 1
99 99 blocksize += cspc.size
100 100
101 101 self.dataOut.data_spc = spc
102 102 self.dataOut.data_cspc = cspc
103 103 self.dataOut.data_dc = dc
104 104 self.dataOut.blockSize = blocksize
105 105 self.dataOut.flagShiftFFT = True
106 106
107 107 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None):
108 108
109 109 self.dataOut.flagNoData = True
110 110
111 111 if self.dataIn.type == "Spectra":
112 112 self.dataOut.copy(self.dataIn)
113 113 return True
114 114
115 115 if self.dataIn.type == "Voltage":
116 116
117 117 if nFFTPoints == None:
118 118 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
119 119
120 120 if nProfiles == None:
121 121 nProfiles = nFFTPoints
122 122
123 123 if ippFactor == None:
124 124 ippFactor = 1
125 125
126 126 self.dataOut.ippFactor = ippFactor
127 127
128 128 self.dataOut.nFFTPoints = nFFTPoints
129 129 self.dataOut.pairsList = pairsList
130 130
131 131 if self.buffer is None:
132 132 self.buffer = numpy.zeros( (self.dataIn.nChannels,
133 133 nProfiles,
134 134 self.dataIn.nHeights),
135 135 dtype='complex')
136 136
137 137 if self.dataIn.flagDataAsBlock:
138 #data dimension: [nChannels, nProfiles, nSamples]
139 nVoltProfiles = self.dataIn.data.shape[1]
140 nVoltProfiles = self.dataIn.nProfiles
138 141
139 if self.dataIn.nProfiles == nProfiles:
142 if nVoltProfiles == nProfiles:
140 143 self.buffer = self.dataIn.data.copy()
141 self.profIndex = nProfiles
144 self.profIndex = nVoltProfiles
142 145
143 elif self.dataIn.nProfiles < nProfiles:
146 elif nVoltProfiles < nProfiles:
144 147
145 148 if self.profIndex == 0:
146 149 self.id_min = 0
147 self.id_max = self.dataIn.nProfiles
150 self.id_max = nVoltProfiles
148 151
149 152 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
150 self.profIndex += self.dataIn.nProfiles
151 self.id_min += self.dataIn.data.shape[1]
152 self.id_max += self.dataIn.data.shape[1]
153 self.profIndex += nVoltProfiles
154 self.id_min += nVoltProfiles
155 self.id_max += nVoltProfiles
153 156 else:
154 157 raise ValueError, "The type object %s has %d profiles, it should be equal to %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
155 158 self.dataOut.flagNoData = True
156 159 return 0
157 160 else:
158 161 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
159 162 self.profIndex += 1
160 163
161 164 if self.firstdatatime == None:
162 165 self.firstdatatime = self.dataIn.utctime
163 166
164 167 if self.profIndex == nProfiles:
165 168 self.__updateSpecFromVoltage()
166 169 self.__getFft()
167 170
168 171 self.dataOut.flagNoData = False
169 172 self.firstdatatime = None
170 173 self.profIndex = 0
171 174
172 175 return True
173 176
174 177 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
175 178
176 179 def __selectPairs(self, channelList=None):
177 180
178 181 if channelList == None:
179 182 return
180 183
181 184 pairsIndexListSelected = []
182 185 for pairIndex in self.dataOut.pairsIndexList:
183 186 #First pair
184 187 if self.dataOut.pairsList[pairIndex][0] not in channelList:
185 188 continue
186 189 #Second pair
187 190 if self.dataOut.pairsList[pairIndex][1] not in channelList:
188 191 continue
189 192
190 193 pairsIndexListSelected.append(pairIndex)
191 194
192 195 if not pairsIndexListSelected:
193 196 self.dataOut.data_cspc = None
194 197 self.dataOut.pairsList = []
195 198 return
196 199
197 200 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
198 201 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
199 202
200 203 return
201 204
202 205 def selectChannels(self, channelList):
203 206
204 207 channelIndexList = []
205 208
206 209 for channel in channelList:
207 210 if channel not in self.dataOut.channelList:
208 211 raise ValueError, "Error selecting channels: The value %d in channelList is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
209 212
210 213 index = self.dataOut.channelList.index(channel)
211 214 channelIndexList.append(index)
212 215
213 216 self.selectChannelsByIndex(channelIndexList)
214 217
215 218 def selectChannelsByIndex(self, channelIndexList):
216 219 """
217 220 Selecciona un bloque de datos en base a canales segun el channelIndexList
218 221
219 222 Input:
220 223 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
221 224
222 225 Affected:
223 226 self.dataOut.data_spc
224 227 self.dataOut.channelIndexList
225 228 self.dataOut.nChannels
226 229
227 230 Return:
228 231 None
229 232 """
230 233
231 234 for channelIndex in channelIndexList:
232 235 if channelIndex not in self.dataOut.channelIndexList:
233 236 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
234 237
235 238 # nChannels = len(channelIndexList)
236 239
237 240 data_spc = self.dataOut.data_spc[channelIndexList,:]
238 241 data_dc = self.dataOut.data_dc[channelIndexList,:]
239 242
240 243 self.dataOut.data_spc = data_spc
241 244 self.dataOut.data_dc = data_dc
242 245
243 246 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
244 247 # self.dataOut.nChannels = nChannels
245 248
246 249 self.__selectPairs(self.dataOut.channelList)
247 250
248 251 return 1
249 252
250 253 def selectHeights(self, minHei, maxHei):
251 254 """
252 255 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
253 256 minHei <= height <= maxHei
254 257
255 258 Input:
256 259 minHei : valor minimo de altura a considerar
257 260 maxHei : valor maximo de altura a considerar
258 261
259 262 Affected:
260 263 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
261 264
262 265 Return:
263 266 1 si el metodo se ejecuto con exito caso contrario devuelve 0
264 267 """
265 268
266 269 if (minHei > maxHei):
267 270 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
268 271
269 272 if (minHei < self.dataOut.heightList[0]):
270 273 minHei = self.dataOut.heightList[0]
271 274
272 275 if (maxHei > self.dataOut.heightList[-1]):
273 276 maxHei = self.dataOut.heightList[-1]
274 277
275 278 minIndex = 0
276 279 maxIndex = 0
277 280 heights = self.dataOut.heightList
278 281
279 282 inda = numpy.where(heights >= minHei)
280 283 indb = numpy.where(heights <= maxHei)
281 284
282 285 try:
283 286 minIndex = inda[0][0]
284 287 except:
285 288 minIndex = 0
286 289
287 290 try:
288 291 maxIndex = indb[0][-1]
289 292 except:
290 293 maxIndex = len(heights)
291 294
292 295 self.selectHeightsByIndex(minIndex, maxIndex)
293 296
294 297 return 1
295 298
296 299 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
297 300 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
298 301
299 302 if hei_ref != None:
300 303 newheis = numpy.where(self.dataOut.heightList>hei_ref)
301 304
302 305 minIndex = min(newheis[0])
303 306 maxIndex = max(newheis[0])
304 307 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
305 308 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
306 309
307 310 # determina indices
308 311 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
309 312 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
310 313 beacon_dB = numpy.sort(avg_dB)[-nheis:]
311 314 beacon_heiIndexList = []
312 315 for val in avg_dB.tolist():
313 316 if val >= beacon_dB[0]:
314 317 beacon_heiIndexList.append(avg_dB.tolist().index(val))
315 318
316 319 #data_spc = data_spc[:,:,beacon_heiIndexList]
317 320 data_cspc = None
318 321 if self.dataOut.data_cspc is not None:
319 322 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
320 323 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
321 324
322 325 data_dc = None
323 326 if self.dataOut.data_dc is not None:
324 327 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
325 328 #data_dc = data_dc[:,beacon_heiIndexList]
326 329
327 330 self.dataOut.data_spc = data_spc
328 331 self.dataOut.data_cspc = data_cspc
329 332 self.dataOut.data_dc = data_dc
330 333 self.dataOut.heightList = heightList
331 334 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
332 335
333 336 return 1
334 337
335 338
336 339 def selectHeightsByIndex(self, minIndex, maxIndex):
337 340 """
338 341 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
339 342 minIndex <= index <= maxIndex
340 343
341 344 Input:
342 345 minIndex : valor de indice minimo de altura a considerar
343 346 maxIndex : valor de indice maximo de altura a considerar
344 347
345 348 Affected:
346 349 self.dataOut.data_spc
347 350 self.dataOut.data_cspc
348 351 self.dataOut.data_dc
349 352 self.dataOut.heightList
350 353
351 354 Return:
352 355 1 si el metodo se ejecuto con exito caso contrario devuelve 0
353 356 """
354 357
355 358 if (minIndex < 0) or (minIndex > maxIndex):
356 359 raise ValueError, "Error selecting heights by index: Index range in (%d,%d) is not valid" % (minIndex, maxIndex)
357 360
358 361 if (maxIndex >= self.dataOut.nHeights):
359 362 maxIndex = self.dataOut.nHeights-1
360 363
361 364 #Spectra
362 365 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
363 366
364 367 data_cspc = None
365 368 if self.dataOut.data_cspc is not None:
366 369 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
367 370
368 371 data_dc = None
369 372 if self.dataOut.data_dc is not None:
370 373 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
371 374
372 375 self.dataOut.data_spc = data_spc
373 376 self.dataOut.data_cspc = data_cspc
374 377 self.dataOut.data_dc = data_dc
375 378
376 379 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
377 380
378 381 return 1
379 382
380 383 def removeDC(self, mode = 2):
381 384 jspectra = self.dataOut.data_spc
382 385 jcspectra = self.dataOut.data_cspc
383 386
384 387
385 388 num_chan = jspectra.shape[0]
386 389 num_hei = jspectra.shape[2]
387 390
388 391 if jcspectra is not None:
389 392 jcspectraExist = True
390 393 num_pairs = jcspectra.shape[0]
391 394 else: jcspectraExist = False
392 395
393 396 freq_dc = jspectra.shape[1]/2
394 397 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
395 398
396 399 if ind_vel[0]<0:
397 400 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
398 401
399 402 if mode == 1:
400 403 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
401 404
402 405 if jcspectraExist:
403 406 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
404 407
405 408 if mode == 2:
406 409
407 410 vel = numpy.array([-2,-1,1,2])
408 411 xx = numpy.zeros([4,4])
409 412
410 413 for fil in range(4):
411 414 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
412 415
413 416 xx_inv = numpy.linalg.inv(xx)
414 417 xx_aux = xx_inv[0,:]
415 418
416 419 for ich in range(num_chan):
417 420 yy = jspectra[ich,ind_vel,:]
418 421 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
419 422
420 423 junkid = jspectra[ich,freq_dc,:]<=0
421 424 cjunkid = sum(junkid)
422 425
423 426 if cjunkid.any():
424 427 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
425 428
426 429 if jcspectraExist:
427 430 for ip in range(num_pairs):
428 431 yy = jcspectra[ip,ind_vel,:]
429 432 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
430 433
431 434
432 435 self.dataOut.data_spc = jspectra
433 436 self.dataOut.data_cspc = jcspectra
434 437
435 438 return 1
436 439
437 440 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
438 441
439 442 jspectra = self.dataOut.data_spc
440 443 jcspectra = self.dataOut.data_cspc
441 444 jnoise = self.dataOut.getNoise()
442 445 num_incoh = self.dataOut.nIncohInt
443 446
444 447 num_channel = jspectra.shape[0]
445 448 num_prof = jspectra.shape[1]
446 449 num_hei = jspectra.shape[2]
447 450
448 451 #hei_interf
449 452 if hei_interf is None:
450 453 count_hei = num_hei/2 #Como es entero no importa
451 454 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
452 455 hei_interf = numpy.asarray(hei_interf)[0]
453 456 #nhei_interf
454 457 if (nhei_interf == None):
455 458 nhei_interf = 5
456 459 if (nhei_interf < 1):
457 460 nhei_interf = 1
458 461 if (nhei_interf > count_hei):
459 462 nhei_interf = count_hei
460 463 if (offhei_interf == None):
461 464 offhei_interf = 0
462 465
463 466 ind_hei = range(num_hei)
464 467 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
465 468 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
466 469 mask_prof = numpy.asarray(range(num_prof))
467 470 num_mask_prof = mask_prof.size
468 471 comp_mask_prof = [0, num_prof/2]
469 472
470 473
471 474 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
472 475 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
473 476 jnoise = numpy.nan
474 477 noise_exist = jnoise[0] < numpy.Inf
475 478
476 479 #Subrutina de Remocion de la Interferencia
477 480 for ich in range(num_channel):
478 481 #Se ordena los espectros segun su potencia (menor a mayor)
479 482 power = jspectra[ich,mask_prof,:]
480 483 power = power[:,hei_interf]
481 484 power = power.sum(axis = 0)
482 485 psort = power.ravel().argsort()
483 486
484 487 #Se estima la interferencia promedio en los Espectros de Potencia empleando
485 488 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
486 489
487 490 if noise_exist:
488 491 # tmp_noise = jnoise[ich] / num_prof
489 492 tmp_noise = jnoise[ich]
490 493 junkspc_interf = junkspc_interf - tmp_noise
491 494 #junkspc_interf[:,comp_mask_prof] = 0
492 495
493 496 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
494 497 jspc_interf = jspc_interf.transpose()
495 498 #Calculando el espectro de interferencia promedio
496 499 noiseid = numpy.where(jspc_interf <= tmp_noise/ math.sqrt(num_incoh))
497 500 noiseid = noiseid[0]
498 501 cnoiseid = noiseid.size
499 502 interfid = numpy.where(jspc_interf > tmp_noise/ math.sqrt(num_incoh))
500 503 interfid = interfid[0]
501 504 cinterfid = interfid.size
502 505
503 506 if (cnoiseid > 0): jspc_interf[noiseid] = 0
504 507
505 508 #Expandiendo los perfiles a limpiar
506 509 if (cinterfid > 0):
507 510 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
508 511 new_interfid = numpy.asarray(new_interfid)
509 512 new_interfid = {x for x in new_interfid}
510 513 new_interfid = numpy.array(list(new_interfid))
511 514 new_cinterfid = new_interfid.size
512 515 else: new_cinterfid = 0
513 516
514 517 for ip in range(new_cinterfid):
515 518 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
516 519 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
517 520
518 521
519 522 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
520 523
521 524 #Removiendo la interferencia del punto de mayor interferencia
522 525 ListAux = jspc_interf[mask_prof].tolist()
523 526 maxid = ListAux.index(max(ListAux))
524 527
525 528
526 529 if cinterfid > 0:
527 530 for ip in range(cinterfid*(interf == 2) - 1):
528 531 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/math.sqrt(num_incoh))).nonzero()
529 532 cind = len(ind)
530 533
531 534 if (cind > 0):
532 535 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/math.sqrt(num_incoh))
533 536
534 537 ind = numpy.array([-2,-1,1,2])
535 538 xx = numpy.zeros([4,4])
536 539
537 540 for id1 in range(4):
538 541 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
539 542
540 543 xx_inv = numpy.linalg.inv(xx)
541 544 xx = xx_inv[:,0]
542 545 ind = (ind + maxid + num_mask_prof)%num_mask_prof
543 546 yy = jspectra[ich,mask_prof[ind],:]
544 547 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
545 548
546 549
547 550 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/math.sqrt(num_incoh))).nonzero()
548 551 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/math.sqrt(num_incoh))
549 552
550 553 #Remocion de Interferencia en el Cross Spectra
551 554 if jcspectra is None: return jspectra, jcspectra
552 555 num_pairs = jcspectra.size/(num_prof*num_hei)
553 556 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
554 557
555 558 for ip in range(num_pairs):
556 559
557 560 #-------------------------------------------
558 561
559 562 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
560 563 cspower = cspower[:,hei_interf]
561 564 cspower = cspower.sum(axis = 0)
562 565
563 566 cspsort = cspower.ravel().argsort()
564 567 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
565 568 junkcspc_interf = junkcspc_interf.transpose()
566 569 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
567 570
568 571 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
569 572
570 573 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
571 574 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
572 575 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
573 576
574 577 for iprof in range(num_prof):
575 578 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
576 579 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
577 580
578 581 #Removiendo la Interferencia
579 582 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
580 583
581 584 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
582 585 maxid = ListAux.index(max(ListAux))
583 586
584 587 ind = numpy.array([-2,-1,1,2])
585 588 xx = numpy.zeros([4,4])
586 589
587 590 for id1 in range(4):
588 591 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
589 592
590 593 xx_inv = numpy.linalg.inv(xx)
591 594 xx = xx_inv[:,0]
592 595
593 596 ind = (ind + maxid + num_mask_prof)%num_mask_prof
594 597 yy = jcspectra[ip,mask_prof[ind],:]
595 598 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
596 599
597 600 #Guardar Resultados
598 601 self.dataOut.data_spc = jspectra
599 602 self.dataOut.data_cspc = jcspectra
600 603
601 604 return 1
602 605
603 606 def setRadarFrequency(self, frequency=None):
604 607
605 608 if frequency != None:
606 609 self.dataOut.frequency = frequency
607 610
608 611 return 1
609 612
610 613 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
611 614 #validacion de rango
612 615 if minHei == None:
613 616 minHei = self.dataOut.heightList[0]
614 617
615 618 if maxHei == None:
616 619 maxHei = self.dataOut.heightList[-1]
617 620
618 621 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
619 622 print 'minHei: %.2f is out of the heights range'%(minHei)
620 623 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
621 624 minHei = self.dataOut.heightList[0]
622 625
623 626 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
624 627 print 'maxHei: %.2f is out of the heights range'%(maxHei)
625 628 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
626 629 maxHei = self.dataOut.heightList[-1]
627 630
628 631 # validacion de velocidades
629 632 velrange = self.dataOut.getVelRange(1)
630 633
631 634 if minVel == None:
632 635 minVel = velrange[0]
633 636
634 637 if maxVel == None:
635 638 maxVel = velrange[-1]
636 639
637 640 if (minVel < velrange[0]) or (minVel > maxVel):
638 641 print 'minVel: %.2f is out of the velocity range'%(minVel)
639 642 print 'minVel is setting to %.2f'%(velrange[0])
640 643 minVel = velrange[0]
641 644
642 645 if (maxVel > velrange[-1]) or (maxVel < minVel):
643 646 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
644 647 print 'maxVel is setting to %.2f'%(velrange[-1])
645 648 maxVel = velrange[-1]
646 649
647 650 # seleccion de indices para rango
648 651 minIndex = 0
649 652 maxIndex = 0
650 653 heights = self.dataOut.heightList
651 654
652 655 inda = numpy.where(heights >= minHei)
653 656 indb = numpy.where(heights <= maxHei)
654 657
655 658 try:
656 659 minIndex = inda[0][0]
657 660 except:
658 661 minIndex = 0
659 662
660 663 try:
661 664 maxIndex = indb[0][-1]
662 665 except:
663 666 maxIndex = len(heights)
664 667
665 668 if (minIndex < 0) or (minIndex > maxIndex):
666 669 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
667 670
668 671 if (maxIndex >= self.dataOut.nHeights):
669 672 maxIndex = self.dataOut.nHeights-1
670 673
671 674 # seleccion de indices para velocidades
672 675 indminvel = numpy.where(velrange >= minVel)
673 676 indmaxvel = numpy.where(velrange <= maxVel)
674 677 try:
675 678 minIndexVel = indminvel[0][0]
676 679 except:
677 680 minIndexVel = 0
678 681
679 682 try:
680 683 maxIndexVel = indmaxvel[0][-1]
681 684 except:
682 685 maxIndexVel = len(velrange)
683 686
684 687 #seleccion del espectro
685 688 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
686 689 #estimacion de ruido
687 690 noise = numpy.zeros(self.dataOut.nChannels)
688 691
689 692 for channel in range(self.dataOut.nChannels):
690 693 daux = data_spc[channel,:,:]
691 694 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
692 695
693 696 self.dataOut.noise_estimation = noise.copy()
694 697
695 698 return 1
696 699
697 700 class IncohInt(Operation):
698 701
699 702
700 703 __profIndex = 0
701 704 __withOverapping = False
702 705
703 706 __byTime = False
704 707 __initime = None
705 708 __lastdatatime = None
706 709 __integrationtime = None
707 710
708 711 __buffer_spc = None
709 712 __buffer_cspc = None
710 713 __buffer_dc = None
711 714
712 715 __dataReady = False
713 716
714 717 __timeInterval = None
715 718
716 719 n = None
717 720
718 721
719 722
720 723 def __init__(self):
721 724
722 725 Operation.__init__(self)
723 726 # self.isConfig = False
724 727
725 728 def setup(self, n=None, timeInterval=None, overlapping=False):
726 729 """
727 730 Set the parameters of the integration class.
728 731
729 732 Inputs:
730 733
731 734 n : Number of coherent integrations
732 735 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
733 736 overlapping :
734 737
735 738 """
736 739
737 740 self.__initime = None
738 741 self.__lastdatatime = 0
739 742
740 743 self.__buffer_spc = 0
741 744 self.__buffer_cspc = 0
742 745 self.__buffer_dc = 0
743 746
744 747 self.__profIndex = 0
745 748 self.__dataReady = False
746 749 self.__byTime = False
747 750
748 751 if n is None and timeInterval is None:
749 752 raise ValueError, "n or timeInterval should be specified ..."
750 753
751 754 if n is not None:
752 755 self.n = int(n)
753 756 else:
754 757 self.__integrationtime = int(timeInterval) #if (type(timeInterval)!=integer) -> change this line
755 758 self.n = None
756 759 self.__byTime = True
757 760
758 761 def putData(self, data_spc, data_cspc, data_dc):
759 762
760 763 """
761 764 Add a profile to the __buffer_spc and increase in one the __profileIndex
762 765
763 766 """
764 767
765 768 self.__buffer_spc += data_spc
766 769
767 770 if data_cspc is None:
768 771 self.__buffer_cspc = None
769 772 else:
770 773 self.__buffer_cspc += data_cspc
771 774
772 775 if data_dc is None:
773 776 self.__buffer_dc = None
774 777 else:
775 778 self.__buffer_dc += data_dc
776 779
777 780 self.__profIndex += 1
778 781
779 782 return
780 783
781 784 def pushData(self):
782 785 """
783 786 Return the sum of the last profiles and the profiles used in the sum.
784 787
785 788 Affected:
786 789
787 790 self.__profileIndex
788 791
789 792 """
790 793
791 794 data_spc = self.__buffer_spc
792 795 data_cspc = self.__buffer_cspc
793 796 data_dc = self.__buffer_dc
794 797 n = self.__profIndex
795 798
796 799 self.__buffer_spc = 0
797 800 self.__buffer_cspc = 0
798 801 self.__buffer_dc = 0
799 802 self.__profIndex = 0
800 803
801 804 return data_spc, data_cspc, data_dc, n
802 805
803 806 def byProfiles(self, *args):
804 807
805 808 self.__dataReady = False
806 809 avgdata_spc = None
807 810 avgdata_cspc = None
808 811 avgdata_dc = None
809 812
810 813 self.putData(*args)
811 814
812 815 if self.__profIndex == self.n:
813 816
814 817 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
815 818 self.n = n
816 819 self.__dataReady = True
817 820
818 821 return avgdata_spc, avgdata_cspc, avgdata_dc
819 822
820 823 def byTime(self, datatime, *args):
821 824
822 825 self.__dataReady = False
823 826 avgdata_spc = None
824 827 avgdata_cspc = None
825 828 avgdata_dc = None
826 829
827 830 self.putData(*args)
828 831
829 832 if (datatime - self.__initime) >= self.__integrationtime:
830 833 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
831 834 self.n = n
832 835 self.__dataReady = True
833 836
834 837 return avgdata_spc, avgdata_cspc, avgdata_dc
835 838
836 839 def integrate(self, datatime, *args):
837 840
838 841 if self.__profIndex == 0:
839 842 self.__initime = datatime
840 843
841 844 if self.__byTime:
842 845 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
843 846 else:
844 847 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
845 848
846 849 if not self.__dataReady:
847 850 return None, None, None, None
848 851
849 852 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
850 853
851 854 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
852 855
853 856 if n==1:
854 857 return
855 858
856 859 dataOut.flagNoData = True
857 860
858 861 if not self.isConfig:
859 862 self.setup(n, timeInterval, overlapping)
860 863 self.isConfig = True
861 864
862 865 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
863 866 dataOut.data_spc,
864 867 dataOut.data_cspc,
865 868 dataOut.data_dc)
866 869
867 870 if self.__dataReady:
868 871
869 872 dataOut.data_spc = avgdata_spc
870 873 dataOut.data_cspc = avgdata_cspc
871 874 dataOut.data_dc = avgdata_dc
872 875
873 876 dataOut.nIncohInt *= self.n
874 877 dataOut.utctime = avgdatatime
875 878 dataOut.flagNoData = False
Show file before | Show file after | Hide comments
@@ -1,1071 +1,1072
1 1 import numpy
2 2
3 3 from jroproc_base import ProcessingUnit, Operation
4 4 from schainpy.model.data.jrodata import Voltage
5 5
6 6 class VoltageProc(ProcessingUnit):
7 7
8 8
9 9 def __init__(self):
10 10
11 11 ProcessingUnit.__init__(self)
12 12
13 13 # self.objectDict = {}
14 14 self.dataOut = Voltage()
15 15 self.flip = 1
16 16
17 17 def run(self):
18 18 if self.dataIn.type == 'AMISR':
19 19 self.__updateObjFromAmisrInput()
20 20
21 21 if self.dataIn.type == 'Voltage':
22 22 self.dataOut.copy(self.dataIn)
23 23
24 24 # self.dataOut.copy(self.dataIn)
25 25
26 26 def __updateObjFromAmisrInput(self):
27 27
28 28 self.dataOut.timeZone = self.dataIn.timeZone
29 29 self.dataOut.dstFlag = self.dataIn.dstFlag
30 30 self.dataOut.errorCount = self.dataIn.errorCount
31 31 self.dataOut.useLocalTime = self.dataIn.useLocalTime
32 32
33 33 self.dataOut.flagNoData = self.dataIn.flagNoData
34 34 self.dataOut.data = self.dataIn.data
35 35 self.dataOut.utctime = self.dataIn.utctime
36 36 self.dataOut.channelList = self.dataIn.channelList
37 37 # self.dataOut.timeInterval = self.dataIn.timeInterval
38 38 self.dataOut.heightList = self.dataIn.heightList
39 39 self.dataOut.nProfiles = self.dataIn.nProfiles
40 40
41 41 self.dataOut.nCohInt = self.dataIn.nCohInt
42 42 self.dataOut.ippSeconds = self.dataIn.ippSeconds
43 43 self.dataOut.frequency = self.dataIn.frequency
44 44
45 45 self.dataOut.azimuth = self.dataIn.azimuth
46 46 self.dataOut.zenith = self.dataIn.zenith
47 47
48 48 self.dataOut.beam.codeList = self.dataIn.beam.codeList
49 49 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
50 50 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
51 51 #
52 52 # pass#
53 53 #
54 54 # def init(self):
55 55 #
56 56 #
57 57 # if self.dataIn.type == 'AMISR':
58 58 # self.__updateObjFromAmisrInput()
59 59 #
60 60 # if self.dataIn.type == 'Voltage':
61 61 # self.dataOut.copy(self.dataIn)
62 62 # # No necesita copiar en cada init() los atributos de dataIn
63 63 # # la copia deberia hacerse por cada nuevo bloque de datos
64 64
65 65 def selectChannels(self, channelList):
66 66
67 67 channelIndexList = []
68 68
69 69 for channel in channelList:
70 70 if channel not in self.dataOut.channelList:
71 71 raise ValueError, "Channel %d is not in %s" %(channel, str(self.dataOut.channelList))
72 72
73 73 index = self.dataOut.channelList.index(channel)
74 74 channelIndexList.append(index)
75 75
76 76 self.selectChannelsByIndex(channelIndexList)
77 77
78 78 def selectChannelsByIndex(self, channelIndexList):
79 79 """
80 80 Selecciona un bloque de datos en base a canales segun el channelIndexList
81 81
82 82 Input:
83 83 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
84 84
85 85 Affected:
86 86 self.dataOut.data
87 87 self.dataOut.channelIndexList
88 88 self.dataOut.nChannels
89 89 self.dataOut.m_ProcessingHeader.totalSpectra
90 90 self.dataOut.systemHeaderObj.numChannels
91 91 self.dataOut.m_ProcessingHeader.blockSize
92 92
93 93 Return:
94 94 None
95 95 """
96 96
97 97 for channelIndex in channelIndexList:
98 98 if channelIndex not in self.dataOut.channelIndexList:
99 99 print channelIndexList
100 100 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
101 101
102 102 if self.dataOut.flagDataAsBlock:
103 103 """
104 104 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
105 105 """
106 106 data = self.dataOut.data[channelIndexList,:,:]
107 107 else:
108 108 data = self.dataOut.data[channelIndexList,:]
109 109
110 110 self.dataOut.data = data
111 111 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
112 112 # self.dataOut.nChannels = nChannels
113 113
114 114 return 1
115 115
116 116 def selectHeights(self, minHei=None, maxHei=None):
117 117 """
118 118 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
119 119 minHei <= height <= maxHei
120 120
121 121 Input:
122 122 minHei : valor minimo de altura a considerar
123 123 maxHei : valor maximo de altura a considerar
124 124
125 125 Affected:
126 126 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
127 127
128 128 Return:
129 129 1 si el metodo se ejecuto con exito caso contrario devuelve 0
130 130 """
131 131
132 132 if minHei == None:
133 133 minHei = self.dataOut.heightList[0]
134 134
135 135 if maxHei == None:
136 136 maxHei = self.dataOut.heightList[-1]
137 137
138 138 if (minHei < self.dataOut.heightList[0]):
139 139 minHei = self.dataOut.heightList[0]
140 140
141 141 if (maxHei > self.dataOut.heightList[-1]):
142 142 maxHei = self.dataOut.heightList[-1]
143 143
144 144 minIndex = 0
145 145 maxIndex = 0
146 146 heights = self.dataOut.heightList
147 147
148 148 inda = numpy.where(heights >= minHei)
149 149 indb = numpy.where(heights <= maxHei)
150 150
151 151 try:
152 152 minIndex = inda[0][0]
153 153 except:
154 154 minIndex = 0
155 155
156 156 try:
157 157 maxIndex = indb[0][-1]
158 158 except:
159 159 maxIndex = len(heights)
160 160
161 161 self.selectHeightsByIndex(minIndex, maxIndex)
162 162
163 163 return 1
164 164
165 165
166 166 def selectHeightsByIndex(self, minIndex, maxIndex):
167 167 """
168 168 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
169 169 minIndex <= index <= maxIndex
170 170
171 171 Input:
172 172 minIndex : valor de indice minimo de altura a considerar
173 173 maxIndex : valor de indice maximo de altura a considerar
174 174
175 175 Affected:
176 176 self.dataOut.data
177 177 self.dataOut.heightList
178 178
179 179 Return:
180 180 1 si el metodo se ejecuto con exito caso contrario devuelve 0
181 181 """
182 182
183 183 if (minIndex < 0) or (minIndex > maxIndex):
184 184 raise ValueError, "Height index range (%d,%d) is not valid" % (minIndex, maxIndex)
185 185
186 186 if (maxIndex >= self.dataOut.nHeights):
187 187 maxIndex = self.dataOut.nHeights
188 188
189 189 #voltage
190 190 if self.dataOut.flagDataAsBlock:
191 191 """
192 192 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
193 193 """
194 194 data = self.dataOut.data[:,:, minIndex:maxIndex]
195 195 else:
196 196 data = self.dataOut.data[:, minIndex:maxIndex]
197 197
198 198 # firstHeight = self.dataOut.heightList[minIndex]
199 199
200 200 self.dataOut.data = data
201 201 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
202 202
203 203 if self.dataOut.nHeights <= 1:
204 204 raise ValueError, "selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights)
205 205
206 206 return 1
207 207
208 208
209 209 def filterByHeights(self, window):
210 210
211 211 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
212 212
213 213 if window == None:
214 214 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
215 215
216 216 newdelta = deltaHeight * window
217 217 r = self.dataOut.nHeights % window
218 218 newheights = (self.dataOut.nHeights-r)/window
219 219
220 220 if newheights <= 1:
221 221 raise ValueError, "filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window)
222 222
223 223 if self.dataOut.flagDataAsBlock:
224 224 """
225 225 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
226 226 """
227 227 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
228 228 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
229 229 buffer = numpy.sum(buffer,3)
230 230
231 231 else:
232 232 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
233 233 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
234 234 buffer = numpy.sum(buffer,2)
235 235
236 236 self.dataOut.data = buffer
237 237 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
238 238 self.dataOut.windowOfFilter = window
239 239
240 240 def setH0(self, h0, deltaHeight = None):
241 241
242 242 if not deltaHeight:
243 243 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
244 244
245 245 nHeights = self.dataOut.nHeights
246 246
247 247 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
248 248
249 249 self.dataOut.heightList = newHeiRange
250 250
251 251 def deFlip(self, channelList = []):
252 252
253 253 data = self.dataOut.data.copy()
254 254
255 255 if self.dataOut.flagDataAsBlock:
256 256 flip = self.flip
257 257 profileList = range(self.dataOut.nProfiles)
258 258
259 259 if not channelList:
260 260 for thisProfile in profileList:
261 261 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
262 262 flip *= -1.0
263 263 else:
264 264 for thisChannel in channelList:
265 265 if thisChannel not in self.dataOut.channelList:
266 266 continue
267 267
268 268 for thisProfile in profileList:
269 269 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
270 270 flip *= -1.0
271 271
272 272 self.flip = flip
273 273
274 274 else:
275 275 if not channelList:
276 276 data[:,:] = data[:,:]*self.flip
277 277 else:
278 278 for thisChannel in channelList:
279 279 if thisChannel not in self.dataOut.channelList:
280 280 continue
281 281
282 282 data[thisChannel,:] = data[thisChannel,:]*self.flip
283 283
284 284 self.flip *= -1.
285 285
286 286 self.dataOut.data = data
287 287
288 288 def setRadarFrequency(self, frequency=None):
289 289
290 290 if frequency != None:
291 291 self.dataOut.frequency = frequency
292 292
293 293 return 1
294 294
295 295 class CohInt(Operation):
296 296
297 297 isConfig = False
298 298
299 299 __profIndex = 0
300 300 __withOverapping = False
301 301
302 302 __byTime = False
303 303 __initime = None
304 304 __lastdatatime = None
305 305 __integrationtime = None
306 306
307 307 __buffer = None
308 308
309 309 __dataReady = False
310 310
311 311 n = None
312 312
313 313
314 314 def __init__(self):
315 315
316 316 Operation.__init__(self)
317 317
318 318 # self.isConfig = False
319 319
320 320 def setup(self, n=None, timeInterval=None, overlapping=False, byblock=False):
321 321 """
322 322 Set the parameters of the integration class.
323 323
324 324 Inputs:
325 325
326 326 n : Number of coherent integrations
327 327 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
328 328 overlapping :
329 329
330 330 """
331 331
332 332 self.__initime = None
333 333 self.__lastdatatime = 0
334 334 self.__buffer = None
335 335 self.__dataReady = False
336 336 self.byblock = byblock
337 337
338 338 if n == None and timeInterval == None:
339 339 raise ValueError, "n or timeInterval should be specified ..."
340 340
341 341 if n != None:
342 342 self.n = n
343 343 self.__byTime = False
344 344 else:
345 345 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
346 346 self.n = 9999
347 347 self.__byTime = True
348 348
349 349 if overlapping:
350 350 self.__withOverapping = True
351 351 self.__buffer = None
352 352 else:
353 353 self.__withOverapping = False
354 354 self.__buffer = 0
355 355
356 356 self.__profIndex = 0
357 357
358 358 def putData(self, data):
359 359
360 360 """
361 361 Add a profile to the __buffer and increase in one the __profileIndex
362 362
363 363 """
364 364
365 365 if not self.__withOverapping:
366 366 self.__buffer += data.copy()
367 367 self.__profIndex += 1
368 368 return
369 369
370 370 #Overlapping data
371 371 nChannels, nHeis = data.shape
372 372 data = numpy.reshape(data, (1, nChannels, nHeis))
373 373
374 374 #If the buffer is empty then it takes the data value
375 375 if self.__buffer is None:
376 376 self.__buffer = data
377 377 self.__profIndex += 1
378 378 return
379 379
380 380 #If the buffer length is lower than n then stakcing the data value
381 381 if self.__profIndex < self.n:
382 382 self.__buffer = numpy.vstack((self.__buffer, data))
383 383 self.__profIndex += 1
384 384 return
385 385
386 386 #If the buffer length is equal to n then replacing the last buffer value with the data value
387 387 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
388 388 self.__buffer[self.n-1] = data
389 389 self.__profIndex = self.n
390 390 return
391 391
392 392
393 393 def pushData(self):
394 394 """
395 395 Return the sum of the last profiles and the profiles used in the sum.
396 396
397 397 Affected:
398 398
399 399 self.__profileIndex
400 400
401 401 """
402 402
403 403 if not self.__withOverapping:
404 404 data = self.__buffer
405 405 n = self.__profIndex
406 406
407 407 self.__buffer = 0
408 408 self.__profIndex = 0
409 409
410 410 return data, n
411 411
412 412 #Integration with Overlapping
413 413 data = numpy.sum(self.__buffer, axis=0)
414 414 n = self.__profIndex
415 415
416 416 return data, n
417 417
418 418 def byProfiles(self, data):
419 419
420 420 self.__dataReady = False
421 421 avgdata = None
422 422 # n = None
423 423
424 424 self.putData(data)
425 425
426 426 if self.__profIndex == self.n:
427 427
428 428 avgdata, n = self.pushData()
429 429 self.__dataReady = True
430 430
431 431 return avgdata
432 432
433 433 def byTime(self, data, datatime):
434 434
435 435 self.__dataReady = False
436 436 avgdata = None
437 437 n = None
438 438
439 439 self.putData(data)
440 440
441 441 if (datatime - self.__initime) >= self.__integrationtime:
442 442 avgdata, n = self.pushData()
443 443 self.n = n
444 444 self.__dataReady = True
445 445
446 446 return avgdata
447 447
448 448 def integrate(self, data, datatime=None):
449 449
450 450 if self.__initime == None:
451 451 self.__initime = datatime
452 452
453 453 if self.__byTime:
454 454 avgdata = self.byTime(data, datatime)
455 455 else:
456 456 avgdata = self.byProfiles(data)
457 457
458 458
459 459 self.__lastdatatime = datatime
460 460
461 461 if avgdata is None:
462 462 return None, None
463 463
464 464 avgdatatime = self.__initime
465 465
466 466 deltatime = datatime -self.__lastdatatime
467 467
468 468 if not self.__withOverapping:
469 469 self.__initime = datatime
470 470 else:
471 471 self.__initime += deltatime
472 472
473 473 return avgdata, avgdatatime
474 474
475 475 def integrateByBlock(self, dataOut):
476 476
477 477 times = int(dataOut.data.shape[1]/self.n)
478 478 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
479 479
480 480 id_min = 0
481 481 id_max = self.n
482 482
483 483 for i in range(times):
484 484 junk = dataOut.data[:,id_min:id_max,:]
485 485 avgdata[:,i,:] = junk.sum(axis=1)
486 486 id_min += self.n
487 487 id_max += self.n
488 488
489 489 timeInterval = dataOut.ippSeconds*self.n
490 490 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
491 491 self.__dataReady = True
492 492 return avgdata, avgdatatime
493 493
494 494 def run(self, dataOut, **kwargs):
495 495
496 496 if not self.isConfig:
497 497 self.setup(**kwargs)
498 498 self.isConfig = True
499 499
500 500 if dataOut.flagDataAsBlock:
501 501 """
502 502 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
503 503 """
504 504 avgdata, avgdatatime = self.integrateByBlock(dataOut)
505 dataOut.nProfiles /= self.n
505 506 else:
506 507 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
507 508
508 509 # dataOut.timeInterval *= n
509 510 dataOut.flagNoData = True
510 511
511 512 if self.__dataReady:
512 513 dataOut.data = avgdata
513 514 dataOut.nCohInt *= self.n
514 515 dataOut.utctime = avgdatatime
515 516 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
516 517 dataOut.flagNoData = False
517 518
518 519 class Decoder(Operation):
519 520
520 521 isConfig = False
521 522 __profIndex = 0
522 523
523 524 code = None
524 525
525 526 nCode = None
526 527 nBaud = None
527 528
528 529
529 530 def __init__(self):
530 531
531 532 Operation.__init__(self)
532 533
533 534 self.times = None
534 535 self.osamp = None
535 536 # self.__setValues = False
536 537 self.isConfig = False
537 538
538 539 def setup(self, code, osamp, dataOut):
539 540
540 541 self.__profIndex = 0
541 542
542 543 self.code = code
543 544
544 545 self.nCode = len(code)
545 546 self.nBaud = len(code[0])
546 547
547 548 if (osamp != None) and (osamp >1):
548 549 self.osamp = osamp
549 550 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
550 551 self.nBaud = self.nBaud*self.osamp
551 552
552 553 self.__nChannels = dataOut.nChannels
553 554 self.__nProfiles = dataOut.nProfiles
554 555 self.__nHeis = dataOut.nHeights
555 556
556 557 if self.__nHeis < self.nBaud:
557 558 raise ValueError, 'Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)
558 559
559 560 #Frequency
560 561 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
561 562
562 563 __codeBuffer[:,0:self.nBaud] = self.code
563 564
564 565 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
565 566
566 567 if dataOut.flagDataAsBlock:
567 568
568 569 self.ndatadec = self.__nHeis #- self.nBaud + 1
569 570
570 571 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
571 572
572 573 else:
573 574
574 575 #Time
575 576 self.ndatadec = self.__nHeis #- self.nBaud + 1
576 577
577 578 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
578 579
579 580 def __convolutionInFreq(self, data):
580 581
581 582 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
582 583
583 584 fft_data = numpy.fft.fft(data, axis=1)
584 585
585 586 conv = fft_data*fft_code
586 587
587 588 data = numpy.fft.ifft(conv,axis=1)
588 589
589 590 return data
590 591
591 592 def __convolutionInFreqOpt(self, data):
592 593
593 594 raise NotImplementedError
594 595
595 596 def __convolutionInTime(self, data):
596 597
597 598 code = self.code[self.__profIndex]
598 599
599 600 for i in range(self.__nChannels):
600 601 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
601 602
602 603 return self.datadecTime
603 604
604 605 def __convolutionByBlockInTime(self, data):
605 606
606 607 repetitions = self.__nProfiles / self.nCode
607 608
608 609 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
609 610 junk = junk.flatten()
610 611 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
611 612
612 613 for i in range(self.__nChannels):
613 614 for j in range(self.__nProfiles):
614 615 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
615 616
616 617 return self.datadecTime
617 618
618 619 def __convolutionByBlockInFreq(self, data):
619 620
620 621 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
621 622
622 623 fft_data = numpy.fft.fft(data, axis=2)
623 624
624 625 conv = fft_data*fft_code
625 626
626 627 data = numpy.fft.ifft(conv,axis=2)
627 628
628 629 return data
629 630
630 631 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
631 632
632 633 if dataOut.flagDecodeData:
633 634 print "This data is already decoded, recoding again ..."
634 635
635 636 if not self.isConfig:
636 637
637 638 if code is None:
638 639 if dataOut.code is None:
639 640 raise ValueError, "Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type
640 641
641 642 code = dataOut.code
642 643 else:
643 644 code = numpy.array(code).reshape(nCode,nBaud)
644 645
645 646 self.setup(code, osamp, dataOut)
646 647
647 648 self.isConfig = True
648 649
649 650 if self.code is None:
650 651 print "Fail decoding: Code is not defined."
651 652 return
652 653
653 654 datadec = None
654 655
655 656 if dataOut.flagDataAsBlock:
656 657 """
657 658 Decoding when data have been read as block,
658 659 """
659 660 if mode == 0:
660 661 datadec = self.__convolutionByBlockInTime(dataOut.data)
661 662 if mode == 1:
662 663 datadec = self.__convolutionByBlockInFreq(dataOut.data)
663 664 else:
664 665 """
665 666 Decoding when data have been read profile by profile
666 667 """
667 668 if mode == 0:
668 669 datadec = self.__convolutionInTime(dataOut.data)
669 670
670 671 if mode == 1:
671 672 datadec = self.__convolutionInFreq(dataOut.data)
672 673
673 674 if mode == 2:
674 675 datadec = self.__convolutionInFreqOpt(dataOut.data)
675 676
676 677 if datadec is None:
677 678 raise ValueError, "Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode
678 679
679 680 dataOut.code = self.code
680 681 dataOut.nCode = self.nCode
681 682 dataOut.nBaud = self.nBaud
682 683
683 684 dataOut.data = datadec
684 685
685 686 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
686 687
687 688 dataOut.flagDecodeData = True #asumo q la data esta decodificada
688 689
689 690 if self.__profIndex == self.nCode-1:
690 691 self.__profIndex = 0
691 692 return 1
692 693
693 694 self.__profIndex += 1
694 695
695 696 return 1
696 697 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
697 698
698 699
699 700 class ProfileConcat(Operation):
700 701
701 702 isConfig = False
702 703 buffer = None
703 704
704 705 def __init__(self):
705 706
706 707 Operation.__init__(self)
707 708 self.profileIndex = 0
708 709
709 710 def reset(self):
710 711 self.buffer = numpy.zeros_like(self.buffer)
711 712 self.start_index = 0
712 713 self.times = 1
713 714
714 715 def setup(self, data, m, n=1):
715 716 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
716 717 self.nHeights = data.nHeights
717 718 self.start_index = 0
718 719 self.times = 1
719 720
720 721 def concat(self, data):
721 722
722 723 self.buffer[:,self.start_index:self.profiles*self.times] = data.copy()
723 724 self.start_index = self.start_index + self.nHeights
724 725
725 726 def run(self, dataOut, m):
726 727
727 728 dataOut.flagNoData = True
728 729
729 730 if not self.isConfig:
730 731 self.setup(dataOut.data, m, 1)
731 732 self.isConfig = True
732 733
733 734 if dataOut.flagDataAsBlock:
734 735 raise ValueError, "ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False"
735 736
736 737 else:
737 738 self.concat(dataOut.data)
738 739 self.times += 1
739 740 if self.times > m:
740 741 dataOut.data = self.buffer
741 742 self.reset()
742 743 dataOut.flagNoData = False
743 744 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
744 745 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
745 746 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
746 747 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
747 748 dataOut.ippSeconds *= m
748 749
749 750 class ProfileSelector(Operation):
750 751
751 752 profileIndex = None
752 753 # Tamanho total de los perfiles
753 754 nProfiles = None
754 755
755 756 def __init__(self):
756 757
757 758 Operation.__init__(self)
758 759 self.profileIndex = 0
759 760
760 761 def incIndex(self):
761 762
762 763 self.profileIndex += 1
763 764
764 765 if self.profileIndex >= self.nProfiles:
765 766 self.profileIndex = 0
766 767
767 768 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
768 769
769 770 if profileIndex < minIndex:
770 771 return False
771 772
772 773 if profileIndex > maxIndex:
773 774 return False
774 775
775 776 return True
776 777
777 778 def isThisProfileInList(self, profileIndex, profileList):
778 779
779 780 if profileIndex not in profileList:
780 781 return False
781 782
782 783 return True
783 784
784 785 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
785 786
786 787 """
787 788 ProfileSelector:
788 789
789 790 Inputs:
790 791 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
791 792
792 793 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
793 794
794 795 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
795 796
796 797 """
797 798
798 799 dataOut.flagNoData = True
799 800
800 801 if dataOut.flagDataAsBlock:
801 802 """
802 803 data dimension = [nChannels, nProfiles, nHeis]
803 804 """
804 805 if profileList != None:
805 806 dataOut.data = dataOut.data[:,profileList,:]
806 807 dataOut.nProfiles = len(profileList)
807 808 dataOut.profileIndex = dataOut.nProfiles - 1
808 809
809 810 if profileRangeList != None:
810 811 minIndex = profileRangeList[0]
811 812 maxIndex = profileRangeList[1]
812 813
813 814 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
814 815 dataOut.nProfiles = maxIndex - minIndex + 1
815 816 dataOut.profileIndex = dataOut.nProfiles - 1
816 817
817 818 if rangeList != None:
818 819 raise ValueError, "Profile Selector: Invalid argument rangeList. Not implemented for getByBlock yet"
819 820
820 821 dataOut.flagNoData = False
821 822
822 823 return True
823 824
824 825 """
825 826 data dimension = [nChannels, nHeis]
826 827 """
827 828
828 829 if nProfiles:
829 830 self.nProfiles = nProfiles
830 831 else:
831 832 self.nProfiles = dataOut.nProfiles
832 833
833 834 if profileList != None:
834 835
835 836 dataOut.nProfiles = len(profileList)
836 837
837 838 if self.isThisProfileInList(dataOut.profileIndex, profileList):
838 839 dataOut.flagNoData = False
839 840 dataOut.profileIndex = self.profileIndex
840 841
841 842 self.incIndex()
842 843 return True
843 844
844 845 if profileRangeList != None:
845 846
846 847 minIndex = profileRangeList[0]
847 848 maxIndex = profileRangeList[1]
848 849
849 850 dataOut.nProfiles = maxIndex - minIndex + 1
850 851
851 852 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
852 853 dataOut.flagNoData = False
853 854 dataOut.profileIndex = self.profileIndex
854 855
855 856 self.incIndex()
856 857 return True
857 858
858 859 if rangeList != None:
859 860
860 861 nProfiles = 0
861 862
862 863 for thisRange in rangeList:
863 864 minIndex = thisRange[0]
864 865 maxIndex = thisRange[1]
865 866
866 867 nProfiles += maxIndex - minIndex + 1
867 868
868 869 dataOut.nProfiles = nProfiles
869 870
870 871 for thisRange in rangeList:
871 872
872 873 minIndex = thisRange[0]
873 874 maxIndex = thisRange[1]
874 875
875 876 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
876 877
877 878 # print "profileIndex = ", dataOut.profileIndex
878 879
879 880 dataOut.flagNoData = False
880 881 dataOut.profileIndex = self.profileIndex
881 882
882 883 self.incIndex()
883 884 break
884 885 return True
885 886
886 887
887 888 if beam != None: #beam is only for AMISR data
888 889 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
889 890 dataOut.flagNoData = False
890 891 dataOut.profileIndex = self.profileIndex
891 892
892 893 self.incIndex()
893 894
894 895 return True
895 896
896 897 raise ValueError, "ProfileSelector needs profileList, profileRangeList or rangeList parameter"
897 898
898 899 return False
899 900
900 901
901 902
902 903 class Reshaper(Operation):
903 904
904 905 def __init__(self):
905 906
906 907 Operation.__init__(self)
907 908 self.updateNewHeights = True
908 909
909 910 def run(self, dataOut, shape):
910 911
911 912 if not dataOut.flagDataAsBlock:
912 913 raise ValueError, "Reshaper can only be used when voltage have been read as Block, getBlock = True"
913 914
914 915 if len(shape) != 3:
915 916 raise ValueError, "shape len should be equal to 3, (nChannels, nProfiles, nHeis)"
916 917
917 918 shape_tuple = tuple(shape)
918 919 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
919 920 dataOut.flagNoData = False
920 921
921 922 if self.updateNewHeights:
922 923
923 924 old_nheights = dataOut.nHeights
924 925 new_nheights = dataOut.data.shape[2]
925 926 factor = 1.0*new_nheights / old_nheights
926 927
927 928 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
928 929
929 930 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * factor
930 931
931 932 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
932 933
933 934 dataOut.nProfiles = dataOut.data.shape[1]
934 935
935 936 dataOut.ippSeconds *= factor
936 937 #
937 938 # import collections
938 939 # from scipy.stats import mode
939 940 #
940 941 # class Synchronize(Operation):
941 942 #
942 943 # isConfig = False
943 944 # __profIndex = 0
944 945 #
945 946 # def __init__(self):
946 947 #
947 948 # Operation.__init__(self)
948 949 # # self.isConfig = False
949 950 # self.__powBuffer = None
950 951 # self.__startIndex = 0
951 952 # self.__pulseFound = False
952 953 #
953 954 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
954 955 #
955 956 # #Read data
956 957 #
957 958 # powerdB = dataOut.getPower(channel = channel)
958 959 # noisedB = dataOut.getNoise(channel = channel)[0]
959 960 #
960 961 # self.__powBuffer.extend(powerdB.flatten())
961 962 #
962 963 # dataArray = numpy.array(self.__powBuffer)
963 964 #
964 965 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
965 966 #
966 967 # maxValue = numpy.nanmax(filteredPower)
967 968 #
968 969 # if maxValue < noisedB + 10:
969 970 # #No se encuentra ningun pulso de transmision
970 971 # return None
971 972 #
972 973 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
973 974 #
974 975 # if len(maxValuesIndex) < 2:
975 976 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
976 977 # return None
977 978 #
978 979 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
979 980 #
980 981 # #Seleccionar solo valores con un espaciamiento de nSamples
981 982 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
982 983 #
983 984 # if len(pulseIndex) < 2:
984 985 # #Solo se encontro un pulso de transmision con ancho mayor a 1
985 986 # return None
986 987 #
987 988 # spacing = pulseIndex[1:] - pulseIndex[:-1]
988 989 #
989 990 # #remover senales que se distancien menos de 10 unidades o muestras
990 991 # #(No deberian existir IPP menor a 10 unidades)
991 992 #
992 993 # realIndex = numpy.where(spacing > 10 )[0]
993 994 #
994 995 # if len(realIndex) < 2:
995 996 # #Solo se encontro un pulso de transmision con ancho mayor a 1
996 997 # return None
997 998 #
998 999 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
999 1000 # realPulseIndex = pulseIndex[realIndex]
1000 1001 #
1001 1002 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1002 1003 #
1003 1004 # print "IPP = %d samples" %period
1004 1005 #
1005 1006 # self.__newNSamples = dataOut.nHeights #int(period)
1006 1007 # self.__startIndex = int(realPulseIndex[0])
1007 1008 #
1008 1009 # return 1
1009 1010 #
1010 1011 #
1011 1012 # def setup(self, nSamples, nChannels, buffer_size = 4):
1012 1013 #
1013 1014 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1014 1015 # maxlen = buffer_size*nSamples)
1015 1016 #
1016 1017 # bufferList = []
1017 1018 #
1018 1019 # for i in range(nChannels):
1019 1020 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1020 1021 # maxlen = buffer_size*nSamples)
1021 1022 #
1022 1023 # bufferList.append(bufferByChannel)
1023 1024 #
1024 1025 # self.__nSamples = nSamples
1025 1026 # self.__nChannels = nChannels
1026 1027 # self.__bufferList = bufferList
1027 1028 #
1028 1029 # def run(self, dataOut, channel = 0):
1029 1030 #
1030 1031 # if not self.isConfig:
1031 1032 # nSamples = dataOut.nHeights
1032 1033 # nChannels = dataOut.nChannels
1033 1034 # self.setup(nSamples, nChannels)
1034 1035 # self.isConfig = True
1035 1036 #
1036 1037 # #Append new data to internal buffer
1037 1038 # for thisChannel in range(self.__nChannels):
1038 1039 # bufferByChannel = self.__bufferList[thisChannel]
1039 1040 # bufferByChannel.extend(dataOut.data[thisChannel])
1040 1041 #
1041 1042 # if self.__pulseFound:
1042 1043 # self.__startIndex -= self.__nSamples
1043 1044 #
1044 1045 # #Finding Tx Pulse
1045 1046 # if not self.__pulseFound:
1046 1047 # indexFound = self.__findTxPulse(dataOut, channel)
1047 1048 #
1048 1049 # if indexFound == None:
1049 1050 # dataOut.flagNoData = True
1050 1051 # return
1051 1052 #
1052 1053 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1053 1054 # self.__pulseFound = True
1054 1055 # self.__startIndex = indexFound
1055 1056 #
1056 1057 # #If pulse was found ...
1057 1058 # for thisChannel in range(self.__nChannels):
1058 1059 # bufferByChannel = self.__bufferList[thisChannel]
1059 1060 # #print self.__startIndex
1060 1061 # x = numpy.array(bufferByChannel)
1061 1062 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1062 1063 #
1063 1064 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1064 1065 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1065 1066 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1066 1067 #
1067 1068 # dataOut.data = self.__arrayBuffer
1068 1069 #
1069 1070 # self.__startIndex += self.__newNSamples
1070 1071 #
1071 1072 # return No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now