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