##// END OF EJS Templates
schain
RSS
Update ACF first step
avaldez -
r1240:ce4b76019d00
parent child
Show More
Side by Side Unified
Context file:
r1240:ce4b76019d00 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1251 +1,1254
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 from schainpy import cSchain
13 13
14 14
15 15 def getNumpyDtype(dataTypeCode):
16 16
17 17 if dataTypeCode == 0:
18 18 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
19 19 elif dataTypeCode == 1:
20 20 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
21 21 elif dataTypeCode == 2:
22 22 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
23 23 elif dataTypeCode == 3:
24 24 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
25 25 elif dataTypeCode == 4:
26 26 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
27 27 elif dataTypeCode == 5:
28 28 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
29 29 else:
30 30 raise ValueError, 'dataTypeCode was not defined'
31 31
32 32 return numpyDtype
33 33
34 34
35 35 def getDataTypeCode(numpyDtype):
36 36
37 37 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
38 38 datatype = 0
39 39 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
40 40 datatype = 1
41 41 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
42 42 datatype = 2
43 43 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
44 44 datatype = 3
45 45 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
46 46 datatype = 4
47 47 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
48 48 datatype = 5
49 49 else:
50 50 datatype = None
51 51
52 52 return datatype
53 53
54 54
55 55 def hildebrand_sekhon(data, navg):
56 56 """
57 57 This method is for the objective determination of the noise level in Doppler spectra. This
58 58 implementation technique is based on the fact that the standard deviation of the spectral
59 59 densities is equal to the mean spectral density for white Gaussian noise
60 60
61 61 Inputs:
62 62 Data : heights
63 63 navg : numbers of averages
64 64
65 65 Return:
66 66 -1 : any error
67 67 anoise : noise's level
68 68 """
69 69
70 70 sortdata = numpy.sort(data, axis=None)
71 71 # lenOfData = len(sortdata)
72 72 # nums_min = lenOfData*0.2
73 73 #
74 74 # if nums_min <= 5:
75 75 # nums_min = 5
76 76 #
77 77 # sump = 0.
78 78 #
79 79 # sumq = 0.
80 80 #
81 81 # j = 0
82 82 #
83 83 # cont = 1
84 84 #
85 85 # while((cont==1)and(j<lenOfData)):
86 86 #
87 87 # sump += sortdata[j]
88 88 #
89 89 # sumq += sortdata[j]**2
90 90 #
91 91 # if j > nums_min:
92 92 # rtest = float(j)/(j-1) + 1.0/navg
93 93 # if ((sumq*j) > (rtest*sump**2)):
94 94 # j = j - 1
95 95 # sump = sump - sortdata[j]
96 96 # sumq = sumq - sortdata[j]**2
97 97 # cont = 0
98 98 #
99 99 # j += 1
100 100 #
101 101 # lnoise = sump /j
102 102 #
103 103 # return lnoise
104 104
105 105 return cSchain.hildebrand_sekhon(sortdata, navg)
106 106
107 107
108 108 class Beam:
109 109
110 110 def __init__(self):
111 111 self.codeList = []
112 112 self.azimuthList = []
113 113 self.zenithList = []
114 114
115 115
116 116 class GenericData(object):
117 117
118 118 flagNoData = True
119 119
120 120 def copy(self, inputObj=None):
121 121
122 122 if inputObj == None:
123 123 return copy.deepcopy(self)
124 124
125 125 for key in inputObj.__dict__.keys():
126 126
127 127 attribute = inputObj.__dict__[key]
128 128
129 129 # If this attribute is a tuple or list
130 130 if type(inputObj.__dict__[key]) in (tuple, list):
131 131 self.__dict__[key] = attribute[:]
132 132 continue
133 133
134 134 # If this attribute is another object or instance
135 135 if hasattr(attribute, '__dict__'):
136 136 self.__dict__[key] = attribute.copy()
137 137 continue
138 138
139 139 self.__dict__[key] = inputObj.__dict__[key]
140 140
141 141 def deepcopy(self):
142 142
143 143 return copy.deepcopy(self)
144 144
145 145 def isEmpty(self):
146 146
147 147 return self.flagNoData
148 148
149 149
150 150 class JROData(GenericData):
151 151
152 152 # m_BasicHeader = BasicHeader()
153 153 # m_ProcessingHeader = ProcessingHeader()
154 154
155 155 systemHeaderObj = SystemHeader()
156 156
157 157 radarControllerHeaderObj = RadarControllerHeader()
158 158
159 159 # data = None
160 160
161 161 type = None
162 162
163 163 datatype = None # dtype but in string
164 164
165 165 # dtype = None
166 166
167 167 # nChannels = None
168 168
169 169 # nHeights = None
170 170
171 171 nProfiles = None
172 172
173 173 heightList = None
174 174
175 175 channelList = None
176 176
177 177 flagDiscontinuousBlock = False
178 178
179 179 useLocalTime = False
180 180
181 181 utctime = None
182 182
183 183 timeZone = None
184 184
185 185 dstFlag = None
186 186
187 187 errorCount = None
188 188
189 189 blocksize = None
190 190
191 191 # nCode = None
192 192 #
193 193 # nBaud = None
194 194 #
195 195 # code = None
196 196
197 197 flagDecodeData = False # asumo q la data no esta decodificada
198 198
199 199 flagDeflipData = False # asumo q la data no esta sin flip
200 200
201 201 flagShiftFFT = False
202 202
203 203 # ippSeconds = None
204 204
205 205 # timeInterval = None
206 206
207 207 nCohInt = None
208 208
209 209 # noise = None
210 210
211 211 windowOfFilter = 1
212 212
213 213 # Speed of ligth
214 214 C = 3e8
215 215
216 216 frequency = 49.92e6
217 217
218 218 realtime = False
219 219
220 220 beacon_heiIndexList = None
221 221
222 222 last_block = None
223 223
224 224 blocknow = None
225 225
226 226 azimuth = None
227 227
228 228 zenith = None
229 229
230 230 beam = Beam()
231 231
232 232 profileIndex = None
233 233
234 234 def getNoise(self):
235 235
236 236 raise NotImplementedError
237 237
238 238 def getNChannels(self):
239 239
240 240 return len(self.channelList)
241 241
242 242 def getChannelIndexList(self):
243 243
244 244 return range(self.nChannels)
245 245
246 246 def getNHeights(self):
247 247
248 248 return len(self.heightList)
249 249
250 250 def getHeiRange(self, extrapoints=0):
251 251
252 252 heis = self.heightList
253 253 # deltah = self.heightList[1] - self.heightList[0]
254 254 #
255 255 # heis.append(self.heightList[-1])
256 256
257 257 return heis
258 258
259 259 def getDeltaH(self):
260 260
261 261 delta = self.heightList[1] - self.heightList[0]
262 262
263 263 return delta
264 264
265 265 def getltctime(self):
266 266
267 267 if self.useLocalTime:
268 268 return self.utctime - self.timeZone * 60
269 269
270 270 return self.utctime
271 271
272 272 def getDatatime(self):
273 273
274 274 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
275 275 return datatimeValue
276 276
277 277 def getTimeRange(self):
278 278
279 279 datatime = []
280 280
281 281 datatime.append(self.ltctime)
282 282 datatime.append(self.ltctime + self.timeInterval + 1)
283 283
284 284 datatime = numpy.array(datatime)
285 285
286 286 return datatime
287 287
288 288 def getFmaxTimeResponse(self):
289 289
290 290 period = (10**-6) * self.getDeltaH() / (0.15)
291 291
292 292 PRF = 1. / (period * self.nCohInt)
293 293
294 294 fmax = PRF
295 295
296 296 return fmax
297 297
298 298 def getFmax(self):
299 299 PRF = 1. / (self.ippSeconds * self.nCohInt)
300 300
301 301 fmax = PRF
302 302 return fmax
303 303
304 304 def getVmax(self):
305 305
306 306 _lambda = self.C / self.frequency
307 307
308 308 vmax = self.getFmax() * _lambda / 2
309 309
310 310 return vmax
311 311
312 312 def get_ippSeconds(self):
313 313 '''
314 314 '''
315 315 return self.radarControllerHeaderObj.ippSeconds
316 316
317 317 def set_ippSeconds(self, ippSeconds):
318 318 '''
319 319 '''
320 320
321 321 self.radarControllerHeaderObj.ippSeconds = ippSeconds
322 322
323 323 return
324 324
325 325 def get_dtype(self):
326 326 '''
327 327 '''
328 328 return getNumpyDtype(self.datatype)
329 329
330 330 def set_dtype(self, numpyDtype):
331 331 '''
332 332 '''
333 333
334 334 self.datatype = getDataTypeCode(numpyDtype)
335 335
336 336 def get_code(self):
337 337 '''
338 338 '''
339 339 return self.radarControllerHeaderObj.code
340 340
341 341 def set_code(self, code):
342 342 '''
343 343 '''
344 344 self.radarControllerHeaderObj.code = code
345 345
346 346 return
347 347
348 348 def get_ncode(self):
349 349 '''
350 350 '''
351 351 return self.radarControllerHeaderObj.nCode
352 352
353 353 def set_ncode(self, nCode):
354 354 '''
355 355 '''
356 356 self.radarControllerHeaderObj.nCode = nCode
357 357
358 358 return
359 359
360 360 def get_nbaud(self):
361 361 '''
362 362 '''
363 363 return self.radarControllerHeaderObj.nBaud
364 364
365 365 def set_nbaud(self, nBaud):
366 366 '''
367 367 '''
368 368 self.radarControllerHeaderObj.nBaud = nBaud
369 369
370 370 return
371 371
372 372 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
373 373 channelIndexList = property(
374 374 getChannelIndexList, "I'm the 'channelIndexList' property.")
375 375 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
376 376 #noise = property(getNoise, "I'm the 'nHeights' property.")
377 377 datatime = property(getDatatime, "I'm the 'datatime' property")
378 378 ltctime = property(getltctime, "I'm the 'ltctime' property")
379 379 ippSeconds = property(get_ippSeconds, set_ippSeconds)
380 380 dtype = property(get_dtype, set_dtype)
381 381 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
382 382 code = property(get_code, set_code)
383 383 nCode = property(get_ncode, set_ncode)
384 384 nBaud = property(get_nbaud, set_nbaud)
385 385
386 386
387 387 class Voltage(JROData):
388 388
389 389 # data es un numpy array de 2 dmensiones (canales, alturas)
390 390 data = None
391 391
392 392 def __init__(self):
393 393 '''
394 394 Constructor
395 395 '''
396 396
397 397 self.useLocalTime = True
398 398
399 399 self.radarControllerHeaderObj = RadarControllerHeader()
400 400
401 401 self.systemHeaderObj = SystemHeader()
402 402
403 403 self.type = "Voltage"
404 404
405 405 self.data = None
406 406
407 407 # self.dtype = None
408 408
409 409 # self.nChannels = 0
410 410
411 411 # self.nHeights = 0
412 412
413 413 self.nProfiles = None
414 414
415 415 self.heightList = None
416 416
417 417 self.channelList = None
418 418
419 419 # self.channelIndexList = None
420 420
421 421 self.flagNoData = True
422 422
423 423 self.flagDiscontinuousBlock = False
424 424
425 425 self.utctime = None
426 426
427 427 self.timeZone = None
428 428
429 429 self.dstFlag = None
430 430
431 431 self.errorCount = None
432 432
433 433 self.nCohInt = None
434 434
435 435 self.blocksize = None
436 436
437 437 self.flagDecodeData = False # asumo q la data no esta decodificada
438 438
439 439 self.flagDeflipData = False # asumo q la data no esta sin flip
440 440
441 441 self.flagShiftFFT = False
442 442
443 443 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
444 444
445 445 self.profileIndex = 0
446 446
447 447 def getNoisebyHildebrand(self, channel=None):
448 448 """
449 449 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
450 450
451 451 Return:
452 452 noiselevel
453 453 """
454 454
455 455 if channel != None:
456 456 data = self.data[channel]
457 457 nChannels = 1
458 458 else:
459 459 data = self.data
460 460 nChannels = self.nChannels
461 461
462 462 noise = numpy.zeros(nChannels)
463 463 power = data * numpy.conjugate(data)
464 464
465 465 for thisChannel in range(nChannels):
466 466 if nChannels == 1:
467 467 daux = power[:].real
468 468 else:
469 469 daux = power[thisChannel, :].real
470 470 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
471 471
472 472 return noise
473 473
474 474 def getNoise(self, type=1, channel=None):
475 475
476 476 if type == 1:
477 477 noise = self.getNoisebyHildebrand(channel)
478 478
479 479 return noise
480 480
481 481 def getPower(self, channel=None):
482 482
483 483 if channel != None:
484 484 data = self.data[channel]
485 485 else:
486 486 data = self.data
487 487
488 488 power = data * numpy.conjugate(data)
489 489 powerdB = 10 * numpy.log10(power.real)
490 490 powerdB = numpy.squeeze(powerdB)
491 491
492 492 return powerdB
493 493
494 494 def getTimeInterval(self):
495 495
496 496 timeInterval = self.ippSeconds * self.nCohInt
497 497
498 498 return timeInterval
499 499
500 500 noise = property(getNoise, "I'm the 'nHeights' property.")
501 501 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
502 502
503 503
504 504 class Spectra(JROData):
505 505
506 506 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
507 507 data_spc = None
508 508
509 509 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
510 510 data_cspc = None
511 511
512 512 # data dc es un numpy array de 2 dmensiones (canales, alturas)
513 513 data_dc = None
514 514
515 515 # data power
516 516 data_pwr = None
517 517
518 518 nFFTPoints = None
519 519
520 520 # nPairs = None
521 521
522 522 pairsList = None
523 523
524 524 nIncohInt = None
525 525
526 526 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
527 527
528 528 nCohInt = None # se requiere para determinar el valor de timeInterval
529 529
530 530 ippFactor = None
531 531
532 532 profileIndex = 0
533 533
534 534 plotting = "spectra"
535 535
536 536 def __init__(self):
537 537 '''
538 538 Constructor
539 539 '''
540 540
541 541 self.useLocalTime = True
542 542
543 543 self.radarControllerHeaderObj = RadarControllerHeader()
544 544
545 545 self.systemHeaderObj = SystemHeader()
546 546
547 547 self.type = "Spectra"
548 548
549 549 # self.data = None
550 550
551 551 # self.dtype = None
552 552
553 553 # self.nChannels = 0
554 554
555 555 # self.nHeights = 0
556 556
557 557 self.nProfiles = None
558 558
559 559 self.heightList = None
560 560
561 561 self.channelList = None
562 562
563 563 # self.channelIndexList = None
564 564
565 565 self.pairsList = None
566 566
567 567 self.flagNoData = True
568 568
569 569 self.flagDiscontinuousBlock = False
570 570
571 571 self.utctime = None
572 572
573 573 self.nCohInt = None
574 574
575 575 self.nIncohInt = None
576 576
577 577 self.blocksize = None
578 578
579 579 self.nFFTPoints = None
580 580
581 581 self.wavelength = None
582 582
583 583 self.flagDecodeData = False # asumo q la data no esta decodificada
584 584
585 585 self.flagDeflipData = False # asumo q la data no esta sin flip
586 586
587 587 self.flagShiftFFT = False
588 588
589 589 self.ippFactor = 1
590 590
591 591 #self.noise = None
592 592
593 593 self.beacon_heiIndexList = []
594 594
595 595 self.noise_estimation = None
596 596
597 597 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
598 598 """
599 599 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
600 600
601 601 Return:
602 602 noiselevel
603 603 """
604 604
605 605 noise = numpy.zeros(self.nChannels)
606 606
607 607 for channel in range(self.nChannels):
608 #print "confuse",self.data_spc.dtype
608 609 daux = self.data_spc[channel,
609 610 xmin_index:xmax_index, ymin_index:ymax_index]
611
612 #print "HI3.0",(daux.dtype),daux.shape
610 613 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
611 614
612 615 return noise
613 616
614 617 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
615 618
616 619 if self.noise_estimation is not None:
617 620 # this was estimated by getNoise Operation defined in jroproc_spectra.py
618 621 return self.noise_estimation
619 622 else:
620 623 noise = self.getNoisebyHildebrand(
621 624 xmin_index, xmax_index, ymin_index, ymax_index)
622 625 return noise
623 626
624 627 def getFreqRangeTimeResponse(self, extrapoints=0):
625 628
626 629 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
627 630 freqrange = deltafreq * \
628 631 (numpy.arange(self.nFFTPoints + extrapoints) -
629 632 self.nFFTPoints / 2.) - deltafreq / 2
630 633
631 634 return freqrange
632 635
633 636 def getAcfRange(self, extrapoints=0):
634 637
635 638 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
636 639 freqrange = deltafreq * \
637 640 (numpy.arange(self.nFFTPoints + extrapoints) -
638 641 self.nFFTPoints / 2.) - deltafreq / 2
639 642
640 643 return freqrange
641 644
642 645 def getFreqRange(self, extrapoints=0):
643 646
644 647 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
645 648 freqrange = deltafreq * \
646 649 (numpy.arange(self.nFFTPoints + extrapoints) -
647 650 self.nFFTPoints / 2.) - deltafreq / 2
648 651
649 652 return freqrange
650 653
651 654 def getVelRange(self, extrapoints=0):
652 655
653 656 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
654 657 velrange = deltav * (numpy.arange(self.nFFTPoints +
655 658 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
656 659
657 660 return velrange
658 661
659 662 def getNPairs(self):
660 663
661 664 return len(self.pairsList)
662 665
663 666 def getPairsIndexList(self):
664 667
665 668 return range(self.nPairs)
666 669
667 670 def getNormFactor(self):
668 671
669 672 pwcode = 1
670 673
671 674 if self.flagDecodeData:
672 675 pwcode = numpy.sum(self.code[0]**2)
673 676 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
674 677 normFactor = self.nProfiles * self.nIncohInt * \
675 678 self.nCohInt * pwcode * self.windowOfFilter
676 679
677 680 return normFactor
678 681
679 682 def getFlagCspc(self):
680 683
681 684 if self.data_cspc is None:
682 685 return True
683 686
684 687 return False
685 688
686 689 def getFlagDc(self):
687 690
688 691 if self.data_dc is None:
689 692 return True
690 693
691 694 return False
692 695
693 696 def getTimeInterval(self):
694 697
695 698 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
696 699
697 700 return timeInterval
698 701
699 702 def getPower(self):
700 703
701 704 factor = self.normFactor
702 705 z = self.data_spc / factor
703 706 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
704 707 avg = numpy.average(z, axis=1)
705 708
706 709 return 10 * numpy.log10(avg)
707 710
708 711 def getCoherence(self, pairsList=None, phase=False):
709 712
710 713 z = []
711 714 if pairsList is None:
712 715 pairsIndexList = self.pairsIndexList
713 716 else:
714 717 pairsIndexList = []
715 718 for pair in pairsList:
716 719 if pair not in self.pairsList:
717 720 raise ValueError, "Pair %s is not in dataOut.pairsList" % (
718 721 pair)
719 722 pairsIndexList.append(self.pairsList.index(pair))
720 723 for i in range(len(pairsIndexList)):
721 724 pair = self.pairsList[pairsIndexList[i]]
722 725 ccf = numpy.average(
723 726 self.data_cspc[pairsIndexList[i], :, :], axis=0)
724 727 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
725 728 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
726 729 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
727 730 if phase:
728 731 data = numpy.arctan2(avgcoherenceComplex.imag,
729 732 avgcoherenceComplex.real) * 180 / numpy.pi
730 733 else:
731 734 data = numpy.abs(avgcoherenceComplex)
732 735
733 736 z.append(data)
734 737
735 738 return numpy.array(z)
736 739
737 740 def setValue(self, value):
738 741
739 742 print "This property should not be initialized"
740 743
741 744 return
742 745
743 746 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
744 747 pairsIndexList = property(
745 748 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
746 749 normFactor = property(getNormFactor, setValue,
747 750 "I'm the 'getNormFactor' property.")
748 751 flag_cspc = property(getFlagCspc, setValue)
749 752 flag_dc = property(getFlagDc, setValue)
750 753 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
751 754 timeInterval = property(getTimeInterval, setValue,
752 755 "I'm the 'timeInterval' property")
753 756
754 757
755 758 class SpectraHeis(Spectra):
756 759
757 760 data_spc = None
758 761
759 762 data_cspc = None
760 763
761 764 data_dc = None
762 765
763 766 nFFTPoints = None
764 767
765 768 # nPairs = None
766 769
767 770 pairsList = None
768 771
769 772 nCohInt = None
770 773
771 774 nIncohInt = None
772 775
773 776 def __init__(self):
774 777
775 778 self.radarControllerHeaderObj = RadarControllerHeader()
776 779
777 780 self.systemHeaderObj = SystemHeader()
778 781
779 782 self.type = "SpectraHeis"
780 783
781 784 # self.dtype = None
782 785
783 786 # self.nChannels = 0
784 787
785 788 # self.nHeights = 0
786 789
787 790 self.nProfiles = None
788 791
789 792 self.heightList = None
790 793
791 794 self.channelList = None
792 795
793 796 # self.channelIndexList = None
794 797
795 798 self.flagNoData = True
796 799
797 800 self.flagDiscontinuousBlock = False
798 801
799 802 # self.nPairs = 0
800 803
801 804 self.utctime = None
802 805
803 806 self.blocksize = None
804 807
805 808 self.profileIndex = 0
806 809
807 810 self.nCohInt = 1
808 811
809 812 self.nIncohInt = 1
810 813
811 814 def getNormFactor(self):
812 815 pwcode = 1
813 816 if self.flagDecodeData:
814 817 pwcode = numpy.sum(self.code[0]**2)
815 818
816 819 normFactor = self.nIncohInt * self.nCohInt * pwcode
817 820
818 821 return normFactor
819 822
820 823 def getTimeInterval(self):
821 824
822 825 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
823 826
824 827 return timeInterval
825 828
826 829 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
827 830 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
828 831
829 832
830 833 class Fits(JROData):
831 834
832 835 heightList = None
833 836
834 837 channelList = None
835 838
836 839 flagNoData = True
837 840
838 841 flagDiscontinuousBlock = False
839 842
840 843 useLocalTime = False
841 844
842 845 utctime = None
843 846
844 847 timeZone = None
845 848
846 849 # ippSeconds = None
847 850
848 851 # timeInterval = None
849 852
850 853 nCohInt = None
851 854
852 855 nIncohInt = None
853 856
854 857 noise = None
855 858
856 859 windowOfFilter = 1
857 860
858 861 # Speed of ligth
859 862 C = 3e8
860 863
861 864 frequency = 49.92e6
862 865
863 866 realtime = False
864 867
865 868 def __init__(self):
866 869
867 870 self.type = "Fits"
868 871
869 872 self.nProfiles = None
870 873
871 874 self.heightList = None
872 875
873 876 self.channelList = None
874 877
875 878 # self.channelIndexList = None
876 879
877 880 self.flagNoData = True
878 881
879 882 self.utctime = None
880 883
881 884 self.nCohInt = 1
882 885
883 886 self.nIncohInt = 1
884 887
885 888 self.useLocalTime = True
886 889
887 890 self.profileIndex = 0
888 891
889 892 # self.utctime = None
890 893 # self.timeZone = None
891 894 # self.ltctime = None
892 895 # self.timeInterval = None
893 896 # self.header = None
894 897 # self.data_header = None
895 898 # self.data = None
896 899 # self.datatime = None
897 900 # self.flagNoData = False
898 901 # self.expName = ''
899 902 # self.nChannels = None
900 903 # self.nSamples = None
901 904 # self.dataBlocksPerFile = None
902 905 # self.comments = ''
903 906 #
904 907
905 908 def getltctime(self):
906 909
907 910 if self.useLocalTime:
908 911 return self.utctime - self.timeZone * 60
909 912
910 913 return self.utctime
911 914
912 915 def getDatatime(self):
913 916
914 917 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
915 918 return datatime
916 919
917 920 def getTimeRange(self):
918 921
919 922 datatime = []
920 923
921 924 datatime.append(self.ltctime)
922 925 datatime.append(self.ltctime + self.timeInterval)
923 926
924 927 datatime = numpy.array(datatime)
925 928
926 929 return datatime
927 930
928 931 def getHeiRange(self):
929 932
930 933 heis = self.heightList
931 934
932 935 return heis
933 936
934 937 def getNHeights(self):
935 938
936 939 return len(self.heightList)
937 940
938 941 def getNChannels(self):
939 942
940 943 return len(self.channelList)
941 944
942 945 def getChannelIndexList(self):
943 946
944 947 return range(self.nChannels)
945 948
946 949 def getNoise(self, type=1):
947 950
948 951 #noise = numpy.zeros(self.nChannels)
949 952
950 953 if type == 1:
951 954 noise = self.getNoisebyHildebrand()
952 955
953 956 if type == 2:
954 957 noise = self.getNoisebySort()
955 958
956 959 if type == 3:
957 960 noise = self.getNoisebyWindow()
958 961
959 962 return noise
960 963
961 964 def getTimeInterval(self):
962 965
963 966 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
964 967
965 968 return timeInterval
966 969
967 970 datatime = property(getDatatime, "I'm the 'datatime' property")
968 971 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
969 972 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
970 973 channelIndexList = property(
971 974 getChannelIndexList, "I'm the 'channelIndexList' property.")
972 975 noise = property(getNoise, "I'm the 'nHeights' property.")
973 976
974 977 ltctime = property(getltctime, "I'm the 'ltctime' property")
975 978 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
976 979
977 980
978 981 class Correlation(JROData):
979 982
980 983 noise = None
981 984
982 985 SNR = None
983 986
984 987 #--------------------------------------------------
985 988
986 989 mode = None
987 990
988 991 split = False
989 992
990 993 data_cf = None
991 994
992 995 lags = None
993 996
994 997 lagRange = None
995 998
996 999 pairsList = None
997 1000
998 1001 normFactor = None
999 1002
1000 1003 #--------------------------------------------------
1001 1004
1002 1005 # calculateVelocity = None
1003 1006
1004 1007 nLags = None
1005 1008
1006 1009 nPairs = None
1007 1010
1008 1011 nAvg = None
1009 1012
1010 1013 def __init__(self):
1011 1014 '''
1012 1015 Constructor
1013 1016 '''
1014 1017 self.radarControllerHeaderObj = RadarControllerHeader()
1015 1018
1016 1019 self.systemHeaderObj = SystemHeader()
1017 1020
1018 1021 self.type = "Correlation"
1019 1022
1020 1023 self.data = None
1021 1024
1022 1025 self.dtype = None
1023 1026
1024 1027 self.nProfiles = None
1025 1028
1026 1029 self.heightList = None
1027 1030
1028 1031 self.channelList = None
1029 1032
1030 1033 self.flagNoData = True
1031 1034
1032 1035 self.flagDiscontinuousBlock = False
1033 1036
1034 1037 self.utctime = None
1035 1038
1036 1039 self.timeZone = None
1037 1040
1038 1041 self.dstFlag = None
1039 1042
1040 1043 self.errorCount = None
1041 1044
1042 1045 self.blocksize = None
1043 1046
1044 1047 self.flagDecodeData = False # asumo q la data no esta decodificada
1045 1048
1046 1049 self.flagDeflipData = False # asumo q la data no esta sin flip
1047 1050
1048 1051 self.pairsList = None
1049 1052
1050 1053 self.nPoints = None
1051 1054
1052 1055 def getPairsList(self):
1053 1056
1054 1057 return self.pairsList
1055 1058
1056 1059 def getNoise(self, mode=2):
1057 1060
1058 1061 indR = numpy.where(self.lagR == 0)[0][0]
1059 1062 indT = numpy.where(self.lagT == 0)[0][0]
1060 1063
1061 1064 jspectra0 = self.data_corr[:, :, indR, :]
1062 1065 jspectra = copy.copy(jspectra0)
1063 1066
1064 1067 num_chan = jspectra.shape[0]
1065 1068 num_hei = jspectra.shape[2]
1066 1069
1067 1070 freq_dc = jspectra.shape[1] / 2
1068 1071 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
1069 1072
1070 1073 if ind_vel[0] < 0:
1071 1074 ind_vel[range(0, 1)] = ind_vel[range(0, 1)] + self.num_prof
1072 1075
1073 1076 if mode == 1:
1074 1077 jspectra[:, freq_dc, :] = (
1075 1078 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
1076 1079
1077 1080 if mode == 2:
1078 1081
1079 1082 vel = numpy.array([-2, -1, 1, 2])
1080 1083 xx = numpy.zeros([4, 4])
1081 1084
1082 1085 for fil in range(4):
1083 1086 xx[fil, :] = vel[fil]**numpy.asarray(range(4))
1084 1087
1085 1088 xx_inv = numpy.linalg.inv(xx)
1086 1089 xx_aux = xx_inv[0, :]
1087 1090
1088 1091 for ich in range(num_chan):
1089 1092 yy = jspectra[ich, ind_vel, :]
1090 1093 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
1091 1094
1092 1095 junkid = jspectra[ich, freq_dc, :] <= 0
1093 1096 cjunkid = sum(junkid)
1094 1097
1095 1098 if cjunkid.any():
1096 1099 jspectra[ich, freq_dc, junkid.nonzero()] = (
1097 1100 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
1098 1101
1099 1102 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
1100 1103
1101 1104 return noise
1102 1105
1103 1106 def getTimeInterval(self):
1104 1107
1105 1108 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1106 1109
1107 1110 return timeInterval
1108 1111
1109 1112 def splitFunctions(self):
1110 1113
1111 1114 pairsList = self.pairsList
1112 1115 ccf_pairs = []
1113 1116 acf_pairs = []
1114 1117 ccf_ind = []
1115 1118 acf_ind = []
1116 1119 for l in range(len(pairsList)):
1117 1120 chan0 = pairsList[l][0]
1118 1121 chan1 = pairsList[l][1]
1119 1122
1120 1123 # Obteniendo pares de Autocorrelacion
1121 1124 if chan0 == chan1:
1122 1125 acf_pairs.append(chan0)
1123 1126 acf_ind.append(l)
1124 1127 else:
1125 1128 ccf_pairs.append(pairsList[l])
1126 1129 ccf_ind.append(l)
1127 1130
1128 1131 data_acf = self.data_cf[acf_ind]
1129 1132 data_ccf = self.data_cf[ccf_ind]
1130 1133
1131 1134 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1132 1135
1133 1136 def getNormFactor(self):
1134 1137 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1135 1138 acf_pairs = numpy.array(acf_pairs)
1136 1139 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1137 1140
1138 1141 for p in range(self.nPairs):
1139 1142 pair = self.pairsList[p]
1140 1143
1141 1144 ch0 = pair[0]
1142 1145 ch1 = pair[1]
1143 1146
1144 1147 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1145 1148 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1146 1149 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1147 1150
1148 1151 return normFactor
1149 1152
1150 1153 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1151 1154 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1152 1155
1153 1156
1154 1157 class Parameters(Spectra):
1155 1158
1156 1159 experimentInfo = None # Information about the experiment
1157 1160
1158 1161 # Information from previous data
1159 1162
1160 1163 inputUnit = None # Type of data to be processed
1161 1164
1162 1165 operation = None # Type of operation to parametrize
1163 1166
1164 1167 # normFactor = None #Normalization Factor
1165 1168
1166 1169 groupList = None # List of Pairs, Groups, etc
1167 1170
1168 1171 # Parameters
1169 1172
1170 1173 data_param = None # Parameters obtained
1171 1174
1172 1175 data_pre = None # Data Pre Parametrization
1173 1176
1174 1177 data_SNR = None # Signal to Noise Ratio
1175 1178
1176 1179 # heightRange = None #Heights
1177 1180
1178 1181 abscissaList = None # Abscissa, can be velocities, lags or time
1179 1182
1180 1183 # noise = None #Noise Potency
1181 1184
1182 1185 utctimeInit = None # Initial UTC time
1183 1186
1184 1187 paramInterval = None # Time interval to calculate Parameters in seconds
1185 1188
1186 1189 useLocalTime = True
1187 1190
1188 1191 # Fitting
1189 1192
1190 1193 data_error = None # Error of the estimation
1191 1194
1192 1195 constants = None
1193 1196
1194 1197 library = None
1195 1198
1196 1199 # Output signal
1197 1200
1198 1201 outputInterval = None # Time interval to calculate output signal in seconds
1199 1202
1200 1203 data_output = None # Out signal
1201 1204
1202 1205 nAvg = None
1203 1206
1204 1207 noise_estimation = None
1205 1208
1206 1209 GauSPC = None # Fit gaussian SPC
1207 1210
1208 1211 def __init__(self):
1209 1212 '''
1210 1213 Constructor
1211 1214 '''
1212 1215 self.radarControllerHeaderObj = RadarControllerHeader()
1213 1216
1214 1217 self.systemHeaderObj = SystemHeader()
1215 1218
1216 1219 self.type = "Parameters"
1217 1220
1218 1221 def getTimeRange1(self, interval):
1219 1222
1220 1223 datatime = []
1221 1224
1222 1225 if self.useLocalTime:
1223 1226 time1 = self.utctimeInit - self.timeZone * 60
1224 1227 else:
1225 1228 time1 = self.utctimeInit
1226 1229
1227 1230 datatime.append(time1)
1228 1231 datatime.append(time1 + interval)
1229 1232 datatime = numpy.array(datatime)
1230 1233
1231 1234 return datatime
1232 1235
1233 1236 def getTimeInterval(self):
1234 1237
1235 1238 if hasattr(self, 'timeInterval1'):
1236 1239 return self.timeInterval1
1237 1240 else:
1238 1241 return self.paramInterval
1239 1242
1240 1243 def setValue(self, value):
1241 1244
1242 1245 print "This property should not be initialized"
1243 1246
1244 1247 return
1245 1248
1246 1249 def getNoise(self):
1247 1250
1248 1251 return self.spc_noise
1249 1252
1250 1253 timeInterval = property(getTimeInterval)
1251 1254 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
Show file before | Show file after | Hide comments
@@ -1,1542 +1,1543
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9
10 10 from figure import Figure, isRealtime, isTimeInHourRange
11 11 from plotting_codes import *
12 12
13 13
14 14 class SpectraPlot(Figure):
15 15
16 16 isConfig = None
17 17 __nsubplots = None
18 18
19 19 WIDTHPROF = None
20 20 HEIGHTPROF = None
21 21 PREFIX = 'spc'
22 22
23 23 def __init__(self, **kwargs):
24 24 Figure.__init__(self, **kwargs)
25 25 self.isConfig = False
26 26 self.__nsubplots = 1
27 27
28 28 self.WIDTH = 250
29 29 self.HEIGHT = 250
30 30 self.WIDTHPROF = 120
31 31 self.HEIGHTPROF = 0
32 32 self.counter_imagwr = 0
33 33
34 34 self.PLOT_CODE = SPEC_CODE
35 35
36 36 self.FTP_WEI = None
37 37 self.EXP_CODE = None
38 38 self.SUB_EXP_CODE = None
39 39 self.PLOT_POS = None
40 40
41 41 self.__xfilter_ena = False
42 42 self.__yfilter_ena = False
43 43
44 44 def getSubplots(self):
45 45
46 46 ncol = int(numpy.sqrt(self.nplots)+0.9)
47 47 nrow = int(self.nplots*1./ncol + 0.9)
48 48
49 49 return nrow, ncol
50 50
51 51 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
52 52
53 53 self.__showprofile = showprofile
54 54 self.nplots = nplots
55 55
56 56 ncolspan = 1
57 57 colspan = 1
58 58 if showprofile:
59 59 ncolspan = 3
60 60 colspan = 2
61 61 self.__nsubplots = 2
62 62
63 63 self.createFigure(id = id,
64 64 wintitle = wintitle,
65 65 widthplot = self.WIDTH + self.WIDTHPROF,
66 66 heightplot = self.HEIGHT + self.HEIGHTPROF,
67 67 show=show)
68 68
69 69 nrow, ncol = self.getSubplots()
70 70
71 71 counter = 0
72 72 for y in range(nrow):
73 73 for x in range(ncol):
74 74
75 75 if counter >= self.nplots:
76 76 break
77 77
78 78 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
79 79
80 80 if showprofile:
81 81 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
82 82
83 83 counter += 1
84 84
85 85 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
86 86 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
87 87 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
88 88 server=None, folder=None, username=None, password=None,
89 89 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
90 90 xaxis="frequency", colormap='jet', normFactor=None):
91 91
92 92 """
93 93
94 94 Input:
95 95 dataOut :
96 96 id :
97 97 wintitle :
98 98 channelList :
99 99 showProfile :
100 100 xmin : None,
101 101 xmax : None,
102 102 ymin : None,
103 103 ymax : None,
104 104 zmin : None,
105 105 zmax : None
106 106 """
107 107 if realtime:
108 108 if not(isRealtime(utcdatatime = dataOut.utctime)):
109 109 print 'Skipping this plot function'
110 110 return
111 111
112 112 if channelList == None:
113 113 channelIndexList = dataOut.channelIndexList
114 114 else:
115 115 channelIndexList = []
116 116 for channel in channelList:
117 117 if channel not in dataOut.channelList:
118 118 raise ValueError, "Channel %d is not in dataOut.channelList" %channel
119 119 channelIndexList.append(dataOut.channelList.index(channel))
120 120
121 121 if normFactor is None:
122 122 factor = dataOut.normFactor
123 123 else:
124 124 factor = normFactor
125 125 if xaxis == "frequency":
126 126 x = dataOut.getFreqRange(1)/1000.
127 127 xlabel = "Frequency (kHz)"
128 128
129 129 elif xaxis == "time":
130 130 x = dataOut.getAcfRange(1)
131 131 xlabel = "Time (ms)"
132 132
133 133 else:
134 134 x = dataOut.getVelRange(1)
135 135 xlabel = "Velocity (m/s)"
136 136
137 137 ylabel = "Range (Km)"
138 138
139 139 y = dataOut.getHeiRange()
140 140
141 141 z = dataOut.data_spc/factor
142 142 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
143 143 zdB = 10*numpy.log10(z)
144 144
145 #print "a000",dataOut.data_spc.dtype
145 146 avg = numpy.average(z, axis=1)
146 147 avgdB = 10*numpy.log10(avg)
147
148 #print "before plot"
148 149 noise = dataOut.getNoise()/factor
149 150 noisedB = 10*numpy.log10(noise)
150 151
151 152 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
152 153 title = wintitle + " Spectra"
153 154 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
154 155 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
155 156
156 157 if not self.isConfig:
157 158
158 159 nplots = len(channelIndexList)
159 160
160 161 self.setup(id=id,
161 162 nplots=nplots,
162 163 wintitle=wintitle,
163 164 showprofile=showprofile,
164 165 show=show)
165 166
166 167 if xmin == None: xmin = numpy.nanmin(x)
167 168 if xmax == None: xmax = numpy.nanmax(x)
168 169 if ymin == None: ymin = numpy.nanmin(y)
169 170 if ymax == None: ymax = numpy.nanmax(y)
170 171 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
171 172 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
172 173
173 174 self.FTP_WEI = ftp_wei
174 175 self.EXP_CODE = exp_code
175 176 self.SUB_EXP_CODE = sub_exp_code
176 177 self.PLOT_POS = plot_pos
177 178
178 179 self.isConfig = True
179 180
180 181 self.setWinTitle(title)
181 182
182 183 for i in range(self.nplots):
183 184 index = channelIndexList[i]
184 185 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
185 186 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
186 187 if len(dataOut.beam.codeList) != 0:
187 188 title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime)
188 189
189 190 axes = self.axesList[i*self.__nsubplots]
190 191 axes.pcolor(x, y, zdB[index,:,:],
191 192 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
192 193 xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap,
193 194 ticksize=9, cblabel='')
194 195
195 196 if self.__showprofile:
196 197 axes = self.axesList[i*self.__nsubplots +1]
197 198 axes.pline(avgdB[index,:], y,
198 199 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
199 200 xlabel='dB', ylabel='', title='',
200 201 ytick_visible=False,
201 202 grid='x')
202 203
203 204 noiseline = numpy.repeat(noisedB[index], len(y))
204 205 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
205 206
206 207 self.draw()
207 208
208 209 if figfile == None:
209 210 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
210 211 name = str_datetime
211 212 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
212 213 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
213 214 figfile = self.getFilename(name)
214 215
215 216 self.save(figpath=figpath,
216 217 figfile=figfile,
217 218 save=save,
218 219 ftp=ftp,
219 220 wr_period=wr_period,
220 221 thisDatetime=thisDatetime)
221 222
222 223 class CrossSpectraPlot(Figure):
223 224
224 225 isConfig = None
225 226 __nsubplots = None
226 227
227 228 WIDTH = None
228 229 HEIGHT = None
229 230 WIDTHPROF = None
230 231 HEIGHTPROF = None
231 232 PREFIX = 'cspc'
232 233
233 234 def __init__(self, **kwargs):
234 235 Figure.__init__(self, **kwargs)
235 236 self.isConfig = False
236 237 self.__nsubplots = 4
237 238 self.counter_imagwr = 0
238 239 self.WIDTH = 250
239 240 self.HEIGHT = 250
240 241 self.WIDTHPROF = 0
241 242 self.HEIGHTPROF = 0
242 243
243 244 self.PLOT_CODE = CROSS_CODE
244 245 self.FTP_WEI = None
245 246 self.EXP_CODE = None
246 247 self.SUB_EXP_CODE = None
247 248 self.PLOT_POS = None
248 249
249 250 def getSubplots(self):
250 251
251 252 ncol = 4
252 253 nrow = self.nplots
253 254
254 255 return nrow, ncol
255 256
256 257 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
257 258
258 259 self.__showprofile = showprofile
259 260 self.nplots = nplots
260 261
261 262 ncolspan = 1
262 263 colspan = 1
263 264
264 265 self.createFigure(id = id,
265 266 wintitle = wintitle,
266 267 widthplot = self.WIDTH + self.WIDTHPROF,
267 268 heightplot = self.HEIGHT + self.HEIGHTPROF,
268 269 show=True)
269 270
270 271 nrow, ncol = self.getSubplots()
271 272
272 273 counter = 0
273 274 for y in range(nrow):
274 275 for x in range(ncol):
275 276 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
276 277
277 278 counter += 1
278 279
279 280 def run(self, dataOut, id, wintitle="", pairsList=None,
280 281 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
281 282 coh_min=None, coh_max=None, phase_min=None, phase_max=None,
282 283 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
283 284 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
284 285 server=None, folder=None, username=None, password=None,
285 286 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None,
286 287 xaxis='frequency'):
287 288
288 289 """
289 290
290 291 Input:
291 292 dataOut :
292 293 id :
293 294 wintitle :
294 295 channelList :
295 296 showProfile :
296 297 xmin : None,
297 298 xmax : None,
298 299 ymin : None,
299 300 ymax : None,
300 301 zmin : None,
301 302 zmax : None
302 303 """
303 304
304 305 if pairsList == None:
305 306 pairsIndexList = dataOut.pairsIndexList
306 307 else:
307 308 pairsIndexList = []
308 309 for pair in pairsList:
309 310 if pair not in dataOut.pairsList:
310 311 raise ValueError, "Pair %s is not in dataOut.pairsList" %str(pair)
311 312 pairsIndexList.append(dataOut.pairsList.index(pair))
312 313
313 314 if not pairsIndexList:
314 315 return
315 316
316 317 if len(pairsIndexList) > 4:
317 318 pairsIndexList = pairsIndexList[0:4]
318
319
319 320 if normFactor is None:
320 321 factor = dataOut.normFactor
321 322 else:
322 323 factor = normFactor
323 324 x = dataOut.getVelRange(1)
324 325 y = dataOut.getHeiRange()
325 326 z = dataOut.data_spc[:,:,:]/factor
326 327 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
327 328
328 329 noise = dataOut.noise/factor
329 330
330 331 zdB = 10*numpy.log10(z)
331 332 noisedB = 10*numpy.log10(noise)
332 333
333 334 if coh_min == None:
334 335 coh_min = 0.0
335 336 if coh_max == None:
336 337 coh_max = 1.0
337 338
338 339 if phase_min == None:
339 340 phase_min = -180
340 341 if phase_max == None:
341 342 phase_max = 180
342 343
343 344 #thisDatetime = dataOut.datatime
344 345 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
345 346 title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
346 347 # xlabel = "Velocity (m/s)"
347 348 ylabel = "Range (Km)"
348 349
349 350 if xaxis == "frequency":
350 351 x = dataOut.getFreqRange(1)/1000.
351 352 xlabel = "Frequency (kHz)"
352 353
353 354 elif xaxis == "time":
354 355 x = dataOut.getAcfRange(1)
355 356 xlabel = "Time (ms)"
356 357
357 358 else:
358 359 x = dataOut.getVelRange(1)
359 360 xlabel = "Velocity (m/s)"
360 361
361 362 if not self.isConfig:
362 363
363 364 nplots = len(pairsIndexList)
364 365
365 366 self.setup(id=id,
366 367 nplots=nplots,
367 368 wintitle=wintitle,
368 369 showprofile=False,
369 370 show=show)
370 371
371 372 avg = numpy.abs(numpy.average(z, axis=1))
372 373 avgdB = 10*numpy.log10(avg)
373 374
374 375 if xmin == None: xmin = numpy.nanmin(x)
375 376 if xmax == None: xmax = numpy.nanmax(x)
376 377 if ymin == None: ymin = numpy.nanmin(y)
377 378 if ymax == None: ymax = numpy.nanmax(y)
378 379 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
379 380 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
380 381
381 382 self.FTP_WEI = ftp_wei
382 383 self.EXP_CODE = exp_code
383 384 self.SUB_EXP_CODE = sub_exp_code
384 385 self.PLOT_POS = plot_pos
385 386
386 387 self.isConfig = True
387 388
388 389 self.setWinTitle(title)
389 390
390 391 for i in range(self.nplots):
391 392 pair = dataOut.pairsList[pairsIndexList[i]]
392 393
393 394 chan_index0 = dataOut.channelList.index(pair[0])
394 395 chan_index1 = dataOut.channelList.index(pair[1])
395 396
396 397 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
397 398 title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[chan_index0], str_datetime)
398 399 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index0,:,:]/factor)
399 400 axes0 = self.axesList[i*self.__nsubplots]
400 401 axes0.pcolor(x, y, zdB,
401 402 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
402 403 xlabel=xlabel, ylabel=ylabel, title=title,
403 404 ticksize=9, colormap=power_cmap, cblabel='')
404 405
405 406 title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[chan_index1], str_datetime)
406 407 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index1,:,:]/factor)
407 408 axes0 = self.axesList[i*self.__nsubplots+1]
408 409 axes0.pcolor(x, y, zdB,
409 410 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
410 411 xlabel=xlabel, ylabel=ylabel, title=title,
411 412 ticksize=9, colormap=power_cmap, cblabel='')
412 413
413 414 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:])
414 415 coherence = numpy.abs(coherenceComplex)
415 416 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
416 417 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
417 418
418 419 title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1])
419 420 axes0 = self.axesList[i*self.__nsubplots+2]
420 421 axes0.pcolor(x, y, coherence,
421 422 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=coh_min, zmax=coh_max,
422 423 xlabel=xlabel, ylabel=ylabel, title=title,
423 424 ticksize=9, colormap=coherence_cmap, cblabel='')
424 425
425 426 title = "Phase Ch%d * Ch%d" %(pair[0], pair[1])
426 427 axes0 = self.axesList[i*self.__nsubplots+3]
427 428 axes0.pcolor(x, y, phase,
428 429 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
429 430 xlabel=xlabel, ylabel=ylabel, title=title,
430 431 ticksize=9, colormap=phase_cmap, cblabel='')
431 432
432 433
433 434
434 435 self.draw()
435 436
436 437 self.save(figpath=figpath,
437 438 figfile=figfile,
438 439 save=save,
439 440 ftp=ftp,
440 441 wr_period=wr_period,
441 442 thisDatetime=thisDatetime)
442 443
443 444
444 445 class RTIPlot(Figure):
445 446
446 447 __isConfig = None
447 448 __nsubplots = None
448 449
449 450 WIDTHPROF = None
450 451 HEIGHTPROF = None
451 452 PREFIX = 'rti'
452 453
453 454 def __init__(self, **kwargs):
454 455
455 456 Figure.__init__(self, **kwargs)
456 457 self.timerange = None
457 458 self.isConfig = False
458 459 self.__nsubplots = 1
459 460
460 461 self.WIDTH = 800
461 462 self.HEIGHT = 180
462 463 self.WIDTHPROF = 120
463 464 self.HEIGHTPROF = 0
464 465 self.counter_imagwr = 0
465 466
466 467 self.PLOT_CODE = RTI_CODE
467 468
468 469 self.FTP_WEI = None
469 470 self.EXP_CODE = None
470 471 self.SUB_EXP_CODE = None
471 472 self.PLOT_POS = None
472 473 self.tmin = None
473 474 self.tmax = None
474 475
475 476 self.xmin = None
476 477 self.xmax = None
477 478
478 479 self.figfile = None
479 480
480 481 def getSubplots(self):
481 482
482 483 ncol = 1
483 484 nrow = self.nplots
484 485
485 486 return nrow, ncol
486 487
487 488 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
488 489
489 490 self.__showprofile = showprofile
490 491 self.nplots = nplots
491 492
492 493 ncolspan = 1
493 494 colspan = 1
494 495 if showprofile:
495 496 ncolspan = 7
496 497 colspan = 6
497 498 self.__nsubplots = 2
498 499
499 500 self.createFigure(id = id,
500 501 wintitle = wintitle,
501 502 widthplot = self.WIDTH + self.WIDTHPROF,
502 503 heightplot = self.HEIGHT + self.HEIGHTPROF,
503 504 show=show)
504 505
505 506 nrow, ncol = self.getSubplots()
506 507
507 508 counter = 0
508 509 for y in range(nrow):
509 510 for x in range(ncol):
510 511
511 512 if counter >= self.nplots:
512 513 break
513 514
514 515 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
515 516
516 517 if showprofile:
517 518 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
518 519
519 520 counter += 1
520 521
521 522 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
522 523 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
523 524 timerange=None, colormap='jet',
524 525 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
525 526 server=None, folder=None, username=None, password=None,
526 527 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, HEIGHT=None):
527 528
528 529 """
529 530
530 531 Input:
531 532 dataOut :
532 533 id :
533 534 wintitle :
534 535 channelList :
535 536 showProfile :
536 537 xmin : None,
537 538 xmax : None,
538 539 ymin : None,
539 540 ymax : None,
540 541 zmin : None,
541 542 zmax : None
542 543 """
543 544
544 545 #colormap = kwargs.get('colormap', 'jet')
545 546 if HEIGHT is not None:
546 547 self.HEIGHT = HEIGHT
547
548
548 549 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
549 550 return
550 551
551 552 if channelList == None:
552 553 channelIndexList = dataOut.channelIndexList
553 554 else:
554 555 channelIndexList = []
555 556 for channel in channelList:
556 557 if channel not in dataOut.channelList:
557 558 raise ValueError, "Channel %d is not in dataOut.channelList"
558 559 channelIndexList.append(dataOut.channelList.index(channel))
559 560
560 561 if normFactor is None:
561 562 factor = dataOut.normFactor
562 563 else:
563 564 factor = normFactor
564 565
565 566 # factor = dataOut.normFactor
566 567 x = dataOut.getTimeRange()
567 568 y = dataOut.getHeiRange()
568 569
569 570 z = dataOut.data_spc/factor
570 571 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
571 572 avg = numpy.average(z, axis=1)
572 573 avgdB = 10.*numpy.log10(avg)
573 574 # avgdB = dataOut.getPower()
574 575
575 576
576 577 thisDatetime = dataOut.datatime
577 578 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
578 579 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
579 580 xlabel = ""
580 581 ylabel = "Range (Km)"
581 582
582 583 update_figfile = False
583 584
584 585 if dataOut.ltctime >= self.xmax:
585 586 self.counter_imagwr = wr_period
586 587 self.isConfig = False
587 588 update_figfile = True
588 589
589 590 if not self.isConfig:
590 591
591 592 nplots = len(channelIndexList)
592 593
593 594 self.setup(id=id,
594 595 nplots=nplots,
595 596 wintitle=wintitle,
596 597 showprofile=showprofile,
597 598 show=show)
598 599
599 600 if timerange != None:
600 601 self.timerange = timerange
601 602
602 603 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
603 604
604 605 noise = dataOut.noise/factor
605 606 noisedB = 10*numpy.log10(noise)
606 607
607 608 if ymin == None: ymin = numpy.nanmin(y)
608 609 if ymax == None: ymax = numpy.nanmax(y)
609 610 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
610 611 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
611 612
612 613 self.FTP_WEI = ftp_wei
613 614 self.EXP_CODE = exp_code
614 615 self.SUB_EXP_CODE = sub_exp_code
615 616 self.PLOT_POS = plot_pos
616 617
617 618 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
618 619 self.isConfig = True
619 620 self.figfile = figfile
620 621 update_figfile = True
621 622
622 623 self.setWinTitle(title)
623 624
624 625 for i in range(self.nplots):
625 626 index = channelIndexList[i]
626 627 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
627 628 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
628 629 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
629 630 axes = self.axesList[i*self.__nsubplots]
630 631 zdB = avgdB[index].reshape((1,-1))
631 632 axes.pcolorbuffer(x, y, zdB,
632 633 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
633 634 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
634 635 ticksize=9, cblabel='', cbsize="1%", colormap=colormap)
635 636
636 637 if self.__showprofile:
637 638 axes = self.axesList[i*self.__nsubplots +1]
638 639 axes.pline(avgdB[index], y,
639 640 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
640 641 xlabel='dB', ylabel='', title='',
641 642 ytick_visible=False,
642 643 grid='x')
643 644
644 645 self.draw()
645 646
646 647 self.save(figpath=figpath,
647 648 figfile=figfile,
648 649 save=save,
649 650 ftp=ftp,
650 651 wr_period=wr_period,
651 652 thisDatetime=thisDatetime,
652 653 update_figfile=update_figfile)
653 654
654 655 class CoherenceMap(Figure):
655 656 isConfig = None
656 657 __nsubplots = None
657 658
658 659 WIDTHPROF = None
659 660 HEIGHTPROF = None
660 661 PREFIX = 'cmap'
661 662
662 663 def __init__(self, **kwargs):
663 664 Figure.__init__(self, **kwargs)
664 665 self.timerange = 2*60*60
665 666 self.isConfig = False
666 667 self.__nsubplots = 1
667 668
668 669 self.WIDTH = 800
669 670 self.HEIGHT = 180
670 671 self.WIDTHPROF = 120
671 672 self.HEIGHTPROF = 0
672 673 self.counter_imagwr = 0
673 674
674 675 self.PLOT_CODE = COH_CODE
675 676
676 677 self.FTP_WEI = None
677 678 self.EXP_CODE = None
678 679 self.SUB_EXP_CODE = None
679 680 self.PLOT_POS = None
680 681 self.counter_imagwr = 0
681 682
682 683 self.xmin = None
683 684 self.xmax = None
684 685
685 686 def getSubplots(self):
686 687 ncol = 1
687 688 nrow = self.nplots*2
688 689
689 690 return nrow, ncol
690 691
691 692 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
692 693 self.__showprofile = showprofile
693 694 self.nplots = nplots
694 695
695 696 ncolspan = 1
696 697 colspan = 1
697 698 if showprofile:
698 699 ncolspan = 7
699 700 colspan = 6
700 701 self.__nsubplots = 2
701 702
702 703 self.createFigure(id = id,
703 704 wintitle = wintitle,
704 705 widthplot = self.WIDTH + self.WIDTHPROF,
705 706 heightplot = self.HEIGHT + self.HEIGHTPROF,
706 707 show=True)
707 708
708 709 nrow, ncol = self.getSubplots()
709 710
710 711 for y in range(nrow):
711 712 for x in range(ncol):
712 713
713 714 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
714 715
715 716 if showprofile:
716 717 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
717 718
718 719 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
719 720 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
720 721 timerange=None, phase_min=None, phase_max=None,
721 722 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
722 723 coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
723 724 server=None, folder=None, username=None, password=None,
724 725 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
725 726
726 727 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
727 728 return
728 729
729 730 if pairsList == None:
730 731 pairsIndexList = dataOut.pairsIndexList
731 732 else:
732 733 pairsIndexList = []
733 734 for pair in pairsList:
734 735 if pair not in dataOut.pairsList:
735 736 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
736 737 pairsIndexList.append(dataOut.pairsList.index(pair))
737 738
738 739 if pairsIndexList == []:
739 740 return
740 741
741 742 if len(pairsIndexList) > 4:
742 743 pairsIndexList = pairsIndexList[0:4]
743 744
744 745 if phase_min == None:
745 746 phase_min = -180
746 747 if phase_max == None:
747 748 phase_max = 180
748 749
749 750 x = dataOut.getTimeRange()
750 751 y = dataOut.getHeiRange()
751 752
752 753 thisDatetime = dataOut.datatime
753 754
754 755 title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
755 756 xlabel = ""
756 757 ylabel = "Range (Km)"
757 758 update_figfile = False
758 759
759 760 if not self.isConfig:
760 761 nplots = len(pairsIndexList)
761 762 self.setup(id=id,
762 763 nplots=nplots,
763 764 wintitle=wintitle,
764 765 showprofile=showprofile,
765 766 show=show)
766 767
767 768 if timerange != None:
768 769 self.timerange = timerange
769 770
770 771 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
771 772
772 773 if ymin == None: ymin = numpy.nanmin(y)
773 774 if ymax == None: ymax = numpy.nanmax(y)
774 775 if zmin == None: zmin = 0.
775 776 if zmax == None: zmax = 1.
776 777
777 778 self.FTP_WEI = ftp_wei
778 779 self.EXP_CODE = exp_code
779 780 self.SUB_EXP_CODE = sub_exp_code
780 781 self.PLOT_POS = plot_pos
781 782
782 783 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
783 784
784 785 self.isConfig = True
785 786 update_figfile = True
786 787
787 788 self.setWinTitle(title)
788 789
789 790 for i in range(self.nplots):
790 791
791 792 pair = dataOut.pairsList[pairsIndexList[i]]
792 793
793 794 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
794 795 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
795 796 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
796 797
797 798
798 799 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
799 800 coherence = numpy.abs(avgcoherenceComplex)
800 801
801 802 z = coherence.reshape((1,-1))
802 803
803 804 counter = 0
804 805
805 806 title = "Coherence Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
806 807 axes = self.axesList[i*self.__nsubplots*2]
807 808 axes.pcolorbuffer(x, y, z,
808 809 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
809 810 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
810 811 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
811 812
812 813 if self.__showprofile:
813 814 counter += 1
814 815 axes = self.axesList[i*self.__nsubplots*2 + counter]
815 816 axes.pline(coherence, y,
816 817 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
817 818 xlabel='', ylabel='', title='', ticksize=7,
818 819 ytick_visible=False, nxticks=5,
819 820 grid='x')
820 821
821 822 counter += 1
822 823
823 824 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
824 825
825 826 z = phase.reshape((1,-1))
826 827
827 828 title = "Phase Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
828 829 axes = self.axesList[i*self.__nsubplots*2 + counter]
829 830 axes.pcolorbuffer(x, y, z,
830 831 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
831 832 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
832 833 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
833 834
834 835 if self.__showprofile:
835 836 counter += 1
836 837 axes = self.axesList[i*self.__nsubplots*2 + counter]
837 838 axes.pline(phase, y,
838 839 xmin=phase_min, xmax=phase_max, ymin=ymin, ymax=ymax,
839 840 xlabel='', ylabel='', title='', ticksize=7,
840 841 ytick_visible=False, nxticks=4,
841 842 grid='x')
842 843
843 844 self.draw()
844 845
845 846 if dataOut.ltctime >= self.xmax:
846 847 self.counter_imagwr = wr_period
847 848 self.isConfig = False
848 849 update_figfile = True
849 850
850 851 self.save(figpath=figpath,
851 852 figfile=figfile,
852 853 save=save,
853 854 ftp=ftp,
854 855 wr_period=wr_period,
855 856 thisDatetime=thisDatetime,
856 857 update_figfile=update_figfile)
857 858
858 859 class PowerProfilePlot(Figure):
859 860
860 861 isConfig = None
861 862 __nsubplots = None
862 863
863 864 WIDTHPROF = None
864 865 HEIGHTPROF = None
865 866 PREFIX = 'spcprofile'
866 867
867 868 def __init__(self, **kwargs):
868 869 Figure.__init__(self, **kwargs)
869 870 self.isConfig = False
870 871 self.__nsubplots = 1
871 872
872 873 self.PLOT_CODE = POWER_CODE
873 874
874 875 self.WIDTH = 300
875 876 self.HEIGHT = 500
876 877 self.counter_imagwr = 0
877 878
878 879 def getSubplots(self):
879 880 ncol = 1
880 881 nrow = 1
881 882
882 883 return nrow, ncol
883 884
884 885 def setup(self, id, nplots, wintitle, show):
885 886
886 887 self.nplots = nplots
887 888
888 889 ncolspan = 1
889 890 colspan = 1
890 891
891 892 self.createFigure(id = id,
892 893 wintitle = wintitle,
893 894 widthplot = self.WIDTH,
894 895 heightplot = self.HEIGHT,
895 896 show=show)
896 897
897 898 nrow, ncol = self.getSubplots()
898 899
899 900 counter = 0
900 901 for y in range(nrow):
901 902 for x in range(ncol):
902 903 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
903 904
904 905 def run(self, dataOut, id, wintitle="", channelList=None,
905 906 xmin=None, xmax=None, ymin=None, ymax=None,
906 907 save=False, figpath='./', figfile=None, show=True,
907 908 ftp=False, wr_period=1, server=None,
908 909 folder=None, username=None, password=None):
909 910
910 911
911 912 if channelList == None:
912 913 channelIndexList = dataOut.channelIndexList
913 914 channelList = dataOut.channelList
914 915 else:
915 916 channelIndexList = []
916 917 for channel in channelList:
917 918 if channel not in dataOut.channelList:
918 919 raise ValueError, "Channel %d is not in dataOut.channelList"
919 920 channelIndexList.append(dataOut.channelList.index(channel))
920 921
921 922 factor = dataOut.normFactor
922 923
923 924 y = dataOut.getHeiRange()
924 925
925 926 #for voltage
926 927 if dataOut.type == 'Voltage':
927 928 x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
928 929 x = x.real
929 930 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
930 931
931 932 #for spectra
932 933 if dataOut.type == 'Spectra':
933 934 x = dataOut.data_spc[channelIndexList,:,:]/factor
934 935 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
935 936 x = numpy.average(x, axis=1)
936 937
937 938
938 939 xdB = 10*numpy.log10(x)
939 940
940 941 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
941 942 title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
942 943 xlabel = "dB"
943 944 ylabel = "Range (Km)"
944 945
945 946 if not self.isConfig:
946 947
947 948 nplots = 1
948 949
949 950 self.setup(id=id,
950 951 nplots=nplots,
951 952 wintitle=wintitle,
952 953 show=show)
953 954
954 955 if ymin == None: ymin = numpy.nanmin(y)
955 956 if ymax == None: ymax = numpy.nanmax(y)
956 957 if xmin == None: xmin = numpy.nanmin(xdB)*0.9
957 958 if xmax == None: xmax = numpy.nanmax(xdB)*1.1
958 959
959 960 self.isConfig = True
960 961
961 962 self.setWinTitle(title)
962 963
963 964 title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
964 965 axes = self.axesList[0]
965 966
966 967 legendlabels = ["channel %d"%x for x in channelList]
967 968 axes.pmultiline(xdB, y,
968 969 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
969 970 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
970 971 ytick_visible=True, nxticks=5,
971 972 grid='x')
972 973
973 974 self.draw()
974 975
975 976 self.save(figpath=figpath,
976 977 figfile=figfile,
977 978 save=save,
978 979 ftp=ftp,
979 980 wr_period=wr_period,
980 981 thisDatetime=thisDatetime)
981 982
982 983 class SpectraCutPlot(Figure):
983 984
984 985 isConfig = None
985 986 __nsubplots = None
986 987
987 988 WIDTHPROF = None
988 989 HEIGHTPROF = None
989 990 PREFIX = 'spc_cut'
990 991
991 992 def __init__(self, **kwargs):
992 993 Figure.__init__(self, **kwargs)
993 994 self.isConfig = False
994 995 self.__nsubplots = 1
995 996
996 997 self.PLOT_CODE = POWER_CODE
997 998
998 999 self.WIDTH = 700
999 1000 self.HEIGHT = 500
1000 1001 self.counter_imagwr = 0
1001 1002
1002 1003 def getSubplots(self):
1003 1004 ncol = 1
1004 1005 nrow = 1
1005 1006
1006 1007 return nrow, ncol
1007 1008
1008 1009 def setup(self, id, nplots, wintitle, show):
1009 1010
1010 1011 self.nplots = nplots
1011 1012
1012 1013 ncolspan = 1
1013 1014 colspan = 1
1014 1015
1015 1016 self.createFigure(id = id,
1016 1017 wintitle = wintitle,
1017 1018 widthplot = self.WIDTH,
1018 1019 heightplot = self.HEIGHT,
1019 1020 show=show)
1020 1021
1021 1022 nrow, ncol = self.getSubplots()
1022 1023
1023 1024 counter = 0
1024 1025 for y in range(nrow):
1025 1026 for x in range(ncol):
1026 1027 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1027 1028
1028 1029 def run(self, dataOut, id, wintitle="", channelList=None,
1029 1030 xmin=None, xmax=None, ymin=None, ymax=None,
1030 1031 save=False, figpath='./', figfile=None, show=True,
1031 1032 ftp=False, wr_period=1, server=None,
1032 1033 folder=None, username=None, password=None,
1033 1034 xaxis="frequency"):
1034 1035
1035 1036
1036 1037 if channelList == None:
1037 1038 channelIndexList = dataOut.channelIndexList
1038 1039 channelList = dataOut.channelList
1039 1040 else:
1040 1041 channelIndexList = []
1041 1042 for channel in channelList:
1042 1043 if channel not in dataOut.channelList:
1043 1044 raise ValueError, "Channel %d is not in dataOut.channelList"
1044 1045 channelIndexList.append(dataOut.channelList.index(channel))
1045 1046
1046 1047 factor = dataOut.normFactor
1047 1048
1048 1049 y = dataOut.getHeiRange()
1049 1050
1050 1051 z = dataOut.data_spc/factor
1051 1052 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1052 1053
1053 1054 hei_index = numpy.arange(25)*3 + 20
1054 1055
1055 1056 if xaxis == "frequency":
1056 1057 x = dataOut.getFreqRange()/1000.
1057 1058 zdB = 10*numpy.log10(z[0,:,hei_index])
1058 1059 xlabel = "Frequency (kHz)"
1059 1060 ylabel = "Power (dB)"
1060 1061
1061 1062 elif xaxis == "time":
1062 1063 x = dataOut.getAcfRange()
1063 1064 zdB = z[0,:,hei_index]
1064 1065 xlabel = "Time (ms)"
1065 1066 ylabel = "ACF"
1066 1067
1067 1068 else:
1068 1069 x = dataOut.getVelRange()
1069 1070 zdB = 10*numpy.log10(z[0,:,hei_index])
1070 1071 xlabel = "Velocity (m/s)"
1071 1072 ylabel = "Power (dB)"
1072 1073
1073 1074 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1074 1075 title = wintitle + " Range Cuts %s" %(thisDatetime.strftime("%d-%b-%Y"))
1075 1076
1076 1077 if not self.isConfig:
1077 1078
1078 1079 nplots = 1
1079 1080
1080 1081 self.setup(id=id,
1081 1082 nplots=nplots,
1082 1083 wintitle=wintitle,
1083 1084 show=show)
1084 1085
1085 1086 if xmin == None: xmin = numpy.nanmin(x)*0.9
1086 1087 if xmax == None: xmax = numpy.nanmax(x)*1.1
1087 1088 if ymin == None: ymin = numpy.nanmin(zdB)
1088 1089 if ymax == None: ymax = numpy.nanmax(zdB)
1089 1090
1090 1091 self.isConfig = True
1091 1092
1092 1093 self.setWinTitle(title)
1093 1094
1094 1095 title = "Spectra Cuts: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1095 1096 axes = self.axesList[0]
1096 1097
1097 1098 legendlabels = ["Range = %dKm" %y[i] for i in hei_index]
1098 1099
1099 1100 axes.pmultilineyaxis( x, zdB,
1100 1101 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1101 1102 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1102 1103 ytick_visible=True, nxticks=5,
1103 1104 grid='x')
1104 1105
1105 1106 self.draw()
1106 1107
1107 1108 self.save(figpath=figpath,
1108 1109 figfile=figfile,
1109 1110 save=save,
1110 1111 ftp=ftp,
1111 1112 wr_period=wr_period,
1112 1113 thisDatetime=thisDatetime)
1113 1114
1114 1115 class Noise(Figure):
1115 1116
1116 1117 isConfig = None
1117 1118 __nsubplots = None
1118 1119
1119 1120 PREFIX = 'noise'
1120 1121
1121 1122
1122 1123 def __init__(self, **kwargs):
1123 1124 Figure.__init__(self, **kwargs)
1124 1125 self.timerange = 24*60*60
1125 1126 self.isConfig = False
1126 1127 self.__nsubplots = 1
1127 1128 self.counter_imagwr = 0
1128 1129 self.WIDTH = 800
1129 1130 self.HEIGHT = 400
1130 1131 self.WIDTHPROF = 120
1131 1132 self.HEIGHTPROF = 0
1132 1133 self.xdata = None
1133 1134 self.ydata = None
1134 1135
1135 1136 self.PLOT_CODE = NOISE_CODE
1136 1137
1137 1138 self.FTP_WEI = None
1138 1139 self.EXP_CODE = None
1139 1140 self.SUB_EXP_CODE = None
1140 1141 self.PLOT_POS = None
1141 1142 self.figfile = None
1142 1143
1143 1144 self.xmin = None
1144 1145 self.xmax = None
1145 1146
1146 1147 def getSubplots(self):
1147 1148
1148 1149 ncol = 1
1149 1150 nrow = 1
1150 1151
1151 1152 return nrow, ncol
1152 1153
1153 1154 def openfile(self, filename):
1154 1155 dirname = os.path.dirname(filename)
1155 1156
1156 1157 if not os.path.exists(dirname):
1157 1158 os.mkdir(dirname)
1158 1159
1159 1160 f = open(filename,'w+')
1160 1161 f.write('\n\n')
1161 1162 f.write('JICAMARCA RADIO OBSERVATORY - Noise \n')
1162 1163 f.write('DD MM YYYY HH MM SS Channel0 Channel1 Channel2 Channel3\n\n' )
1163 1164 f.close()
1164 1165
1165 1166 def save_data(self, filename_phase, data, data_datetime):
1166 1167
1167 1168 f=open(filename_phase,'a')
1168 1169
1169 1170 timetuple_data = data_datetime.timetuple()
1170 1171 day = str(timetuple_data.tm_mday)
1171 1172 month = str(timetuple_data.tm_mon)
1172 1173 year = str(timetuple_data.tm_year)
1173 1174 hour = str(timetuple_data.tm_hour)
1174 1175 minute = str(timetuple_data.tm_min)
1175 1176 second = str(timetuple_data.tm_sec)
1176 1177
1177 1178 data_msg = ''
1178 1179 for i in range(len(data)):
1179 1180 data_msg += str(data[i]) + ' '
1180 1181
1181 1182 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' ' + data_msg + '\n')
1182 1183 f.close()
1183 1184
1184 1185
1185 1186 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1186 1187
1187 1188 self.__showprofile = showprofile
1188 1189 self.nplots = nplots
1189 1190
1190 1191 ncolspan = 7
1191 1192 colspan = 6
1192 1193 self.__nsubplots = 2
1193 1194
1194 1195 self.createFigure(id = id,
1195 1196 wintitle = wintitle,
1196 1197 widthplot = self.WIDTH+self.WIDTHPROF,
1197 1198 heightplot = self.HEIGHT+self.HEIGHTPROF,
1198 1199 show=show)
1199 1200
1200 1201 nrow, ncol = self.getSubplots()
1201 1202
1202 1203 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1203 1204
1204 1205
1205 1206 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1206 1207 xmin=None, xmax=None, ymin=None, ymax=None,
1207 1208 timerange=None,
1208 1209 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1209 1210 server=None, folder=None, username=None, password=None,
1210 1211 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1211 1212
1212 1213 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1213 1214 return
1214 1215
1215 1216 if channelList == None:
1216 1217 channelIndexList = dataOut.channelIndexList
1217 1218 channelList = dataOut.channelList
1218 1219 else:
1219 1220 channelIndexList = []
1220 1221 for channel in channelList:
1221 1222 if channel not in dataOut.channelList:
1222 1223 raise ValueError, "Channel %d is not in dataOut.channelList"
1223 1224 channelIndexList.append(dataOut.channelList.index(channel))
1224 1225
1225 1226 x = dataOut.getTimeRange()
1226 1227 #y = dataOut.getHeiRange()
1227 1228 factor = dataOut.normFactor
1228 1229 noise = dataOut.noise[channelIndexList]/factor
1229 1230 noisedB = 10*numpy.log10(noise)
1230 1231
1231 1232 thisDatetime = dataOut.datatime
1232 1233
1233 1234 title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1234 1235 xlabel = ""
1235 1236 ylabel = "Intensity (dB)"
1236 1237 update_figfile = False
1237 1238
1238 1239 if not self.isConfig:
1239 1240
1240 1241 nplots = 1
1241 1242
1242 1243 self.setup(id=id,
1243 1244 nplots=nplots,
1244 1245 wintitle=wintitle,
1245 1246 showprofile=showprofile,
1246 1247 show=show)
1247 1248
1248 1249 if timerange != None:
1249 1250 self.timerange = timerange
1250 1251
1251 1252 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1252 1253
1253 1254 if ymin == None: ymin = numpy.floor(numpy.nanmin(noisedB)) - 10.0
1254 1255 if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
1255 1256
1256 1257 self.FTP_WEI = ftp_wei
1257 1258 self.EXP_CODE = exp_code
1258 1259 self.SUB_EXP_CODE = sub_exp_code
1259 1260 self.PLOT_POS = plot_pos
1260 1261
1261 1262
1262 1263 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1263 1264 self.isConfig = True
1264 1265 self.figfile = figfile
1265 1266 self.xdata = numpy.array([])
1266 1267 self.ydata = numpy.array([])
1267 1268
1268 1269 update_figfile = True
1269 1270
1270 1271 #open file beacon phase
1271 1272 path = '%s%03d' %(self.PREFIX, self.id)
1272 1273 noise_file = os.path.join(path,'%s.txt'%self.name)
1273 1274 self.filename_noise = os.path.join(figpath,noise_file)
1274 1275
1275 1276 self.setWinTitle(title)
1276 1277
1277 1278 title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1278 1279
1279 1280 legendlabels = ["channel %d"%(idchannel) for idchannel in channelList]
1280 1281 axes = self.axesList[0]
1281 1282
1282 1283 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1283 1284
1284 1285 if len(self.ydata)==0:
1285 1286 self.ydata = noisedB.reshape(-1,1)
1286 1287 else:
1287 1288 self.ydata = numpy.hstack((self.ydata, noisedB.reshape(-1,1)))
1288 1289
1289 1290
1290 1291 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1291 1292 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1292 1293 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1293 1294 XAxisAsTime=True, grid='both'
1294 1295 )
1295 1296
1296 1297 self.draw()
1297 1298
1298 1299 if dataOut.ltctime >= self.xmax:
1299 1300 self.counter_imagwr = wr_period
1300 1301 self.isConfig = False
1301 1302 update_figfile = True
1302 1303
1303 1304 self.save(figpath=figpath,
1304 1305 figfile=figfile,
1305 1306 save=save,
1306 1307 ftp=ftp,
1307 1308 wr_period=wr_period,
1308 1309 thisDatetime=thisDatetime,
1309 1310 update_figfile=update_figfile)
1310 1311
1311 1312 #store data beacon phase
1312 1313 if save:
1313 1314 self.save_data(self.filename_noise, noisedB, thisDatetime)
1314 1315
1315 1316 class BeaconPhase(Figure):
1316 1317
1317 1318 __isConfig = None
1318 1319 __nsubplots = None
1319 1320
1320 1321 PREFIX = 'beacon_phase'
1321 1322
1322 1323 def __init__(self, **kwargs):
1323 1324 Figure.__init__(self, **kwargs)
1324 1325 self.timerange = 24*60*60
1325 1326 self.isConfig = False
1326 1327 self.__nsubplots = 1
1327 1328 self.counter_imagwr = 0
1328 1329 self.WIDTH = 800
1329 1330 self.HEIGHT = 400
1330 1331 self.WIDTHPROF = 120
1331 1332 self.HEIGHTPROF = 0
1332 1333 self.xdata = None
1333 1334 self.ydata = None
1334 1335
1335 1336 self.PLOT_CODE = BEACON_CODE
1336 1337
1337 1338 self.FTP_WEI = None
1338 1339 self.EXP_CODE = None
1339 1340 self.SUB_EXP_CODE = None
1340 1341 self.PLOT_POS = None
1341 1342
1342 1343 self.filename_phase = None
1343 1344
1344 1345 self.figfile = None
1345 1346
1346 1347 self.xmin = None
1347 1348 self.xmax = None
1348 1349
1349 1350 def getSubplots(self):
1350 1351
1351 1352 ncol = 1
1352 1353 nrow = 1
1353 1354
1354 1355 return nrow, ncol
1355 1356
1356 1357 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1357 1358
1358 1359 self.__showprofile = showprofile
1359 1360 self.nplots = nplots
1360 1361
1361 1362 ncolspan = 7
1362 1363 colspan = 6
1363 1364 self.__nsubplots = 2
1364 1365
1365 1366 self.createFigure(id = id,
1366 1367 wintitle = wintitle,
1367 1368 widthplot = self.WIDTH+self.WIDTHPROF,
1368 1369 heightplot = self.HEIGHT+self.HEIGHTPROF,
1369 1370 show=show)
1370 1371
1371 1372 nrow, ncol = self.getSubplots()
1372 1373
1373 1374 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1374 1375
1375 1376 def save_phase(self, filename_phase):
1376 1377 f = open(filename_phase,'w+')
1377 1378 f.write('\n\n')
1378 1379 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1379 1380 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1380 1381 f.close()
1381 1382
1382 1383 def save_data(self, filename_phase, data, data_datetime):
1383 1384 f=open(filename_phase,'a')
1384 1385 timetuple_data = data_datetime.timetuple()
1385 1386 day = str(timetuple_data.tm_mday)
1386 1387 month = str(timetuple_data.tm_mon)
1387 1388 year = str(timetuple_data.tm_year)
1388 1389 hour = str(timetuple_data.tm_hour)
1389 1390 minute = str(timetuple_data.tm_min)
1390 1391 second = str(timetuple_data.tm_sec)
1391 1392 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1392 1393 f.close()
1393 1394
1394 1395
1395 1396 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1396 1397 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1397 1398 timerange=None,
1398 1399 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1399 1400 server=None, folder=None, username=None, password=None,
1400 1401 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1401 1402
1402 1403 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1403 1404 return
1404 1405
1405 1406 if pairsList == None:
1406 1407 pairsIndexList = dataOut.pairsIndexList[:10]
1407 1408 else:
1408 1409 pairsIndexList = []
1409 1410 for pair in pairsList:
1410 1411 if pair not in dataOut.pairsList:
1411 1412 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
1412 1413 pairsIndexList.append(dataOut.pairsList.index(pair))
1413 1414
1414 1415 if pairsIndexList == []:
1415 1416 return
1416 1417
1417 1418 # if len(pairsIndexList) > 4:
1418 1419 # pairsIndexList = pairsIndexList[0:4]
1419 1420
1420 1421 hmin_index = None
1421 1422 hmax_index = None
1422 1423
1423 1424 if hmin != None and hmax != None:
1424 1425 indexes = numpy.arange(dataOut.nHeights)
1425 1426 hmin_list = indexes[dataOut.heightList >= hmin]
1426 1427 hmax_list = indexes[dataOut.heightList <= hmax]
1427 1428
1428 1429 if hmin_list.any():
1429 1430 hmin_index = hmin_list[0]
1430 1431
1431 1432 if hmax_list.any():
1432 1433 hmax_index = hmax_list[-1]+1
1433 1434
1434 1435 x = dataOut.getTimeRange()
1435 1436 #y = dataOut.getHeiRange()
1436 1437
1437 1438
1438 1439 thisDatetime = dataOut.datatime
1439 1440
1440 1441 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1441 1442 xlabel = "Local Time"
1442 1443 ylabel = "Phase (degrees)"
1443 1444
1444 1445 update_figfile = False
1445 1446
1446 1447 nplots = len(pairsIndexList)
1447 1448 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1448 1449 phase_beacon = numpy.zeros(len(pairsIndexList))
1449 1450 for i in range(nplots):
1450 1451 pair = dataOut.pairsList[pairsIndexList[i]]
1451 1452 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1452 1453 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1453 1454 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1454 1455 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1455 1456 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1456 1457
1457 1458 #print "Phase %d%d" %(pair[0], pair[1])
1458 1459 #print phase[dataOut.beacon_heiIndexList]
1459 1460
1460 1461 if dataOut.beacon_heiIndexList:
1461 1462 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1462 1463 else:
1463 1464 phase_beacon[i] = numpy.average(phase)
1464 1465
1465 1466 if not self.isConfig:
1466 1467
1467 1468 nplots = len(pairsIndexList)
1468 1469
1469 1470 self.setup(id=id,
1470 1471 nplots=nplots,
1471 1472 wintitle=wintitle,
1472 1473 showprofile=showprofile,
1473 1474 show=show)
1474 1475
1475 1476 if timerange != None:
1476 1477 self.timerange = timerange
1477 1478
1478 1479 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1479 1480
1480 1481 if ymin == None: ymin = 0
1481 1482 if ymax == None: ymax = 360
1482 1483
1483 1484 self.FTP_WEI = ftp_wei
1484 1485 self.EXP_CODE = exp_code
1485 1486 self.SUB_EXP_CODE = sub_exp_code
1486 1487 self.PLOT_POS = plot_pos
1487 1488
1488 1489 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1489 1490 self.isConfig = True
1490 1491 self.figfile = figfile
1491 1492 self.xdata = numpy.array([])
1492 1493 self.ydata = numpy.array([])
1493 1494
1494 1495 update_figfile = True
1495 1496
1496 1497 #open file beacon phase
1497 1498 path = '%s%03d' %(self.PREFIX, self.id)
1498 1499 beacon_file = os.path.join(path,'%s.txt'%self.name)
1499 1500 self.filename_phase = os.path.join(figpath,beacon_file)
1500 1501 #self.save_phase(self.filename_phase)
1501 1502
1502 1503
1503 1504 #store data beacon phase
1504 1505 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1505 1506
1506 1507 self.setWinTitle(title)
1507 1508
1508 1509
1509 1510 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1510 1511
1511 1512 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1512 1513
1513 1514 axes = self.axesList[0]
1514 1515
1515 1516 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1516 1517
1517 1518 if len(self.ydata)==0:
1518 1519 self.ydata = phase_beacon.reshape(-1,1)
1519 1520 else:
1520 1521 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1521 1522
1522 1523
1523 1524 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1524 1525 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1525 1526 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1526 1527 XAxisAsTime=True, grid='both'
1527 1528 )
1528 1529
1529 1530 self.draw()
1530 1531
1531 1532 if dataOut.ltctime >= self.xmax:
1532 1533 self.counter_imagwr = wr_period
1533 1534 self.isConfig = False
1534 1535 update_figfile = True
1535 1536
1536 1537 self.save(figpath=figpath,
1537 1538 figfile=figfile,
1538 1539 save=save,
1539 1540 ftp=ftp,
1540 1541 wr_period=wr_period,
1541 1542 thisDatetime=thisDatetime,
1542 1543 update_figfile=update_figfile)
Show file before | Show file after | Hide comments
@@ -1,958 +1,960
1 1 import itertools
2 2
3 3 import numpy
4 4
5 5 from jroproc_base import ProcessingUnit, Operation
6 6 from schainpy.model.data.jrodata import Spectra
7 7 from schainpy.model.data.jrodata import hildebrand_sekhon
8 8
9 9
10 10 class SpectraProc(ProcessingUnit):
11 11
12 12 def __init__(self, **kwargs):
13 13
14 14 ProcessingUnit.__init__(self, **kwargs)
15 15
16 16 self.buffer = None
17 17 self.firstdatatime = None
18 18 self.profIndex = 0
19 19 self.dataOut = Spectra()
20 20 self.id_min = None
21 21 self.id_max = None
22 22
23 23 def __updateSpecFromVoltage(self):
24 24
25 25 self.dataOut.timeZone = self.dataIn.timeZone
26 26 self.dataOut.dstFlag = self.dataIn.dstFlag
27 27 self.dataOut.errorCount = self.dataIn.errorCount
28 28 self.dataOut.useLocalTime = self.dataIn.useLocalTime
29 29 try:
30 30 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
31 31 except:
32 32 pass
33 33 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
34 34 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
35 35 self.dataOut.channelList = self.dataIn.channelList
36 36 self.dataOut.heightList = self.dataIn.heightList
37 37 #print self.dataOut.heightList.shape,"spec4"
38 38 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
39 39
40 40 self.dataOut.nBaud = self.dataIn.nBaud
41 41 self.dataOut.nCode = self.dataIn.nCode
42 42 self.dataOut.code = self.dataIn.code
43 43 self.dataOut.nProfiles = self.dataOut.nFFTPoints
44 44
45 45 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
46 46 self.dataOut.utctime = self.firstdatatime
47 47 # asumo q la data esta decodificada
48 48 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
49 49 # asumo q la data esta sin flip
50 50 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
51 51 self.dataOut.flagShiftFFT = False
52 52
53 53 self.dataOut.nCohInt = self.dataIn.nCohInt
54 54 self.dataOut.nIncohInt = 1
55 55
56 56 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
57 57
58 58 self.dataOut.frequency = self.dataIn.frequency
59 59 self.dataOut.realtime = self.dataIn.realtime
60 60
61 61 self.dataOut.azimuth = self.dataIn.azimuth
62 62 self.dataOut.zenith = self.dataIn.zenith
63 63
64 64 self.dataOut.beam.codeList = self.dataIn.beam.codeList
65 65 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
66 66 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
67 67
68 68 def __getFft(self):
69 69 """
70 70 Convierte valores de Voltaje a Spectra
71 71
72 72 Affected:
73 73 self.dataOut.data_spc
74 74 self.dataOut.data_cspc
75 75 self.dataOut.data_dc
76 76 self.dataOut.heightList
77 77 self.profIndex
78 78 self.buffer
79 79 self.dataOut.flagNoData
80 80 """
81 81 fft_volt = numpy.fft.fft(
82 82 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
83 83 fft_volt = fft_volt.astype(numpy.dtype('complex'))
84 84 dc = fft_volt[:, 0, :]
85 85
86 86 # calculo de self-spectra
87 87 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
88 #print "spec dtype 0",fft_volt.dtype
88 89 spc = fft_volt * numpy.conjugate(fft_volt)
89 90 spc = spc.real
91 #print "spec dtype 1",spc.dtype
90 92
91 93 blocksize = 0
92 94 blocksize += dc.size
93 95 blocksize += spc.size
94 96
95 97 cspc = None
96 98 pairIndex = 0
97 99 if self.dataOut.pairsList != None:
98 100 # calculo de cross-spectra
99 101 cspc = numpy.zeros(
100 102 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
101 103 for pair in self.dataOut.pairsList:
102 104 if pair[0] not in self.dataOut.channelList:
103 105 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
104 106 str(pair), str(self.dataOut.channelList))
105 107 if pair[1] not in self.dataOut.channelList:
106 108 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
107 109 str(pair), str(self.dataOut.channelList))
108 110
109 111 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
110 112 numpy.conjugate(fft_volt[pair[1], :, :])
111 113 pairIndex += 1
112 114 blocksize += cspc.size
113 115
114 116 self.dataOut.data_spc = spc
115 117 self.dataOut.data_cspc = cspc
116 118 self.dataOut.data_dc = dc
117 119 self.dataOut.blockSize = blocksize
118 120 self.dataOut.flagShiftFFT = True
119 121
120 122 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None, shift_fft=False):
121 123
122 124 self.dataOut.flagNoData = True
123 125
124 126 if self.dataIn.type == "Spectra":
125 127 self.dataOut.copy(self.dataIn)
126 128 # if not pairsList:
127 129 # pairsList = itertools.combinations(self.dataOut.channelList, 2)
128 130 # if self.dataOut.data_cspc is not None:
129 131 # self.__selectPairs(pairsList)
130 132 if shift_fft:
131 133 #desplaza a la derecha en el eje 2 determinadas posiciones
132 134 shift = int(self.dataOut.nFFTPoints/2)
133 135 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
134 136
135 137 if self.dataOut.data_cspc is not None:
136 138 #desplaza a la derecha en el eje 2 determinadas posiciones
137 139 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
138 140
139 141 return True
140 142
141 143 if self.dataIn.type == "Voltage":
142 144
143 145 if nFFTPoints == None:
144 146 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
145 147
146 148 if nProfiles == None:
147 149 nProfiles = nFFTPoints
148 150
149 151 if ippFactor == None:
150 152 ippFactor = 1
151 153
152 154 self.dataOut.ippFactor = ippFactor
153 155
154 156 self.dataOut.nFFTPoints = nFFTPoints
155 157 self.dataOut.pairsList = pairsList
156 158
157 159 if self.buffer is None:
158 160 self.buffer = numpy.zeros((self.dataIn.nChannels,
159 161 nProfiles,
160 162 self.dataIn.heightList.shape[0]),
161 163 dtype='complex')
162 164
163 165 #print self.buffer.shape,"spec2"
164 166 #print self.dataIn.heightList.shape[0],"spec3"
165 167
166 168 if self.dataIn.flagDataAsBlock:
167 169 # data dimension: [nChannels, nProfiles, nSamples]
168 170 nVoltProfiles = self.dataIn.data.shape[1]
169 171 # nVoltProfiles = self.dataIn.nProfiles
170 172
171 173 #print nVoltProfiles,"spec1"
172 174 #print nProfiles
173 175 if nVoltProfiles == nProfiles:
174 176 self.buffer = self.dataIn.data.copy()
175 177 self.profIndex = nVoltProfiles
176 178
177 179 elif nVoltProfiles < nProfiles:
178 180
179 181 if self.profIndex == 0:
180 182 self.id_min = 0
181 183 self.id_max = nVoltProfiles
182 184
183 185 self.buffer[:, self.id_min:self.id_max,:] = self.dataIn.data
184 186 self.profIndex += nVoltProfiles
185 187 self.id_min += nVoltProfiles
186 188 self.id_max += nVoltProfiles
187 189 else:
188 190 raise ValueError, "The type object %s has %d profiles, it should just has %d profiles" % (
189 191 self.dataIn.type, self.dataIn.data.shape[1], nProfiles)
190 192 self.dataOut.flagNoData = True
191 193 return 0
192 194 else:
193 195 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
194 196 self.profIndex += 1
195 197 #print self.profIndex,"spectra D"
196 198
197 199 if self.firstdatatime == None:
198 200 self.firstdatatime = self.dataIn.utctime
199 201
200 202 if self.profIndex == nProfiles:
201 203 self.__updateSpecFromVoltage()
202 204 self.__getFft()
203 205
204 206 self.dataOut.flagNoData = False
205 207 self.firstdatatime = None
206 208 self.profIndex = 0
207 209
208 210 return True
209 211
210 212 raise ValueError, "The type of input object '%s' is not valid" % (
211 213 self.dataIn.type)
212 214
213 215 def __selectPairs(self, pairsList):
214 216
215 217 if not pairsList:
216 218 return
217 219
218 220 pairs = []
219 221 pairsIndex = []
220 222
221 223 for pair in pairsList:
222 224 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
223 225 continue
224 226 pairs.append(pair)
225 227 pairsIndex.append(pairs.index(pair))
226 228
227 229 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
228 230 self.dataOut.pairsList = pairs
229 231
230 232 return
231 233
232 234 def __selectPairsByChannel(self, channelList=None):
233 235
234 236 if channelList == None:
235 237 return
236 238
237 239 pairsIndexListSelected = []
238 240 for pairIndex in self.dataOut.pairsIndexList:
239 241 # First pair
240 242 if self.dataOut.pairsList[pairIndex][0] not in channelList:
241 243 continue
242 244 # Second pair
243 245 if self.dataOut.pairsList[pairIndex][1] not in channelList:
244 246 continue
245 247
246 248 pairsIndexListSelected.append(pairIndex)
247 249
248 250 if not pairsIndexListSelected:
249 251 self.dataOut.data_cspc = None
250 252 self.dataOut.pairsList = []
251 253 return
252 254
253 255 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
254 256 self.dataOut.pairsList = [self.dataOut.pairsList[i]
255 257 for i in pairsIndexListSelected]
256 258
257 259 return
258 260
259 261 def selectChannels(self, channelList):
260 262
261 263 channelIndexList = []
262 264
263 265 for channel in channelList:
264 266 if channel not in self.dataOut.channelList:
265 267 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" % (
266 268 channel, str(self.dataOut.channelList))
267 269
268 270 index = self.dataOut.channelList.index(channel)
269 271 channelIndexList.append(index)
270 272
271 273 self.selectChannelsByIndex(channelIndexList)
272 274
273 275 def selectChannelsByIndex(self, channelIndexList):
274 276 """
275 277 Selecciona un bloque de datos en base a canales segun el channelIndexList
276 278
277 279 Input:
278 280 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
279 281
280 282 Affected:
281 283 self.dataOut.data_spc
282 284 self.dataOut.channelIndexList
283 285 self.dataOut.nChannels
284 286
285 287 Return:
286 288 None
287 289 """
288 290
289 291 for channelIndex in channelIndexList:
290 292 if channelIndex not in self.dataOut.channelIndexList:
291 293 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " % (
292 294 channelIndex, self.dataOut.channelIndexList)
293 295
294 296 # nChannels = len(channelIndexList)
295 297
296 298 data_spc = self.dataOut.data_spc[channelIndexList, :]
297 299 data_dc = self.dataOut.data_dc[channelIndexList, :]
298 300
299 301 self.dataOut.data_spc = data_spc
300 302 self.dataOut.data_dc = data_dc
301 303
302 304 self.dataOut.channelList = [
303 305 self.dataOut.channelList[i] for i in channelIndexList]
304 306 # self.dataOut.nChannels = nChannels
305 307
306 308 self.__selectPairsByChannel(self.dataOut.channelList)
307 309
308 310 return 1
309 311
310 312 def selectHeights(self, minHei, maxHei):
311 313 """
312 314 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
313 315 minHei <= height <= maxHei
314 316
315 317 Input:
316 318 minHei : valor minimo de altura a considerar
317 319 maxHei : valor maximo de altura a considerar
318 320
319 321 Affected:
320 322 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
321 323
322 324 Return:
323 325 1 si el metodo se ejecuto con exito caso contrario devuelve 0
324 326 """
325 327
326 328 if (minHei > maxHei):
327 329 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (
328 330 minHei, maxHei)
329 331
330 332 if (minHei < self.dataOut.heightList[0]):
331 333 minHei = self.dataOut.heightList[0]
332 334
333 335 if (maxHei > self.dataOut.heightList[-1]):
334 336 maxHei = self.dataOut.heightList[-1]
335 337
336 338 minIndex = 0
337 339 maxIndex = 0
338 340 heights = self.dataOut.heightList
339 341
340 342 inda = numpy.where(heights >= minHei)
341 343 indb = numpy.where(heights <= maxHei)
342 344
343 345 try:
344 346 minIndex = inda[0][0]
345 347 except:
346 348 minIndex = 0
347 349
348 350 try:
349 351 maxIndex = indb[0][-1]
350 352 except:
351 353 maxIndex = len(heights)
352 354
353 355 self.selectHeightsByIndex(minIndex, maxIndex)
354 356
355 357 return 1
356 358
357 359 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
358 360 newheis = numpy.where(
359 361 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
360 362
361 363 if hei_ref != None:
362 364 newheis = numpy.where(self.dataOut.heightList > hei_ref)
363 365
364 366 minIndex = min(newheis[0])
365 367 maxIndex = max(newheis[0])
366 368 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
367 369 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
368 370
369 371 # determina indices
370 372 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
371 373 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
372 374 avg_dB = 10 * \
373 375 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
374 376 beacon_dB = numpy.sort(avg_dB)[-nheis:]
375 377 beacon_heiIndexList = []
376 378 for val in avg_dB.tolist():
377 379 if val >= beacon_dB[0]:
378 380 beacon_heiIndexList.append(avg_dB.tolist().index(val))
379 381
380 382 #data_spc = data_spc[:,:,beacon_heiIndexList]
381 383 data_cspc = None
382 384 if self.dataOut.data_cspc is not None:
383 385 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
384 386 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
385 387
386 388 data_dc = None
387 389 if self.dataOut.data_dc is not None:
388 390 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
389 391 #data_dc = data_dc[:,beacon_heiIndexList]
390 392
391 393 self.dataOut.data_spc = data_spc
392 394 self.dataOut.data_cspc = data_cspc
393 395 self.dataOut.data_dc = data_dc
394 396 self.dataOut.heightList = heightList
395 397 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
396 398
397 399 return 1
398 400
399 401 def selectHeightsByIndex(self, minIndex, maxIndex):
400 402 """
401 403 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
402 404 minIndex <= index <= maxIndex
403 405
404 406 Input:
405 407 minIndex : valor de indice minimo de altura a considerar
406 408 maxIndex : valor de indice maximo de altura a considerar
407 409
408 410 Affected:
409 411 self.dataOut.data_spc
410 412 self.dataOut.data_cspc
411 413 self.dataOut.data_dc
412 414 self.dataOut.heightList
413 415
414 416 Return:
415 417 1 si el metodo se ejecuto con exito caso contrario devuelve 0
416 418 """
417 419
418 420 if (minIndex < 0) or (minIndex > maxIndex):
419 421 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (
420 422 minIndex, maxIndex)
421 423
422 424 if (maxIndex >= self.dataOut.nHeights):
423 425 maxIndex = self.dataOut.nHeights - 1
424 426
425 427 # Spectra
426 428 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
427 429
428 430 data_cspc = None
429 431 if self.dataOut.data_cspc is not None:
430 432 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
431 433
432 434 data_dc = None
433 435 if self.dataOut.data_dc is not None:
434 436 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
435 437
436 438 self.dataOut.data_spc = data_spc
437 439 self.dataOut.data_cspc = data_cspc
438 440 self.dataOut.data_dc = data_dc
439 441
440 442 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
441 443
442 444 return 1
443 445
444 446 def removeDC(self, mode=2):
445 447 jspectra = self.dataOut.data_spc
446 448 jcspectra = self.dataOut.data_cspc
447 449
448 450 num_chan = jspectra.shape[0]
449 451 num_hei = jspectra.shape[2]
450 452
451 453 if jcspectra is not None:
452 454 jcspectraExist = True
453 455 num_pairs = jcspectra.shape[0]
454 456 else:
455 457 jcspectraExist = False
456 458
457 459 freq_dc = jspectra.shape[1] / 2
458 460 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
459 461
460 462 if ind_vel[0] < 0:
461 463 ind_vel[range(0, 1)] = ind_vel[range(0, 1)] + self.num_prof
462 464
463 465 if mode == 1:
464 466 jspectra[:, freq_dc, :] = (
465 467 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
466 468
467 469 if jcspectraExist:
468 470 jcspectra[:, freq_dc, :] = (
469 471 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
470 472
471 473 if mode == 2:
472 474
473 475 vel = numpy.array([-2, -1, 1, 2])
474 476 xx = numpy.zeros([4, 4])
475 477
476 478 for fil in range(4):
477 479 xx[fil, :] = vel[fil]**numpy.asarray(range(4))
478 480
479 481 xx_inv = numpy.linalg.inv(xx)
480 482 xx_aux = xx_inv[0, :]
481 483
482 484 for ich in range(num_chan):
483 485 yy = jspectra[ich, ind_vel, :]
484 486 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
485 487
486 488 junkid = jspectra[ich, freq_dc, :] <= 0
487 489 cjunkid = sum(junkid)
488 490
489 491 if cjunkid.any():
490 492 jspectra[ich, freq_dc, junkid.nonzero()] = (
491 493 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
492 494
493 495 if jcspectraExist:
494 496 for ip in range(num_pairs):
495 497 yy = jcspectra[ip, ind_vel, :]
496 498 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
497 499
498 500 self.dataOut.data_spc = jspectra
499 501 self.dataOut.data_cspc = jcspectra
500 502
501 503 return 1
502 504
503 505 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
504 506
505 507 jspectra = self.dataOut.data_spc
506 508 jcspectra = self.dataOut.data_cspc
507 509 jnoise = self.dataOut.getNoise()
508 510 num_incoh = self.dataOut.nIncohInt
509 511
510 512 num_channel = jspectra.shape[0]
511 513 num_prof = jspectra.shape[1]
512 514 num_hei = jspectra.shape[2]
513 515
514 516 # hei_interf
515 517 if hei_interf is None:
516 518 count_hei = num_hei / 2 # Como es entero no importa
517 519 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
518 520 hei_interf = numpy.asarray(hei_interf)[0]
519 521 # nhei_interf
520 522 if (nhei_interf == None):
521 523 nhei_interf = 5
522 524 if (nhei_interf < 1):
523 525 nhei_interf = 1
524 526 if (nhei_interf > count_hei):
525 527 nhei_interf = count_hei
526 528 if (offhei_interf == None):
527 529 offhei_interf = 0
528 530
529 531 ind_hei = range(num_hei)
530 532 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
531 533 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
532 534 mask_prof = numpy.asarray(range(num_prof))
533 535 num_mask_prof = mask_prof.size
534 536 comp_mask_prof = [0, num_prof / 2]
535 537
536 538 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
537 539 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
538 540 jnoise = numpy.nan
539 541 noise_exist = jnoise[0] < numpy.Inf
540 542
541 543 # Subrutina de Remocion de la Interferencia
542 544 for ich in range(num_channel):
543 545 # Se ordena los espectros segun su potencia (menor a mayor)
544 546 power = jspectra[ich, mask_prof, :]
545 547 power = power[:, hei_interf]
546 548 power = power.sum(axis=0)
547 549 psort = power.ravel().argsort()
548 550
549 551 # Se estima la interferencia promedio en los Espectros de Potencia empleando
550 552 junkspc_interf = jspectra[ich, :, hei_interf[psort[range(
551 553 offhei_interf, nhei_interf + offhei_interf)]]]
552 554
553 555 if noise_exist:
554 556 # tmp_noise = jnoise[ich] / num_prof
555 557 tmp_noise = jnoise[ich]
556 558 junkspc_interf = junkspc_interf - tmp_noise
557 559 #junkspc_interf[:,comp_mask_prof] = 0
558 560
559 561 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
560 562 jspc_interf = jspc_interf.transpose()
561 563 # Calculando el espectro de interferencia promedio
562 564 noiseid = numpy.where(
563 565 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
564 566 noiseid = noiseid[0]
565 567 cnoiseid = noiseid.size
566 568 interfid = numpy.where(
567 569 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
568 570 interfid = interfid[0]
569 571 cinterfid = interfid.size
570 572
571 573 if (cnoiseid > 0):
572 574 jspc_interf[noiseid] = 0
573 575
574 576 # Expandiendo los perfiles a limpiar
575 577 if (cinterfid > 0):
576 578 new_interfid = (
577 579 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
578 580 new_interfid = numpy.asarray(new_interfid)
579 581 new_interfid = {x for x in new_interfid}
580 582 new_interfid = numpy.array(list(new_interfid))
581 583 new_cinterfid = new_interfid.size
582 584 else:
583 585 new_cinterfid = 0
584 586
585 587 for ip in range(new_cinterfid):
586 588 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
587 589 jspc_interf[new_interfid[ip]
588 590 ] = junkspc_interf[ind[nhei_interf / 2], new_interfid[ip]]
589 591
590 592 jspectra[ich, :, ind_hei] = jspectra[ich, :,
591 593 ind_hei] - jspc_interf # Corregir indices
592 594
593 595 # Removiendo la interferencia del punto de mayor interferencia
594 596 ListAux = jspc_interf[mask_prof].tolist()
595 597 maxid = ListAux.index(max(ListAux))
596 598
597 599 if cinterfid > 0:
598 600 for ip in range(cinterfid * (interf == 2) - 1):
599 601 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
600 602 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
601 603 cind = len(ind)
602 604
603 605 if (cind > 0):
604 606 jspectra[ich, interfid[ip], ind] = tmp_noise * \
605 607 (1 + (numpy.random.uniform(cind) - 0.5) /
606 608 numpy.sqrt(num_incoh))
607 609
608 610 ind = numpy.array([-2, -1, 1, 2])
609 611 xx = numpy.zeros([4, 4])
610 612
611 613 for id1 in range(4):
612 614 xx[:, id1] = ind[id1]**numpy.asarray(range(4))
613 615
614 616 xx_inv = numpy.linalg.inv(xx)
615 617 xx = xx_inv[:, 0]
616 618 ind = (ind + maxid + num_mask_prof) % num_mask_prof
617 619 yy = jspectra[ich, mask_prof[ind], :]
618 620 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
619 621 yy.transpose(), xx)
620 622
621 623 indAux = (jspectra[ich, :, :] < tmp_noise *
622 624 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
623 625 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
624 626 (1 - 1 / numpy.sqrt(num_incoh))
625 627
626 628 # Remocion de Interferencia en el Cross Spectra
627 629 if jcspectra is None:
628 630 return jspectra, jcspectra
629 631 num_pairs = jcspectra.size / (num_prof * num_hei)
630 632 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
631 633
632 634 for ip in range(num_pairs):
633 635
634 636 #-------------------------------------------
635 637
636 638 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
637 639 cspower = cspower[:, hei_interf]
638 640 cspower = cspower.sum(axis=0)
639 641
640 642 cspsort = cspower.ravel().argsort()
641 643 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[range(
642 644 offhei_interf, nhei_interf + offhei_interf)]]]
643 645 junkcspc_interf = junkcspc_interf.transpose()
644 646 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
645 647
646 648 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
647 649
648 650 median_real = numpy.median(numpy.real(
649 651 junkcspc_interf[mask_prof[ind[range(3 * num_prof / 4)]], :]))
650 652 median_imag = numpy.median(numpy.imag(
651 653 junkcspc_interf[mask_prof[ind[range(3 * num_prof / 4)]], :]))
652 654 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
653 655 median_real, median_imag)
654 656
655 657 for iprof in range(num_prof):
656 658 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
657 659 jcspc_interf[iprof] = junkcspc_interf[iprof,
658 660 ind[nhei_interf / 2]]
659 661
660 662 # Removiendo la Interferencia
661 663 jcspectra[ip, :, ind_hei] = jcspectra[ip,
662 664 :, ind_hei] - jcspc_interf
663 665
664 666 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
665 667 maxid = ListAux.index(max(ListAux))
666 668
667 669 ind = numpy.array([-2, -1, 1, 2])
668 670 xx = numpy.zeros([4, 4])
669 671
670 672 for id1 in range(4):
671 673 xx[:, id1] = ind[id1]**numpy.asarray(range(4))
672 674
673 675 xx_inv = numpy.linalg.inv(xx)
674 676 xx = xx_inv[:, 0]
675 677
676 678 ind = (ind + maxid + num_mask_prof) % num_mask_prof
677 679 yy = jcspectra[ip, mask_prof[ind], :]
678 680 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
679 681
680 682 # Guardar Resultados
681 683 self.dataOut.data_spc = jspectra
682 684 self.dataOut.data_cspc = jcspectra
683 685
684 686 return 1
685 687
686 688 def setRadarFrequency(self, frequency=None):
687 689
688 690 if frequency != None:
689 691 self.dataOut.frequency = frequency
690 692
691 693 return 1
692 694
693 695 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
694 696 # validacion de rango
695 697 if minHei == None:
696 698 minHei = self.dataOut.heightList[0]
697 699
698 700 if maxHei == None:
699 701 maxHei = self.dataOut.heightList[-1]
700 702
701 703 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
702 704 print 'minHei: %.2f is out of the heights range' % (minHei)
703 705 print 'minHei is setting to %.2f' % (self.dataOut.heightList[0])
704 706 minHei = self.dataOut.heightList[0]
705 707
706 708 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
707 709 print 'maxHei: %.2f is out of the heights range' % (maxHei)
708 710 print 'maxHei is setting to %.2f' % (self.dataOut.heightList[-1])
709 711 maxHei = self.dataOut.heightList[-1]
710 712
711 713 # validacion de velocidades
712 714 velrange = self.dataOut.getVelRange(1)
713 715
714 716 if minVel == None:
715 717 minVel = velrange[0]
716 718
717 719 if maxVel == None:
718 720 maxVel = velrange[-1]
719 721
720 722 if (minVel < velrange[0]) or (minVel > maxVel):
721 723 print 'minVel: %.2f is out of the velocity range' % (minVel)
722 724 print 'minVel is setting to %.2f' % (velrange[0])
723 725 minVel = velrange[0]
724 726
725 727 if (maxVel > velrange[-1]) or (maxVel < minVel):
726 728 print 'maxVel: %.2f is out of the velocity range' % (maxVel)
727 729 print 'maxVel is setting to %.2f' % (velrange[-1])
728 730 maxVel = velrange[-1]
729 731
730 732 # seleccion de indices para rango
731 733 minIndex = 0
732 734 maxIndex = 0
733 735 heights = self.dataOut.heightList
734 736
735 737 inda = numpy.where(heights >= minHei)
736 738 indb = numpy.where(heights <= maxHei)
737 739
738 740 try:
739 741 minIndex = inda[0][0]
740 742 except:
741 743 minIndex = 0
742 744
743 745 try:
744 746 maxIndex = indb[0][-1]
745 747 except:
746 748 maxIndex = len(heights)
747 749
748 750 if (minIndex < 0) or (minIndex > maxIndex):
749 751 raise ValueError, "some value in (%d,%d) is not valid" % (
750 752 minIndex, maxIndex)
751 753
752 754 if (maxIndex >= self.dataOut.nHeights):
753 755 maxIndex = self.dataOut.nHeights - 1
754 756
755 757 # seleccion de indices para velocidades
756 758 indminvel = numpy.where(velrange >= minVel)
757 759 indmaxvel = numpy.where(velrange <= maxVel)
758 760 try:
759 761 minIndexVel = indminvel[0][0]
760 762 except:
761 763 minIndexVel = 0
762 764
763 765 try:
764 766 maxIndexVel = indmaxvel[0][-1]
765 767 except:
766 768 maxIndexVel = len(velrange)
767 769
768 770 # seleccion del espectro
769 771 data_spc = self.dataOut.data_spc[:,
770 772 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
771 773 # estimacion de ruido
772 774 noise = numpy.zeros(self.dataOut.nChannels)
773 775
774 776 for channel in range(self.dataOut.nChannels):
775 777 daux = data_spc[channel, :, :]
776 778 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
777 779
778 780 self.dataOut.noise_estimation = noise.copy()
779 781
780 782 return 1
781 783
782 784
783 785 class IncohInt(Operation):
784 786
785 787 __profIndex = 0
786 788 __withOverapping = False
787 789
788 790 __byTime = False
789 791 __initime = None
790 792 __lastdatatime = None
791 793 __integrationtime = None
792 794
793 795 __buffer_spc = None
794 796 __buffer_cspc = None
795 797 __buffer_dc = None
796 798
797 799 __dataReady = False
798 800
799 801 __timeInterval = None
800 802
801 803 n = None
802 804
803 805 def __init__(self, **kwargs):
804 806
805 807 Operation.__init__(self, **kwargs)
806 808 # self.isConfig = False
807 809
808 810 def setup(self, n=None, timeInterval=None, overlapping=False):
809 811 """
810 812 Set the parameters of the integration class.
811 813
812 814 Inputs:
813 815
814 816 n : Number of coherent integrations
815 817 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
816 818 overlapping :
817 819
818 820 """
819 821
820 822 self.__initime = None
821 823 self.__lastdatatime = 0
822 824
823 825 self.__buffer_spc = 0
824 826 self.__buffer_cspc = 0
825 827 self.__buffer_dc = 0
826 828
827 829 self.__profIndex = 0
828 830 self.__dataReady = False
829 831 self.__byTime = False
830 832
831 833 if n is None and timeInterval is None:
832 834 raise ValueError, "n or timeInterval should be specified ..."
833 835
834 836 if n is not None:
835 837 self.n = int(n)
836 838 else:
837 839 # if (type(timeInterval)!=integer) -> change this line
838 840 self.__integrationtime = int(timeInterval)
839 841 self.n = None
840 842 self.__byTime = True
841 843
842 844 def putData(self, data_spc, data_cspc, data_dc):
843 845 """
844 846 Add a profile to the __buffer_spc and increase in one the __profileIndex
845 847
846 848 """
847 849
848 850 self.__buffer_spc += data_spc
849 851
850 852 if data_cspc is None:
851 853 self.__buffer_cspc = None
852 854 else:
853 855 self.__buffer_cspc += data_cspc
854 856
855 857 if data_dc is None:
856 858 self.__buffer_dc = None
857 859 else:
858 860 self.__buffer_dc += data_dc
859 861
860 862 self.__profIndex += 1
861 863
862 864 return
863 865
864 866 def pushData(self):
865 867 """
866 868 Return the sum of the last profiles and the profiles used in the sum.
867 869
868 870 Affected:
869 871
870 872 self.__profileIndex
871 873
872 874 """
873 875
874 876 data_spc = self.__buffer_spc
875 877 data_cspc = self.__buffer_cspc
876 878 data_dc = self.__buffer_dc
877 879 n = self.__profIndex
878 880
879 881 self.__buffer_spc = 0
880 882 self.__buffer_cspc = 0
881 883 self.__buffer_dc = 0
882 884 self.__profIndex = 0
883 885
884 886 return data_spc, data_cspc, data_dc, n
885 887
886 888 def byProfiles(self, *args):
887 889
888 890 self.__dataReady = False
889 891 avgdata_spc = None
890 892 avgdata_cspc = None
891 893 avgdata_dc = None
892 894
893 895 self.putData(*args)
894 896
895 897 if self.__profIndex == self.n:
896 898
897 899 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
898 900 self.n = n
899 901 self.__dataReady = True
900 902
901 903 return avgdata_spc, avgdata_cspc, avgdata_dc
902 904
903 905 def byTime(self, datatime, *args):
904 906
905 907 self.__dataReady = False
906 908 avgdata_spc = None
907 909 avgdata_cspc = None
908 910 avgdata_dc = None
909 911
910 912 self.putData(*args)
911 913
912 914 if (datatime - self.__initime) >= self.__integrationtime:
913 915 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
914 916 self.n = n
915 917 self.__dataReady = True
916 918
917 919 return avgdata_spc, avgdata_cspc, avgdata_dc
918 920
919 921 def integrate(self, datatime, *args):
920 922
921 923 if self.__profIndex == 0:
922 924 self.__initime = datatime
923 925
924 926 if self.__byTime:
925 927 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
926 928 datatime, *args)
927 929 else:
928 930 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
929 931
930 932 if not self.__dataReady:
931 933 return None, None, None, None
932 934
933 935 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
934 936
935 937 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
936 938 if n == 1:
937 939 return
938 940
939 941 dataOut.flagNoData = True
940 942
941 943 if not self.isConfig:
942 944 self.setup(n, timeInterval, overlapping)
943 945 self.isConfig = True
944 946
945 947 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
946 948 dataOut.data_spc,
947 949 dataOut.data_cspc,
948 950 dataOut.data_dc)
949 951
950 952 if self.__dataReady:
951 953
952 954 dataOut.data_spc = avgdata_spc
953 955 dataOut.data_cspc = avgdata_cspc
954 956 dataOut.data_dc = avgdata_dc
955 957
956 958 dataOut.nIncohInt *= self.n
957 959 dataOut.utctime = avgdatatime
958 960 dataOut.flagNoData = False
Show file before | Show file after | Hide comments
@@ -1,736 +1,757
1 1 import numpy
2 2
3 3 from jroproc_base import ProcessingUnit, Operation
4 4 from schainpy.model.data.jrodata import Spectra
5 5 from schainpy.model.data.jrodata import hildebrand_sekhon
6 6
7 7 class SpectraAFCProc(ProcessingUnit):
8 8
9 9 def __init__(self, **kwargs):
10 10
11 11 ProcessingUnit.__init__(self, **kwargs)
12 12
13 13 self.buffer = None
14 14 self.firstdatatime = None
15 15 self.profIndex = 0
16 16 self.dataOut = Spectra()
17 17 self.id_min = None
18 18 self.id_max = None
19 19
20 20 def __updateSpecFromVoltage(self):
21 21
22 22 self.dataOut.plotting = "spectra_acf"
23 23
24 24 self.dataOut.timeZone = self.dataIn.timeZone
25 25 self.dataOut.dstFlag = self.dataIn.dstFlag
26 26 self.dataOut.errorCount = self.dataIn.errorCount
27 27 self.dataOut.useLocalTime = self.dataIn.useLocalTime
28 28
29 29 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
30 30 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
31 31 self.dataOut.ippSeconds = self.dataIn.getDeltaH()*(10**-6)/0.15
32 32
33 33 self.dataOut.channelList = self.dataIn.channelList
34 34 self.dataOut.heightList = self.dataIn.heightList
35 35 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
36 36
37 37 self.dataOut.nBaud = self.dataIn.nBaud
38 38 self.dataOut.nCode = self.dataIn.nCode
39 39 self.dataOut.code = self.dataIn.code
40 40 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
41 41
42 42 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
43 43 self.dataOut.utctime = self.firstdatatime
44 44 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
45 45 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
46 46 self.dataOut.flagShiftFFT = False
47 47
48 48 self.dataOut.nCohInt = self.dataIn.nCohInt
49 49 self.dataOut.nIncohInt = 1
50 50
51 51 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
52 52
53 53 self.dataOut.frequency = self.dataIn.frequency
54 54 self.dataOut.realtime = self.dataIn.realtime
55 55
56 56 self.dataOut.azimuth = self.dataIn.azimuth
57 57 self.dataOut.zenith = self.dataIn.zenith
58 58
59 59 self.dataOut.beam.codeList = self.dataIn.beam.codeList
60 60 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
61 61 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
62 62
63 63 def __decodeData(self, nProfiles, code):
64 64
65 65 if code is None:
66 66 return
67 67
68 68 for i in range(nProfiles):
69 69 self.buffer[:,i,:] = self.buffer[:,i,:]*code[0][i]
70 70
71 71 def __getFft(self):
72 72 """
73 73 Convierte valores de Voltaje a Spectra
74 74
75 75 Affected:
76 76 self.dataOut.data_spc
77 77 self.dataOut.data_cspc
78 78 self.dataOut.data_dc
79 79 self.dataOut.heightList
80 80 self.profIndex
81 81 self.buffer
82 82 self.dataOut.flagNoData
83 83 """
84 84 nsegments = self.dataOut.nHeights
85 85
86 86 _fft_buffer = numpy.zeros((self.dataOut.nChannels, self.dataOut.nProfiles, nsegments), dtype='complex')
87 87
88 88 for i in range(nsegments):
89 89 try:
90 90 _fft_buffer[:,:,i] = self.buffer[:,i:i+self.dataOut.nProfiles]
91 91
92 92 if self.code is not None:
93 93 _fft_buffer[:,:,i] = _fft_buffer[:,:,i]*self.code[0]
94 94 except:
95 95 pass
96 96
97 97 fft_volt = numpy.fft.fft(_fft_buffer, n=self.dataOut.nFFTPoints, axis=1)
98 98 fft_volt = fft_volt.astype(numpy.dtype('complex'))
99 99 dc = fft_volt[:,0,:]
100 100
101 101 #calculo de self-spectra
102 102 # fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
103 103 spc = fft_volt * numpy.conjugate(fft_volt)
104 104
105
105 106 data = numpy.fft.ifft(spc, axis=1)
106 107 data = numpy.fft.fftshift(data, axes=(1,))
107 108
108 spc = data
109 spc = data.real
110
111
109 112
110 113 blocksize = 0
111 114 blocksize += dc.size
112 115 blocksize += spc.size
113 116
114 117 cspc = None
115 118 pairIndex = 0
116 119
117 120 if self.dataOut.pairsList != None:
118 121 #calculo de cross-spectra
119 122 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
120 123 for pair in self.dataOut.pairsList:
121 124 if pair[0] not in self.dataOut.channelList:
122 125 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
123 126 if pair[1] not in self.dataOut.channelList:
124 127 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
125 128
126 129 chan_index0 = self.dataOut.channelList.index(pair[0])
127 130 chan_index1 = self.dataOut.channelList.index(pair[1])
128 131
129 132 cspc[pairIndex,:,:] = fft_volt[chan_index0,:,:] * numpy.conjugate(fft_volt[chan_index1,:,:])
130 133 pairIndex += 1
131 134 blocksize += cspc.size
132 135
133 136 self.dataOut.data_spc = spc
134 137 self.dataOut.data_cspc = cspc
135 138 self.dataOut.data_dc = dc
136 139 self.dataOut.blockSize = blocksize
137 140 self.dataOut.flagShiftFFT = True
138 141
139 142 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], code=None, nCode=1, nBaud=1):
140 143
141 144 self.dataOut.flagNoData = True
142 145
146 if self.dataIn.type == "Spectra":
147 self.dataOut.copy(self.dataIn)
148 spc= self.dataOut.data_spc
149 data = numpy.fft.ifft(spc, axis=1)
150 data = numpy.fft.fftshift(data, axes=(1,))
151 spc = data.real
152 spc = spc[0,:,0] / numpy.max(numpy.abs(spc[0,:,0]))
153 print spc
154 import matplotlib.pyplot as plt
155 #plt.plot(spc[10:])
156 plt.show()
157
158
159 self.dataOut.data_spc = spc
160
161 return True
162
163
143 164 if code is not None:
144 165 self.code = numpy.array(code).reshape(nCode,nBaud)
145 166 else:
146 167 self.code = None
147 168
148 169 if self.dataIn.type == "Voltage":
149 170
150 171 if nFFTPoints == None:
151 172 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
152 173
153 174 if nProfiles == None:
154 175 nProfiles = nFFTPoints
155 176
156 177 self.dataOut.ippFactor = 1
157 178
158 179 self.dataOut.nFFTPoints = nFFTPoints
159 180 self.dataOut.nProfiles = nProfiles
160 181 self.dataOut.pairsList = pairsList
161 182
162 183 # if self.buffer is None:
163 184 # self.buffer = numpy.zeros( (self.dataIn.nChannels, nProfiles, self.dataIn.nHeights),
164 185 # dtype='complex')
165 186
166 187 if not self.dataIn.flagDataAsBlock:
167 188 self.buffer = self.dataIn.data.copy()
168 189
169 190 # for i in range(self.dataIn.nHeights):
170 191 # self.buffer[:, self.profIndex, self.profIndex:] = voltage_data[:,:self.dataIn.nHeights - self.profIndex]
171 192 #
172 193 # self.profIndex += 1
173 194
174 195 else:
175 196 raise ValueError, ""
176 197
177 198 self.firstdatatime = self.dataIn.utctime
178 199
179 200 self.profIndex == nProfiles
180 201
181 202 self.__updateSpecFromVoltage()
182 203
183 204 self.__getFft()
184 205
185 206 self.dataOut.flagNoData = False
186 207
187 208 return True
188 209
189 210 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
190 211
191 212 def __selectPairs(self, pairsList):
192 213
193 214 if channelList == None:
194 215 return
195 216
196 217 pairsIndexListSelected = []
197 218
198 219 for thisPair in pairsList:
199 220
200 221 if thisPair not in self.dataOut.pairsList:
201 222 continue
202 223
203 224 pairIndex = self.dataOut.pairsList.index(thisPair)
204 225
205 226 pairsIndexListSelected.append(pairIndex)
206 227
207 228 if not pairsIndexListSelected:
208 229 self.dataOut.data_cspc = None
209 230 self.dataOut.pairsList = []
210 231 return
211 232
212 233 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
213 234 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
214 235
215 236 return
216 237
217 238 def __selectPairsByChannel(self, channelList=None):
218 239
219 240 if channelList == None:
220 241 return
221 242
222 243 pairsIndexListSelected = []
223 244 for pairIndex in self.dataOut.pairsIndexList:
224 245 #First pair
225 246 if self.dataOut.pairsList[pairIndex][0] not in channelList:
226 247 continue
227 248 #Second pair
228 249 if self.dataOut.pairsList[pairIndex][1] not in channelList:
229 250 continue
230 251
231 252 pairsIndexListSelected.append(pairIndex)
232 253
233 254 if not pairsIndexListSelected:
234 255 self.dataOut.data_cspc = None
235 256 self.dataOut.pairsList = []
236 257 return
237 258
238 259 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
239 260 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
240 261
241 262 return
242 263
243 264 def selectChannels(self, channelList):
244 265
245 266 channelIndexList = []
246 267
247 268 for channel in channelList:
248 269 if channel not in self.dataOut.channelList:
249 270 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
250 271
251 272 index = self.dataOut.channelList.index(channel)
252 273 channelIndexList.append(index)
253 274
254 275 self.selectChannelsByIndex(channelIndexList)
255 276
256 277 def selectChannelsByIndex(self, channelIndexList):
257 278 """
258 279 Selecciona un bloque de datos en base a canales segun el channelIndexList
259 280
260 281 Input:
261 282 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
262 283
263 284 Affected:
264 285 self.dataOut.data_spc
265 286 self.dataOut.channelIndexList
266 287 self.dataOut.nChannels
267 288
268 289 Return:
269 290 None
270 291 """
271 292
272 293 for channelIndex in channelIndexList:
273 294 if channelIndex not in self.dataOut.channelIndexList:
274 295 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
275 296
276 297 # nChannels = len(channelIndexList)
277 298
278 299 data_spc = self.dataOut.data_spc[channelIndexList,:]
279 300 data_dc = self.dataOut.data_dc[channelIndexList,:]
280 301
281 302 self.dataOut.data_spc = data_spc
282 303 self.dataOut.data_dc = data_dc
283 304
284 305 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
285 306 # self.dataOut.nChannels = nChannels
286 307
287 308 self.__selectPairsByChannel(self.dataOut.channelList)
288 309
289 310 return 1
290 311
291 312 def selectHeights(self, minHei, maxHei):
292 313 """
293 314 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
294 315 minHei <= height <= maxHei
295 316
296 317 Input:
297 318 minHei : valor minimo de altura a considerar
298 319 maxHei : valor maximo de altura a considerar
299 320
300 321 Affected:
301 322 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
302 323
303 324 Return:
304 325 1 si el metodo se ejecuto con exito caso contrario devuelve 0
305 326 """
306 327
307 328 if (minHei > maxHei):
308 329 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
309 330
310 331 if (minHei < self.dataOut.heightList[0]):
311 332 minHei = self.dataOut.heightList[0]
312 333
313 334 if (maxHei > self.dataOut.heightList[-1]):
314 335 maxHei = self.dataOut.heightList[-1]
315 336
316 337 minIndex = 0
317 338 maxIndex = 0
318 339 heights = self.dataOut.heightList
319 340
320 341 inda = numpy.where(heights >= minHei)
321 342 indb = numpy.where(heights <= maxHei)
322 343
323 344 try:
324 345 minIndex = inda[0][0]
325 346 except:
326 347 minIndex = 0
327 348
328 349 try:
329 350 maxIndex = indb[0][-1]
330 351 except:
331 352 maxIndex = len(heights)
332 353
333 354 self.selectHeightsByIndex(minIndex, maxIndex)
334 355
335 356 return 1
336 357
337 358 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
338 359 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
339 360
340 361 if hei_ref != None:
341 362 newheis = numpy.where(self.dataOut.heightList>hei_ref)
342 363
343 364 minIndex = min(newheis[0])
344 365 maxIndex = max(newheis[0])
345 366 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
346 367 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
347 368
348 369 # determina indices
349 370 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
350 371 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
351 372 beacon_dB = numpy.sort(avg_dB)[-nheis:]
352 373 beacon_heiIndexList = []
353 374 for val in avg_dB.tolist():
354 375 if val >= beacon_dB[0]:
355 376 beacon_heiIndexList.append(avg_dB.tolist().index(val))
356 377
357 378 #data_spc = data_spc[:,:,beacon_heiIndexList]
358 379 data_cspc = None
359 380 if self.dataOut.data_cspc is not None:
360 381 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
361 382 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
362 383
363 384 data_dc = None
364 385 if self.dataOut.data_dc is not None:
365 386 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
366 387 #data_dc = data_dc[:,beacon_heiIndexList]
367 388
368 389 self.dataOut.data_spc = data_spc
369 390 self.dataOut.data_cspc = data_cspc
370 391 self.dataOut.data_dc = data_dc
371 392 self.dataOut.heightList = heightList
372 393 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
373 394
374 395 return 1
375 396
376 397
377 398 def selectHeightsByIndex(self, minIndex, maxIndex):
378 399 """
379 400 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
380 401 minIndex <= index <= maxIndex
381 402
382 403 Input:
383 404 minIndex : valor de indice minimo de altura a considerar
384 405 maxIndex : valor de indice maximo de altura a considerar
385 406
386 407 Affected:
387 408 self.dataOut.data_spc
388 409 self.dataOut.data_cspc
389 410 self.dataOut.data_dc
390 411 self.dataOut.heightList
391 412
392 413 Return:
393 414 1 si el metodo se ejecuto con exito caso contrario devuelve 0
394 415 """
395 416
396 417 if (minIndex < 0) or (minIndex > maxIndex):
397 418 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
398 419
399 420 if (maxIndex >= self.dataOut.nHeights):
400 421 maxIndex = self.dataOut.nHeights-1
401 422
402 423 #Spectra
403 424 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
404 425
405 426 data_cspc = None
406 427 if self.dataOut.data_cspc is not None:
407 428 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
408 429
409 430 data_dc = None
410 431 if self.dataOut.data_dc is not None:
411 432 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
412 433
413 434 self.dataOut.data_spc = data_spc
414 435 self.dataOut.data_cspc = data_cspc
415 436 self.dataOut.data_dc = data_dc
416 437
417 438 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
418 439
419 440 return 1
420 441
421 442 def removeDC(self, mode = 2):
422 443 jspectra = self.dataOut.data_spc
423 444 jcspectra = self.dataOut.data_cspc
424 445
425 446
426 447 num_chan = jspectra.shape[0]
427 448 num_hei = jspectra.shape[2]
428 449
429 450 if jcspectra is not None:
430 451 jcspectraExist = True
431 452 num_pairs = jcspectra.shape[0]
432 453 else: jcspectraExist = False
433 454
434 455 freq_dc = jspectra.shape[1]/2
435 456 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
436 457
437 458 if ind_vel[0]<0:
438 459 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
439 460
440 461 if mode == 1:
441 462 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
442 463
443 464 if jcspectraExist:
444 465 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
445 466
446 467 if mode == 2:
447 468
448 469 vel = numpy.array([-2,-1,1,2])
449 470 xx = numpy.zeros([4,4])
450 471
451 472 for fil in range(4):
452 473 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
453 474
454 475 xx_inv = numpy.linalg.inv(xx)
455 476 xx_aux = xx_inv[0,:]
456 477
457 478 for ich in range(num_chan):
458 479 yy = jspectra[ich,ind_vel,:]
459 480 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
460 481
461 482 junkid = jspectra[ich,freq_dc,:]<=0
462 483 cjunkid = sum(junkid)
463 484
464 485 if cjunkid.any():
465 486 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
466 487
467 488 if jcspectraExist:
468 489 for ip in range(num_pairs):
469 490 yy = jcspectra[ip,ind_vel,:]
470 491 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
471 492
472 493
473 494 self.dataOut.data_spc = jspectra
474 495 self.dataOut.data_cspc = jcspectra
475 496
476 497 return 1
477 498
478 499 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
479 500
480 501 jspectra = self.dataOut.data_spc
481 502 jcspectra = self.dataOut.data_cspc
482 503 jnoise = self.dataOut.getNoise()
483 504 num_incoh = self.dataOut.nIncohInt
484 505
485 506 num_channel = jspectra.shape[0]
486 507 num_prof = jspectra.shape[1]
487 508 num_hei = jspectra.shape[2]
488 509
489 510 #hei_interf
490 511 if hei_interf is None:
491 512 count_hei = num_hei/2 #Como es entero no importa
492 513 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
493 514 hei_interf = numpy.asarray(hei_interf)[0]
494 515 #nhei_interf
495 516 if (nhei_interf == None):
496 517 nhei_interf = 5
497 518 if (nhei_interf < 1):
498 519 nhei_interf = 1
499 520 if (nhei_interf > count_hei):
500 521 nhei_interf = count_hei
501 522 if (offhei_interf == None):
502 523 offhei_interf = 0
503 524
504 525 ind_hei = range(num_hei)
505 526 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
506 527 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
507 528 mask_prof = numpy.asarray(range(num_prof))
508 529 num_mask_prof = mask_prof.size
509 530 comp_mask_prof = [0, num_prof/2]
510 531
511 532
512 533 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
513 534 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
514 535 jnoise = numpy.nan
515 536 noise_exist = jnoise[0] < numpy.Inf
516 537
517 538 #Subrutina de Remocion de la Interferencia
518 539 for ich in range(num_channel):
519 540 #Se ordena los espectros segun su potencia (menor a mayor)
520 541 power = jspectra[ich,mask_prof,:]
521 542 power = power[:,hei_interf]
522 543 power = power.sum(axis = 0)
523 544 psort = power.ravel().argsort()
524 545
525 546 #Se estima la interferencia promedio en los Espectros de Potencia empleando
526 547 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
527 548
528 549 if noise_exist:
529 550 # tmp_noise = jnoise[ich] / num_prof
530 551 tmp_noise = jnoise[ich]
531 552 junkspc_interf = junkspc_interf - tmp_noise
532 553 #junkspc_interf[:,comp_mask_prof] = 0
533 554
534 555 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
535 556 jspc_interf = jspc_interf.transpose()
536 557 #Calculando el espectro de interferencia promedio
537 558 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
538 559 noiseid = noiseid[0]
539 560 cnoiseid = noiseid.size
540 561 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
541 562 interfid = interfid[0]
542 563 cinterfid = interfid.size
543 564
544 565 if (cnoiseid > 0): jspc_interf[noiseid] = 0
545 566
546 567 #Expandiendo los perfiles a limpiar
547 568 if (cinterfid > 0):
548 569 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
549 570 new_interfid = numpy.asarray(new_interfid)
550 571 new_interfid = {x for x in new_interfid}
551 572 new_interfid = numpy.array(list(new_interfid))
552 573 new_cinterfid = new_interfid.size
553 574 else: new_cinterfid = 0
554 575
555 576 for ip in range(new_cinterfid):
556 577 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
557 578 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
558 579
559 580
560 581 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
561 582
562 583 #Removiendo la interferencia del punto de mayor interferencia
563 584 ListAux = jspc_interf[mask_prof].tolist()
564 585 maxid = ListAux.index(max(ListAux))
565 586
566 587
567 588 if cinterfid > 0:
568 589 for ip in range(cinterfid*(interf == 2) - 1):
569 590 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
570 591 cind = len(ind)
571 592
572 593 if (cind > 0):
573 594 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
574 595
575 596 ind = numpy.array([-2,-1,1,2])
576 597 xx = numpy.zeros([4,4])
577 598
578 599 for id1 in range(4):
579 600 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
580 601
581 602 xx_inv = numpy.linalg.inv(xx)
582 603 xx = xx_inv[:,0]
583 604 ind = (ind + maxid + num_mask_prof)%num_mask_prof
584 605 yy = jspectra[ich,mask_prof[ind],:]
585 606 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
586 607
587 608
588 609 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
589 610 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
590 611
591 612 #Remocion de Interferencia en el Cross Spectra
592 613 if jcspectra is None: return jspectra, jcspectra
593 614 num_pairs = jcspectra.size/(num_prof*num_hei)
594 615 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
595 616
596 617 for ip in range(num_pairs):
597 618
598 619 #-------------------------------------------
599 620
600 621 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
601 622 cspower = cspower[:,hei_interf]
602 623 cspower = cspower.sum(axis = 0)
603 624
604 625 cspsort = cspower.ravel().argsort()
605 626 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
606 627 junkcspc_interf = junkcspc_interf.transpose()
607 628 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
608 629
609 630 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
610 631
611 632 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
612 633 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
613 634 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
614 635
615 636 for iprof in range(num_prof):
616 637 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
617 638 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
618 639
619 640 #Removiendo la Interferencia
620 641 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
621 642
622 643 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
623 644 maxid = ListAux.index(max(ListAux))
624 645
625 646 ind = numpy.array([-2,-1,1,2])
626 647 xx = numpy.zeros([4,4])
627 648
628 649 for id1 in range(4):
629 650 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
630 651
631 652 xx_inv = numpy.linalg.inv(xx)
632 653 xx = xx_inv[:,0]
633 654
634 655 ind = (ind + maxid + num_mask_prof)%num_mask_prof
635 656 yy = jcspectra[ip,mask_prof[ind],:]
636 657 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
637 658
638 659 #Guardar Resultados
639 660 self.dataOut.data_spc = jspectra
640 661 self.dataOut.data_cspc = jcspectra
641 662
642 663 return 1
643 664
644 665 def setRadarFrequency(self, frequency=None):
645 666
646 667 if frequency != None:
647 668 self.dataOut.frequency = frequency
648 669
649 670 return 1
650 671
651 672 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
652 673 #validacion de rango
653 674 if minHei == None:
654 675 minHei = self.dataOut.heightList[0]
655 676
656 677 if maxHei == None:
657 678 maxHei = self.dataOut.heightList[-1]
658 679
659 680 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
660 681 print 'minHei: %.2f is out of the heights range'%(minHei)
661 682 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
662 683 minHei = self.dataOut.heightList[0]
663 684
664 685 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
665 686 print 'maxHei: %.2f is out of the heights range'%(maxHei)
666 687 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
667 688 maxHei = self.dataOut.heightList[-1]
668 689
669 690 # validacion de velocidades
670 691 velrange = self.dataOut.getVelRange(1)
671 692
672 693 if minVel == None:
673 694 minVel = velrange[0]
674 695
675 696 if maxVel == None:
676 697 maxVel = velrange[-1]
677 698
678 699 if (minVel < velrange[0]) or (minVel > maxVel):
679 700 print 'minVel: %.2f is out of the velocity range'%(minVel)
680 701 print 'minVel is setting to %.2f'%(velrange[0])
681 702 minVel = velrange[0]
682 703
683 704 if (maxVel > velrange[-1]) or (maxVel < minVel):
684 705 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
685 706 print 'maxVel is setting to %.2f'%(velrange[-1])
686 707 maxVel = velrange[-1]
687 708
688 709 # seleccion de indices para rango
689 710 minIndex = 0
690 711 maxIndex = 0
691 712 heights = self.dataOut.heightList
692 713
693 714 inda = numpy.where(heights >= minHei)
694 715 indb = numpy.where(heights <= maxHei)
695 716
696 717 try:
697 718 minIndex = inda[0][0]
698 719 except:
699 720 minIndex = 0
700 721
701 722 try:
702 723 maxIndex = indb[0][-1]
703 724 except:
704 725 maxIndex = len(heights)
705 726
706 727 if (minIndex < 0) or (minIndex > maxIndex):
707 728 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
708 729
709 730 if (maxIndex >= self.dataOut.nHeights):
710 731 maxIndex = self.dataOut.nHeights-1
711 732
712 733 # seleccion de indices para velocidades
713 734 indminvel = numpy.where(velrange >= minVel)
714 735 indmaxvel = numpy.where(velrange <= maxVel)
715 736 try:
716 737 minIndexVel = indminvel[0][0]
717 738 except:
718 739 minIndexVel = 0
719 740
720 741 try:
721 742 maxIndexVel = indmaxvel[0][-1]
722 743 except:
723 744 maxIndexVel = len(velrange)
724 745
725 746 #seleccion del espectro
726 747 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
727 748 #estimacion de ruido
728 749 noise = numpy.zeros(self.dataOut.nChannels)
729 750
730 751 for channel in range(self.dataOut.nChannels):
731 752 daux = data_spc[channel,:,:]
732 753 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
733 754
734 755 self.dataOut.noise_estimation = noise.copy()
735 756
736 757 return 1
Show file before | Show file after | Hide comments
@@ -1,1395 +1,1395
1 1 import sys
2 2 import numpy
3 3 from scipy import interpolate
4 4 from schainpy import cSchain
5 5 from jroproc_base import ProcessingUnit, Operation
6 6 from schainpy.model.data.jrodata import Voltage
7 7 from time import time
8 8
9 9 class VoltageProc(ProcessingUnit):
10 10
11 11
12 12 def __init__(self, **kwargs):
13 13
14 14 ProcessingUnit.__init__(self, **kwargs)
15 15
16 16 # self.objectDict = {}
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19
20 20 def run(self):
21 21 if self.dataIn.type == 'AMISR':
22 22 self.__updateObjFromAmisrInput()
23 23
24 24 if self.dataIn.type == 'Voltage':
25 25 self.dataOut.copy(self.dataIn)
26 26
27 27 # self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 # self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54 #
55 55 # pass#
56 56 #
57 57 # def init(self):
58 58 #
59 59 #
60 60 # if self.dataIn.type == 'AMISR':
61 61 # self.__updateObjFromAmisrInput()
62 62 #
63 63 # if self.dataIn.type == 'Voltage':
64 64 # self.dataOut.copy(self.dataIn)
65 65 # # No necesita copiar en cada init() los atributos de dataIn
66 66 # # la copia deberia hacerse por cada nuevo bloque de datos
67 67
68 68 def selectChannels(self, channelList):
69 69
70 70 channelIndexList = []
71 71
72 72 for channel in channelList:
73 73 if channel not in self.dataOut.channelList:
74 74 raise ValueError, "Channel %d is not in %s" %(channel, str(self.dataOut.channelList))
75 75
76 76 index = self.dataOut.channelList.index(channel)
77 77 channelIndexList.append(index)
78 78
79 79 self.selectChannelsByIndex(channelIndexList)
80 80
81 81 def selectChannelsByIndex(self, channelIndexList):
82 82 """
83 83 Selecciona un bloque de datos en base a canales segun el channelIndexList
84 84
85 85 Input:
86 86 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
87 87
88 88 Affected:
89 89 self.dataOut.data
90 90 self.dataOut.channelIndexList
91 91 self.dataOut.nChannels
92 92 self.dataOut.m_ProcessingHeader.totalSpectra
93 93 self.dataOut.systemHeaderObj.numChannels
94 94 self.dataOut.m_ProcessingHeader.blockSize
95 95
96 96 Return:
97 97 None
98 98 """
99 99
100 100 for channelIndex in channelIndexList:
101 101 if channelIndex not in self.dataOut.channelIndexList:
102 102 print channelIndexList
103 103 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
104 104
105 105 if self.dataOut.flagDataAsBlock:
106 106 """
107 107 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
108 108 """
109 109 data = self.dataOut.data[channelIndexList,:,:]
110 110 else:
111 111 data = self.dataOut.data[channelIndexList,:]
112 112
113 113 self.dataOut.data = data
114 114 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 # self.dataOut.nChannels = nChannels
116 116
117 117 return 1
118 118
119 119 def selectHeights(self, minHei=None, maxHei=None):
120 120 """
121 121 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
122 122 minHei <= height <= maxHei
123 123
124 124 Input:
125 125 minHei : valor minimo de altura a considerar
126 126 maxHei : valor maximo de altura a considerar
127 127
128 128 Affected:
129 129 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
130 130
131 131 Return:
132 132 1 si el metodo se ejecuto con exito caso contrario devuelve 0
133 133 """
134 134
135 135 if minHei == None:
136 136 minHei = self.dataOut.heightList[0]
137 137
138 138 if maxHei == None:
139 139 maxHei = self.dataOut.heightList[-1]
140 140
141 141 if (minHei < self.dataOut.heightList[0]):
142 142 minHei = self.dataOut.heightList[0]
143 143
144 144 if (maxHei > self.dataOut.heightList[-1]):
145 145 maxHei = self.dataOut.heightList[-1]
146 146
147 147 minIndex = 0
148 148 maxIndex = 0
149 149 heights = self.dataOut.heightList
150 150
151 151 inda = numpy.where(heights >= minHei)
152 152 indb = numpy.where(heights <= maxHei)
153 153
154 154 try:
155 155 minIndex = inda[0][0]
156 156 except:
157 157 minIndex = 0
158 158
159 159 try:
160 160 maxIndex = indb[0][-1]
161 161 except:
162 162 maxIndex = len(heights)
163 163
164 164 self.selectHeightsByIndex(minIndex, maxIndex)
165 165
166 166 return 1
167 167
168 168
169 169 def selectHeightsByIndex(self, minIndex, maxIndex):
170 170 """
171 171 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
172 172 minIndex <= index <= maxIndex
173 173
174 174 Input:
175 175 minIndex : valor de indice minimo de altura a considerar
176 176 maxIndex : valor de indice maximo de altura a considerar
177 177
178 178 Affected:
179 179 self.dataOut.data
180 180 self.dataOut.heightList
181 181
182 182 Return:
183 183 1 si el metodo se ejecuto con exito caso contrario devuelve 0
184 184 """
185 185
186 186 if (minIndex < 0) or (minIndex > maxIndex):
187 187 raise ValueError, "Height index range (%d,%d) is not valid" % (minIndex, maxIndex)
188 188
189 189 if (maxIndex >= self.dataOut.nHeights):
190 190 maxIndex = self.dataOut.nHeights
191 191
192 192 #voltage
193 193 if self.dataOut.flagDataAsBlock:
194 194 """
195 195 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
196 196 """
197 197 data = self.dataOut.data[:,:, minIndex:maxIndex]
198 198 else:
199 199 data = self.dataOut.data[:, minIndex:maxIndex]
200 200
201 201 # firstHeight = self.dataOut.heightList[minIndex]
202 202
203 203 self.dataOut.data = data
204 204 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
205 205
206 206 if self.dataOut.nHeights <= 1:
207 207 raise ValueError, "selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights)
208 208
209 209 return 1
210 210
211 211
212 212 def filterByHeights(self, window):
213 213
214 214 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
215 215
216 216 if window == None:
217 217 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
218 218
219 219 newdelta = deltaHeight * window
220 220 r = self.dataOut.nHeights % window
221 221 newheights = (self.dataOut.nHeights-r)/window
222 222
223 223 if newheights <= 1:
224 224 raise ValueError, "filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window)
225 225
226 226 if self.dataOut.flagDataAsBlock:
227 227 """
228 228 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
229 229 """
230 230 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
231 231 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
232 232 buffer = numpy.sum(buffer,3)
233 233
234 234 else:
235 235 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
236 236 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
237 237 buffer = numpy.sum(buffer,2)
238 238
239 239 self.dataOut.data = buffer
240 240 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
241 241 self.dataOut.windowOfFilter = window
242 242
243 243 def setH0(self, h0, deltaHeight = None):
244 244
245 245 if not deltaHeight:
246 246 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
247 247
248 248 nHeights = self.dataOut.nHeights
249 249
250 250 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
251 251
252 252 self.dataOut.heightList = newHeiRange
253 253
254 254 def deFlip(self, channelList = []):
255 255
256 256 data = self.dataOut.data.copy()
257 257
258 258 if self.dataOut.flagDataAsBlock:
259 259 flip = self.flip
260 260 profileList = range(self.dataOut.nProfiles)
261 261
262 262 if not channelList:
263 263 for thisProfile in profileList:
264 264 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
265 265 flip *= -1.0
266 266 else:
267 267 for thisChannel in channelList:
268 268 if thisChannel not in self.dataOut.channelList:
269 269 continue
270 270
271 271 for thisProfile in profileList:
272 272 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
273 273 flip *= -1.0
274 274
275 275 self.flip = flip
276 276
277 277 else:
278 278 if not channelList:
279 279 data[:,:] = data[:,:]*self.flip
280 280 else:
281 281 for thisChannel in channelList:
282 282 if thisChannel not in self.dataOut.channelList:
283 283 continue
284 284
285 285 data[thisChannel,:] = data[thisChannel,:]*self.flip
286 286
287 287 self.flip *= -1.
288 288
289 289 self.dataOut.data = data
290 290
291 291 def setRadarFrequency(self, frequency=None):
292 292
293 293 if frequency != None:
294 294 self.dataOut.frequency = frequency
295 295
296 296 return 1
297 297
298 298 def interpolateHeights(self, topLim, botLim):
299 299 #69 al 72 para julia
300 300 #82-84 para meteoros
301 301 if len(numpy.shape(self.dataOut.data))==2:
302 302 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
303 303 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
304 304 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
305 305 self.dataOut.data[:,botLim:topLim+1] = sampInterp
306 306 else:
307 307 nHeights = self.dataOut.data.shape[2]
308 308 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
309 309 y = self.dataOut.data[:,:,range(botLim)+range(topLim+1,nHeights)]
310 310 f = interpolate.interp1d(x, y, axis = 2)
311 311 xnew = numpy.arange(botLim,topLim+1)
312 312 ynew = f(xnew)
313 313
314 314 self.dataOut.data[:,:,botLim:topLim+1] = ynew
315 315
316 316 # import collections
317 317
318 318 class CohInt(Operation):
319 319
320 320 isConfig = False
321 321 __profIndex = 0
322 322 __byTime = False
323 323 __initime = None
324 324 __lastdatatime = None
325 325 __integrationtime = None
326 326 __buffer = None
327 327 __bufferStride = []
328 328 __dataReady = False
329 329 __profIndexStride = 0
330 330 __dataToPutStride = False
331 331 n = None
332 332
333 333 def __init__(self, **kwargs):
334 334
335 335 Operation.__init__(self, **kwargs)
336 336
337 337 # self.isConfig = False
338 338
339 339 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
340 340 """
341 341 Set the parameters of the integration class.
342 342
343 343 Inputs:
344 344
345 345 n : Number of coherent integrations
346 346 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
347 347 overlapping :
348 348 """
349 349
350 350 self.__initime = None
351 351 self.__lastdatatime = 0
352 352 self.__buffer = None
353 353 self.__dataReady = False
354 354 self.byblock = byblock
355 355 self.stride = stride
356 356
357 357 if n == None and timeInterval == None:
358 358 raise ValueError, "n or timeInterval should be specified ..."
359 359
360 360 if n != None:
361 361 self.n = n
362 362 self.__byTime = False
363 363 else:
364 364 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
365 365 self.n = 9999
366 366 self.__byTime = True
367 367
368 368 if overlapping:
369 369 self.__withOverlapping = True
370 370 self.__buffer = None
371 371 else:
372 372 self.__withOverlapping = False
373 373 self.__buffer = 0
374 374
375 375 self.__profIndex = 0
376 376
377 377 def putData(self, data):
378 378
379 379 """
380 380 Add a profile to the __buffer and increase in one the __profileIndex
381 381
382 382 """
383 383
384 384 if not self.__withOverlapping:
385 385 self.__buffer += data.copy()
386 386 self.__profIndex += 1
387 387 return
388 388
389 389 #Overlapping data
390 390 nChannels, nHeis = data.shape
391 391 data = numpy.reshape(data, (1, nChannels, nHeis))
392 392
393 393 #If the buffer is empty then it takes the data value
394 394 if self.__buffer is None:
395 395 self.__buffer = data
396 396 self.__profIndex += 1
397 397 return
398 398
399 399 #If the buffer length is lower than n then stakcing the data value
400 400 if self.__profIndex < self.n:
401 401 self.__buffer = numpy.vstack((self.__buffer, data))
402 402 self.__profIndex += 1
403 403 return
404 404
405 405 #If the buffer length is equal to n then replacing the last buffer value with the data value
406 406 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
407 407 self.__buffer[self.n-1] = data
408 408 self.__profIndex = self.n
409 409 return
410 410
411 411
412 412 def pushData(self):
413 413 """
414 414 Return the sum of the last profiles and the profiles used in the sum.
415 415
416 416 Affected:
417 417
418 418 self.__profileIndex
419 419
420 420 """
421 421
422 422 if not self.__withOverlapping:
423 423 data = self.__buffer
424 424 n = self.__profIndex
425 425
426 426 self.__buffer = 0
427 427 self.__profIndex = 0
428 428
429 429 return data, n
430 430
431 431 #Integration with Overlapping
432 432 data = numpy.sum(self.__buffer, axis=0)
433 433 # print data
434 434 # raise
435 435 n = self.__profIndex
436 436
437 437 return data, n
438 438
439 439 def byProfiles(self, data):
440 440
441 441 self.__dataReady = False
442 442 avgdata = None
443 443 # n = None
444 444 # print data
445 445 # raise
446 446 self.putData(data)
447 447
448 448 if self.__profIndex == self.n:
449 449 avgdata, n = self.pushData()
450 450 self.__dataReady = True
451 451
452 452 return avgdata
453 453
454 454 def byTime(self, data, datatime):
455 455
456 456 self.__dataReady = False
457 457 avgdata = None
458 458 n = None
459 459
460 460 self.putData(data)
461 461
462 462 if (datatime - self.__initime) >= self.__integrationtime:
463 463 avgdata, n = self.pushData()
464 464 self.n = n
465 465 self.__dataReady = True
466 466
467 467 return avgdata
468 468
469 469 def integrateByStride(self, data, datatime):
470 470 # print data
471 471 if self.__profIndex == 0:
472 472 self.__buffer = [[data.copy(), datatime]]
473 473 else:
474 474 self.__buffer.append([data.copy(),datatime])
475 475 self.__profIndex += 1
476 476 self.__dataReady = False
477 477
478 478 if self.__profIndex == self.n * self.stride :
479 479 self.__dataToPutStride = True
480 480 self.__profIndexStride = 0
481 481 self.__profIndex = 0
482 482 self.__bufferStride = []
483 483 for i in range(self.stride):
484 484 current = self.__buffer[i::self.stride]
485 485 data = numpy.sum([t[0] for t in current], axis=0)
486 486 avgdatatime = numpy.average([t[1] for t in current])
487 487 # print data
488 488 self.__bufferStride.append((data, avgdatatime))
489 489
490 490 if self.__dataToPutStride:
491 491 self.__dataReady = True
492 492 self.__profIndexStride += 1
493 493 if self.__profIndexStride == self.stride:
494 494 self.__dataToPutStride = False
495 495 # print self.__bufferStride[self.__profIndexStride - 1]
496 496 # raise
497 497 return self.__bufferStride[self.__profIndexStride - 1]
498 498
499 499
500 500 return None, None
501 501
502 502 def integrate(self, data, datatime=None):
503 503
504 504 if self.__initime == None:
505 505 self.__initime = datatime
506 506
507 507 if self.__byTime:
508 508 avgdata = self.byTime(data, datatime)
509 509 else:
510 510 avgdata = self.byProfiles(data)
511 511
512 512
513 513 self.__lastdatatime = datatime
514 514
515 515 if avgdata is None:
516 516 return None, None
517 517
518 518 avgdatatime = self.__initime
519 519
520 520 deltatime = datatime - self.__lastdatatime
521 521
522 522 if not self.__withOverlapping:
523 523 self.__initime = datatime
524 524 else:
525 525 self.__initime += deltatime
526 526
527 527 return avgdata, avgdatatime
528 528
529 529 def integrateByBlock(self, dataOut):
530 530
531 531 times = int(dataOut.data.shape[1]/self.n)
532 532 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
533 533
534 534 id_min = 0
535 535 id_max = self.n
536 536
537 537 for i in range(times):
538 538 junk = dataOut.data[:,id_min:id_max,:]
539 539 avgdata[:,i,:] = junk.sum(axis=1)
540 540 id_min += self.n
541 541 id_max += self.n
542 542
543 543 timeInterval = dataOut.ippSeconds*self.n
544 544 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
545 545 self.__dataReady = True
546 546 return avgdata, avgdatatime
547 547
548 548 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
549 549 if not self.isConfig:
550 550 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
551 551 self.isConfig = True
552 552
553 553 if dataOut.flagDataAsBlock:
554 554 """
555 555 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
556 556 """
557 557 avgdata, avgdatatime = self.integrateByBlock(dataOut)
558 558 dataOut.nProfiles /= self.n
559 559 else:
560 560 if stride is None:
561 561 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
562 562 else:
563 563 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
564 564
565 565
566 566 # dataOut.timeInterval *= n
567 567 dataOut.flagNoData = True
568 568
569 569 if self.__dataReady:
570 570 dataOut.data = avgdata
571 571 dataOut.nCohInt *= self.n
572 572 dataOut.utctime = avgdatatime
573 573 # print avgdata, avgdatatime
574 574 # raise
575 575 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
576 576 dataOut.flagNoData = False
577 577
578 578 class Decoder(Operation):
579 579
580 580 isConfig = False
581 581 __profIndex = 0
582 582
583 583 code = None
584 584
585 585 nCode = None
586 586 nBaud = None
587 587
588 588 def __init__(self, **kwargs):
589 589
590 590 Operation.__init__(self, **kwargs)
591 591
592 592 self.times = None
593 593 self.osamp = None
594 594 # self.__setValues = False
595 595 self.isConfig = False
596 596
597 597 def setup(self, code, osamp, dataOut):
598 598
599 599 self.__profIndex = 0
600 600
601 601 self.code = code
602 602
603 603 self.nCode = len(code)
604 604 self.nBaud = len(code[0])
605 605
606 606 if (osamp != None) and (osamp >1):
607 607 self.osamp = osamp
608 608 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
609 609 self.nBaud = self.nBaud*self.osamp
610 610
611 611 self.__nChannels = dataOut.nChannels
612 612 self.__nProfiles = dataOut.nProfiles
613 613 self.__nHeis = dataOut.nHeights
614 614
615 615 if self.__nHeis < self.nBaud:
616 616 raise ValueError, 'Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)
617 617
618 618 #Frequency
619 619 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
620 620
621 621 __codeBuffer[:,0:self.nBaud] = self.code
622 622
623 623 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
624 624
625 625 if dataOut.flagDataAsBlock:
626 626
627 627 self.ndatadec = self.__nHeis #- self.nBaud + 1
628 628
629 629 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
630 630
631 631 else:
632 632
633 633 #Time
634 634 self.ndatadec = self.__nHeis #- self.nBaud + 1
635 635
636 636 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
637 637
638 638 def __convolutionInFreq(self, data):
639 639
640 640 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
641 641
642 642 fft_data = numpy.fft.fft(data, axis=1)
643 643
644 644 conv = fft_data*fft_code
645 645
646 646 data = numpy.fft.ifft(conv,axis=1)
647 647
648 648 return data
649 649
650 650 def __convolutionInFreqOpt(self, data):
651 651
652 652 raise NotImplementedError
653 653
654 654 def __convolutionInTime(self, data):
655 655
656 656 code = self.code[self.__profIndex]
657 657 for i in range(self.__nChannels):
658 658 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
659 659
660 660 return self.datadecTime
661 661
662 662 def __convolutionByBlockInTime(self, data):
663 663
664 664 repetitions = self.__nProfiles / self.nCode
665 665
666 666 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
667 667 junk = junk.flatten()
668 668 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
669 669 profilesList = xrange(self.__nProfiles)
670 670
671 671 for i in range(self.__nChannels):
672 672 for j in profilesList:
673 673 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
674 674 return self.datadecTime
675 675
676 676 def __convolutionByBlockInFreq(self, data):
677 677
678 678 raise NotImplementedError, "Decoder by frequency fro Blocks not implemented"
679 679
680 680
681 681 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
682 682
683 683 fft_data = numpy.fft.fft(data, axis=2)
684 684
685 685 conv = fft_data*fft_code
686 686
687 687 data = numpy.fft.ifft(conv,axis=2)
688 688
689 689 return data
690 690
691 691
692 692 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
693 693
694 694 if dataOut.flagDecodeData:
695 695 print "This data is already decoded, recoding again ..."
696 696
697 697 if not self.isConfig:
698 698
699 699 if code is None:
700 700 if dataOut.code is None:
701 701 raise ValueError, "Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type
702 702
703 703 code = dataOut.code
704 704 else:
705 705 code = numpy.array(code).reshape(nCode,nBaud)
706 706 self.setup(code, osamp, dataOut)
707 707
708 708 self.isConfig = True
709 709
710 710 if mode == 3:
711 711 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
712 712
713 713 if times != None:
714 714 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
715 715
716 716 if self.code is None:
717 717 print "Fail decoding: Code is not defined."
718 718 return
719 719
720 720 self.__nProfiles = dataOut.nProfiles
721 721 datadec = None
722 722
723 723 if mode == 3:
724 724 mode = 0
725 725
726 726 if dataOut.flagDataAsBlock:
727 727 """
728 728 Decoding when data have been read as block,
729 729 """
730 730
731 731 if mode == 0:
732 732 datadec = self.__convolutionByBlockInTime(dataOut.data)
733 733 if mode == 1:
734 734 datadec = self.__convolutionByBlockInFreq(dataOut.data)
735 735 else:
736 736 """
737 737 Decoding when data have been read profile by profile
738 738 """
739 739 if mode == 0:
740 740 datadec = self.__convolutionInTime(dataOut.data)
741 741
742 742 if mode == 1:
743 743 datadec = self.__convolutionInFreq(dataOut.data)
744 744
745 745 if mode == 2:
746 746 datadec = self.__convolutionInFreqOpt(dataOut.data)
747 747
748 748 if datadec is None:
749 749 raise ValueError, "Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode
750 750
751 751 dataOut.code = self.code
752 752 dataOut.nCode = self.nCode
753 753 dataOut.nBaud = self.nBaud
754 754
755 755 dataOut.data = datadec
756 756
757 757 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
758 758
759 759 dataOut.flagDecodeData = True #asumo q la data esta decodificada
760 760
761 761 if self.__profIndex == self.nCode-1:
762 762 self.__profIndex = 0
763 763 return 1
764 764
765 765 self.__profIndex += 1
766 766
767 767 return 1
768 768 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
769 769
770 770
771 771 class ProfileConcat(Operation):
772 772
773 773 isConfig = False
774 774 buffer = None
775 775
776 776 def __init__(self, **kwargs):
777 777
778 778 Operation.__init__(self, **kwargs)
779 779 self.profileIndex = 0
780 780
781 781 def reset(self):
782 782 self.buffer = numpy.zeros_like(self.buffer)
783 783 self.start_index = 0
784 784 self.times = 1
785 785
786 786 def setup(self, data, m, n=1):
787 787 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
788 788 self.nHeights = data.shape[1]#.nHeights
789 789 self.start_index = 0
790 790 self.times = 1
791 791
792 792 def concat(self, data):
793 793
794 794 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
795 795 self.start_index = self.start_index + self.nHeights
796 796
797 797 def run(self, dataOut, m):
798 798
799 799 dataOut.flagNoData = True
800 800
801 801 if not self.isConfig:
802 802 self.setup(dataOut.data, m, 1)
803 803 self.isConfig = True
804 804
805 805 if dataOut.flagDataAsBlock:
806 806 raise ValueError, "ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False"
807 807
808 808 else:
809 809 self.concat(dataOut.data)
810 810 self.times += 1
811 811 if self.times > m:
812 812 dataOut.data = self.buffer
813 813 self.reset()
814 814 dataOut.flagNoData = False
815 815 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
816 816 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
817 817 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
818 818 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
819 819 dataOut.ippSeconds *= m
820 820
821 821 class ProfileSelector(Operation):
822 822
823 823 profileIndex = None
824 824 # Tamanho total de los perfiles
825 825 nProfiles = None
826 826
827 827 def __init__(self, **kwargs):
828 828
829 829 Operation.__init__(self, **kwargs)
830 830 self.profileIndex = 0
831 831
832 832 def incProfileIndex(self):
833 833
834 834 self.profileIndex += 1
835 835
836 836 if self.profileIndex >= self.nProfiles:
837 837 self.profileIndex = 0
838 838
839 839 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
840 840
841 841 if profileIndex < minIndex:
842 842 return False
843 843
844 844 if profileIndex > maxIndex:
845 845 return False
846 846
847 847 return True
848 848
849 849 def isThisProfileInList(self, profileIndex, profileList):
850 850
851 851 if profileIndex not in profileList:
852 852 return False
853 853
854 854 return True
855 855
856 856 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
857 857
858 858 """
859 859 ProfileSelector:
860 860
861 861 Inputs:
862 862 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
863 863
864 864 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
865 865
866 866 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
867 867
868 868 """
869 869
870 870 if rangeList is not None:
871 871 if type(rangeList[0]) not in (tuple, list):
872 872 rangeList = [rangeList]
873 873
874 874 dataOut.flagNoData = True
875 875
876 876 if dataOut.flagDataAsBlock:
877 877 """
878 878 data dimension = [nChannels, nProfiles, nHeis]
879 879 """
880 880 if profileList != None:
881 881 dataOut.data = dataOut.data[:,profileList,:]
882 882
883 883 if profileRangeList != None:
884 884 minIndex = profileRangeList[0]
885 885 maxIndex = profileRangeList[1]
886 886 profileList = range(minIndex, maxIndex+1)
887 887
888 888 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
889 889
890 890 if rangeList != None:
891 891
892 892 profileList = []
893 893
894 894 for thisRange in rangeList:
895 895 minIndex = thisRange[0]
896 896 maxIndex = thisRange[1]
897 897
898 898 profileList.extend(range(minIndex, maxIndex+1))
899 899
900 900 dataOut.data = dataOut.data[:,profileList,:]
901 901
902 902 dataOut.nProfiles = len(profileList)
903 903 dataOut.profileIndex = dataOut.nProfiles - 1
904 904 dataOut.flagNoData = False
905 905
906 906 return True
907 907
908 908 """
909 909 data dimension = [nChannels, nHeis]
910 910 """
911 911
912 912 if profileList != None:
913 913
914 914 if self.isThisProfileInList(dataOut.profileIndex, profileList):
915 915
916 916 self.nProfiles = len(profileList)
917 917 dataOut.nProfiles = self.nProfiles
918 918 dataOut.profileIndex = self.profileIndex
919 919 dataOut.flagNoData = False
920 920
921 921 self.incProfileIndex()
922 922 return True
923 923
924 924 if profileRangeList != None:
925 925
926 926 minIndex = profileRangeList[0]
927 927 maxIndex = profileRangeList[1]
928 928
929 929 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
930 930
931 931 self.nProfiles = maxIndex - minIndex + 1
932 932 dataOut.nProfiles = self.nProfiles
933 933 dataOut.profileIndex = self.profileIndex
934 934 dataOut.flagNoData = False
935 935
936 936 self.incProfileIndex()
937 937 return True
938 938
939 939 if rangeList != None:
940 940
941 941 nProfiles = 0
942 942
943 943 for thisRange in rangeList:
944 944 minIndex = thisRange[0]
945 945 maxIndex = thisRange[1]
946 946
947 947 nProfiles += maxIndex - minIndex + 1
948 948
949 949 for thisRange in rangeList:
950 950
951 951 minIndex = thisRange[0]
952 952 maxIndex = thisRange[1]
953 953
954 954 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
955 955
956 956 self.nProfiles = nProfiles
957 957 dataOut.nProfiles = self.nProfiles
958 958 dataOut.profileIndex = self.profileIndex
959 959 dataOut.flagNoData = False
960 960
961 961 self.incProfileIndex()
962 962
963 963 break
964 964
965 965 return True
966 966
967 967
968 968 if beam != None: #beam is only for AMISR data
969 969 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
970 970 dataOut.flagNoData = False
971 971 dataOut.profileIndex = self.profileIndex
972 972
973 973 self.incProfileIndex()
974 974
975 975 return True
976 976
977 977 raise ValueError, "ProfileSelector needs profileList, profileRangeList or rangeList parameter"
978 978
979 979 return False
980 980
981 981 class Reshaper(Operation):
982 982
983 983 def __init__(self, **kwargs):
984 984
985 985 Operation.__init__(self, **kwargs)
986 986
987 987 self.__buffer = None
988 988 self.__nitems = 0
989 989
990 990 def __appendProfile(self, dataOut, nTxs):
991 991
992 992 if self.__buffer is None:
993 993 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
994 994 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
995 995
996 996 ini = dataOut.nHeights * self.__nitems
997 997 end = ini + dataOut.nHeights
998 998
999 999 self.__buffer[:, ini:end] = dataOut.data
1000 1000
1001 1001 self.__nitems += 1
1002 1002
1003 1003 return int(self.__nitems*nTxs)
1004 1004
1005 1005 def __getBuffer(self):
1006 1006
1007 1007 if self.__nitems == int(1./self.__nTxs):
1008 1008
1009 1009 self.__nitems = 0
1010 1010
1011 1011 return self.__buffer.copy()
1012 1012
1013 1013 return None
1014 1014
1015 1015 def __checkInputs(self, dataOut, shape, nTxs):
1016 1016
1017 1017 if shape is None and nTxs is None:
1018 1018 raise ValueError, "Reshaper: shape of factor should be defined"
1019 1019
1020 1020 if nTxs:
1021 1021 if nTxs < 0:
1022 1022 raise ValueError, "nTxs should be greater than 0"
1023 1023
1024 1024 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1025 1025 raise ValueError, "nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs))
1026 1026
1027 1027 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1028 1028
1029 1029 return shape, nTxs
1030 1030
1031 1031 if len(shape) != 2 and len(shape) != 3:
1032 1032 raise ValueError, "shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights)
1033 1033
1034 1034 if len(shape) == 2:
1035 1035 shape_tuple = [dataOut.nChannels]
1036 1036 shape_tuple.extend(shape)
1037 1037 else:
1038 1038 shape_tuple = list(shape)
1039 1039
1040 1040 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1041 1041
1042 1042 return shape_tuple, nTxs
1043 1043
1044 1044 def run(self, dataOut, shape=None, nTxs=None):
1045 1045
1046 1046 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1047 1047
1048 1048 dataOut.flagNoData = True
1049 1049 profileIndex = None
1050 1050
1051 1051 if dataOut.flagDataAsBlock:
1052 1052
1053 1053 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1054 1054 dataOut.flagNoData = False
1055 1055
1056 1056 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1057 1057
1058 1058 else:
1059 1059
1060 1060 if self.__nTxs < 1:
1061 1061
1062 1062 self.__appendProfile(dataOut, self.__nTxs)
1063 1063 new_data = self.__getBuffer()
1064 1064
1065 1065 if new_data is not None:
1066 1066 dataOut.data = new_data
1067 1067 dataOut.flagNoData = False
1068 1068
1069 1069 profileIndex = dataOut.profileIndex*nTxs
1070 1070
1071 1071 else:
1072 1072 raise ValueError, "nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)"
1073 1073
1074 1074 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1075 1075
1076 1076 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1077 1077
1078 1078 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1079 1079
1080 1080 dataOut.profileIndex = profileIndex
1081 1081
1082 1082 dataOut.ippSeconds /= self.__nTxs
1083 1083
1084 1084 class SplitProfiles(Operation):
1085 1085
1086 1086 def __init__(self, **kwargs):
1087 1087
1088 1088 Operation.__init__(self, **kwargs)
1089 1089
1090 1090 def run(self, dataOut, n):
1091 1091
1092 1092 dataOut.flagNoData = True
1093 1093 profileIndex = None
1094 1094
1095 1095 if dataOut.flagDataAsBlock:
1096 1096
1097 1097 #nchannels, nprofiles, nsamples
1098 1098 shape = dataOut.data.shape
1099 1099
1100 1100 if shape[2] % n != 0:
1101 1101 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[2])
1102 1102
1103 1103 new_shape = shape[0], shape[1]*n, shape[2]/n
1104 1104
1105 1105 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1106 1106 dataOut.flagNoData = False
1107 1107
1108 1108 profileIndex = int(dataOut.nProfiles/n) - 1
1109 1109
1110 1110 else:
1111 1111
1112 1112 raise ValueError, "Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)"
1113 1113
1114 1114 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1115 1115
1116 1116 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1117 1117
1118 1118 dataOut.nProfiles = int(dataOut.nProfiles*n)
1119 1119
1120 1120 dataOut.profileIndex = profileIndex
1121 1121
1122 1122 dataOut.ippSeconds /= n
1123 1123
1124 1124 class CombineProfiles(Operation):
1125 1125
1126 1126 def __init__(self, **kwargs):
1127 1127
1128 1128 Operation.__init__(self, **kwargs)
1129 1129
1130 1130 self.__remData = None
1131 1131 self.__profileIndex = 0
1132 1132
1133 1133 def run(self, dataOut, n):
1134 1134
1135 1135 dataOut.flagNoData = True
1136 1136 profileIndex = None
1137 1137
1138 1138 if dataOut.flagDataAsBlock:
1139 1139
1140 1140 #nchannels, nprofiles, nsamples
1141 1141 shape = dataOut.data.shape
1142 1142 new_shape = shape[0], shape[1]/n, shape[2]*n
1143 1143
1144 1144 if shape[1] % n != 0:
1145 1145 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[1])
1146 1146
1147 1147 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1148 1148 dataOut.flagNoData = False
1149 1149
1150 1150 profileIndex = int(dataOut.nProfiles*n) - 1
1151 1151
1152 1152 else:
1153 1153
1154 1154 #nchannels, nsamples
1155 1155 if self.__remData is None:
1156 1156 newData = dataOut.data
1157 1157 else:
1158 1158 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1159 1159
1160 1160 self.__profileIndex += 1
1161 1161
1162 1162 if self.__profileIndex < n:
1163 1163 self.__remData = newData
1164 1164 #continue
1165 1165 return
1166 1166
1167 1167 self.__profileIndex = 0
1168 1168 self.__remData = None
1169 1169
1170 1170 dataOut.data = newData
1171 1171 dataOut.flagNoData = False
1172 1172
1173 1173 profileIndex = dataOut.profileIndex/n
1174 1174
1175 1175
1176 1176 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1177 1177
1178 1178 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1179 1179
1180 1180 dataOut.nProfiles = int(dataOut.nProfiles/n)
1181 1181
1182 1182 dataOut.profileIndex = profileIndex
1183 1183
1184 1184 dataOut.ippSeconds *= n
1185 1185
1186 1186
1187 1187 class SSheightProfiles(Operation):
1188 1188
1189 1189 step = None
1190 1190 nsamples = None
1191 1191 bufferShape = None
1192 1192 profileShape= None
1193 1193 sshProfiles = None
1194 1194 profileIndex= None
1195 1195
1196 1196 def __init__(self, **kwargs):
1197 1197
1198 1198 Operation.__init__(self, **kwargs)
1199 1199 self.isConfig = False
1200 1200
1201 1201 def setup(self,dataOut ,step = None , nsamples = None):
1202 1202
1203 1203 if step == None and nsamples == None:
1204 1204 raise ValueError, "step or nheights should be specified ..."
1205 1205
1206 1206 self.step = step
1207 1207 self.nsamples = nsamples
1208 1208 self.__nChannels = dataOut.nChannels
1209 1209 self.__nProfiles = dataOut.nProfiles
1210 1210 self.__nHeis = dataOut.nHeights
1211 1211 shape = dataOut.data.shape #nchannels, nprofiles, nsamples
1212
1212 print "shape",shape
1213 1213 #last test
1214 1214 residue = (shape[1] - self.nsamples) % self.step
1215 1215 if residue != 0:
1216 1216 print "The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue)
1217 1217
1218 1218 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1219 1219 numberProfile = self.nsamples
1220 1220 numberSamples = (shape[1] - self.nsamples)/self.step
1221 1221
1222 print "New number of profile: %d, number of height: %d, Resolution %d"%(numberProfile,numberSamples,deltaHeight*self.step)
1222 print "New number of profile: %d, number of height: %d, Resolution %d Km"%(numberProfile,numberSamples,deltaHeight*self.step)
1223 1223
1224 1224 self.bufferShape = shape[0], numberSamples, numberProfile # nchannels, nsamples , nprofiles
1225 1225 self.profileShape = shape[0], numberProfile, numberSamples # nchannels, nprofiles, nsamples
1226 1226
1227 1227 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1228 1228 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
1229 1229
1230 1230 def run(self, dataOut, step, nsamples):
1231 1231
1232 1232 dataOut.flagNoData = True
1233 1233 dataOut.flagDataAsBlock =False
1234 1234 profileIndex = None
1235 1235
1236 1236 if not self.isConfig:
1237 1237 self.setup(dataOut, step=step , nsamples=nsamples)
1238 1238 self.isConfig = True
1239 1239
1240 1240 for i in range(self.buffer.shape[1]):
1241 1241 self.buffer[:,i] = numpy.flip(dataOut.data[:,i*self.step:i*self.step + self.nsamples])
1242 1242 #self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights])
1243 1243
1244 1244 for j in range(self.buffer.shape[0]):
1245 1245 self.sshProfiles[j] = numpy.transpose(self.buffer[j])
1246 1246
1247 1247 profileIndex = self.nsamples
1248 1248 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1249 1249 ippSeconds = (deltaHeight*1.0e-6)/(0.15)
1250 1250
1251 1251 dataOut.data = self.sshProfiles
1252 1252 dataOut.flagNoData = False
1253 1253 dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0]
1254 1254 dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
1255 1255 dataOut.profileIndex = profileIndex
1256 1256 dataOut.flagDataAsBlock = True
1257 1257 dataOut.ippSeconds = ippSeconds
1258 1258
1259 1259
1260 1260
1261 1261 # import collections
1262 1262 # from scipy.stats import mode
1263 1263 #
1264 1264 # class Synchronize(Operation):
1265 1265 #
1266 1266 # isConfig = False
1267 1267 # __profIndex = 0
1268 1268 #
1269 1269 # def __init__(self, **kwargs):
1270 1270 #
1271 1271 # Operation.__init__(self, **kwargs)
1272 1272 # # self.isConfig = False
1273 1273 # self.__powBuffer = None
1274 1274 # self.__startIndex = 0
1275 1275 # self.__pulseFound = False
1276 1276 #
1277 1277 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1278 1278 #
1279 1279 # #Read data
1280 1280 #
1281 1281 # powerdB = dataOut.getPower(channel = channel)
1282 1282 # noisedB = dataOut.getNoise(channel = channel)[0]
1283 1283 #
1284 1284 # self.__powBuffer.extend(powerdB.flatten())
1285 1285 #
1286 1286 # dataArray = numpy.array(self.__powBuffer)
1287 1287 #
1288 1288 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1289 1289 #
1290 1290 # maxValue = numpy.nanmax(filteredPower)
1291 1291 #
1292 1292 # if maxValue < noisedB + 10:
1293 1293 # #No se encuentra ningun pulso de transmision
1294 1294 # return None
1295 1295 #
1296 1296 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1297 1297 #
1298 1298 # if len(maxValuesIndex) < 2:
1299 1299 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1300 1300 # return None
1301 1301 #
1302 1302 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1303 1303 #
1304 1304 # #Seleccionar solo valores con un espaciamiento de nSamples
1305 1305 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1306 1306 #
1307 1307 # if len(pulseIndex) < 2:
1308 1308 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1309 1309 # return None
1310 1310 #
1311 1311 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1312 1312 #
1313 1313 # #remover senales que se distancien menos de 10 unidades o muestras
1314 1314 # #(No deberian existir IPP menor a 10 unidades)
1315 1315 #
1316 1316 # realIndex = numpy.where(spacing > 10 )[0]
1317 1317 #
1318 1318 # if len(realIndex) < 2:
1319 1319 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1320 1320 # return None
1321 1321 #
1322 1322 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1323 1323 # realPulseIndex = pulseIndex[realIndex]
1324 1324 #
1325 1325 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1326 1326 #
1327 1327 # print "IPP = %d samples" %period
1328 1328 #
1329 1329 # self.__newNSamples = dataOut.nHeights #int(period)
1330 1330 # self.__startIndex = int(realPulseIndex[0])
1331 1331 #
1332 1332 # return 1
1333 1333 #
1334 1334 #
1335 1335 # def setup(self, nSamples, nChannels, buffer_size = 4):
1336 1336 #
1337 1337 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1338 1338 # maxlen = buffer_size*nSamples)
1339 1339 #
1340 1340 # bufferList = []
1341 1341 #
1342 1342 # for i in range(nChannels):
1343 1343 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1344 1344 # maxlen = buffer_size*nSamples)
1345 1345 #
1346 1346 # bufferList.append(bufferByChannel)
1347 1347 #
1348 1348 # self.__nSamples = nSamples
1349 1349 # self.__nChannels = nChannels
1350 1350 # self.__bufferList = bufferList
1351 1351 #
1352 1352 # def run(self, dataOut, channel = 0):
1353 1353 #
1354 1354 # if not self.isConfig:
1355 1355 # nSamples = dataOut.nHeights
1356 1356 # nChannels = dataOut.nChannels
1357 1357 # self.setup(nSamples, nChannels)
1358 1358 # self.isConfig = True
1359 1359 #
1360 1360 # #Append new data to internal buffer
1361 1361 # for thisChannel in range(self.__nChannels):
1362 1362 # bufferByChannel = self.__bufferList[thisChannel]
1363 1363 # bufferByChannel.extend(dataOut.data[thisChannel])
1364 1364 #
1365 1365 # if self.__pulseFound:
1366 1366 # self.__startIndex -= self.__nSamples
1367 1367 #
1368 1368 # #Finding Tx Pulse
1369 1369 # if not self.__pulseFound:
1370 1370 # indexFound = self.__findTxPulse(dataOut, channel)
1371 1371 #
1372 1372 # if indexFound == None:
1373 1373 # dataOut.flagNoData = True
1374 1374 # return
1375 1375 #
1376 1376 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1377 1377 # self.__pulseFound = True
1378 1378 # self.__startIndex = indexFound
1379 1379 #
1380 1380 # #If pulse was found ...
1381 1381 # for thisChannel in range(self.__nChannels):
1382 1382 # bufferByChannel = self.__bufferList[thisChannel]
1383 1383 # #print self.__startIndex
1384 1384 # x = numpy.array(bufferByChannel)
1385 1385 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1386 1386 #
1387 1387 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1388 1388 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1389 1389 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1390 1390 #
1391 1391 # dataOut.data = self.__arrayBuffer
1392 1392 #
1393 1393 # self.__startIndex += self.__newNSamples
1394 1394 #
1395 1395 # return
General Comments 0
You need to be logged in to leave comments. Login now