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