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