##// END OF EJS Templates
schain
RSS
jrodata.py: el atributo ippFactor se inicia en 1...
Daniel Valdez -
r448:318977a743a9
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,723 +1,725
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 os, sys
8 8 import copy
9 9 import numpy
10 10 import datetime
11 11
12 12 from jroheaderIO import SystemHeader, RadarControllerHeader
13 13
14 14 def hildebrand_sekhon(data, navg):
15 15 """
16 16 This method is for the objective determination of de noise level in Doppler spectra. This
17 17 implementation technique is based on the fact that the standard deviation of the spectral
18 18 densities is equal to the mean spectral density for white Gaussian noise
19 19
20 20 Inputs:
21 21 Data : heights
22 22 navg : numbers of averages
23 23
24 24 Return:
25 25 -1 : any error
26 26 anoise : noise's level
27 27 """
28 28
29 29 dataflat = data.copy().reshape(-1)
30 30 dataflat.sort()
31 31 npts = dataflat.size #numbers of points of the data
32 32 npts_noise = 0.2*npts
33 33
34 34 if npts < 32:
35 35 print "error in noise - requires at least 32 points"
36 36 return -1.0
37 37
38 38 dataflat2 = numpy.power(dataflat,2)
39 39
40 40 cs = numpy.cumsum(dataflat)
41 41 cs2 = numpy.cumsum(dataflat2)
42 42
43 43 # data sorted in ascending order
44 44 nmin = int((npts + 7.)/8)
45 45
46 46 for i in range(nmin, npts):
47 47 s = cs[i]
48 48 s2 = cs2[i]
49 49 p = s / float(i);
50 50 p2 = p**2;
51 51 q = s2 / float(i) - p2;
52 52 leftc = p2;
53 53 rightc = q * float(navg);
54 54 R2 = leftc/rightc
55 55
56 56 # Signal detect: R2 < 1 (R2 = leftc/rightc)
57 57 if R2 < 1:
58 58 npts_noise = i
59 59 break
60 60
61 61
62 62 anoise = numpy.average(dataflat[0:npts_noise])
63 63
64 64 return anoise;
65 65
66 66 def sorting_bruce(data, navg):
67 67
68 68 data = data.copy()
69 69
70 70 sortdata = numpy.sort(data)
71 71 lenOfData = len(data)
72 72 nums_min = lenOfData/10
73 73
74 74 if (lenOfData/10) > 0:
75 75 nums_min = lenOfData/10
76 76 else:
77 77 nums_min = 0
78 78
79 79 rtest = 1.0 + 1.0/navg
80 80
81 81 sump = 0.
82 82
83 83 sumq = 0.
84 84
85 85 j = 0
86 86
87 87 cont = 1
88 88
89 89 while((cont==1)and(j<lenOfData)):
90 90
91 91 sump += sortdata[j]
92 92
93 93 sumq += sortdata[j]**2
94 94
95 95 j += 1
96 96
97 97 if j > nums_min:
98 98 if ((sumq*j) <= (rtest*sump**2)):
99 99 lnoise = sump / j
100 100 else:
101 101 j = j - 1
102 102 sump = sump - sortdata[j]
103 103 sumq = sumq - sortdata[j]**2
104 104 cont = 0
105 105
106 106 if j == nums_min:
107 107 lnoise = sump /j
108 108
109 109 return lnoise
110 110
111 111 class JROData:
112 112
113 113 # m_BasicHeader = BasicHeader()
114 114 # m_ProcessingHeader = ProcessingHeader()
115 115
116 116 systemHeaderObj = SystemHeader()
117 117
118 118 radarControllerHeaderObj = RadarControllerHeader()
119 119
120 120 # data = None
121 121
122 122 type = None
123 123
124 124 dtype = None
125 125
126 126 # nChannels = None
127 127
128 128 # nHeights = None
129 129
130 130 nProfiles = None
131 131
132 132 heightList = None
133 133
134 134 channelList = None
135 135
136 136 flagNoData = True
137 137
138 138 flagTimeBlock = False
139 139
140 140 useLocalTime = False
141 141
142 142 utctime = None
143 143
144 144 timeZone = None
145 145
146 146 dstFlag = None
147 147
148 148 errorCount = None
149 149
150 150 blocksize = None
151 151
152 152 nCode = None
153 153
154 154 nBaud = None
155 155
156 156 code = None
157 157
158 158 flagDecodeData = False #asumo q la data no esta decodificada
159 159
160 160 flagDeflipData = False #asumo q la data no esta sin flip
161 161
162 162 flagShiftFFT = False
163 163
164 164 ippSeconds = None
165 165
166 166 timeInterval = None
167 167
168 168 nCohInt = None
169 169
170 170 noise = None
171 171
172 172 windowOfFilter = 1
173 173
174 174 #Speed of ligth
175 175 C = 3e8
176 176
177 177 frequency = 49.92e6
178 178
179 179 realtime = False
180 180
181 181 def __init__(self):
182 182
183 183 raise ValueError, "This class has not been implemented"
184 184
185 185 def copy(self, inputObj=None):
186 186
187 187 if inputObj == None:
188 188 return copy.deepcopy(self)
189 189
190 190 for key in inputObj.__dict__.keys():
191 191 self.__dict__[key] = inputObj.__dict__[key]
192 192
193 193 def deepcopy(self):
194 194
195 195 return copy.deepcopy(self)
196 196
197 197 def isEmpty(self):
198 198
199 199 return self.flagNoData
200 200
201 201 def getNoise(self):
202 202
203 203 raise ValueError, "Not implemented"
204 204
205 205 def getNChannels(self):
206 206
207 207 return len(self.channelList)
208 208
209 209 def getChannelIndexList(self):
210 210
211 211 return range(self.nChannels)
212 212
213 213 def getNHeights(self):
214 214
215 215 return len(self.heightList)
216 216
217 217 def getHeiRange(self, extrapoints=0):
218 218
219 219 heis = self.heightList
220 220 # deltah = self.heightList[1] - self.heightList[0]
221 221 #
222 222 # heis.append(self.heightList[-1])
223 223
224 224 return heis
225 225
226 226 def getltctime(self):
227 227
228 228 if self.useLocalTime:
229 229 return self.utctime - self.timeZone*60
230 230
231 231 return self.utctime
232 232
233 233 def getDatatime(self):
234 234
235 235 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
236 236 return datatime
237 237
238 238 def getTimeRange(self):
239 239
240 240 datatime = []
241 241
242 242 datatime.append(self.ltctime)
243 243 datatime.append(self.ltctime + self.timeInterval)
244 244
245 245 datatime = numpy.array(datatime)
246 246
247 247 return datatime
248 248
249 249 def getFmax(self):
250 250
251 251 PRF = 1./(self.ippSeconds * self.nCohInt)
252 252
253 253 fmax = PRF/2.
254 254
255 255 return fmax
256 256
257 257 def getVmax(self):
258 258
259 259 _lambda = self.C/self.frequency
260 260
261 261 vmax = self.getFmax() * _lambda
262 262
263 263 return vmax
264 264
265 265 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
266 266 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
267 267 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
268 268 noise = property(getNoise, "I'm the 'nHeights' property.")
269 269 datatime = property(getDatatime, "I'm the 'datatime' property")
270 270 ltctime = property(getltctime, "I'm the 'ltctime' property")
271 271
272 272 class Voltage(JROData):
273 273
274 274 #data es un numpy array de 2 dmensiones (canales, alturas)
275 275 data = None
276 276
277 277 def __init__(self):
278 278 '''
279 279 Constructor
280 280 '''
281 281
282 282 self.radarControllerHeaderObj = RadarControllerHeader()
283 283
284 284 self.systemHeaderObj = SystemHeader()
285 285
286 286 self.type = "Voltage"
287 287
288 288 self.data = None
289 289
290 290 self.dtype = None
291 291
292 292 # self.nChannels = 0
293 293
294 294 # self.nHeights = 0
295 295
296 296 self.nProfiles = None
297 297
298 298 self.heightList = None
299 299
300 300 self.channelList = None
301 301
302 302 # self.channelIndexList = None
303 303
304 304 self.flagNoData = True
305 305
306 306 self.flagTimeBlock = False
307 307
308 308 self.utctime = None
309 309
310 310 self.timeZone = None
311 311
312 312 self.dstFlag = None
313 313
314 314 self.errorCount = None
315 315
316 316 self.nCohInt = None
317 317
318 318 self.blocksize = None
319 319
320 320 self.flagDecodeData = False #asumo q la data no esta decodificada
321 321
322 322 self.flagDeflipData = False #asumo q la data no esta sin flip
323 323
324 324 self.flagShiftFFT = False
325 325
326 326
327 327 def getNoisebyHildebrand(self):
328 328 """
329 329 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
330 330
331 331 Return:
332 332 noiselevel
333 333 """
334 334
335 335 for channel in range(self.nChannels):
336 336 daux = self.data_spc[channel,:,:]
337 337 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
338 338
339 339 return self.noise
340 340
341 341 def getNoise(self, type = 1):
342 342
343 343 self.noise = numpy.zeros(self.nChannels)
344 344
345 345 if type == 1:
346 346 noise = self.getNoisebyHildebrand()
347 347
348 348 return 10*numpy.log10(noise)
349 349
350 350 class Spectra(JROData):
351 351
352 352 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
353 353 data_spc = None
354 354
355 355 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
356 356 data_cspc = None
357 357
358 358 #data es un numpy array de 2 dmensiones (canales, alturas)
359 359 data_dc = None
360 360
361 361 nFFTPoints = None
362 362
363 363 nPairs = None
364 364
365 365 pairsList = None
366 366
367 367 nIncohInt = None
368 368
369 369 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
370 370
371 371 nCohInt = None #se requiere para determinar el valor de timeInterval
372 372
373 373 ippFactor = None
374 374
375 375 def __init__(self):
376 376 '''
377 377 Constructor
378 378 '''
379 379
380 380 self.radarControllerHeaderObj = RadarControllerHeader()
381 381
382 382 self.systemHeaderObj = SystemHeader()
383 383
384 384 self.type = "Spectra"
385 385
386 386 # self.data = None
387 387
388 388 self.dtype = None
389 389
390 390 # self.nChannels = 0
391 391
392 392 # self.nHeights = 0
393 393
394 394 self.nProfiles = None
395 395
396 396 self.heightList = None
397 397
398 398 self.channelList = None
399 399
400 400 # self.channelIndexList = None
401 401
402 402 self.flagNoData = True
403 403
404 404 self.flagTimeBlock = False
405 405
406 406 self.utctime = None
407 407
408 408 self.nCohInt = None
409 409
410 410 self.nIncohInt = None
411 411
412 412 self.blocksize = None
413 413
414 414 self.nFFTPoints = None
415 415
416 416 self.wavelength = None
417 417
418 418 self.flagDecodeData = False #asumo q la data no esta decodificada
419 419
420 420 self.flagDeflipData = False #asumo q la data no esta sin flip
421 421
422 422 self.flagShiftFFT = False
423
424 self.ippFactor = 1
423 425
424 426 def getNoisebyHildebrand(self):
425 427 """
426 428 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
427 429
428 430 Return:
429 431 noiselevel
430 432 """
431 433
432 434 for channel in range(self.nChannels):
433 435 daux = self.data_spc[channel,:,:]
434 436 self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
435 437
436 438 return self.noise
437 439
438 440 def getNoisebyWindow(self, heiIndexMin=0, heiIndexMax=-1, freqIndexMin=0, freqIndexMax=-1):
439 441 """
440 442 Determina el ruido del canal utilizando la ventana indicada con las coordenadas:
441 443 (heiIndexMIn, freqIndexMin) hasta (heiIndexMax, freqIndexMAx)
442 444
443 445 Inputs:
444 446 heiIndexMin: Limite inferior del eje de alturas
445 447 heiIndexMax: Limite superior del eje de alturas
446 448 freqIndexMin: Limite inferior del eje de frecuencia
447 449 freqIndexMax: Limite supoerior del eje de frecuencia
448 450 """
449 451
450 452 data = self.data_spc[:, heiIndexMin:heiIndexMax, freqIndexMin:freqIndexMax]
451 453
452 454 for channel in range(self.nChannels):
453 455 daux = data[channel,:,:]
454 456 self.noise[channel] = numpy.average(daux)
455 457
456 458 return self.noise
457 459
458 460 def getNoisebySort(self):
459 461
460 462 for channel in range(self.nChannels):
461 463 daux = self.data_spc[channel,:,:]
462 464 self.noise[channel] = sorting_bruce(daux, self.nIncohInt)
463 465
464 466 return self.noise
465 467
466 468 def getNoise(self, type = 1):
469 if self.noise == None:
470 self.noise = numpy.zeros(self.nChannels)
471
472 if type == 1:
473 self.noise = self.getNoisebyHildebrand()
474
475 if type == 2:
476 self.noise = self.getNoisebySort()
477
478 if type == 3:
479 self.noise = self.getNoisebyWindow()
467 480
468 self.noise = numpy.zeros(self.nChannels)
469
470 if type == 1:
471 noise = self.getNoisebyHildebrand()
472
473 if type == 2:
474 noise = self.getNoisebySort()
475
476 if type == 3:
477 noise = self.getNoisebyWindow()
478
479 return noise
481 return self.noise
480 482
481 483
482 484 def getFreqRange(self, extrapoints=0):
483 485
484 486 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
485 487 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
486 488
487 489 return freqrange
488 490
489 491 def getVelRange(self, extrapoints=0):
490 492
491 493 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
492 494 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
493 495
494 496 return velrange
495 497
496 498 def getNPairs(self):
497 499
498 500 return len(self.pairsList)
499 501
500 502 def getPairsIndexList(self):
501 503
502 504 return range(self.nPairs)
503 505
504 506 def getNormFactor(self):
505 507 pwcode = 1
506 508 if self.flagDecodeData:
507 509 pwcode = numpy.sum(self.code[0]**2)
508 510 normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode
509 511
510 512 return normFactor
511 513
512 514 def getFlagCspc(self):
513 515
514 516 if self.data_cspc == None:
515 517 return True
516 518
517 519 return False
518 520
519 521 def getFlagDc(self):
520 522
521 523 if self.data_dc == None:
522 524 return True
523 525
524 526 return False
525 527
526 528 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
527 529 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
528 530 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
529 531 flag_cspc = property(getFlagCspc)
530 532 flag_dc = property(getFlagDc)
531 533
532 534 class SpectraHeis(JROData):
533 535
534 536 data_spc = None
535 537
536 538 data_cspc = None
537 539
538 540 data_dc = None
539 541
540 542 nFFTPoints = None
541 543
542 544 nPairs = None
543 545
544 546 pairsList = None
545 547
546 548 nIncohInt = None
547 549
548 550 def __init__(self):
549 551
550 552 self.radarControllerHeaderObj = RadarControllerHeader()
551 553
552 554 self.systemHeaderObj = SystemHeader()
553 555
554 556 self.type = "SpectraHeis"
555 557
556 558 self.dtype = None
557 559
558 560 # self.nChannels = 0
559 561
560 562 # self.nHeights = 0
561 563
562 564 self.nProfiles = None
563 565
564 566 self.heightList = None
565 567
566 568 self.channelList = None
567 569
568 570 # self.channelIndexList = None
569 571
570 572 self.flagNoData = True
571 573
572 574 self.flagTimeBlock = False
573 575
574 576 self.nPairs = 0
575 577
576 578 self.utctime = None
577 579
578 580 self.blocksize = None
579 581
580 582 class Fits:
581 583
582 584 heightList = None
583 585
584 586 channelList = None
585 587
586 588 flagNoData = True
587 589
588 590 flagTimeBlock = False
589 591
590 592 useLocalTime = False
591 593
592 594 utctime = None
593 595
594 596 timeZone = None
595 597
596 598 ippSeconds = None
597 599
598 600 timeInterval = None
599 601
600 602 nCohInt = None
601 603
602 604 nIncohInt = None
603 605
604 606 noise = None
605 607
606 608 windowOfFilter = 1
607 609
608 610 #Speed of ligth
609 611 C = 3e8
610 612
611 613 frequency = 49.92e6
612 614
613 615 realtime = False
614 616
615 617
616 618 def __init__(self):
617 619
618 620 self.type = "Fits"
619 621
620 622 self.nProfiles = None
621 623
622 624 self.heightList = None
623 625
624 626 self.channelList = None
625 627
626 628 # self.channelIndexList = None
627 629
628 630 self.flagNoData = True
629 631
630 632 self.utctime = None
631 633
632 634 self.nCohInt = None
633 635
634 636 self.nIncohInt = None
635 637
636 638 self.useLocalTime = True
637 639
638 640 # self.utctime = None
639 641 # self.timeZone = None
640 642 # self.ltctime = None
641 643 # self.timeInterval = None
642 644 # self.header = None
643 645 # self.data_header = None
644 646 # self.data = None
645 647 # self.datatime = None
646 648 # self.flagNoData = False
647 649 # self.expName = ''
648 650 # self.nChannels = None
649 651 # self.nSamples = None
650 652 # self.dataBlocksPerFile = None
651 653 # self.comments = ''
652 654 #
653 655
654 656
655 657 def getltctime(self):
656 658
657 659 if self.useLocalTime:
658 660 return self.utctime - self.timeZone*60
659 661
660 662 return self.utctime
661 663
662 664 def getDatatime(self):
663 665
664 666 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
665 667 return datatime
666 668
667 669 def getTimeRange(self):
668 670
669 671 datatime = []
670 672
671 673 datatime.append(self.ltctime)
672 674 datatime.append(self.ltctime + self.timeInterval)
673 675
674 676 datatime = numpy.array(datatime)
675 677
676 678 return datatime
677 679
678 680 def getHeiRange(self):
679 681
680 682 heis = self.heightList
681 683
682 684 return heis
683 685
684 686 def isEmpty(self):
685 687
686 688 return self.flagNoData
687 689
688 690 def getNHeights(self):
689 691
690 692 return len(self.heightList)
691 693
692 694 def getNChannels(self):
693 695
694 696 return len(self.channelList)
695 697
696 698 def getChannelIndexList(self):
697 699
698 700 return range(self.nChannels)
699 701
700 702 def getNoise(self, type = 1):
701 703
702 704 self.noise = numpy.zeros(self.nChannels)
703 705
704 706 if type == 1:
705 707 noise = self.getNoisebyHildebrand()
706 708
707 709 if type == 2:
708 710 noise = self.getNoisebySort()
709 711
710 712 if type == 3:
711 713 noise = self.getNoisebyWindow()
712 714
713 715 return noise
714 716
715 717 datatime = property(getDatatime, "I'm the 'datatime' property")
716 718 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
717 719 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
718 720 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
719 721 noise = property(getNoise, "I'm the 'nHeights' property.")
720 722 datatime = property(getDatatime, "I'm the 'datatime' property")
721 723 ltctime = property(getltctime, "I'm the 'ltctime' property")
722 724
723 725 ltctime = property(getltctime, "I'm the 'ltctime' property") No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1515 +1,1515
1 1 import numpy
2 2 import time, datetime, os
3 3 from graphics.figure import *
4 4 def isRealtime(utcdatatime):
5 5 utcnow = time.mktime(time.localtime())
6 6 delta = abs(utcnow - utcdatatime) # abs
7 7 if delta >= 30.:
8 8 return False
9 9 return True
10 10
11 11 class CrossSpectraPlot(Figure):
12 12
13 13 __isConfig = None
14 14 __nsubplots = None
15 15
16 16 WIDTH = None
17 17 HEIGHT = None
18 18 WIDTHPROF = None
19 19 HEIGHTPROF = None
20 20 PREFIX = 'cspc'
21 21
22 22 def __init__(self):
23 23
24 24 self.__isConfig = False
25 25 self.__nsubplots = 4
26 26 self.counter_imagwr = 0
27 27 self.WIDTH = 250
28 28 self.HEIGHT = 250
29 29 self.WIDTHPROF = 0
30 30 self.HEIGHTPROF = 0
31 31
32 32 self.PLOT_CODE = 1
33 33 self.FTP_WEI = None
34 34 self.EXP_CODE = None
35 35 self.SUB_EXP_CODE = None
36 36 self.PLOT_POS = None
37 37
38 38 def getSubplots(self):
39 39
40 40 ncol = 4
41 41 nrow = self.nplots
42 42
43 43 return nrow, ncol
44 44
45 45 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
46 46
47 47 self.__showprofile = showprofile
48 48 self.nplots = nplots
49 49
50 50 ncolspan = 1
51 51 colspan = 1
52 52
53 53 self.createFigure(id = id,
54 54 wintitle = wintitle,
55 55 widthplot = self.WIDTH + self.WIDTHPROF,
56 56 heightplot = self.HEIGHT + self.HEIGHTPROF,
57 57 show=True)
58 58
59 59 nrow, ncol = self.getSubplots()
60 60
61 61 counter = 0
62 62 for y in range(nrow):
63 63 for x in range(ncol):
64 64 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
65 65
66 66 counter += 1
67 67
68 68 def run(self, dataOut, id, wintitle="", pairsList=None,
69 69 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
70 70 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
71 71 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
72 72 server=None, folder=None, username=None, password=None,
73 73 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
74 74
75 75 """
76 76
77 77 Input:
78 78 dataOut :
79 79 id :
80 80 wintitle :
81 81 channelList :
82 82 showProfile :
83 83 xmin : None,
84 84 xmax : None,
85 85 ymin : None,
86 86 ymax : None,
87 87 zmin : None,
88 88 zmax : None
89 89 """
90 90
91 91 if pairsList == None:
92 92 pairsIndexList = dataOut.pairsIndexList
93 93 else:
94 94 pairsIndexList = []
95 95 for pair in pairsList:
96 96 if pair not in dataOut.pairsList:
97 97 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
98 98 pairsIndexList.append(dataOut.pairsList.index(pair))
99 99
100 100 if pairsIndexList == []:
101 101 return
102 102
103 103 if len(pairsIndexList) > 4:
104 104 pairsIndexList = pairsIndexList[0:4]
105 factor = dataOut.normFactor
105 #factor = dataOut.normFactor
106 106 x = dataOut.getVelRange(1)
107 107 y = dataOut.getHeiRange()
108 z = dataOut.data_spc[:,:,:]/factor
108 z = dataOut.data_spc[:,:,:]#/factor
109 109 # z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
110 110 avg = numpy.abs(numpy.average(z, axis=1))
111 noise = dataOut.getNoise()/factor
111 noise = dataOut.getNoise()#/factor
112 112
113 113 zdB = 10*numpy.log10(z)
114 114 avgdB = 10*numpy.log10(avg)
115 115 noisedB = 10*numpy.log10(noise)
116 116
117 117
118 118 #thisDatetime = dataOut.datatime
119 119 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
120 120 title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
121 121 xlabel = "Velocity (m/s)"
122 122 ylabel = "Range (Km)"
123 123
124 124 if not self.__isConfig:
125 125
126 126 nplots = len(pairsIndexList)
127 127
128 128 self.setup(id=id,
129 129 nplots=nplots,
130 130 wintitle=wintitle,
131 131 showprofile=False,
132 132 show=show)
133 133
134 134 if xmin == None: xmin = numpy.nanmin(x)
135 135 if xmax == None: xmax = numpy.nanmax(x)
136 136 if ymin == None: ymin = numpy.nanmin(y)
137 137 if ymax == None: ymax = numpy.nanmax(y)
138 138 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
139 139 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
140 140
141 141 self.FTP_WEI = ftp_wei
142 142 self.EXP_CODE = exp_code
143 143 self.SUB_EXP_CODE = sub_exp_code
144 144 self.PLOT_POS = plot_pos
145 145
146 146 self.__isConfig = True
147 147
148 148 self.setWinTitle(title)
149 149
150 150 for i in range(self.nplots):
151 151 pair = dataOut.pairsList[pairsIndexList[i]]
152 152 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
153 153 title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[pair[0]], str_datetime)
154 zdB = 10.*numpy.log10(dataOut.data_spc[pair[0],:,:]/factor)
154 zdB = 10.*numpy.log10(dataOut.data_spc[pair[0],:,:])
155 155 axes0 = self.axesList[i*self.__nsubplots]
156 156 axes0.pcolor(x, y, zdB,
157 157 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
158 158 xlabel=xlabel, ylabel=ylabel, title=title,
159 159 ticksize=9, colormap=power_cmap, cblabel='')
160 160
161 161 title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[pair[1]], str_datetime)
162 zdB = 10.*numpy.log10(dataOut.data_spc[pair[1],:,:]/factor)
162 zdB = 10.*numpy.log10(dataOut.data_spc[pair[1],:,:])
163 163 axes0 = self.axesList[i*self.__nsubplots+1]
164 164 axes0.pcolor(x, y, zdB,
165 165 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
166 166 xlabel=xlabel, ylabel=ylabel, title=title,
167 167 ticksize=9, colormap=power_cmap, cblabel='')
168 168
169 169 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:])
170 170 coherence = numpy.abs(coherenceComplex)
171 171 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
172 172 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
173 173
174 174 title = "Coherence %d%d" %(pair[0], pair[1])
175 175 axes0 = self.axesList[i*self.__nsubplots+2]
176 176 axes0.pcolor(x, y, coherence,
177 177 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=1,
178 178 xlabel=xlabel, ylabel=ylabel, title=title,
179 179 ticksize=9, colormap=coherence_cmap, cblabel='')
180 180
181 181 title = "Phase %d%d" %(pair[0], pair[1])
182 182 axes0 = self.axesList[i*self.__nsubplots+3]
183 183 axes0.pcolor(x, y, phase,
184 184 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
185 185 xlabel=xlabel, ylabel=ylabel, title=title,
186 186 ticksize=9, colormap=phase_cmap, cblabel='')
187 187
188 188
189 189
190 190 self.draw()
191 191
192 192 if save:
193 193
194 194 self.counter_imagwr += 1
195 195 if (self.counter_imagwr==wr_period):
196 196 if figfile == None:
197 197 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
198 198 figfile = self.getFilename(name = str_datetime)
199 199
200 200 self.saveFigure(figpath, figfile)
201 201
202 202 if ftp:
203 203 #provisionalmente envia archivos en el formato de la web en tiempo real
204 204 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
205 205 path = '%s%03d' %(self.PREFIX, self.id)
206 206 ftp_file = os.path.join(path,'ftp','%s.png'%name)
207 207 self.saveFigure(figpath, ftp_file)
208 208 ftp_filename = os.path.join(figpath,ftp_file)
209 209
210 210 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
211 211 self.counter_imagwr = 0
212 212
213 213 self.counter_imagwr = 0
214 214
215 215
216 216 class RTIPlot(Figure):
217 217
218 218 __isConfig = None
219 219 __nsubplots = None
220 220
221 221 WIDTHPROF = None
222 222 HEIGHTPROF = None
223 223 PREFIX = 'rti'
224 224
225 225 def __init__(self):
226 226
227 227 self.timerange = 2*60*60
228 228 self.__isConfig = False
229 229 self.__nsubplots = 1
230 230
231 231 self.WIDTH = 800
232 232 self.HEIGHT = 150
233 233 self.WIDTHPROF = 120
234 234 self.HEIGHTPROF = 0
235 235 self.counter_imagwr = 0
236 236
237 237 self.PLOT_CODE = 0
238 238 self.FTP_WEI = None
239 239 self.EXP_CODE = None
240 240 self.SUB_EXP_CODE = None
241 241 self.PLOT_POS = None
242 242
243 243 def getSubplots(self):
244 244
245 245 ncol = 1
246 246 nrow = self.nplots
247 247
248 248 return nrow, ncol
249 249
250 250 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
251 251
252 252 self.__showprofile = showprofile
253 253 self.nplots = nplots
254 254
255 255 ncolspan = 1
256 256 colspan = 1
257 257 if showprofile:
258 258 ncolspan = 7
259 259 colspan = 6
260 260 self.__nsubplots = 2
261 261
262 262 self.createFigure(id = id,
263 263 wintitle = wintitle,
264 264 widthplot = self.WIDTH + self.WIDTHPROF,
265 265 heightplot = self.HEIGHT + self.HEIGHTPROF,
266 266 show=show)
267 267
268 268 nrow, ncol = self.getSubplots()
269 269
270 270 counter = 0
271 271 for y in range(nrow):
272 272 for x in range(ncol):
273 273
274 274 if counter >= self.nplots:
275 275 break
276 276
277 277 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
278 278
279 279 if showprofile:
280 280 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
281 281
282 282 counter += 1
283 283
284 284 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
285 285 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
286 286 timerange=None,
287 287 save=False, figpath='./', figfile=None, ftp=False, wr_period=1, show=True,
288 288 server=None, folder=None, username=None, password=None,
289 289 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
290 290
291 291 """
292 292
293 293 Input:
294 294 dataOut :
295 295 id :
296 296 wintitle :
297 297 channelList :
298 298 showProfile :
299 299 xmin : None,
300 300 xmax : None,
301 301 ymin : None,
302 302 ymax : None,
303 303 zmin : None,
304 304 zmax : None
305 305 """
306 306
307 307 if channelList == None:
308 308 channelIndexList = dataOut.channelIndexList
309 309 else:
310 310 channelIndexList = []
311 311 for channel in channelList:
312 312 if channel not in dataOut.channelList:
313 313 raise ValueError, "Channel %d is not in dataOut.channelList"
314 314 channelIndexList.append(dataOut.channelList.index(channel))
315 315
316 316 if timerange != None:
317 317 self.timerange = timerange
318 318
319 319 tmin = None
320 320 tmax = None
321 factor = dataOut.normFactor
321 #factor = dataOut.normFactor
322 322 x = dataOut.getTimeRange()
323 323 y = dataOut.getHeiRange()
324 324
325 z = dataOut.data_spc[channelIndexList,:,:]/factor
325 z = dataOut.data_spc[channelIndexList,:,:]#/factor
326 326 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
327 327 avg = numpy.average(z, axis=1)
328 328
329 329 avgdB = 10.*numpy.log10(avg)
330 330
331 331
332 332 # thisDatetime = dataOut.datatime
333 333 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
334 334 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
335 335 xlabel = ""
336 336 ylabel = "Range (Km)"
337 337
338 338 if not self.__isConfig:
339 339
340 340 nplots = len(channelIndexList)
341 341
342 342 self.setup(id=id,
343 343 nplots=nplots,
344 344 wintitle=wintitle,
345 345 showprofile=showprofile,
346 346 show=show)
347 347
348 348 tmin, tmax = self.getTimeLim(x, xmin, xmax)
349 349 if ymin == None: ymin = numpy.nanmin(y)
350 350 if ymax == None: ymax = numpy.nanmax(y)
351 351 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
352 352 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
353 353
354 354 self.FTP_WEI = ftp_wei
355 355 self.EXP_CODE = exp_code
356 356 self.SUB_EXP_CODE = sub_exp_code
357 357 self.PLOT_POS = plot_pos
358 358
359 359 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
360 360 self.__isConfig = True
361 361
362 362
363 363 self.setWinTitle(title)
364 364
365 365 for i in range(self.nplots):
366 366 title = "Channel %d: %s" %(dataOut.channelList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
367 367 axes = self.axesList[i*self.__nsubplots]
368 368 zdB = avgdB[i].reshape((1,-1))
369 369 axes.pcolorbuffer(x, y, zdB,
370 370 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
371 371 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
372 372 ticksize=9, cblabel='', cbsize="1%")
373 373
374 374 if self.__showprofile:
375 375 axes = self.axesList[i*self.__nsubplots +1]
376 376 axes.pline(avgdB[i], y,
377 377 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
378 378 xlabel='dB', ylabel='', title='',
379 379 ytick_visible=False,
380 380 grid='x')
381 381
382 382 self.draw()
383 383
384 384 if save:
385 385
386 386 self.counter_imagwr += 1
387 387 if (self.counter_imagwr==wr_period):
388 388 if figfile == None:
389 389 figfile = self.getFilename(name = self.name)
390 390 self.saveFigure(figpath, figfile)
391 391
392 392 if ftp:
393 393 #provisionalmente envia archivos en el formato de la web en tiempo real
394 394 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
395 395 path = '%s%03d' %(self.PREFIX, self.id)
396 396 ftp_file = os.path.join(path,'ftp','%s.png'%name)
397 397 self.saveFigure(figpath, ftp_file)
398 398 ftp_filename = os.path.join(figpath,ftp_file)
399 399 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
400 400 self.counter_imagwr = 0
401 401
402 402 self.counter_imagwr = 0
403 403
404 404 if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax:
405 405 self.__isConfig = False
406 406
407 407 class SpectraPlot(Figure):
408 408
409 409 __isConfig = None
410 410 __nsubplots = None
411 411
412 412 WIDTHPROF = None
413 413 HEIGHTPROF = None
414 414 PREFIX = 'spc'
415 415
416 416 def __init__(self):
417 417
418 418 self.__isConfig = False
419 419 self.__nsubplots = 1
420 420
421 421 self.WIDTH = 280
422 422 self.HEIGHT = 250
423 423 self.WIDTHPROF = 120
424 424 self.HEIGHTPROF = 0
425 425 self.counter_imagwr = 0
426 426
427 427 self.PLOT_CODE = 1
428 428 self.FTP_WEI = None
429 429 self.EXP_CODE = None
430 430 self.SUB_EXP_CODE = None
431 431 self.PLOT_POS = None
432 432
433 433 def getSubplots(self):
434 434
435 435 ncol = int(numpy.sqrt(self.nplots)+0.9)
436 436 nrow = int(self.nplots*1./ncol + 0.9)
437 437
438 438 return nrow, ncol
439 439
440 440 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
441 441
442 442 self.__showprofile = showprofile
443 443 self.nplots = nplots
444 444
445 445 ncolspan = 1
446 446 colspan = 1
447 447 if showprofile:
448 448 ncolspan = 3
449 449 colspan = 2
450 450 self.__nsubplots = 2
451 451
452 452 self.createFigure(id = id,
453 453 wintitle = wintitle,
454 454 widthplot = self.WIDTH + self.WIDTHPROF,
455 455 heightplot = self.HEIGHT + self.HEIGHTPROF,
456 456 show=show)
457 457
458 458 nrow, ncol = self.getSubplots()
459 459
460 460 counter = 0
461 461 for y in range(nrow):
462 462 for x in range(ncol):
463 463
464 464 if counter >= self.nplots:
465 465 break
466 466
467 467 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
468 468
469 469 if showprofile:
470 470 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
471 471
472 472 counter += 1
473 473
474 474 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
475 475 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
476 476 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
477 477 server=None, folder=None, username=None, password=None,
478 478 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
479 479
480 480 """
481 481
482 482 Input:
483 483 dataOut :
484 484 id :
485 485 wintitle :
486 486 channelList :
487 487 showProfile :
488 488 xmin : None,
489 489 xmax : None,
490 490 ymin : None,
491 491 ymax : None,
492 492 zmin : None,
493 493 zmax : None
494 494 """
495 495
496 496 if realtime:
497 497 if not(isRealtime(utcdatatime = dataOut.utctime)):
498 498 print 'Skipping this plot function'
499 499 return
500 500
501 501 if channelList == None:
502 502 channelIndexList = dataOut.channelIndexList
503 503 else:
504 504 channelIndexList = []
505 505 for channel in channelList:
506 506 if channel not in dataOut.channelList:
507 507 raise ValueError, "Channel %d is not in dataOut.channelList"
508 508 channelIndexList.append(dataOut.channelList.index(channel))
509 factor = dataOut.normFactor
509 #factor = dataOut.normFactor
510 510 x = dataOut.getVelRange(1)
511 511 y = dataOut.getHeiRange()
512 512
513 z = dataOut.data_spc[channelIndexList,:,:]/factor
513 z = dataOut.data_spc[channelIndexList,:,:]#/factor
514 514 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
515 515 avg = numpy.average(z, axis=1)
516 noise = dataOut.getNoise()/factor
516 noise = dataOut.getNoise()#/factor
517 517
518 518 zdB = 10*numpy.log10(z)
519 519 avgdB = 10*numpy.log10(avg)
520 520 noisedB = 10*numpy.log10(noise)
521 521
522 522 #thisDatetime = dataOut.datatime
523 523 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
524 524 title = wintitle + " Spectra"
525 525 xlabel = "Velocity (m/s)"
526 526 ylabel = "Range (Km)"
527 527
528 528 if not self.__isConfig:
529 529
530 530 nplots = len(channelIndexList)
531 531
532 532 self.setup(id=id,
533 533 nplots=nplots,
534 534 wintitle=wintitle,
535 535 showprofile=showprofile,
536 536 show=show)
537 537
538 538 if xmin == None: xmin = numpy.nanmin(x)
539 539 if xmax == None: xmax = numpy.nanmax(x)
540 540 if ymin == None: ymin = numpy.nanmin(y)
541 541 if ymax == None: ymax = numpy.nanmax(y)
542 542 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
543 543 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
544 544
545 545 self.FTP_WEI = ftp_wei
546 546 self.EXP_CODE = exp_code
547 547 self.SUB_EXP_CODE = sub_exp_code
548 548 self.PLOT_POS = plot_pos
549 549
550 550 self.__isConfig = True
551 551
552 552 self.setWinTitle(title)
553 553
554 554 for i in range(self.nplots):
555 555 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
556 556 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
557 557 axes = self.axesList[i*self.__nsubplots]
558 558 axes.pcolor(x, y, zdB[i,:,:],
559 559 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
560 560 xlabel=xlabel, ylabel=ylabel, title=title,
561 561 ticksize=9, cblabel='')
562 562
563 563 if self.__showprofile:
564 564 axes = self.axesList[i*self.__nsubplots +1]
565 565 axes.pline(avgdB[i], y,
566 566 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
567 567 xlabel='dB', ylabel='', title='',
568 568 ytick_visible=False,
569 569 grid='x')
570 570
571 571 noiseline = numpy.repeat(noisedB[i], len(y))
572 572 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
573 573
574 574 self.draw()
575 575
576 576 if save:
577 577
578 578 self.counter_imagwr += 1
579 579 if (self.counter_imagwr==wr_period):
580 580 if figfile == None:
581 581 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
582 582 figfile = self.getFilename(name = str_datetime)
583 583
584 584 self.saveFigure(figpath, figfile)
585 585
586 586 if ftp:
587 587 #provisionalmente envia archivos en el formato de la web en tiempo real
588 588 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
589 589 path = '%s%03d' %(self.PREFIX, self.id)
590 590 ftp_file = os.path.join(path,'ftp','%s.png'%name)
591 591 self.saveFigure(figpath, ftp_file)
592 592 ftp_filename = os.path.join(figpath,ftp_file)
593 593 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
594 594 self.counter_imagwr = 0
595 595
596 596
597 597 self.counter_imagwr = 0
598 598
599 599
600 600 class Scope(Figure):
601 601
602 602 __isConfig = None
603 603
604 604 def __init__(self):
605 605
606 606 self.__isConfig = False
607 607 self.WIDTH = 600
608 608 self.HEIGHT = 200
609 609 self.counter_imagwr = 0
610 610
611 611 def getSubplots(self):
612 612
613 613 nrow = self.nplots
614 614 ncol = 3
615 615 return nrow, ncol
616 616
617 617 def setup(self, id, nplots, wintitle, show):
618 618
619 619 self.nplots = nplots
620 620
621 621 self.createFigure(id=id,
622 622 wintitle=wintitle,
623 623 show=show)
624 624
625 625 nrow,ncol = self.getSubplots()
626 626 colspan = 3
627 627 rowspan = 1
628 628
629 629 for i in range(nplots):
630 630 self.addAxes(nrow, ncol, i, 0, colspan, rowspan)
631 631
632 632
633 633
634 634 def run(self, dataOut, id, wintitle="", channelList=None,
635 635 xmin=None, xmax=None, ymin=None, ymax=None, save=False,
636 636 figpath='./', figfile=None, show=True, wr_period=1,
637 637 server=None, folder=None, username=None, password=None):
638 638
639 639 """
640 640
641 641 Input:
642 642 dataOut :
643 643 id :
644 644 wintitle :
645 645 channelList :
646 646 xmin : None,
647 647 xmax : None,
648 648 ymin : None,
649 649 ymax : None,
650 650 """
651 651
652 652 if channelList == None:
653 653 channelIndexList = dataOut.channelIndexList
654 654 else:
655 655 channelIndexList = []
656 656 for channel in channelList:
657 657 if channel not in dataOut.channelList:
658 658 raise ValueError, "Channel %d is not in dataOut.channelList"
659 659 channelIndexList.append(dataOut.channelList.index(channel))
660 660
661 661 x = dataOut.heightList
662 662 y = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
663 663 y = y.real
664 664
665 665 #thisDatetime = dataOut.datatime
666 666 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
667 667 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
668 668 xlabel = "Range (Km)"
669 669 ylabel = "Intensity"
670 670
671 671 if not self.__isConfig:
672 672 nplots = len(channelIndexList)
673 673
674 674 self.setup(id=id,
675 675 nplots=nplots,
676 676 wintitle=wintitle,
677 677 show=show)
678 678
679 679 if xmin == None: xmin = numpy.nanmin(x)
680 680 if xmax == None: xmax = numpy.nanmax(x)
681 681 if ymin == None: ymin = numpy.nanmin(y)
682 682 if ymax == None: ymax = numpy.nanmax(y)
683 683
684 684 self.__isConfig = True
685 685
686 686 self.setWinTitle(title)
687 687
688 688 for i in range(len(self.axesList)):
689 689 title = "Channel %d" %(i)
690 690 axes = self.axesList[i]
691 691 ychannel = y[i,:]
692 692 axes.pline(x, ychannel,
693 693 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
694 694 xlabel=xlabel, ylabel=ylabel, title=title)
695 695
696 696 self.draw()
697 697
698 698 if save:
699 699 date = thisDatetime.strftime("%Y%m%d_%H%M%S")
700 700 if figfile == None:
701 701 figfile = self.getFilename(name = date)
702 702
703 703 self.saveFigure(figpath, figfile)
704 704
705 705 self.counter_imagwr += 1
706 706 if (ftp and (self.counter_imagwr==wr_period)):
707 707 ftp_filename = os.path.join(figpath,figfile)
708 708 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
709 709 self.counter_imagwr = 0
710 710
711 711 class PowerProfilePlot(Figure):
712 712 __isConfig = None
713 713 __nsubplots = None
714 714
715 715 WIDTHPROF = None
716 716 HEIGHTPROF = None
717 717 PREFIX = 'spcprofile'
718 718
719 719 def __init__(self):
720 720 self.__isConfig = False
721 721 self.__nsubplots = 1
722 722
723 723 self.WIDTH = 300
724 724 self.HEIGHT = 500
725 725 self.counter_imagwr = 0
726 726
727 727 def getSubplots(self):
728 728 ncol = 1
729 729 nrow = 1
730 730
731 731 return nrow, ncol
732 732
733 733 def setup(self, id, nplots, wintitle, show):
734 734
735 735 self.nplots = nplots
736 736
737 737 ncolspan = 1
738 738 colspan = 1
739 739
740 740 self.createFigure(id = id,
741 741 wintitle = wintitle,
742 742 widthplot = self.WIDTH,
743 743 heightplot = self.HEIGHT,
744 744 show=show)
745 745
746 746 nrow, ncol = self.getSubplots()
747 747
748 748 counter = 0
749 749 for y in range(nrow):
750 750 for x in range(ncol):
751 751 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
752 752
753 753 def run(self, dataOut, id, wintitle="", channelList=None,
754 754 xmin=None, xmax=None, ymin=None, ymax=None,
755 755 save=False, figpath='./', figfile=None, show=True, wr_period=1,
756 756 server=None, folder=None, username=None, password=None,):
757 757
758 758 if channelList == None:
759 759 channelIndexList = dataOut.channelIndexList
760 760 channelList = dataOut.channelList
761 761 else:
762 762 channelIndexList = []
763 763 for channel in channelList:
764 764 if channel not in dataOut.channelList:
765 765 raise ValueError, "Channel %d is not in dataOut.channelList"
766 766 channelIndexList.append(dataOut.channelList.index(channel))
767 767
768 factor = dataOut.normFactor
768 #factor = dataOut.normFactor
769 769 y = dataOut.getHeiRange()
770 x = dataOut.data_spc[channelIndexList,:,:]/factor
770 x = dataOut.data_spc[channelIndexList,:,:]#/factor
771 771 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
772 772 avg = numpy.average(x, axis=1)
773 773
774 774 avgdB = 10*numpy.log10(avg)
775 775
776 776 #thisDatetime = dataOut.datatime
777 777 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
778 778 title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
779 779 xlabel = "dB"
780 780 ylabel = "Range (Km)"
781 781
782 782 if not self.__isConfig:
783 783
784 784 nplots = 1
785 785
786 786 self.setup(id=id,
787 787 nplots=nplots,
788 788 wintitle=wintitle,
789 789 show=show)
790 790
791 791 if ymin == None: ymin = numpy.nanmin(y)
792 792 if ymax == None: ymax = numpy.nanmax(y)
793 793 if xmin == None: xmin = numpy.nanmin(avgdB)*0.9
794 794 if xmax == None: xmax = numpy.nanmax(avgdB)*0.9
795 795
796 796 self.__isConfig = True
797 797
798 798 self.setWinTitle(title)
799 799
800 800
801 801 title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
802 802 axes = self.axesList[0]
803 803
804 804 legendlabels = ["channel %d"%x for x in channelList]
805 805 axes.pmultiline(avgdB, y,
806 806 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
807 807 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
808 808 ytick_visible=True, nxticks=5,
809 809 grid='x')
810 810
811 811 self.draw()
812 812
813 813 if save:
814 814 date = thisDatetime.strftime("%Y%m%d")
815 815 if figfile == None:
816 816 figfile = self.getFilename(name = date)
817 817
818 818 self.saveFigure(figpath, figfile)
819 819
820 820 self.counter_imagwr += 1
821 821 if (ftp and (self.counter_imagwr==wr_period)):
822 822 ftp_filename = os.path.join(figpath,figfile)
823 823 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
824 824 self.counter_imagwr = 0
825 825
826 826 class CoherenceMap(Figure):
827 827 __isConfig = None
828 828 __nsubplots = None
829 829
830 830 WIDTHPROF = None
831 831 HEIGHTPROF = None
832 832 PREFIX = 'cmap'
833 833
834 834 def __init__(self):
835 835 self.timerange = 2*60*60
836 836 self.__isConfig = False
837 837 self.__nsubplots = 1
838 838
839 839 self.WIDTH = 800
840 840 self.HEIGHT = 150
841 841 self.WIDTHPROF = 120
842 842 self.HEIGHTPROF = 0
843 843 self.counter_imagwr = 0
844 844
845 845 self.PLOT_CODE = 3
846 846 self.FTP_WEI = None
847 847 self.EXP_CODE = None
848 848 self.SUB_EXP_CODE = None
849 849 self.PLOT_POS = None
850 850 self.counter_imagwr = 0
851 851
852 852 def getSubplots(self):
853 853 ncol = 1
854 854 nrow = self.nplots*2
855 855
856 856 return nrow, ncol
857 857
858 858 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
859 859 self.__showprofile = showprofile
860 860 self.nplots = nplots
861 861
862 862 ncolspan = 1
863 863 colspan = 1
864 864 if showprofile:
865 865 ncolspan = 7
866 866 colspan = 6
867 867 self.__nsubplots = 2
868 868
869 869 self.createFigure(id = id,
870 870 wintitle = wintitle,
871 871 widthplot = self.WIDTH + self.WIDTHPROF,
872 872 heightplot = self.HEIGHT + self.HEIGHTPROF,
873 873 show=True)
874 874
875 875 nrow, ncol = self.getSubplots()
876 876
877 877 for y in range(nrow):
878 878 for x in range(ncol):
879 879
880 880 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
881 881
882 882 if showprofile:
883 883 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
884 884
885 885 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
886 886 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
887 887 timerange=None,
888 888 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
889 889 coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
890 890 server=None, folder=None, username=None, password=None,
891 891 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
892 892
893 893 if pairsList == None:
894 894 pairsIndexList = dataOut.pairsIndexList
895 895 else:
896 896 pairsIndexList = []
897 897 for pair in pairsList:
898 898 if pair not in dataOut.pairsList:
899 899 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
900 900 pairsIndexList.append(dataOut.pairsList.index(pair))
901 901
902 902 if timerange != None:
903 903 self.timerange = timerange
904 904
905 905 if pairsIndexList == []:
906 906 return
907 907
908 908 if len(pairsIndexList) > 4:
909 909 pairsIndexList = pairsIndexList[0:4]
910 910
911 911 tmin = None
912 912 tmax = None
913 913 x = dataOut.getTimeRange()
914 914 y = dataOut.getHeiRange()
915 915
916 916 #thisDatetime = dataOut.datatime
917 917 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
918 918 title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
919 919 xlabel = ""
920 920 ylabel = "Range (Km)"
921 921
922 922 if not self.__isConfig:
923 923 nplots = len(pairsIndexList)
924 924 self.setup(id=id,
925 925 nplots=nplots,
926 926 wintitle=wintitle,
927 927 showprofile=showprofile,
928 928 show=show)
929 929
930 930 tmin, tmax = self.getTimeLim(x, xmin, xmax)
931 931 if ymin == None: ymin = numpy.nanmin(y)
932 932 if ymax == None: ymax = numpy.nanmax(y)
933 933 if zmin == None: zmin = 0.
934 934 if zmax == None: zmax = 1.
935 935
936 936 self.FTP_WEI = ftp_wei
937 937 self.EXP_CODE = exp_code
938 938 self.SUB_EXP_CODE = sub_exp_code
939 939 self.PLOT_POS = plot_pos
940 940
941 941 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
942 942
943 943 self.__isConfig = True
944 944
945 945 self.setWinTitle(title)
946 946
947 947 for i in range(self.nplots):
948 948
949 949 pair = dataOut.pairsList[pairsIndexList[i]]
950 950 # coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:])
951 951 # avgcoherenceComplex = numpy.average(coherenceComplex, axis=0)
952 952 # coherence = numpy.abs(avgcoherenceComplex)
953 953
954 954 ## coherence = numpy.abs(coherenceComplex)
955 955 ## avg = numpy.average(coherence, axis=0)
956 956
957 957 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
958 958 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
959 959 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
960 960
961 961
962 962 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
963 963 coherence = numpy.abs(avgcoherenceComplex)
964 964
965 965 z = coherence.reshape((1,-1))
966 966
967 967 counter = 0
968 968
969 969 title = "Coherence %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
970 970 axes = self.axesList[i*self.__nsubplots*2]
971 971 axes.pcolorbuffer(x, y, z,
972 972 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
973 973 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
974 974 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
975 975
976 976 if self.__showprofile:
977 977 counter += 1
978 978 axes = self.axesList[i*self.__nsubplots*2 + counter]
979 979 axes.pline(coherence, y,
980 980 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
981 981 xlabel='', ylabel='', title='', ticksize=7,
982 982 ytick_visible=False, nxticks=5,
983 983 grid='x')
984 984
985 985 counter += 1
986 986 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
987 987 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
988 988 # avg = numpy.average(phase, axis=0)
989 989 z = phase.reshape((1,-1))
990 990
991 991 title = "Phase %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
992 992 axes = self.axesList[i*self.__nsubplots*2 + counter]
993 993 axes.pcolorbuffer(x, y, z,
994 994 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
995 995 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
996 996 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
997 997
998 998 if self.__showprofile:
999 999 counter += 1
1000 1000 axes = self.axesList[i*self.__nsubplots*2 + counter]
1001 1001 axes.pline(phase, y,
1002 1002 xmin=-180, xmax=180, ymin=ymin, ymax=ymax,
1003 1003 xlabel='', ylabel='', title='', ticksize=7,
1004 1004 ytick_visible=False, nxticks=4,
1005 1005 grid='x')
1006 1006
1007 1007 self.draw()
1008 1008
1009 1009 if save:
1010 1010
1011 1011 self.counter_imagwr += 1
1012 1012 if (self.counter_imagwr==wr_period):
1013 1013 if figfile == None:
1014 1014 figfile = self.getFilename(name = self.name)
1015 1015 self.saveFigure(figpath, figfile)
1016 1016
1017 1017 if ftp:
1018 1018 #provisionalmente envia archivos en el formato de la web en tiempo real
1019 1019 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
1020 1020 path = '%s%03d' %(self.PREFIX, self.id)
1021 1021 ftp_file = os.path.join(path,'ftp','%s.png'%name)
1022 1022 self.saveFigure(figpath, ftp_file)
1023 1023 ftp_filename = os.path.join(figpath,ftp_file)
1024 1024 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
1025 1025 self.counter_imagwr = 0
1026 1026
1027 1027 self.counter_imagwr = 0
1028 1028
1029 1029
1030 1030 if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax:
1031 1031 self.__isConfig = False
1032 1032
1033 1033 class Noise(Figure):
1034 1034
1035 1035 __isConfig = None
1036 1036 __nsubplots = None
1037 1037
1038 1038 PREFIX = 'noise'
1039 1039
1040 1040 def __init__(self):
1041 1041
1042 1042 self.timerange = 24*60*60
1043 1043 self.__isConfig = False
1044 1044 self.__nsubplots = 1
1045 1045 self.counter_imagwr = 0
1046 1046 self.WIDTH = 600
1047 1047 self.HEIGHT = 300
1048 1048 self.WIDTHPROF = 120
1049 1049 self.HEIGHTPROF = 0
1050 1050 self.xdata = None
1051 1051 self.ydata = None
1052 1052
1053 1053 self.PLOT_CODE = 77
1054 1054 self.FTP_WEI = None
1055 1055 self.EXP_CODE = None
1056 1056 self.SUB_EXP_CODE = None
1057 1057 self.PLOT_POS = None
1058 1058
1059 1059 def getSubplots(self):
1060 1060
1061 1061 ncol = 1
1062 1062 nrow = 1
1063 1063
1064 1064 return nrow, ncol
1065 1065
1066 1066 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1067 1067
1068 1068 self.__showprofile = showprofile
1069 1069 self.nplots = nplots
1070 1070
1071 1071 ncolspan = 7
1072 1072 colspan = 6
1073 1073 self.__nsubplots = 2
1074 1074
1075 1075 self.createFigure(id = id,
1076 1076 wintitle = wintitle,
1077 1077 widthplot = self.WIDTH+self.WIDTHPROF,
1078 1078 heightplot = self.HEIGHT+self.HEIGHTPROF,
1079 1079 show=show)
1080 1080
1081 1081 nrow, ncol = self.getSubplots()
1082 1082
1083 1083 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1084 1084
1085 1085
1086 1086 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1087 1087 xmin=None, xmax=None, ymin=None, ymax=None,
1088 1088 timerange=None,
1089 1089 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1090 1090 server=None, folder=None, username=None, password=None,
1091 1091 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1092 1092
1093 1093 if channelList == None:
1094 1094 channelIndexList = dataOut.channelIndexList
1095 1095 channelList = dataOut.channelList
1096 1096 else:
1097 1097 channelIndexList = []
1098 1098 for channel in channelList:
1099 1099 if channel not in dataOut.channelList:
1100 1100 raise ValueError, "Channel %d is not in dataOut.channelList"
1101 1101 channelIndexList.append(dataOut.channelList.index(channel))
1102 1102
1103 1103 if timerange != None:
1104 1104 self.timerange = timerange
1105 1105
1106 1106 tmin = None
1107 1107 tmax = None
1108 1108 x = dataOut.getTimeRange()
1109 1109 y = dataOut.getHeiRange()
1110 factor = dataOut.normFactor
1111 noise = dataOut.getNoise()/factor
1110 #factor = dataOut.normFactor
1111 noise = dataOut.getNoise()#/factor
1112 1112 noisedB = 10*numpy.log10(noise)
1113 1113
1114 1114 #thisDatetime = dataOut.datatime
1115 1115 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1116 1116 title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1117 1117 xlabel = ""
1118 1118 ylabel = "Intensity (dB)"
1119 1119
1120 1120 if not self.__isConfig:
1121 1121
1122 1122 nplots = 1
1123 1123
1124 1124 self.setup(id=id,
1125 1125 nplots=nplots,
1126 1126 wintitle=wintitle,
1127 1127 showprofile=showprofile,
1128 1128 show=show)
1129 1129
1130 1130 tmin, tmax = self.getTimeLim(x, xmin, xmax)
1131 1131 if ymin == None: ymin = numpy.nanmin(noisedB) - 10.0
1132 1132 if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
1133 1133
1134 1134 self.FTP_WEI = ftp_wei
1135 1135 self.EXP_CODE = exp_code
1136 1136 self.SUB_EXP_CODE = sub_exp_code
1137 1137 self.PLOT_POS = plot_pos
1138 1138
1139 1139 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1140 1140
1141 1141
1142 1142 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1143 1143 self.__isConfig = True
1144 1144
1145 1145 self.xdata = numpy.array([])
1146 1146 self.ydata = numpy.array([])
1147 1147
1148 1148 self.setWinTitle(title)
1149 1149
1150 1150
1151 1151 title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1152 1152
1153 1153 legendlabels = ["channel %d"%(idchannel+1) for idchannel in channelList]
1154 1154 axes = self.axesList[0]
1155 1155
1156 1156 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1157 1157
1158 1158 if len(self.ydata)==0:
1159 1159 self.ydata = noisedB[channelIndexList].reshape(-1,1)
1160 1160 else:
1161 1161 self.ydata = numpy.hstack((self.ydata, noisedB[channelIndexList].reshape(-1,1)))
1162 1162
1163 1163
1164 1164 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1165 1165 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax,
1166 1166 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1167 1167 XAxisAsTime=True, grid='both'
1168 1168 )
1169 1169
1170 1170 self.draw()
1171 1171
1172 1172 # if save:
1173 1173 #
1174 1174 # if figfile == None:
1175 1175 # figfile = self.getFilename(name = self.name)
1176 1176 #
1177 1177 # self.saveFigure(figpath, figfile)
1178 1178
1179 1179 if save:
1180 1180
1181 1181 self.counter_imagwr += 1
1182 1182 if (self.counter_imagwr==wr_period):
1183 1183 if figfile == None:
1184 1184 figfile = self.getFilename(name = self.name)
1185 1185 self.saveFigure(figpath, figfile)
1186 1186
1187 1187 if ftp:
1188 1188 #provisionalmente envia archivos en el formato de la web en tiempo real
1189 1189 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
1190 1190 path = '%s%03d' %(self.PREFIX, self.id)
1191 1191 ftp_file = os.path.join(path,'ftp','%s.png'%name)
1192 1192 self.saveFigure(figpath, ftp_file)
1193 1193 ftp_filename = os.path.join(figpath,ftp_file)
1194 1194 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
1195 1195 self.counter_imagwr = 0
1196 1196
1197 1197 self.counter_imagwr = 0
1198 1198
1199 1199 if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax:
1200 1200 self.__isConfig = False
1201 1201 del self.xdata
1202 1202 del self.ydata
1203 1203
1204 1204
1205 1205 class SpectraHeisScope(Figure):
1206 1206
1207 1207
1208 1208 __isConfig = None
1209 1209 __nsubplots = None
1210 1210
1211 1211 WIDTHPROF = None
1212 1212 HEIGHTPROF = None
1213 1213 PREFIX = 'spc'
1214 1214
1215 1215 def __init__(self):
1216 1216
1217 1217 self.__isConfig = False
1218 1218 self.__nsubplots = 1
1219 1219
1220 1220 self.WIDTH = 230
1221 1221 self.HEIGHT = 250
1222 1222 self.WIDTHPROF = 120
1223 1223 self.HEIGHTPROF = 0
1224 1224 self.counter_imagwr = 0
1225 1225
1226 1226 def getSubplots(self):
1227 1227
1228 1228 ncol = int(numpy.sqrt(self.nplots)+0.9)
1229 1229 nrow = int(self.nplots*1./ncol + 0.9)
1230 1230
1231 1231 return nrow, ncol
1232 1232
1233 1233 def setup(self, id, nplots, wintitle, show):
1234 1234
1235 1235 showprofile = False
1236 1236 self.__showprofile = showprofile
1237 1237 self.nplots = nplots
1238 1238
1239 1239 ncolspan = 1
1240 1240 colspan = 1
1241 1241 if showprofile:
1242 1242 ncolspan = 3
1243 1243 colspan = 2
1244 1244 self.__nsubplots = 2
1245 1245
1246 1246 self.createFigure(id = id,
1247 1247 wintitle = wintitle,
1248 1248 widthplot = self.WIDTH + self.WIDTHPROF,
1249 1249 heightplot = self.HEIGHT + self.HEIGHTPROF,
1250 1250 show = show)
1251 1251
1252 1252 nrow, ncol = self.getSubplots()
1253 1253
1254 1254 counter = 0
1255 1255 for y in range(nrow):
1256 1256 for x in range(ncol):
1257 1257
1258 1258 if counter >= self.nplots:
1259 1259 break
1260 1260
1261 1261 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1262 1262
1263 1263 if showprofile:
1264 1264 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
1265 1265
1266 1266 counter += 1
1267 1267
1268 1268
1269 1269 def run(self, dataOut, id, wintitle="", channelList=None,
1270 1270 xmin=None, xmax=None, ymin=None, ymax=None, save=False,
1271 1271 figpath='./', figfile=None, ftp=False, wr_period=1, show=True,
1272 1272 server=None, folder=None, username=None, password=None):
1273 1273
1274 1274 """
1275 1275
1276 1276 Input:
1277 1277 dataOut :
1278 1278 id :
1279 1279 wintitle :
1280 1280 channelList :
1281 1281 xmin : None,
1282 1282 xmax : None,
1283 1283 ymin : None,
1284 1284 ymax : None,
1285 1285 """
1286 1286
1287 1287 if dataOut.realtime:
1288 1288 if not(isRealtime(utcdatatime = dataOut.utctime)):
1289 1289 print 'Skipping this plot function'
1290 1290 return
1291 1291
1292 1292 if channelList == None:
1293 1293 channelIndexList = dataOut.channelIndexList
1294 1294 else:
1295 1295 channelIndexList = []
1296 1296 for channel in channelList:
1297 1297 if channel not in dataOut.channelList:
1298 1298 raise ValueError, "Channel %d is not in dataOut.channelList"
1299 1299 channelIndexList.append(dataOut.channelList.index(channel))
1300 1300
1301 1301 # x = dataOut.heightList
1302 1302 c = 3E8
1303 1303 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1304 1304 #deberia cambiar para el caso de 1Mhz y 100KHz
1305 1305 x = numpy.arange(-1*dataOut.nHeights/2.,dataOut.nHeights/2.)*(c/(2*deltaHeight*dataOut.nHeights*1000))
1306 1306 #para 1Mhz descomentar la siguiente linea
1307 1307 #x= x/(10000.0)
1308 1308 # y = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
1309 1309 # y = y.real
1310 1310 datadB = 10.*numpy.log10(dataOut.data_spc)
1311 1311 y = datadB
1312 1312
1313 1313 #thisDatetime = dataOut.datatime
1314 1314 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1315 1315 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1316 1316 xlabel = ""
1317 1317 #para 1Mhz descomentar la siguiente linea
1318 1318 #xlabel = "Frequency x 10000"
1319 1319 ylabel = "Intensity (dB)"
1320 1320
1321 1321 if not self.__isConfig:
1322 1322 nplots = len(channelIndexList)
1323 1323
1324 1324 self.setup(id=id,
1325 1325 nplots=nplots,
1326 1326 wintitle=wintitle,
1327 1327 show=show)
1328 1328
1329 1329 if xmin == None: xmin = numpy.nanmin(x)
1330 1330 if xmax == None: xmax = numpy.nanmax(x)
1331 1331 if ymin == None: ymin = numpy.nanmin(y)
1332 1332 if ymax == None: ymax = numpy.nanmax(y)
1333 1333
1334 1334 self.__isConfig = True
1335 1335
1336 1336 self.setWinTitle(title)
1337 1337
1338 1338 for i in range(len(self.axesList)):
1339 1339 ychannel = y[i,:]
1340 1340 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
1341 1341 title = "Channel %d: %4.2fdB: %s" %(i, numpy.max(ychannel), str_datetime)
1342 1342 axes = self.axesList[i]
1343 1343 axes.pline(x, ychannel,
1344 1344 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1345 1345 xlabel=xlabel, ylabel=ylabel, title=title, grid='both')
1346 1346
1347 1347
1348 1348 self.draw()
1349 1349
1350 1350 if save:
1351 1351 date = thisDatetime.strftime("%Y%m%d_%H%M%S")
1352 1352 if figfile == None:
1353 1353 figfile = self.getFilename(name = date)
1354 1354
1355 1355 self.saveFigure(figpath, figfile)
1356 1356
1357 1357 self.counter_imagwr += 1
1358 1358 if (ftp and (self.counter_imagwr==wr_period)):
1359 1359 ftp_filename = os.path.join(figpath,figfile)
1360 1360 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
1361 1361 self.counter_imagwr = 0
1362 1362
1363 1363
1364 1364 class RTIfromSpectraHeis(Figure):
1365 1365
1366 1366 __isConfig = None
1367 1367 __nsubplots = None
1368 1368
1369 1369 PREFIX = 'rtinoise'
1370 1370
1371 1371 def __init__(self):
1372 1372
1373 1373 self.timerange = 24*60*60
1374 1374 self.__isConfig = False
1375 1375 self.__nsubplots = 1
1376 1376
1377 1377 self.WIDTH = 820
1378 1378 self.HEIGHT = 200
1379 1379 self.WIDTHPROF = 120
1380 1380 self.HEIGHTPROF = 0
1381 1381 self.counter_imagwr = 0
1382 1382 self.xdata = None
1383 1383 self.ydata = None
1384 1384
1385 1385 def getSubplots(self):
1386 1386
1387 1387 ncol = 1
1388 1388 nrow = 1
1389 1389
1390 1390 return nrow, ncol
1391 1391
1392 1392 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1393 1393
1394 1394 self.__showprofile = showprofile
1395 1395 self.nplots = nplots
1396 1396
1397 1397 ncolspan = 7
1398 1398 colspan = 6
1399 1399 self.__nsubplots = 2
1400 1400
1401 1401 self.createFigure(id = id,
1402 1402 wintitle = wintitle,
1403 1403 widthplot = self.WIDTH+self.WIDTHPROF,
1404 1404 heightplot = self.HEIGHT+self.HEIGHTPROF,
1405 1405 show = show)
1406 1406
1407 1407 nrow, ncol = self.getSubplots()
1408 1408
1409 1409 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1410 1410
1411 1411
1412 1412 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1413 1413 xmin=None, xmax=None, ymin=None, ymax=None,
1414 1414 timerange=None,
1415 1415 save=False, figpath='./', figfile=None, ftp=False, wr_period=1, show=True,
1416 1416 server=None, folder=None, username=None, password=None):
1417 1417
1418 1418 if channelList == None:
1419 1419 channelIndexList = dataOut.channelIndexList
1420 1420 channelList = dataOut.channelList
1421 1421 else:
1422 1422 channelIndexList = []
1423 1423 for channel in channelList:
1424 1424 if channel not in dataOut.channelList:
1425 1425 raise ValueError, "Channel %d is not in dataOut.channelList"
1426 1426 channelIndexList.append(dataOut.channelList.index(channel))
1427 1427
1428 1428 if timerange != None:
1429 1429 self.timerange = timerange
1430 1430
1431 1431 tmin = None
1432 1432 tmax = None
1433 1433 x = dataOut.getTimeRange()
1434 1434 y = dataOut.getHeiRange()
1435 1435
1436 factor = 1
1437 data = dataOut.data_spc/factor
1436 #factor = 1
1437 data = dataOut.data_spc#/factor
1438 1438 data = numpy.average(data,axis=1)
1439 1439 datadB = 10*numpy.log10(data)
1440 1440
1441 1441 # factor = dataOut.normFactor
1442 1442 # noise = dataOut.getNoise()/factor
1443 1443 # noisedB = 10*numpy.log10(noise)
1444 1444
1445 1445 #thisDatetime = dataOut.datatime
1446 1446 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1447 1447 title = wintitle + " RTI: %s" %(thisDatetime.strftime("%d-%b-%Y"))
1448 1448 xlabel = "Local Time"
1449 1449 ylabel = "Intensity (dB)"
1450 1450
1451 1451 if not self.__isConfig:
1452 1452
1453 1453 nplots = 1
1454 1454
1455 1455 self.setup(id=id,
1456 1456 nplots=nplots,
1457 1457 wintitle=wintitle,
1458 1458 showprofile=showprofile,
1459 1459 show=show)
1460 1460
1461 1461 tmin, tmax = self.getTimeLim(x, xmin, xmax)
1462 1462 if ymin == None: ymin = numpy.nanmin(datadB)
1463 1463 if ymax == None: ymax = numpy.nanmax(datadB)
1464 1464
1465 1465 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1466 1466 self.__isConfig = True
1467 1467
1468 1468 self.xdata = numpy.array([])
1469 1469 self.ydata = numpy.array([])
1470 1470
1471 1471 self.setWinTitle(title)
1472 1472
1473 1473
1474 1474 # title = "RTI %s" %(thisDatetime.strftime("%d-%b-%Y"))
1475 1475 title = "RTI - %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1476 1476
1477 1477 legendlabels = ["channel %d"%idchannel for idchannel in channelList]
1478 1478 axes = self.axesList[0]
1479 1479
1480 1480 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1481 1481
1482 1482 if len(self.ydata)==0:
1483 1483 self.ydata = datadB[channelIndexList].reshape(-1,1)
1484 1484 else:
1485 1485 self.ydata = numpy.hstack((self.ydata, datadB[channelIndexList].reshape(-1,1)))
1486 1486
1487 1487
1488 1488 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1489 1489 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax,
1490 1490 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='.', markersize=8, linestyle="solid", grid='both',
1491 1491 XAxisAsTime=True
1492 1492 )
1493 1493
1494 1494 self.draw()
1495 1495
1496 1496 if save:
1497 1497
1498 1498 if figfile == None:
1499 1499 figfile = self.getFilename(name = self.name)
1500 1500
1501 1501 self.saveFigure(figpath, figfile)
1502 1502
1503 1503 self.counter_imagwr += 1
1504 1504 if (ftp and (self.counter_imagwr==wr_period)):
1505 1505 ftp_filename = os.path.join(figpath,figfile)
1506 1506 self.sendByFTP_Thread(ftp_filename, server, folder, username, password)
1507 1507 self.counter_imagwr = 0
1508 1508
1509 1509 if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax:
1510 1510 self.__isConfig = False
1511 1511 del self.xdata
1512 1512 del self.ydata
1513 1513
1514 1514
1515 1515 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1917 +1,1960
1 1 '''
2 2
3 3 $Author: dsuarez $
4 4 $Id: Processor.py 1 2012-11-12 18:56:07Z dsuarez $
5 5 '''
6 6 import os
7 7 import numpy
8 8 import datetime
9 9 import time
10 10 import math
11 11 from jrodata import *
12 12 from jrodataIO import *
13 13 from jroplot import *
14 14
15 15 try:
16 16 import cfunctions
17 17 except:
18 18 pass
19 19
20 20 class ProcessingUnit:
21 21
22 22 """
23 23 Esta es la clase base para el procesamiento de datos.
24 24
25 25 Contiene el metodo "call" para llamar operaciones. Las operaciones pueden ser:
26 26 - Metodos internos (callMethod)
27 27 - Objetos del tipo Operation (callObject). Antes de ser llamados, estos objetos
28 28 tienen que ser agreagados con el metodo "add".
29 29
30 30 """
31 31 # objeto de datos de entrada (Voltage, Spectra o Correlation)
32 32 dataIn = None
33 33
34 34 # objeto de datos de entrada (Voltage, Spectra o Correlation)
35 35 dataOut = None
36 36
37 37
38 38 objectDict = None
39 39
40 40 def __init__(self):
41 41
42 42 self.objectDict = {}
43 43
44 44 def init(self):
45 45
46 46 raise ValueError, "Not implemented"
47 47
48 48 def addOperation(self, object, objId):
49 49
50 50 """
51 51 Agrega el objeto "object" a la lista de objetos "self.objectList" y retorna el
52 52 identificador asociado a este objeto.
53 53
54 54 Input:
55 55
56 56 object : objeto de la clase "Operation"
57 57
58 58 Return:
59 59
60 60 objId : identificador del objeto, necesario para ejecutar la operacion
61 61 """
62 62
63 63 self.objectDict[objId] = object
64 64
65 65 return objId
66 66
67 67 def operation(self, **kwargs):
68 68
69 69 """
70 70 Operacion directa sobre la data (dataOut.data). Es necesario actualizar los valores de los
71 71 atributos del objeto dataOut
72 72
73 73 Input:
74 74
75 75 **kwargs : Diccionario de argumentos de la funcion a ejecutar
76 76 """
77 77
78 78 raise ValueError, "ImplementedError"
79 79
80 80 def callMethod(self, name, **kwargs):
81 81
82 82 """
83 83 Ejecuta el metodo con el nombre "name" y con argumentos **kwargs de la propia clase.
84 84
85 85 Input:
86 86 name : nombre del metodo a ejecutar
87 87
88 88 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
89 89
90 90 """
91 91 if name != 'run':
92 92
93 93 if name == 'init' and self.dataIn.isEmpty():
94 94 self.dataOut.flagNoData = True
95 95 return False
96 96
97 97 if name != 'init' and self.dataOut.isEmpty():
98 98 return False
99 99
100 100 methodToCall = getattr(self, name)
101 101
102 102 methodToCall(**kwargs)
103 103
104 104 if name != 'run':
105 105 return True
106 106
107 107 if self.dataOut.isEmpty():
108 108 return False
109 109
110 110 return True
111 111
112 112 def callObject(self, objId, **kwargs):
113 113
114 114 """
115 115 Ejecuta la operacion asociada al identificador del objeto "objId"
116 116
117 117 Input:
118 118
119 119 objId : identificador del objeto a ejecutar
120 120
121 121 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
122 122
123 123 Return:
124 124
125 125 None
126 126 """
127 127
128 128 if self.dataOut.isEmpty():
129 129 return False
130 130
131 131 object = self.objectDict[objId]
132 132
133 133 object.run(self.dataOut, **kwargs)
134 134
135 135 return True
136 136
137 137 def call(self, operationConf, **kwargs):
138 138
139 139 """
140 140 Return True si ejecuta la operacion "operationConf.name" con los
141 141 argumentos "**kwargs". False si la operacion no se ha ejecutado.
142 142 La operacion puede ser de dos tipos:
143 143
144 144 1. Un metodo propio de esta clase:
145 145
146 146 operation.type = "self"
147 147
148 148 2. El metodo "run" de un objeto del tipo Operation o de un derivado de ella:
149 149 operation.type = "other".
150 150
151 151 Este objeto de tipo Operation debe de haber sido agregado antes con el metodo:
152 152 "addOperation" e identificado con el operation.id
153 153
154 154
155 155 con el id de la operacion.
156 156
157 157 Input:
158 158
159 159 Operation : Objeto del tipo operacion con los atributos: name, type y id.
160 160
161 161 """
162 162
163 163 if operationConf.type == 'self':
164 164 sts = self.callMethod(operationConf.name, **kwargs)
165 165
166 166 if operationConf.type == 'other':
167 167 sts = self.callObject(operationConf.id, **kwargs)
168 168
169 169 return sts
170 170
171 171 def setInput(self, dataIn):
172 172
173 173 self.dataIn = dataIn
174 174
175 175 def getOutput(self):
176 176
177 177 return self.dataOut
178 178
179 179 class Operation():
180 180
181 181 """
182 182 Clase base para definir las operaciones adicionales que se pueden agregar a la clase ProcessingUnit
183 183 y necesiten acumular informacion previa de los datos a procesar. De preferencia usar un buffer de
184 184 acumulacion dentro de esta clase
185 185
186 186 Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer)
187 187
188 188 """
189 189
190 190 __buffer = None
191 191 __isConfig = False
192 192
193 193 def __init__(self):
194 194
195 195 pass
196 196
197 197 def run(self, dataIn, **kwargs):
198 198
199 199 """
200 200 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los atributos del objeto dataIn.
201 201
202 202 Input:
203 203
204 204 dataIn : objeto del tipo JROData
205 205
206 206 Return:
207 207
208 208 None
209 209
210 210 Affected:
211 211 __buffer : buffer de recepcion de datos.
212 212
213 213 """
214 214
215 215 raise ValueError, "ImplementedError"
216 216
217 217 class VoltageProc(ProcessingUnit):
218 218
219 219
220 220 def __init__(self):
221 221
222 222 self.objectDict = {}
223 223 self.dataOut = Voltage()
224 224 self.flip = 1
225 225
226 226 def init(self):
227 227
228 228 self.dataOut.copy(self.dataIn)
229 229 # No necesita copiar en cada init() los atributos de dataIn
230 230 # la copia deberia hacerse por cada nuevo bloque de datos
231 231
232 232 def selectChannels(self, channelList):
233 233
234 234 channelIndexList = []
235 235
236 236 for channel in channelList:
237 237 index = self.dataOut.channelList.index(channel)
238 238 channelIndexList.append(index)
239 239
240 240 self.selectChannelsByIndex(channelIndexList)
241 241
242 242 def selectChannelsByIndex(self, channelIndexList):
243 243 """
244 244 Selecciona un bloque de datos en base a canales segun el channelIndexList
245 245
246 246 Input:
247 247 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
248 248
249 249 Affected:
250 250 self.dataOut.data
251 251 self.dataOut.channelIndexList
252 252 self.dataOut.nChannels
253 253 self.dataOut.m_ProcessingHeader.totalSpectra
254 254 self.dataOut.systemHeaderObj.numChannels
255 255 self.dataOut.m_ProcessingHeader.blockSize
256 256
257 257 Return:
258 258 None
259 259 """
260 260
261 261 for channelIndex in channelIndexList:
262 262 if channelIndex not in self.dataOut.channelIndexList:
263 263 print channelIndexList
264 264 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
265 265
266 266 nChannels = len(channelIndexList)
267 267
268 268 data = self.dataOut.data[channelIndexList,:]
269 269
270 270 self.dataOut.data = data
271 271 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
272 272 # self.dataOut.nChannels = nChannels
273 273
274 274 return 1
275 275
276 276 def selectHeights(self, minHei=None, maxHei=None):
277 277 """
278 278 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
279 279 minHei <= height <= maxHei
280 280
281 281 Input:
282 282 minHei : valor minimo de altura a considerar
283 283 maxHei : valor maximo de altura a considerar
284 284
285 285 Affected:
286 286 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
287 287
288 288 Return:
289 289 1 si el metodo se ejecuto con exito caso contrario devuelve 0
290 290 """
291 291
292 292 if minHei == None:
293 293 minHei = self.dataOut.heightList[0]
294 294
295 295 if maxHei == None:
296 296 maxHei = self.dataOut.heightList[-1]
297 297
298 298 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
299 299 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
300 300
301 301
302 302 if (maxHei > self.dataOut.heightList[-1]):
303 303 maxHei = self.dataOut.heightList[-1]
304 304 # raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
305 305
306 306 minIndex = 0
307 307 maxIndex = 0
308 308 heights = self.dataOut.heightList
309 309
310 310 inda = numpy.where(heights >= minHei)
311 311 indb = numpy.where(heights <= maxHei)
312 312
313 313 try:
314 314 minIndex = inda[0][0]
315 315 except:
316 316 minIndex = 0
317 317
318 318 try:
319 319 maxIndex = indb[0][-1]
320 320 except:
321 321 maxIndex = len(heights)
322 322
323 323 self.selectHeightsByIndex(minIndex, maxIndex)
324 324
325 325 return 1
326 326
327 327
328 328 def selectHeightsByIndex(self, minIndex, maxIndex):
329 329 """
330 330 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
331 331 minIndex <= index <= maxIndex
332 332
333 333 Input:
334 334 minIndex : valor de indice minimo de altura a considerar
335 335 maxIndex : valor de indice maximo de altura a considerar
336 336
337 337 Affected:
338 338 self.dataOut.data
339 339 self.dataOut.heightList
340 340
341 341 Return:
342 342 1 si el metodo se ejecuto con exito caso contrario devuelve 0
343 343 """
344 344
345 345 if (minIndex < 0) or (minIndex > maxIndex):
346 346 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
347 347
348 348 if (maxIndex >= self.dataOut.nHeights):
349 349 maxIndex = self.dataOut.nHeights-1
350 350 # raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
351 351
352 352 nHeights = maxIndex - minIndex + 1
353 353
354 354 #voltage
355 355 data = self.dataOut.data[:,minIndex:maxIndex+1]
356 356
357 357 firstHeight = self.dataOut.heightList[minIndex]
358 358
359 359 self.dataOut.data = data
360 360 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
361 361
362 362 return 1
363 363
364 364
365 365 def filterByHeights(self, window):
366 366 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
367 367
368 368 if window == None:
369 369 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
370 370
371 371 newdelta = deltaHeight * window
372 372 r = self.dataOut.data.shape[1] % window
373 373 buffer = self.dataOut.data[:,0:self.dataOut.data.shape[1]-r]
374 374 buffer = buffer.reshape(self.dataOut.data.shape[0],self.dataOut.data.shape[1]/window,window)
375 375 buffer = numpy.sum(buffer,2)
376 376 self.dataOut.data = buffer
377 377 self.dataOut.heightList = numpy.arange(self.dataOut.heightList[0],newdelta*(self.dataOut.nHeights-r)/window,newdelta)
378 378 self.dataOut.windowOfFilter = window
379 379
380 380 def deFlip(self):
381 381 self.dataOut.data *= self.flip
382 382 self.flip *= -1.
383 383
384 384 def setRadarFrequency(self, frequency=None):
385 385 if frequency != None:
386 386 self.dataOut.frequency = frequency
387 387
388 388 return 1
389 389
390 390 class CohInt(Operation):
391 391
392 392 __isConfig = False
393 393
394 394 __profIndex = 0
395 395 __withOverapping = False
396 396
397 397 __byTime = False
398 398 __initime = None
399 399 __lastdatatime = None
400 400 __integrationtime = None
401 401
402 402 __buffer = None
403 403
404 404 __dataReady = False
405 405
406 406 n = None
407 407
408 408
409 409 def __init__(self):
410 410
411 411 self.__isConfig = False
412 412
413 413 def setup(self, n=None, timeInterval=None, overlapping=False):
414 414 """
415 415 Set the parameters of the integration class.
416 416
417 417 Inputs:
418 418
419 419 n : Number of coherent integrations
420 420 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
421 421 overlapping :
422 422
423 423 """
424 424
425 425 self.__initime = None
426 426 self.__lastdatatime = 0
427 427 self.__buffer = None
428 428 self.__dataReady = False
429 429
430 430
431 431 if n == None and timeInterval == None:
432 432 raise ValueError, "n or timeInterval should be specified ..."
433 433
434 434 if n != None:
435 435 self.n = n
436 436 self.__byTime = False
437 437 else:
438 438 self.__integrationtime = timeInterval * 60. #if (type(timeInterval)!=integer) -> change this line
439 439 self.n = 9999
440 440 self.__byTime = True
441 441
442 442 if overlapping:
443 443 self.__withOverapping = True
444 444 self.__buffer = None
445 445 else:
446 446 self.__withOverapping = False
447 447 self.__buffer = 0
448 448
449 449 self.__profIndex = 0
450 450
451 451 def putData(self, data):
452 452
453 453 """
454 454 Add a profile to the __buffer and increase in one the __profileIndex
455 455
456 456 """
457 457
458 458 if not self.__withOverapping:
459 459 self.__buffer += data.copy()
460 460 self.__profIndex += 1
461 461 return
462 462
463 463 #Overlapping data
464 464 nChannels, nHeis = data.shape
465 465 data = numpy.reshape(data, (1, nChannels, nHeis))
466 466
467 467 #If the buffer is empty then it takes the data value
468 468 if self.__buffer == None:
469 469 self.__buffer = data
470 470 self.__profIndex += 1
471 471 return
472 472
473 473 #If the buffer length is lower than n then stakcing the data value
474 474 if self.__profIndex < self.n:
475 475 self.__buffer = numpy.vstack((self.__buffer, data))
476 476 self.__profIndex += 1
477 477 return
478 478
479 479 #If the buffer length is equal to n then replacing the last buffer value with the data value
480 480 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
481 481 self.__buffer[self.n-1] = data
482 482 self.__profIndex = self.n
483 483 return
484 484
485 485
486 486 def pushData(self):
487 487 """
488 488 Return the sum of the last profiles and the profiles used in the sum.
489 489
490 490 Affected:
491 491
492 492 self.__profileIndex
493 493
494 494 """
495 495
496 496 if not self.__withOverapping:
497 497 data = self.__buffer
498 498 n = self.__profIndex
499 499
500 500 self.__buffer = 0
501 501 self.__profIndex = 0
502 502
503 503 return data, n
504 504
505 505 #Integration with Overlapping
506 506 data = numpy.sum(self.__buffer, axis=0)
507 507 n = self.__profIndex
508 508
509 509 return data, n
510 510
511 511 def byProfiles(self, data):
512 512
513 513 self.__dataReady = False
514 514 avgdata = None
515 515 n = None
516 516
517 517 self.putData(data)
518 518
519 519 if self.__profIndex == self.n:
520 520
521 521 avgdata, n = self.pushData()
522 522 self.__dataReady = True
523 523
524 524 return avgdata
525 525
526 526 def byTime(self, data, datatime):
527 527
528 528 self.__dataReady = False
529 529 avgdata = None
530 530 n = None
531 531
532 532 self.putData(data)
533 533
534 534 if (datatime - self.__initime) >= self.__integrationtime:
535 535 avgdata, n = self.pushData()
536 536 self.n = n
537 537 self.__dataReady = True
538 538
539 539 return avgdata
540 540
541 541 def integrate(self, data, datatime=None):
542 542
543 543 if self.__initime == None:
544 544 self.__initime = datatime
545 545
546 546 if self.__byTime:
547 547 avgdata = self.byTime(data, datatime)
548 548 else:
549 549 avgdata = self.byProfiles(data)
550 550
551 551
552 552 self.__lastdatatime = datatime
553 553
554 554 if avgdata == None:
555 555 return None, None
556 556
557 557 avgdatatime = self.__initime
558 558
559 559 deltatime = datatime -self.__lastdatatime
560 560
561 561 if not self.__withOverapping:
562 562 self.__initime = datatime
563 563 else:
564 564 self.__initime += deltatime
565 565
566 566 return avgdata, avgdatatime
567 567
568 568 def run(self, dataOut, **kwargs):
569 569
570 570 if not self.__isConfig:
571 571 self.setup(**kwargs)
572 572 self.__isConfig = True
573 573
574 574 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
575 575
576 576 # dataOut.timeInterval *= n
577 577 dataOut.flagNoData = True
578 578
579 579 if self.__dataReady:
580 580 dataOut.data = avgdata
581 581 dataOut.nCohInt *= self.n
582 582 dataOut.utctime = avgdatatime
583 583 dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
584 584 dataOut.flagNoData = False
585 585
586 586
587 587 class Decoder(Operation):
588 588
589 589 __isConfig = False
590 590 __profIndex = 0
591 591
592 592 code = None
593 593
594 594 nCode = None
595 595 nBaud = None
596 596
597 597 def __init__(self):
598 598
599 599 self.__isConfig = False
600 600
601 601 def setup(self, code, shape):
602 602
603 603 self.__profIndex = 0
604 604
605 605 self.code = code
606 606
607 607 self.nCode = len(code)
608 608 self.nBaud = len(code[0])
609 609
610 610 self.__nChannels, self.__nHeis = shape
611 611
612 612 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
613 613
614 614 __codeBuffer[:,0:self.nBaud] = self.code
615 615
616 616 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
617 617
618 618 self.ndatadec = self.__nHeis - self.nBaud + 1
619 619
620 620 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
621 621
622 622 def convolutionInFreq(self, data):
623 623
624 624 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
625 625
626 626 fft_data = numpy.fft.fft(data, axis=1)
627 627
628 628 conv = fft_data*fft_code
629 629
630 630 data = numpy.fft.ifft(conv,axis=1)
631 631
632 632 datadec = data[:,:-self.nBaud+1]
633 633
634 634 return datadec
635 635
636 636 def convolutionInFreqOpt(self, data):
637 637
638 638 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
639 639
640 640 data = cfunctions.decoder(fft_code, data)
641 641
642 642 datadec = data[:,:-self.nBaud+1]
643 643
644 644 return datadec
645 645
646 646 def convolutionInTime(self, data):
647 647
648 648 code = self.code[self.__profIndex]
649 649
650 650 for i in range(self.__nChannels):
651 651 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='valid')
652 652
653 653 return self.datadecTime
654 654
655 655 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0):
656 656
657 657 if not self.__isConfig:
658 658
659 659 if code == None:
660 660 code = dataOut.code
661 661 else:
662 662 code = numpy.array(code).reshape(nCode,nBaud)
663 663 dataOut.code = code
664 664 dataOut.nCode = nCode
665 665 dataOut.nBaud = nBaud
666 666
667 667 if code == None:
668 668 return 1
669 669
670 670 self.setup(code, dataOut.data.shape)
671 671 self.__isConfig = True
672 672
673 673 if mode == 0:
674 674 datadec = self.convolutionInTime(dataOut.data)
675 675
676 676 if mode == 1:
677 677 datadec = self.convolutionInFreq(dataOut.data)
678 678
679 679 if mode == 2:
680 680 datadec = self.convolutionInFreqOpt(dataOut.data)
681 681
682 682 dataOut.data = datadec
683 683
684 684 dataOut.heightList = dataOut.heightList[0:self.ndatadec]
685 685
686 686 dataOut.flagDecodeData = True #asumo q la data no esta decodificada
687 687
688 688 if self.__profIndex == self.nCode-1:
689 689 self.__profIndex = 0
690 690 return 1
691 691
692 692 self.__profIndex += 1
693 693
694 694 return 1
695 695 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
696 696
697 697
698 698
699 699 class SpectraProc(ProcessingUnit):
700 700
701 701 def __init__(self):
702 702
703 703 self.objectDict = {}
704 704 self.buffer = None
705 705 self.firstdatatime = None
706 706 self.profIndex = 0
707 707 self.dataOut = Spectra()
708 708
709 709 def __updateObjFromInput(self):
710 710
711 711 self.dataOut.timeZone = self.dataIn.timeZone
712 712 self.dataOut.dstFlag = self.dataIn.dstFlag
713 713 self.dataOut.errorCount = self.dataIn.errorCount
714 714 self.dataOut.useLocalTime = self.dataIn.useLocalTime
715 715
716 716 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
717 717 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
718 718 self.dataOut.channelList = self.dataIn.channelList
719 719 self.dataOut.heightList = self.dataIn.heightList
720 720 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
721 721 # self.dataOut.nHeights = self.dataIn.nHeights
722 722 # self.dataOut.nChannels = self.dataIn.nChannels
723 723 self.dataOut.nBaud = self.dataIn.nBaud
724 724 self.dataOut.nCode = self.dataIn.nCode
725 725 self.dataOut.code = self.dataIn.code
726 726 self.dataOut.nProfiles = self.dataOut.nFFTPoints
727 727 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
728 728 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
729 729 self.dataOut.utctime = self.firstdatatime
730 730 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
731 731 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
732 732 # self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
733 733 self.dataOut.nCohInt = self.dataIn.nCohInt
734 734 self.dataOut.nIncohInt = 1
735 735 self.dataOut.ippSeconds = self.dataIn.ippSeconds
736 736 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
737 737
738 738 self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt
739 739 self.dataOut.frequency = self.dataIn.frequency
740 740 self.dataOut.realtime = self.dataIn.realtime
741 741
742 742 def __getFft(self):
743 743 """
744 744 Convierte valores de Voltaje a Spectra
745 745
746 746 Affected:
747 747 self.dataOut.data_spc
748 748 self.dataOut.data_cspc
749 749 self.dataOut.data_dc
750 750 self.dataOut.heightList
751 751 self.profIndex
752 752 self.buffer
753 753 self.dataOut.flagNoData
754 754 """
755 755 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
756 756 fft_volt = fft_volt.astype(numpy.dtype('complex'))
757 757 dc = fft_volt[:,0,:]
758 758
759 759 #calculo de self-spectra
760 760 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
761 761 spc = fft_volt * numpy.conjugate(fft_volt)
762 762 spc = spc.real
763 763
764 764 blocksize = 0
765 765 blocksize += dc.size
766 766 blocksize += spc.size
767 767
768 768 cspc = None
769 769 pairIndex = 0
770 770 if self.dataOut.pairsList != None:
771 771 #calculo de cross-spectra
772 772 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
773 773 for pair in self.dataOut.pairsList:
774 774 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
775 775 pairIndex += 1
776 776 blocksize += cspc.size
777 777
778 778 self.dataOut.data_spc = spc
779 779 self.dataOut.data_cspc = cspc
780 780 self.dataOut.data_dc = dc
781 781 self.dataOut.blockSize = blocksize
782 782 self.dataOut.flagShiftFFT = False
783 783
784 784 def init(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None):
785 785
786 786 self.dataOut.flagNoData = True
787 787
788 788 if self.dataIn.type == "Spectra":
789 789 self.dataOut.copy(self.dataIn)
790 790 return
791 791
792 792 if self.dataIn.type == "Voltage":
793 793
794 794 if nFFTPoints == None:
795 795 raise ValueError, "This SpectraProc.init() need nFFTPoints input variable"
796 796
797 797 if pairsList == None:
798 798 nPairs = 0
799 799 else:
800 800 nPairs = len(pairsList)
801 801
802 802 if ippFactor == None:
803 803 ippFactor = 1
804 804 self.dataOut.ippFactor = ippFactor
805 805
806 806 self.dataOut.nFFTPoints = nFFTPoints
807 807 self.dataOut.pairsList = pairsList
808 808 self.dataOut.nPairs = nPairs
809 809
810 810 if self.buffer == None:
811 811 self.buffer = numpy.zeros((self.dataIn.nChannels,
812 812 nProfiles,
813 813 self.dataIn.nHeights),
814 814 dtype='complex')
815 815
816 816
817 817 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
818 818 self.profIndex += 1
819 819
820 820 if self.firstdatatime == None:
821 821 self.firstdatatime = self.dataIn.utctime
822 822
823 823 if self.profIndex == nProfiles:
824 824 self.__updateObjFromInput()
825 825 self.__getFft()
826 826
827 827 self.dataOut.flagNoData = False
828 828
829 829 self.buffer = None
830 830 self.firstdatatime = None
831 831 self.profIndex = 0
832 832
833 833 return
834 834
835 835 raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
836 836
837 837 def selectChannels(self, channelList):
838 838
839 839 channelIndexList = []
840 840
841 841 for channel in channelList:
842 842 index = self.dataOut.channelList.index(channel)
843 843 channelIndexList.append(index)
844 844
845 845 self.selectChannelsByIndex(channelIndexList)
846 846
847 847 def selectChannelsByIndex(self, channelIndexList):
848 848 """
849 849 Selecciona un bloque de datos en base a canales segun el channelIndexList
850 850
851 851 Input:
852 852 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
853 853
854 854 Affected:
855 855 self.dataOut.data_spc
856 856 self.dataOut.channelIndexList
857 857 self.dataOut.nChannels
858 858
859 859 Return:
860 860 None
861 861 """
862 862
863 863 for channelIndex in channelIndexList:
864 864 if channelIndex not in self.dataOut.channelIndexList:
865 865 print channelIndexList
866 866 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
867 867
868 868 nChannels = len(channelIndexList)
869 869
870 870 data_spc = self.dataOut.data_spc[channelIndexList,:]
871 871
872 872 self.dataOut.data_spc = data_spc
873 873 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
874 874 # self.dataOut.nChannels = nChannels
875 875
876 876 return 1
877 877
878 878 def selectHeights(self, minHei, maxHei):
879 879 """
880 880 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
881 881 minHei <= height <= maxHei
882 882
883 883 Input:
884 884 minHei : valor minimo de altura a considerar
885 885 maxHei : valor maximo de altura a considerar
886 886
887 887 Affected:
888 888 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
889 889
890 890 Return:
891 891 1 si el metodo se ejecuto con exito caso contrario devuelve 0
892 892 """
893 893 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
894 894 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
895 895
896 896 if (maxHei > self.dataOut.heightList[-1]):
897 897 maxHei = self.dataOut.heightList[-1]
898 898 # raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
899 899
900 900 minIndex = 0
901 901 maxIndex = 0
902 902 heights = self.dataOut.heightList
903 903
904 904 inda = numpy.where(heights >= minHei)
905 905 indb = numpy.where(heights <= maxHei)
906 906
907 907 try:
908 908 minIndex = inda[0][0]
909 909 except:
910 910 minIndex = 0
911 911
912 912 try:
913 913 maxIndex = indb[0][-1]
914 914 except:
915 915 maxIndex = len(heights)
916 916
917 917 self.selectHeightsByIndex(minIndex, maxIndex)
918 918
919 919 return 1
920 920
921 921
922 922 def selectHeightsByIndex(self, minIndex, maxIndex):
923 923 """
924 924 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
925 925 minIndex <= index <= maxIndex
926 926
927 927 Input:
928 928 minIndex : valor de indice minimo de altura a considerar
929 929 maxIndex : valor de indice maximo de altura a considerar
930 930
931 931 Affected:
932 932 self.dataOut.data_spc
933 933 self.dataOut.data_cspc
934 934 self.dataOut.data_dc
935 935 self.dataOut.heightList
936 936
937 937 Return:
938 938 1 si el metodo se ejecuto con exito caso contrario devuelve 0
939 939 """
940 940
941 941 if (minIndex < 0) or (minIndex > maxIndex):
942 942 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
943 943
944 944 if (maxIndex >= self.dataOut.nHeights):
945 945 maxIndex = self.dataOut.nHeights-1
946 946 # raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
947 947
948 948 nHeights = maxIndex - minIndex + 1
949 949
950 950 #Spectra
951 951 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
952 952
953 953 data_cspc = None
954 954 if self.dataOut.data_cspc != None:
955 955 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
956 956
957 957 data_dc = None
958 958 if self.dataOut.data_dc != None:
959 959 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
960 960
961 961 self.dataOut.data_spc = data_spc
962 962 self.dataOut.data_cspc = data_cspc
963 963 self.dataOut.data_dc = data_dc
964 964
965 965 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
966 966
967 967 return 1
968 968
969 969 def removeDC(self, mode = 2):
970 970 jspectra = self.dataOut.data_spc
971 971 jcspectra = self.dataOut.data_cspc
972 972
973 973
974 974 num_chan = jspectra.shape[0]
975 975 num_hei = jspectra.shape[2]
976 976
977 977 if jcspectra != None:
978 978 jcspectraExist = True
979 979 num_pairs = jcspectra.shape[0]
980 980 else: jcspectraExist = False
981 981
982 982 freq_dc = jspectra.shape[1]/2
983 983 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
984 984
985 985 if ind_vel[0]<0:
986 986 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
987 987
988 988 if mode == 1:
989 989 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
990 990
991 991 if jcspectraExist:
992 992 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
993 993
994 994 if mode == 2:
995 995
996 996 vel = numpy.array([-2,-1,1,2])
997 997 xx = numpy.zeros([4,4])
998 998
999 999 for fil in range(4):
1000 1000 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1001 1001
1002 1002 xx_inv = numpy.linalg.inv(xx)
1003 1003 xx_aux = xx_inv[0,:]
1004 1004
1005 1005 for ich in range(num_chan):
1006 1006 yy = jspectra[ich,ind_vel,:]
1007 1007 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1008 1008
1009 1009 junkid = jspectra[ich,freq_dc,:]<=0
1010 1010 cjunkid = sum(junkid)
1011 1011
1012 1012 if cjunkid.any():
1013 1013 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1014 1014
1015 1015 if jcspectraExist:
1016 1016 for ip in range(num_pairs):
1017 1017 yy = jcspectra[ip,ind_vel,:]
1018 1018 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
1019 1019
1020 1020
1021 1021 self.dataOut.data_spc = jspectra
1022 1022 self.dataOut.data_cspc = jcspectra
1023 1023
1024 1024 return 1
1025 1025
1026 1026 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
1027 1027
1028 1028 jspectra = self.dataOut.data_spc
1029 1029 jcspectra = self.dataOut.data_cspc
1030 1030 jnoise = self.dataOut.getNoise()
1031 1031 num_incoh = self.dataOut.nIncohInt
1032 1032
1033 1033 num_channel = jspectra.shape[0]
1034 1034 num_prof = jspectra.shape[1]
1035 1035 num_hei = jspectra.shape[2]
1036 1036
1037 1037 #hei_interf
1038 1038 if hei_interf == None:
1039 1039 count_hei = num_hei/2 #Como es entero no importa
1040 1040 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
1041 1041 hei_interf = numpy.asarray(hei_interf)[0]
1042 1042 #nhei_interf
1043 1043 if (nhei_interf == None):
1044 1044 nhei_interf = 5
1045 1045 if (nhei_interf < 1):
1046 1046 nhei_interf = 1
1047 1047 if (nhei_interf > count_hei):
1048 1048 nhei_interf = count_hei
1049 1049 if (offhei_interf == None):
1050 1050 offhei_interf = 0
1051 1051
1052 1052 ind_hei = range(num_hei)
1053 1053 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1054 1054 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1055 1055 mask_prof = numpy.asarray(range(num_prof))
1056 1056 num_mask_prof = mask_prof.size
1057 1057 comp_mask_prof = [0, num_prof/2]
1058 1058
1059 1059
1060 1060 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1061 1061 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1062 1062 jnoise = numpy.nan
1063 1063 noise_exist = jnoise[0] < numpy.Inf
1064 1064
1065 1065 #Subrutina de Remocion de la Interferencia
1066 1066 for ich in range(num_channel):
1067 1067 #Se ordena los espectros segun su potencia (menor a mayor)
1068 1068 power = jspectra[ich,mask_prof,:]
1069 1069 power = power[:,hei_interf]
1070 1070 power = power.sum(axis = 0)
1071 1071 psort = power.ravel().argsort()
1072 1072
1073 1073 #Se estima la interferencia promedio en los Espectros de Potencia empleando
1074 1074 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
1075 1075
1076 1076 if noise_exist:
1077 1077 # tmp_noise = jnoise[ich] / num_prof
1078 1078 tmp_noise = jnoise[ich]
1079 1079 junkspc_interf = junkspc_interf - tmp_noise
1080 1080 #junkspc_interf[:,comp_mask_prof] = 0
1081 1081
1082 1082 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
1083 1083 jspc_interf = jspc_interf.transpose()
1084 1084 #Calculando el espectro de interferencia promedio
1085 1085 noiseid = numpy.where(jspc_interf <= tmp_noise/ math.sqrt(num_incoh))
1086 1086 noiseid = noiseid[0]
1087 1087 cnoiseid = noiseid.size
1088 1088 interfid = numpy.where(jspc_interf > tmp_noise/ math.sqrt(num_incoh))
1089 1089 interfid = interfid[0]
1090 1090 cinterfid = interfid.size
1091 1091
1092 1092 if (cnoiseid > 0): jspc_interf[noiseid] = 0
1093 1093
1094 1094 #Expandiendo los perfiles a limpiar
1095 1095 if (cinterfid > 0):
1096 1096 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
1097 1097 new_interfid = numpy.asarray(new_interfid)
1098 1098 new_interfid = {x for x in new_interfid}
1099 1099 new_interfid = numpy.array(list(new_interfid))
1100 1100 new_cinterfid = new_interfid.size
1101 1101 else: new_cinterfid = 0
1102 1102
1103 1103 for ip in range(new_cinterfid):
1104 1104 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
1105 1105 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
1106 1106
1107 1107
1108 1108 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
1109 1109
1110 1110 #Removiendo la interferencia del punto de mayor interferencia
1111 1111 ListAux = jspc_interf[mask_prof].tolist()
1112 1112 maxid = ListAux.index(max(ListAux))
1113 1113
1114 1114
1115 1115 if cinterfid > 0:
1116 1116 for ip in range(cinterfid*(interf == 2) - 1):
1117 1117 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/math.sqrt(num_incoh))).nonzero()
1118 1118 cind = len(ind)
1119 1119
1120 1120 if (cind > 0):
1121 1121 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/math.sqrt(num_incoh))
1122 1122
1123 1123 ind = numpy.array([-2,-1,1,2])
1124 1124 xx = numpy.zeros([4,4])
1125 1125
1126 1126 for id1 in range(4):
1127 1127 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
1128 1128
1129 1129 xx_inv = numpy.linalg.inv(xx)
1130 1130 xx = xx_inv[:,0]
1131 1131 ind = (ind + maxid + num_mask_prof)%num_mask_prof
1132 1132 yy = jspectra[ich,mask_prof[ind],:]
1133 1133 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
1134 1134
1135 1135
1136 1136 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/math.sqrt(num_incoh))).nonzero()
1137 1137 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/math.sqrt(num_incoh))
1138 1138
1139 1139 #Remocion de Interferencia en el Cross Spectra
1140 1140 if jcspectra == None: return jspectra, jcspectra
1141 1141 num_pairs = jcspectra.size/(num_prof*num_hei)
1142 1142 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1143 1143
1144 1144 for ip in range(num_pairs):
1145 1145
1146 1146 #-------------------------------------------
1147 1147
1148 1148 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
1149 1149 cspower = cspower[:,hei_interf]
1150 1150 cspower = cspower.sum(axis = 0)
1151 1151
1152 1152 cspsort = cspower.ravel().argsort()
1153 1153 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
1154 1154 junkcspc_interf = junkcspc_interf.transpose()
1155 1155 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
1156 1156
1157 1157 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1158 1158
1159 1159 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
1160 1160 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
1161 1161 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
1162 1162
1163 1163 for iprof in range(num_prof):
1164 1164 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
1165 1165 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
1166 1166
1167 1167 #Removiendo la Interferencia
1168 1168 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
1169 1169
1170 1170 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1171 1171 maxid = ListAux.index(max(ListAux))
1172 1172
1173 1173 ind = numpy.array([-2,-1,1,2])
1174 1174 xx = numpy.zeros([4,4])
1175 1175
1176 1176 for id1 in range(4):
1177 1177 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
1178 1178
1179 1179 xx_inv = numpy.linalg.inv(xx)
1180 1180 xx = xx_inv[:,0]
1181 1181
1182 1182 ind = (ind + maxid + num_mask_prof)%num_mask_prof
1183 1183 yy = jcspectra[ip,mask_prof[ind],:]
1184 1184 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
1185 1185
1186 1186 #Guardar Resultados
1187 1187 self.dataOut.data_spc = jspectra
1188 1188 self.dataOut.data_cspc = jcspectra
1189 1189
1190 1190 return 1
1191 1191
1192 1192 def setRadarFrequency(self, frequency=None):
1193 1193 if frequency != None:
1194 1194 self.dataOut.frequency = frequency
1195 1195
1196 1196 return 1
1197
1198 def getNoise(self, minHei, maxHei):
1199
1200 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1201 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
1202
1203 if (maxHei > self.dataOut.heightList[-1]):
1204 maxHei = self.dataOut.heightList[-1]
1205
1206 minIndex = 0
1207 maxIndex = 0
1208 heights = self.dataOut.heightList
1209
1210 inda = numpy.where(heights >= minHei)
1211 indb = numpy.where(heights <= maxHei)
1212
1213 try:
1214 minIndex = inda[0][0]
1215 except:
1216 minIndex = 0
1217
1218 try:
1219 maxIndex = indb[0][-1]
1220 except:
1221 maxIndex = len(heights)
1222
1223 if (minIndex < 0) or (minIndex > maxIndex):
1224 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
1225
1226 if (maxIndex >= self.dataOut.nHeights):
1227 maxIndex = self.dataOut.nHeights-1
1228
1229 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
1230
1231 noise = numpy.zeros(self.dataOut.nChannels)
1232
1233 for channel in range(self.dataOut.nChannels):
1234 daux = data_spc[channel,:,:]
1235 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
1236
1237 self.dataOut.noise = noise.copy()
1238
1239 return 1
1197 1240
1198 1241
1199 1242 class IncohInt(Operation):
1200 1243
1201 1244
1202 1245 __profIndex = 0
1203 1246 __withOverapping = False
1204 1247
1205 1248 __byTime = False
1206 1249 __initime = None
1207 1250 __lastdatatime = None
1208 1251 __integrationtime = None
1209 1252
1210 1253 __buffer_spc = None
1211 1254 __buffer_cspc = None
1212 1255 __buffer_dc = None
1213 1256
1214 1257 __dataReady = False
1215 1258
1216 1259 __timeInterval = None
1217 1260
1218 1261 n = None
1219 1262
1220 1263
1221 1264
1222 1265 def __init__(self):
1223 1266
1224 1267 self.__isConfig = False
1225 1268
1226 1269 def setup(self, n=None, timeInterval=None, overlapping=False):
1227 1270 """
1228 1271 Set the parameters of the integration class.
1229 1272
1230 1273 Inputs:
1231 1274
1232 1275 n : Number of coherent integrations
1233 1276 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1234 1277 overlapping :
1235 1278
1236 1279 """
1237 1280
1238 1281 self.__initime = None
1239 1282 self.__lastdatatime = 0
1240 1283 self.__buffer_spc = None
1241 1284 self.__buffer_cspc = None
1242 1285 self.__buffer_dc = None
1243 1286 self.__dataReady = False
1244 1287
1245 1288
1246 1289 if n == None and timeInterval == None:
1247 1290 raise ValueError, "n or timeInterval should be specified ..."
1248 1291
1249 1292 if n != None:
1250 1293 self.n = n
1251 1294 self.__byTime = False
1252 1295 else:
1253 1296 self.__integrationtime = timeInterval #if (type(timeInterval)!=integer) -> change this line
1254 1297 self.n = 9999
1255 1298 self.__byTime = True
1256 1299
1257 1300 if overlapping:
1258 1301 self.__withOverapping = True
1259 1302 else:
1260 1303 self.__withOverapping = False
1261 1304 self.__buffer_spc = 0
1262 1305 self.__buffer_cspc = 0
1263 1306 self.__buffer_dc = 0
1264 1307
1265 1308 self.__profIndex = 0
1266 1309
1267 1310 def putData(self, data_spc, data_cspc, data_dc):
1268 1311
1269 1312 """
1270 1313 Add a profile to the __buffer_spc and increase in one the __profileIndex
1271 1314
1272 1315 """
1273 1316
1274 1317 if not self.__withOverapping:
1275 1318 self.__buffer_spc += data_spc
1276 1319
1277 1320 if data_cspc == None:
1278 1321 self.__buffer_cspc = None
1279 1322 else:
1280 1323 self.__buffer_cspc += data_cspc
1281 1324
1282 1325 if data_dc == None:
1283 1326 self.__buffer_dc = None
1284 1327 else:
1285 1328 self.__buffer_dc += data_dc
1286 1329
1287 1330 self.__profIndex += 1
1288 1331 return
1289 1332
1290 1333 #Overlapping data
1291 1334 nChannels, nFFTPoints, nHeis = data_spc.shape
1292 1335 data_spc = numpy.reshape(data_spc, (1, nChannels, nFFTPoints, nHeis))
1293 1336 if data_cspc != None:
1294 1337 data_cspc = numpy.reshape(data_cspc, (1, -1, nFFTPoints, nHeis))
1295 1338 if data_dc != None:
1296 1339 data_dc = numpy.reshape(data_dc, (1, -1, nHeis))
1297 1340
1298 1341 #If the buffer is empty then it takes the data value
1299 1342 if self.__buffer_spc == None:
1300 1343 self.__buffer_spc = data_spc
1301 1344
1302 1345 if data_cspc == None:
1303 1346 self.__buffer_cspc = None
1304 1347 else:
1305 1348 self.__buffer_cspc += data_cspc
1306 1349
1307 1350 if data_dc == None:
1308 1351 self.__buffer_dc = None
1309 1352 else:
1310 1353 self.__buffer_dc += data_dc
1311 1354
1312 1355 self.__profIndex += 1
1313 1356 return
1314 1357
1315 1358 #If the buffer length is lower than n then stakcing the data value
1316 1359 if self.__profIndex < self.n:
1317 1360 self.__buffer_spc = numpy.vstack((self.__buffer_spc, data_spc))
1318 1361
1319 1362 if data_cspc != None:
1320 1363 self.__buffer_cspc = numpy.vstack((self.__buffer_cspc, data_cspc))
1321 1364
1322 1365 if data_dc != None:
1323 1366 self.__buffer_dc = numpy.vstack((self.__buffer_dc, data_dc))
1324 1367
1325 1368 self.__profIndex += 1
1326 1369 return
1327 1370
1328 1371 #If the buffer length is equal to n then replacing the last buffer value with the data value
1329 1372 self.__buffer_spc = numpy.roll(self.__buffer_spc, -1, axis=0)
1330 1373 self.__buffer_spc[self.n-1] = data_spc
1331 1374
1332 1375 if data_cspc != None:
1333 1376 self.__buffer_cspc = numpy.roll(self.__buffer_cspc, -1, axis=0)
1334 1377 self.__buffer_cspc[self.n-1] = data_cspc
1335 1378
1336 1379 if data_dc != None:
1337 1380 self.__buffer_dc = numpy.roll(self.__buffer_dc, -1, axis=0)
1338 1381 self.__buffer_dc[self.n-1] = data_dc
1339 1382
1340 1383 self.__profIndex = self.n
1341 1384 return
1342 1385
1343 1386
1344 1387 def pushData(self):
1345 1388 """
1346 1389 Return the sum of the last profiles and the profiles used in the sum.
1347 1390
1348 1391 Affected:
1349 1392
1350 1393 self.__profileIndex
1351 1394
1352 1395 """
1353 1396 data_spc = None
1354 1397 data_cspc = None
1355 1398 data_dc = None
1356 1399
1357 1400 if not self.__withOverapping:
1358 1401 data_spc = self.__buffer_spc
1359 1402 data_cspc = self.__buffer_cspc
1360 1403 data_dc = self.__buffer_dc
1361 1404
1362 1405 n = self.__profIndex
1363 1406
1364 1407 self.__buffer_spc = 0
1365 1408 self.__buffer_cspc = 0
1366 1409 self.__buffer_dc = 0
1367 1410 self.__profIndex = 0
1368 1411
1369 1412 return data_spc, data_cspc, data_dc, n
1370 1413
1371 1414 #Integration with Overlapping
1372 1415 data_spc = numpy.sum(self.__buffer_spc, axis=0)
1373 1416
1374 1417 if self.__buffer_cspc != None:
1375 1418 data_cspc = numpy.sum(self.__buffer_cspc, axis=0)
1376 1419
1377 1420 if self.__buffer_dc != None:
1378 1421 data_dc = numpy.sum(self.__buffer_dc, axis=0)
1379 1422
1380 1423 n = self.__profIndex
1381 1424
1382 1425 return data_spc, data_cspc, data_dc, n
1383 1426
1384 1427 def byProfiles(self, *args):
1385 1428
1386 1429 self.__dataReady = False
1387 1430 avgdata_spc = None
1388 1431 avgdata_cspc = None
1389 1432 avgdata_dc = None
1390 1433 n = None
1391 1434
1392 1435 self.putData(*args)
1393 1436
1394 1437 if self.__profIndex == self.n:
1395 1438
1396 1439 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1397 1440 self.__dataReady = True
1398 1441
1399 1442 return avgdata_spc, avgdata_cspc, avgdata_dc
1400 1443
1401 1444 def byTime(self, datatime, *args):
1402 1445
1403 1446 self.__dataReady = False
1404 1447 avgdata_spc = None
1405 1448 avgdata_cspc = None
1406 1449 avgdata_dc = None
1407 1450 n = None
1408 1451
1409 1452 self.putData(*args)
1410 1453
1411 1454 if (datatime - self.__initime) >= self.__integrationtime:
1412 1455 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1413 1456 self.n = n
1414 1457 self.__dataReady = True
1415 1458
1416 1459 return avgdata_spc, avgdata_cspc, avgdata_dc
1417 1460
1418 1461 def integrate(self, datatime, *args):
1419 1462
1420 1463 if self.__initime == None:
1421 1464 self.__initime = datatime
1422 1465
1423 1466 if self.__byTime:
1424 1467 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
1425 1468 else:
1426 1469 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1427 1470
1428 1471 self.__lastdatatime = datatime
1429 1472
1430 1473 if avgdata_spc == None:
1431 1474 return None, None, None, None
1432 1475
1433 1476 avgdatatime = self.__initime
1434 1477 try:
1435 1478 self.__timeInterval = (self.__lastdatatime - self.__initime)/(self.n - 1)
1436 1479 except:
1437 1480 self.__timeInterval = self.__lastdatatime - self.__initime
1438 1481
1439 1482 deltatime = datatime -self.__lastdatatime
1440 1483
1441 1484 if not self.__withOverapping:
1442 1485 self.__initime = datatime
1443 1486 else:
1444 1487 self.__initime += deltatime
1445 1488
1446 1489 return avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc
1447 1490
1448 1491 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1449 1492
1450 1493 if n==1:
1451 1494 dataOut.flagNoData = False
1452 1495 return
1453 1496
1454 1497 if not self.__isConfig:
1455 1498 self.setup(n, timeInterval, overlapping)
1456 1499 self.__isConfig = True
1457 1500
1458 1501 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1459 1502 dataOut.data_spc,
1460 1503 dataOut.data_cspc,
1461 1504 dataOut.data_dc)
1462 1505
1463 1506 # dataOut.timeInterval *= n
1464 1507 dataOut.flagNoData = True
1465 1508
1466 1509 if self.__dataReady:
1467 1510
1468 1511 dataOut.data_spc = avgdata_spc
1469 1512 dataOut.data_cspc = avgdata_cspc
1470 1513 dataOut.data_dc = avgdata_dc
1471 1514
1472 1515 dataOut.nIncohInt *= self.n
1473 1516 dataOut.utctime = avgdatatime
1474 1517 #dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt * dataOut.nIncohInt * dataOut.nFFTPoints
1475 1518 dataOut.timeInterval = self.__timeInterval*self.n
1476 1519 dataOut.flagNoData = False
1477 1520
1478 1521 class ProfileConcat(Operation):
1479 1522
1480 1523 __isConfig = False
1481 1524 buffer = None
1482 1525
1483 1526 def __init__(self):
1484 1527
1485 1528 self.profileIndex = 0
1486 1529
1487 1530 def reset(self):
1488 1531 self.buffer = numpy.zeros_like(self.buffer)
1489 1532 self.start_index = 0
1490 1533 self.times = 1
1491 1534
1492 1535 def setup(self, data, m, n=1):
1493 1536 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
1494 1537 self.profiles = data.shape[1]
1495 1538 self.start_index = 0
1496 1539 self.times = 1
1497 1540
1498 1541 def concat(self, data):
1499 1542
1500 1543 self.buffer[:,self.start_index:self.profiles*self.times] = data.copy()
1501 1544 self.start_index = self.start_index + self.profiles
1502 1545
1503 1546 def run(self, dataOut, m):
1504 1547
1505 1548 dataOut.flagNoData = True
1506 1549
1507 1550 if not self.__isConfig:
1508 1551 self.setup(dataOut.data, m, 1)
1509 1552 self.__isConfig = True
1510 1553
1511 1554 self.concat(dataOut.data)
1512 1555 self.times += 1
1513 1556 if self.times > m:
1514 1557 dataOut.data = self.buffer
1515 1558 self.reset()
1516 1559 dataOut.flagNoData = False
1517 1560 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
1518 1561 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1519 1562 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * 5
1520 1563 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
1521 1564
1522 1565
1523 1566
1524 1567 class ProfileSelector(Operation):
1525 1568
1526 1569 profileIndex = None
1527 1570 # Tamanho total de los perfiles
1528 1571 nProfiles = None
1529 1572
1530 1573 def __init__(self):
1531 1574
1532 1575 self.profileIndex = 0
1533 1576
1534 1577 def incIndex(self):
1535 1578 self.profileIndex += 1
1536 1579
1537 1580 if self.profileIndex >= self.nProfiles:
1538 1581 self.profileIndex = 0
1539 1582
1540 1583 def isProfileInRange(self, minIndex, maxIndex):
1541 1584
1542 1585 if self.profileIndex < minIndex:
1543 1586 return False
1544 1587
1545 1588 if self.profileIndex > maxIndex:
1546 1589 return False
1547 1590
1548 1591 return True
1549 1592
1550 1593 def isProfileInList(self, profileList):
1551 1594
1552 1595 if self.profileIndex not in profileList:
1553 1596 return False
1554 1597
1555 1598 return True
1556 1599
1557 1600 def run(self, dataOut, profileList=None, profileRangeList=None):
1558 1601
1559 1602 dataOut.flagNoData = True
1560 1603 self.nProfiles = dataOut.nProfiles
1561 1604
1562 1605 if profileList != None:
1563 1606 if self.isProfileInList(profileList):
1564 1607 dataOut.flagNoData = False
1565 1608
1566 1609 self.incIndex()
1567 1610 return 1
1568 1611
1569 1612
1570 1613 elif profileRangeList != None:
1571 1614 minIndex = profileRangeList[0]
1572 1615 maxIndex = profileRangeList[1]
1573 1616 if self.isProfileInRange(minIndex, maxIndex):
1574 1617 dataOut.flagNoData = False
1575 1618
1576 1619 self.incIndex()
1577 1620 return 1
1578 1621
1579 1622 else:
1580 1623 raise ValueError, "ProfileSelector needs profileList or profileRangeList"
1581 1624
1582 1625 return 0
1583 1626
1584 1627 class SpectraHeisProc(ProcessingUnit):
1585 1628 def __init__(self):
1586 1629 self.objectDict = {}
1587 1630 # self.buffer = None
1588 1631 # self.firstdatatime = None
1589 1632 # self.profIndex = 0
1590 1633 self.dataOut = SpectraHeis()
1591 1634
1592 1635 def __updateObjFromInput(self):
1593 1636 self.dataOut.timeZone = self.dataIn.timeZone
1594 1637 self.dataOut.dstFlag = self.dataIn.dstFlag
1595 1638 self.dataOut.errorCount = self.dataIn.errorCount
1596 1639 self.dataOut.useLocalTime = self.dataIn.useLocalTime
1597 1640
1598 1641 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()#
1599 1642 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()#
1600 1643 self.dataOut.channelList = self.dataIn.channelList
1601 1644 self.dataOut.heightList = self.dataIn.heightList
1602 1645 # self.dataOut.dtype = self.dataIn.dtype
1603 1646 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
1604 1647 # self.dataOut.nHeights = self.dataIn.nHeights
1605 1648 # self.dataOut.nChannels = self.dataIn.nChannels
1606 1649 self.dataOut.nBaud = self.dataIn.nBaud
1607 1650 self.dataOut.nCode = self.dataIn.nCode
1608 1651 self.dataOut.code = self.dataIn.code
1609 1652 # self.dataOut.nProfiles = 1
1610 1653 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
1611 1654 self.dataOut.nFFTPoints = self.dataIn.nHeights
1612 1655 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
1613 1656 # self.dataOut.flagNoData = self.dataIn.flagNoData
1614 1657 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
1615 1658 self.dataOut.utctime = self.dataIn.utctime
1616 1659 # self.dataOut.utctime = self.firstdatatime
1617 1660 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
1618 1661 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
1619 1662 # self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
1620 1663 self.dataOut.nCohInt = self.dataIn.nCohInt
1621 1664 self.dataOut.nIncohInt = 1
1622 1665 self.dataOut.ippSeconds= self.dataIn.ippSeconds
1623 1666 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
1624 1667
1625 1668 self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nIncohInt
1626 1669 # self.dataOut.set=self.dataIn.set
1627 1670 # self.dataOut.deltaHeight=self.dataIn.deltaHeight
1628 1671
1629 1672
1630 1673 def __updateObjFromFits(self):
1631 1674 self.dataOut.utctime = self.dataIn.utctime
1632 1675 self.dataOut.channelIndexList = self.dataIn.channelIndexList
1633 1676
1634 1677 self.dataOut.channelList = self.dataIn.channelList
1635 1678 self.dataOut.heightList = self.dataIn.heightList
1636 1679 self.dataOut.data_spc = self.dataIn.data
1637 1680 self.dataOut.timeInterval = self.dataIn.timeInterval
1638 1681 self.dataOut.timeZone = self.dataIn.timeZone
1639 1682 self.dataOut.useLocalTime = True
1640 1683 # self.dataOut.
1641 1684 # self.dataOut.
1642 1685
1643 1686 def __getFft(self):
1644 1687
1645 1688 fft_volt = numpy.fft.fft(self.dataIn.data, axis=1)
1646 1689 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
1647 1690 spc = numpy.abs(fft_volt * numpy.conjugate(fft_volt))/(self.dataOut.nFFTPoints)
1648 1691 self.dataOut.data_spc = spc
1649 1692
1650 1693 def init(self):
1651 1694
1652 1695 self.dataOut.flagNoData = True
1653 1696
1654 1697 if self.dataIn.type == "Fits":
1655 1698 self.__updateObjFromFits()
1656 1699 self.dataOut.flagNoData = False
1657 1700 return
1658 1701
1659 1702 if self.dataIn.type == "SpectraHeis":
1660 1703 self.dataOut.copy(self.dataIn)
1661 1704 return
1662 1705
1663 1706 if self.dataIn.type == "Voltage":
1664 1707 self.__updateObjFromInput()
1665 1708 self.__getFft()
1666 1709 self.dataOut.flagNoData = False
1667 1710
1668 1711 return
1669 1712
1670 1713 raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
1671 1714
1672 1715
1673 1716 def selectChannels(self, channelList):
1674 1717
1675 1718 channelIndexList = []
1676 1719
1677 1720 for channel in channelList:
1678 1721 index = self.dataOut.channelList.index(channel)
1679 1722 channelIndexList.append(index)
1680 1723
1681 1724 self.selectChannelsByIndex(channelIndexList)
1682 1725
1683 1726 def selectChannelsByIndex(self, channelIndexList):
1684 1727 """
1685 1728 Selecciona un bloque de datos en base a canales segun el channelIndexList
1686 1729
1687 1730 Input:
1688 1731 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
1689 1732
1690 1733 Affected:
1691 1734 self.dataOut.data
1692 1735 self.dataOut.channelIndexList
1693 1736 self.dataOut.nChannels
1694 1737 self.dataOut.m_ProcessingHeader.totalSpectra
1695 1738 self.dataOut.systemHeaderObj.numChannels
1696 1739 self.dataOut.m_ProcessingHeader.blockSize
1697 1740
1698 1741 Return:
1699 1742 None
1700 1743 """
1701 1744
1702 1745 for channelIndex in channelIndexList:
1703 1746 if channelIndex not in self.dataOut.channelIndexList:
1704 1747 print channelIndexList
1705 1748 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
1706 1749
1707 1750 nChannels = len(channelIndexList)
1708 1751
1709 1752 data_spc = self.dataOut.data_spc[channelIndexList,:]
1710 1753
1711 1754 self.dataOut.data_spc = data_spc
1712 1755 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
1713 1756
1714 1757 return 1
1715 1758
1716 1759 class IncohInt4SpectraHeis(Operation):
1717 1760
1718 1761 __isConfig = False
1719 1762
1720 1763 __profIndex = 0
1721 1764 __withOverapping = False
1722 1765
1723 1766 __byTime = False
1724 1767 __initime = None
1725 1768 __lastdatatime = None
1726 1769 __integrationtime = None
1727 1770
1728 1771 __buffer = None
1729 1772
1730 1773 __dataReady = False
1731 1774
1732 1775 n = None
1733 1776
1734 1777
1735 1778 def __init__(self):
1736 1779
1737 1780 self.__isConfig = False
1738 1781
1739 1782 def setup(self, n=None, timeInterval=None, overlapping=False):
1740 1783 """
1741 1784 Set the parameters of the integration class.
1742 1785
1743 1786 Inputs:
1744 1787
1745 1788 n : Number of coherent integrations
1746 1789 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1747 1790 overlapping :
1748 1791
1749 1792 """
1750 1793
1751 1794 self.__initime = None
1752 1795 self.__lastdatatime = 0
1753 1796 self.__buffer = None
1754 1797 self.__dataReady = False
1755 1798
1756 1799
1757 1800 if n == None and timeInterval == None:
1758 1801 raise ValueError, "n or timeInterval should be specified ..."
1759 1802
1760 1803 if n != None:
1761 1804 self.n = n
1762 1805 self.__byTime = False
1763 1806 else:
1764 1807 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
1765 1808 self.n = 9999
1766 1809 self.__byTime = True
1767 1810
1768 1811 if overlapping:
1769 1812 self.__withOverapping = True
1770 1813 self.__buffer = None
1771 1814 else:
1772 1815 self.__withOverapping = False
1773 1816 self.__buffer = 0
1774 1817
1775 1818 self.__profIndex = 0
1776 1819
1777 1820 def putData(self, data):
1778 1821
1779 1822 """
1780 1823 Add a profile to the __buffer and increase in one the __profileIndex
1781 1824
1782 1825 """
1783 1826
1784 1827 if not self.__withOverapping:
1785 1828 self.__buffer += data.copy()
1786 1829 self.__profIndex += 1
1787 1830 return
1788 1831
1789 1832 #Overlapping data
1790 1833 nChannels, nHeis = data.shape
1791 1834 data = numpy.reshape(data, (1, nChannels, nHeis))
1792 1835
1793 1836 #If the buffer is empty then it takes the data value
1794 1837 if self.__buffer == None:
1795 1838 self.__buffer = data
1796 1839 self.__profIndex += 1
1797 1840 return
1798 1841
1799 1842 #If the buffer length is lower than n then stakcing the data value
1800 1843 if self.__profIndex < self.n:
1801 1844 self.__buffer = numpy.vstack((self.__buffer, data))
1802 1845 self.__profIndex += 1
1803 1846 return
1804 1847
1805 1848 #If the buffer length is equal to n then replacing the last buffer value with the data value
1806 1849 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
1807 1850 self.__buffer[self.n-1] = data
1808 1851 self.__profIndex = self.n
1809 1852 return
1810 1853
1811 1854
1812 1855 def pushData(self):
1813 1856 """
1814 1857 Return the sum of the last profiles and the profiles used in the sum.
1815 1858
1816 1859 Affected:
1817 1860
1818 1861 self.__profileIndex
1819 1862
1820 1863 """
1821 1864
1822 1865 if not self.__withOverapping:
1823 1866 data = self.__buffer
1824 1867 n = self.__profIndex
1825 1868
1826 1869 self.__buffer = 0
1827 1870 self.__profIndex = 0
1828 1871
1829 1872 return data, n
1830 1873
1831 1874 #Integration with Overlapping
1832 1875 data = numpy.sum(self.__buffer, axis=0)
1833 1876 n = self.__profIndex
1834 1877
1835 1878 return data, n
1836 1879
1837 1880 def byProfiles(self, data):
1838 1881
1839 1882 self.__dataReady = False
1840 1883 avgdata = None
1841 1884 n = None
1842 1885
1843 1886 self.putData(data)
1844 1887
1845 1888 if self.__profIndex == self.n:
1846 1889
1847 1890 avgdata, n = self.pushData()
1848 1891 self.__dataReady = True
1849 1892
1850 1893 return avgdata
1851 1894
1852 1895 def byTime(self, data, datatime):
1853 1896
1854 1897 self.__dataReady = False
1855 1898 avgdata = None
1856 1899 n = None
1857 1900
1858 1901 self.putData(data)
1859 1902
1860 1903 if (datatime - self.__initime) >= self.__integrationtime:
1861 1904 avgdata, n = self.pushData()
1862 1905 self.n = n
1863 1906 self.__dataReady = True
1864 1907
1865 1908 return avgdata
1866 1909
1867 1910 def integrate(self, data, datatime=None):
1868 1911
1869 1912 if self.__initime == None:
1870 1913 self.__initime = datatime
1871 1914
1872 1915 if self.__byTime:
1873 1916 avgdata = self.byTime(data, datatime)
1874 1917 else:
1875 1918 avgdata = self.byProfiles(data)
1876 1919
1877 1920
1878 1921 self.__lastdatatime = datatime
1879 1922
1880 1923 if avgdata == None:
1881 1924 return None, None
1882 1925
1883 1926 avgdatatime = self.__initime
1884 1927
1885 1928 deltatime = datatime -self.__lastdatatime
1886 1929
1887 1930 if not self.__withOverapping:
1888 1931 self.__initime = datatime
1889 1932 else:
1890 1933 self.__initime += deltatime
1891 1934
1892 1935 return avgdata, avgdatatime
1893 1936
1894 1937 def run(self, dataOut, **kwargs):
1895 1938
1896 1939 if not self.__isConfig:
1897 1940 self.setup(**kwargs)
1898 1941 self.__isConfig = True
1899 1942
1900 1943 avgdata, avgdatatime = self.integrate(dataOut.data_spc, dataOut.utctime)
1901 1944
1902 1945 # dataOut.timeInterval *= n
1903 1946 dataOut.flagNoData = True
1904 1947
1905 1948 if self.__dataReady:
1906 1949 dataOut.data_spc = avgdata
1907 1950 dataOut.nIncohInt *= self.n
1908 1951 # dataOut.nCohInt *= self.n
1909 1952 dataOut.utctime = avgdatatime
1910 1953 dataOut.timeInterval = dataOut.ippSeconds * dataOut.nIncohInt
1911 1954 # dataOut.timeInterval = self.__timeInterval*self.n
1912 1955 dataOut.flagNoData = False
1913 1956
1914 1957
1915 1958
1916 1959
1917 1960 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now