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