##// END OF EJS Templates
schain-jc
RSS
7 de marzo 2018
ebocanegra -
r1148:fee9f80824df
parent child
Show More
Side by Side Unified
Context file:
r1148:fee9f80824df Loading diff...:

The requested changes are too big and content was truncated. Show full diff

Show file before | Show file after | Hide comments
@@ -1,1231 +1,1229
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12 from schainpy import cSchain
13 13
14 14
15 15 def getNumpyDtype(dataTypeCode):
16 16
17 17 if dataTypeCode == 0:
18 18 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
19 19 elif dataTypeCode == 1:
20 20 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
21 21 elif dataTypeCode == 2:
22 22 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
23 23 elif dataTypeCode == 3:
24 24 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
25 25 elif dataTypeCode == 4:
26 26 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
27 27 elif dataTypeCode == 5:
28 28 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
29 29 else:
30 30 raise ValueError, 'dataTypeCode was not defined'
31 31
32 32 return numpyDtype
33 33
34 34 def getDataTypeCode(numpyDtype):
35 35
36 36 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
37 37 datatype = 0
38 38 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
39 39 datatype = 1
40 40 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
41 41 datatype = 2
42 42 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
43 43 datatype = 3
44 44 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
45 45 datatype = 4
46 46 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
47 47 datatype = 5
48 48 else:
49 49 datatype = None
50 50
51 51 return datatype
52 52
53 53 def hildebrand_sekhon(data, navg):
54 54 """
55 55 This method is for the objective determination of the noise level in Doppler spectra. This
56 56 implementation technique is based on the fact that the standard deviation of the spectral
57 57 densities is equal to the mean spectral density for white Gaussian noise
58 58
59 59 Inputs:
60 60 Data : heights
61 61 navg : numbers of averages
62 62
63 63 Return:
64 64 -1 : any error
65 65 anoise : noise's level
66 66 """
67 67
68 68 sortdata = numpy.sort(data, axis=None)
69 69 # lenOfData = len(sortdata)
70 70 # nums_min = lenOfData*0.2
71 71 #
72 72 # if nums_min <= 5:
73 73 # nums_min = 5
74 74 #
75 75 # sump = 0.
76 76 #
77 77 # sumq = 0.
78 78 #
79 79 # j = 0
80 80 #
81 81 # cont = 1
82 82 #
83 83 # while((cont==1)and(j<lenOfData)):
84 84 #
85 85 # sump += sortdata[j]
86 86 #
87 87 # sumq += sortdata[j]**2
88 88 #
89 89 # if j > nums_min:
90 90 # rtest = float(j)/(j-1) + 1.0/navg
91 91 # if ((sumq*j) > (rtest*sump**2)):
92 92 # j = j - 1
93 93 # sump = sump - sortdata[j]
94 94 # sumq = sumq - sortdata[j]**2
95 95 # cont = 0
96 96 #
97 97 # j += 1
98 98 #
99 99 # lnoise = sump /j
100 100 #
101 101 # return lnoise
102 102
103 103 return cSchain.hildebrand_sekhon(sortdata, navg)
104 104
105 105
106 106 class Beam:
107 107
108 108 def __init__(self):
109 109 self.codeList = []
110 110 self.azimuthList = []
111 111 self.zenithList = []
112 112
113 113 class GenericData(object):
114 114
115 115 flagNoData = True
116 116
117 117 def copy(self, inputObj=None):
118 118
119 119 if inputObj == None:
120 120 return copy.deepcopy(self)
121 121
122 122 for key in inputObj.__dict__.keys():
123 123
124 124 attribute = inputObj.__dict__[key]
125 125
126 126 #If this attribute is a tuple or list
127 127 if type(inputObj.__dict__[key]) in (tuple, list):
128 128 self.__dict__[key] = attribute[:]
129 129 continue
130 130
131 131 #If this attribute is another object or instance
132 132 if hasattr(attribute, '__dict__'):
133 133 self.__dict__[key] = attribute.copy()
134 134 continue
135 135
136 136 self.__dict__[key] = inputObj.__dict__[key]
137 137
138 138 def deepcopy(self):
139 139
140 140 return copy.deepcopy(self)
141 141
142 142 def isEmpty(self):
143 143
144 144 return self.flagNoData
145 145
146 146 class JROData(GenericData):
147 147
148 148 # m_BasicHeader = BasicHeader()
149 149 # m_ProcessingHeader = ProcessingHeader()
150 150
151 151 systemHeaderObj = SystemHeader()
152 152
153 153 radarControllerHeaderObj = RadarControllerHeader()
154 154
155 155 # data = None
156 156
157 157 type = None
158 158
159 159 datatype = None #dtype but in string
160 160
161 161 # dtype = None
162 162
163 163 # nChannels = None
164 164
165 165 # nHeights = None
166 166
167 167 nProfiles = None
168 168
169 169 heightList = None
170 170
171 171 channelList = None
172 172
173 173 flagDiscontinuousBlock = False
174 174
175 175 useLocalTime = False
176 176
177 177 utctime = None
178 178
179 179 timeZone = None
180 180
181 181 dstFlag = None
182 182
183 183 errorCount = None
184 184
185 185 blocksize = None
186 186
187 187 # nCode = None
188 188 #
189 189 # nBaud = None
190 190 #
191 191 # code = None
192 192
193 193 flagDecodeData = False #asumo q la data no esta decodificada
194 194
195 195 flagDeflipData = False #asumo q la data no esta sin flip
196 196
197 197 flagShiftFFT = False
198 198
199 199 # ippSeconds = None
200 200
201 201 # timeInterval = None
202 202
203 203 nCohInt = None
204 204
205 205 # noise = None
206 206
207 207 windowOfFilter = 1
208 208
209 209 #Speed of ligth
210 210 C = 3e8
211 211
212 212 frequency = 49.92e6
213 213
214 214 realtime = False
215 215
216 216 beacon_heiIndexList = None
217 217
218 218 last_block = None
219 219
220 220 blocknow = None
221 221
222 222 azimuth = None
223 223
224 224 zenith = None
225 225
226 226 beam = Beam()
227 227
228 228 profileIndex = None
229 229
230 230 def getNoise(self):
231 231
232 232 raise NotImplementedError
233 233
234 234 def getNChannels(self):
235 235
236 236 return len(self.channelList)
237 237
238 238 def getChannelIndexList(self):
239 239
240 240 return range(self.nChannels)
241 241
242 242 def getNHeights(self):
243 243
244 244 return len(self.heightList)
245 245
246 246 def getHeiRange(self, extrapoints=0):
247 247
248 248 heis = self.heightList
249 249 # deltah = self.heightList[1] - self.heightList[0]
250 250 #
251 251 # heis.append(self.heightList[-1])
252 252
253 253 return heis
254 254
255 255 def getDeltaH(self):
256 256
257 257 delta = self.heightList[1] - self.heightList[0]
258 258
259 259 return delta
260 260
261 261 def getltctime(self):
262 262
263 263 if self.useLocalTime:
264 264 return self.utctime - self.timeZone*60
265 265
266 266 return self.utctime
267 267
268 268 def getDatatime(self):
269 269
270 270 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
271 271 return datatimeValue
272 272
273 273 def getTimeRange(self):
274 274
275 275 datatime = []
276 276
277 277 datatime.append(self.ltctime)
278 278 datatime.append(self.ltctime + self.timeInterval+1)
279 279
280 280 datatime = numpy.array(datatime)
281 281
282 282 return datatime
283 283
284 284 def getFmaxTimeResponse(self):
285 285
286 286 period = (10**-6)*self.getDeltaH()/(0.15)
287 287
288 288 PRF = 1./(period * self.nCohInt)
289 289
290 290 fmax = PRF
291 291
292 292 return fmax
293 293
294 294 def getFmax(self):
295 295
296 296 PRF = 1./(self.ippSeconds * self.nCohInt)
297 297
298 298 fmax = PRF
299 299
300 300 return fmax
301 301
302 302 def getVmax(self):
303 303
304 304 _lambda = self.C/self.frequency
305 305
306 306 vmax = self.getFmax() * _lambda/2
307 307
308 308 return vmax
309 309
310 310 def get_ippSeconds(self):
311 311 '''
312 312 '''
313 313 return self.radarControllerHeaderObj.ippSeconds
314 314
315 315 def set_ippSeconds(self, ippSeconds):
316 316 '''
317 317 '''
318 318
319 319 self.radarControllerHeaderObj.ippSeconds = ippSeconds
320 320
321 321 return
322 322
323 323 def get_dtype(self):
324 324 '''
325 325 '''
326 326 return getNumpyDtype(self.datatype)
327 327
328 328 def set_dtype(self, numpyDtype):
329 329 '''
330 330 '''
331 331
332 332 self.datatype = getDataTypeCode(numpyDtype)
333 333
334 334 def get_code(self):
335 335 '''
336 336 '''
337 337 return self.radarControllerHeaderObj.code
338 338
339 339 def set_code(self, code):
340 340 '''
341 341 '''
342 342 self.radarControllerHeaderObj.code = code
343 343
344 344 return
345 345
346 346 def get_ncode(self):
347 347 '''
348 348 '''
349 349 return self.radarControllerHeaderObj.nCode
350 350
351 351 def set_ncode(self, nCode):
352 352 '''
353 353 '''
354 354 self.radarControllerHeaderObj.nCode = nCode
355 355
356 356 return
357 357
358 358 def get_nbaud(self):
359 359 '''
360 360 '''
361 361 return self.radarControllerHeaderObj.nBaud
362 362
363 363 def set_nbaud(self, nBaud):
364 364 '''
365 365 '''
366 366 self.radarControllerHeaderObj.nBaud = nBaud
367 367
368 368 return
369 369
370 370 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
371 371 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
372 372 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
373 373 #noise = property(getNoise, "I'm the 'nHeights' property.")
374 374 datatime = property(getDatatime, "I'm the 'datatime' property")
375 375 ltctime = property(getltctime, "I'm the 'ltctime' property")
376 376 ippSeconds = property(get_ippSeconds, set_ippSeconds)
377 377 dtype = property(get_dtype, set_dtype)
378 378 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
379 379 code = property(get_code, set_code)
380 380 nCode = property(get_ncode, set_ncode)
381 381 nBaud = property(get_nbaud, set_nbaud)
382 382
383 383 class Voltage(JROData):
384 384
385 385 #data es un numpy array de 2 dmensiones (canales, alturas)
386 386 data = None
387 387
388 388 def __init__(self):
389 389 '''
390 390 Constructor
391 391 '''
392 392
393 393 self.useLocalTime = True
394 394
395 395 self.radarControllerHeaderObj = RadarControllerHeader()
396 396
397 397 self.systemHeaderObj = SystemHeader()
398 398
399 399 self.type = "Voltage"
400 400
401 401 self.data = None
402 402
403 403 # self.dtype = None
404 404
405 405 # self.nChannels = 0
406 406
407 407 # self.nHeights = 0
408 408
409 409 self.nProfiles = None
410 410
411 411 self.heightList = None
412 412
413 413 self.channelList = None
414 414
415 415 # self.channelIndexList = None
416 416
417 417 self.flagNoData = True
418 418
419 419 self.flagDiscontinuousBlock = False
420 420
421 421 self.utctime = None
422 422
423 423 self.timeZone = None
424 424
425 425 self.dstFlag = None
426 426
427 427 self.errorCount = None
428 428
429 429 self.nCohInt = None
430 430
431 431 self.blocksize = None
432 432
433 433 self.flagDecodeData = False #asumo q la data no esta decodificada
434 434
435 435 self.flagDeflipData = False #asumo q la data no esta sin flip
436 436
437 437 self.flagShiftFFT = False
438 438
439 439 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
440 440
441 441 self.profileIndex = 0
442 442
443 443 def getNoisebyHildebrand(self, channel = None):
444 444 """
445 445 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
446 446
447 447 Return:
448 448 noiselevel
449 449 """
450 450
451 451 if channel != None:
452 452 data = self.data[channel]
453 453 nChannels = 1
454 454 else:
455 455 data = self.data
456 456 nChannels = self.nChannels
457 457
458 458 noise = numpy.zeros(nChannels)
459 459 power = data * numpy.conjugate(data)
460 460
461 461 for thisChannel in range(nChannels):
462 462 if nChannels == 1:
463 463 daux = power[:].real
464 464 else:
465 465 daux = power[thisChannel,:].real
466 466 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
467 467
468 468 return noise
469 469
470 470 def getNoise(self, type = 1, channel = None):
471 471
472 472 if type == 1:
473 473 noise = self.getNoisebyHildebrand(channel)
474 474
475 475 return noise
476 476
477 477 def getPower(self, channel = None):
478 478
479 479 if channel != None:
480 480 data = self.data[channel]
481 481 else:
482 482 data = self.data
483 483
484 484 power = data * numpy.conjugate(data)
485 485 powerdB = 10*numpy.log10(power.real)
486 486 powerdB = numpy.squeeze(powerdB)
487 487
488 488 return powerdB
489 489
490 490 def getTimeInterval(self):
491 491
492 492 timeInterval = self.ippSeconds * self.nCohInt
493 493
494 494 return timeInterval
495 495
496 496 noise = property(getNoise, "I'm the 'nHeights' property.")
497 497 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
498 498
499 499 class Spectra(JROData):
500 500
501 501 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
502 502 data_spc = None
503 503
504 504 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
505 505 data_cspc = None
506 506
507 507 #data dc es un numpy array de 2 dmensiones (canales, alturas)
508 508 data_dc = None
509 509
510 510 #data power
511 511 data_pwr = None
512 512
513 513 nFFTPoints = None
514 514
515 515 # nPairs = None
516 516
517 517 pairsList = None
518 518
519 519 nIncohInt = None
520 520
521 521 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
522 522
523 523 nCohInt = None #se requiere para determinar el valor de timeInterval
524 524
525 525 ippFactor = None
526 526
527 527 profileIndex = 0
528 528
529 529 plotting = "spectra"
530 530
531 531 def __init__(self):
532 532 '''
533 533 Constructor
534 534 '''
535 535
536 536 self.useLocalTime = True
537 537
538 538 self.radarControllerHeaderObj = RadarControllerHeader()
539 539
540 540 self.systemHeaderObj = SystemHeader()
541 541
542 542 self.type = "Spectra"
543 543
544 544 # self.data = None
545 545
546 546 # self.dtype = None
547 547
548 548 # self.nChannels = 0
549 549
550 550 # self.nHeights = 0
551 551
552 552 self.nProfiles = None
553 553
554 554 self.heightList = None
555 555
556 556 self.channelList = None
557 557
558 558 # self.channelIndexList = None
559 559
560 560 self.pairsList = None
561 561
562 562 self.flagNoData = True
563 563
564 564 self.flagDiscontinuousBlock = False
565 565
566 566 self.utctime = None
567 567
568 568 self.nCohInt = None
569 569
570 570 self.nIncohInt = None
571 571
572 572 self.blocksize = None
573 573
574 574 self.nFFTPoints = None
575 575
576 576 self.wavelength = None
577 577
578 578 self.flagDecodeData = False #asumo q la data no esta decodificada
579 579
580 580 self.flagDeflipData = False #asumo q la data no esta sin flip
581 581
582 582 self.flagShiftFFT = False
583 583
584 584 self.ippFactor = 1
585 585
586 586 #self.noise = None
587 587
588 588 self.beacon_heiIndexList = []
589 589
590 590 self.noise_estimation = None
591 591
592 592
593 593 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
594 594 """
595 595 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
596 596
597 597 Return:
598 598 noiselevel
599 599 """
600 600
601 601 noise = numpy.zeros(self.nChannels)
602 602
603 603 for channel in range(self.nChannels):
604 604 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
605 605 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
606 606
607 607 return noise
608 608
609 609 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
610 610
611 611 if self.noise_estimation is not None:
612 612 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
613 613 else:
614 614 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
615 615 return noise
616 616
617 617 def getFreqRangeTimeResponse(self, extrapoints=0):
618 618
619 619 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
620 620 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
621 621
622 622 return freqrange
623 623
624 624 def getAcfRange(self, extrapoints=0):
625 625
626 626 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
627 627 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
628 628
629 629 return freqrange
630 630
631 631 def getFreqRange(self, extrapoints=0):
632 632
633 633 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
634 634 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
635 635
636 636 return freqrange
637 637
638 638 def getVelRange(self, extrapoints=0):
639 639
640 print 'VELMAX', self.getVmax()
641 asdasdasd
642 640 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
643 641 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
644 642
645 643 return velrange
646 644
647 645 def getNPairs(self):
648 646
649 647 return len(self.pairsList)
650 648
651 649 def getPairsIndexList(self):
652 650
653 651 return range(self.nPairs)
654 652
655 653 def getNormFactor(self):
656 654
657 655 pwcode = 1
658 656
659 657 if self.flagDecodeData:
660 658 pwcode = numpy.sum(self.code[0]**2)
661 659 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
662 660 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
663 661
664 662 return normFactor
665 663
666 664 def getFlagCspc(self):
667 665
668 666 if self.data_cspc is None:
669 667 return True
670 668
671 669 return False
672 670
673 671 def getFlagDc(self):
674 672
675 673 if self.data_dc is None:
676 674 return True
677 675
678 676 return False
679 677
680 678 def getTimeInterval(self):
681 679
682 680 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
683 681
684 682 return timeInterval
685 683
686 684 def getPower(self):
687 685
688 686 factor = self.normFactor
689 687 z = self.data_spc/factor
690 688 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
691 689 avg = numpy.average(z, axis=1)
692 690
693 691 return 10*numpy.log10(avg)
694 692
695 693 def getCoherence(self, pairsList=None, phase=False):
696 694
697 695 z = []
698 696 if pairsList is None:
699 697 pairsIndexList = self.pairsIndexList
700 698 else:
701 699 pairsIndexList = []
702 700 for pair in pairsList:
703 701 if pair not in self.pairsList:
704 702 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
705 703 pairsIndexList.append(self.pairsList.index(pair))
706 704 for i in range(len(pairsIndexList)):
707 705 pair = self.pairsList[pairsIndexList[i]]
708 706 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
709 707 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
710 708 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
711 709 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
712 710 if phase:
713 711 data = numpy.arctan2(avgcoherenceComplex.imag,
714 712 avgcoherenceComplex.real)*180/numpy.pi
715 713 else:
716 714 data = numpy.abs(avgcoherenceComplex)
717 715
718 716 z.append(data)
719 717
720 718 return numpy.array(z)
721 719
722 720 def setValue(self, value):
723 721
724 722 print "This property should not be initialized"
725 723
726 724 return
727 725
728 726 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
729 727 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
730 728 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
731 729 flag_cspc = property(getFlagCspc, setValue)
732 730 flag_dc = property(getFlagDc, setValue)
733 731 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
734 732 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
735 733
736 734 class SpectraHeis(Spectra):
737 735
738 736 data_spc = None
739 737
740 738 data_cspc = None
741 739
742 740 data_dc = None
743 741
744 742 nFFTPoints = None
745 743
746 744 # nPairs = None
747 745
748 746 pairsList = None
749 747
750 748 nCohInt = None
751 749
752 750 nIncohInt = None
753 751
754 752 def __init__(self):
755 753
756 754 self.radarControllerHeaderObj = RadarControllerHeader()
757 755
758 756 self.systemHeaderObj = SystemHeader()
759 757
760 758 self.type = "SpectraHeis"
761 759
762 760 # self.dtype = None
763 761
764 762 # self.nChannels = 0
765 763
766 764 # self.nHeights = 0
767 765
768 766 self.nProfiles = None
769 767
770 768 self.heightList = None
771 769
772 770 self.channelList = None
773 771
774 772 # self.channelIndexList = None
775 773
776 774 self.flagNoData = True
777 775
778 776 self.flagDiscontinuousBlock = False
779 777
780 778 # self.nPairs = 0
781 779
782 780 self.utctime = None
783 781
784 782 self.blocksize = None
785 783
786 784 self.profileIndex = 0
787 785
788 786 self.nCohInt = 1
789 787
790 788 self.nIncohInt = 1
791 789
792 790 def getNormFactor(self):
793 791 pwcode = 1
794 792 if self.flagDecodeData:
795 793 pwcode = numpy.sum(self.code[0]**2)
796 794
797 795 normFactor = self.nIncohInt*self.nCohInt*pwcode
798 796
799 797 return normFactor
800 798
801 799 def getTimeInterval(self):
802 800
803 801 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
804 802
805 803 return timeInterval
806 804
807 805 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
808 806 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
809 807
810 808 class Fits(JROData):
811 809
812 810 heightList = None
813 811
814 812 channelList = None
815 813
816 814 flagNoData = True
817 815
818 816 flagDiscontinuousBlock = False
819 817
820 818 useLocalTime = False
821 819
822 820 utctime = None
823 821
824 822 timeZone = None
825 823
826 824 # ippSeconds = None
827 825
828 826 # timeInterval = None
829 827
830 828 nCohInt = None
831 829
832 830 nIncohInt = None
833 831
834 832 noise = None
835 833
836 834 windowOfFilter = 1
837 835
838 836 #Speed of ligth
839 837 C = 3e8
840 838
841 839 frequency = 49.92e6
842 840
843 841 realtime = False
844 842
845 843
846 844 def __init__(self):
847 845
848 846 self.type = "Fits"
849 847
850 848 self.nProfiles = None
851 849
852 850 self.heightList = None
853 851
854 852 self.channelList = None
855 853
856 854 # self.channelIndexList = None
857 855
858 856 self.flagNoData = True
859 857
860 858 self.utctime = None
861 859
862 860 self.nCohInt = 1
863 861
864 862 self.nIncohInt = 1
865 863
866 864 self.useLocalTime = True
867 865
868 866 self.profileIndex = 0
869 867
870 868 # self.utctime = None
871 869 # self.timeZone = None
872 870 # self.ltctime = None
873 871 # self.timeInterval = None
874 872 # self.header = None
875 873 # self.data_header = None
876 874 # self.data = None
877 875 # self.datatime = None
878 876 # self.flagNoData = False
879 877 # self.expName = ''
880 878 # self.nChannels = None
881 879 # self.nSamples = None
882 880 # self.dataBlocksPerFile = None
883 881 # self.comments = ''
884 882 #
885 883
886 884
887 885 def getltctime(self):
888 886
889 887 if self.useLocalTime:
890 888 return self.utctime - self.timeZone*60
891 889
892 890 return self.utctime
893 891
894 892 def getDatatime(self):
895 893
896 894 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
897 895 return datatime
898 896
899 897 def getTimeRange(self):
900 898
901 899 datatime = []
902 900
903 901 datatime.append(self.ltctime)
904 902 datatime.append(self.ltctime + self.timeInterval)
905 903
906 904 datatime = numpy.array(datatime)
907 905
908 906 return datatime
909 907
910 908 def getHeiRange(self):
911 909
912 910 heis = self.heightList
913 911
914 912 return heis
915 913
916 914 def getNHeights(self):
917 915
918 916 return len(self.heightList)
919 917
920 918 def getNChannels(self):
921 919
922 920 return len(self.channelList)
923 921
924 922 def getChannelIndexList(self):
925 923
926 924 return range(self.nChannels)
927 925
928 926 def getNoise(self, type = 1):
929 927
930 928 #noise = numpy.zeros(self.nChannels)
931 929
932 930 if type == 1:
933 931 noise = self.getNoisebyHildebrand()
934 932
935 933 if type == 2:
936 934 noise = self.getNoisebySort()
937 935
938 936 if type == 3:
939 937 noise = self.getNoisebyWindow()
940 938
941 939 return noise
942 940
943 941 def getTimeInterval(self):
944 942
945 943 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
946 944
947 945 return timeInterval
948 946
949 947 datatime = property(getDatatime, "I'm the 'datatime' property")
950 948 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
951 949 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
952 950 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
953 951 noise = property(getNoise, "I'm the 'nHeights' property.")
954 952
955 953 ltctime = property(getltctime, "I'm the 'ltctime' property")
956 954 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
957 955
958 956
959 957 class Correlation(JROData):
960 958
961 959 noise = None
962 960
963 961 SNR = None
964 962
965 963 #--------------------------------------------------
966 964
967 965 mode = None
968 966
969 967 split = False
970 968
971 969 data_cf = None
972 970
973 971 lags = None
974 972
975 973 lagRange = None
976 974
977 975 pairsList = None
978 976
979 977 normFactor = None
980 978
981 979 #--------------------------------------------------
982 980
983 981 # calculateVelocity = None
984 982
985 983 nLags = None
986 984
987 985 nPairs = None
988 986
989 987 nAvg = None
990 988
991 989
992 990 def __init__(self):
993 991 '''
994 992 Constructor
995 993 '''
996 994 self.radarControllerHeaderObj = RadarControllerHeader()
997 995
998 996 self.systemHeaderObj = SystemHeader()
999 997
1000 998 self.type = "Correlation"
1001 999
1002 1000 self.data = None
1003 1001
1004 1002 self.dtype = None
1005 1003
1006 1004 self.nProfiles = None
1007 1005
1008 1006 self.heightList = None
1009 1007
1010 1008 self.channelList = None
1011 1009
1012 1010 self.flagNoData = True
1013 1011
1014 1012 self.flagDiscontinuousBlock = False
1015 1013
1016 1014 self.utctime = None
1017 1015
1018 1016 self.timeZone = None
1019 1017
1020 1018 self.dstFlag = None
1021 1019
1022 1020 self.errorCount = None
1023 1021
1024 1022 self.blocksize = None
1025 1023
1026 1024 self.flagDecodeData = False #asumo q la data no esta decodificada
1027 1025
1028 1026 self.flagDeflipData = False #asumo q la data no esta sin flip
1029 1027
1030 1028 self.pairsList = None
1031 1029
1032 1030 self.nPoints = None
1033 1031
1034 1032 def getPairsList(self):
1035 1033
1036 1034 return self.pairsList
1037 1035
1038 1036 def getNoise(self, mode = 2):
1039 1037
1040 1038 indR = numpy.where(self.lagR == 0)[0][0]
1041 1039 indT = numpy.where(self.lagT == 0)[0][0]
1042 1040
1043 1041 jspectra0 = self.data_corr[:,:,indR,:]
1044 1042 jspectra = copy.copy(jspectra0)
1045 1043
1046 1044 num_chan = jspectra.shape[0]
1047 1045 num_hei = jspectra.shape[2]
1048 1046
1049 1047 freq_dc = jspectra.shape[1]/2
1050 1048 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1051 1049
1052 1050 if ind_vel[0]<0:
1053 1051 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1054 1052
1055 1053 if mode == 1:
1056 1054 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1057 1055
1058 1056 if mode == 2:
1059 1057
1060 1058 vel = numpy.array([-2,-1,1,2])
1061 1059 xx = numpy.zeros([4,4])
1062 1060
1063 1061 for fil in range(4):
1064 1062 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1065 1063
1066 1064 xx_inv = numpy.linalg.inv(xx)
1067 1065 xx_aux = xx_inv[0,:]
1068 1066
1069 1067 for ich in range(num_chan):
1070 1068 yy = jspectra[ich,ind_vel,:]
1071 1069 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1072 1070
1073 1071 junkid = jspectra[ich,freq_dc,:]<=0
1074 1072 cjunkid = sum(junkid)
1075 1073
1076 1074 if cjunkid.any():
1077 1075 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1078 1076
1079 1077 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1080 1078
1081 1079 return noise
1082 1080
1083 1081 def getTimeInterval(self):
1084 1082
1085 1083 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1086 1084
1087 1085 return timeInterval
1088 1086
1089 1087 def splitFunctions(self):
1090 1088
1091 1089 pairsList = self.pairsList
1092 1090 ccf_pairs = []
1093 1091 acf_pairs = []
1094 1092 ccf_ind = []
1095 1093 acf_ind = []
1096 1094 for l in range(len(pairsList)):
1097 1095 chan0 = pairsList[l][0]
1098 1096 chan1 = pairsList[l][1]
1099 1097
1100 1098 #Obteniendo pares de Autocorrelacion
1101 1099 if chan0 == chan1:
1102 1100 acf_pairs.append(chan0)
1103 1101 acf_ind.append(l)
1104 1102 else:
1105 1103 ccf_pairs.append(pairsList[l])
1106 1104 ccf_ind.append(l)
1107 1105
1108 1106 data_acf = self.data_cf[acf_ind]
1109 1107 data_ccf = self.data_cf[ccf_ind]
1110 1108
1111 1109 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1112 1110
1113 1111 def getNormFactor(self):
1114 1112 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1115 1113 acf_pairs = numpy.array(acf_pairs)
1116 1114 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1117 1115
1118 1116 for p in range(self.nPairs):
1119 1117 pair = self.pairsList[p]
1120 1118
1121 1119 ch0 = pair[0]
1122 1120 ch1 = pair[1]
1123 1121
1124 1122 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1125 1123 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1126 1124 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1127 1125
1128 1126 return normFactor
1129 1127
1130 1128 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1131 1129 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1132 1130
1133 1131 class Parameters(Spectra):
1134 1132
1135 1133 experimentInfo = None #Information about the experiment
1136 1134
1137 1135 #Information from previous data
1138 1136
1139 1137 inputUnit = None #Type of data to be processed
1140 1138
1141 1139 operation = None #Type of operation to parametrize
1142 1140
1143 1141 #normFactor = None #Normalization Factor
1144 1142
1145 1143 groupList = None #List of Pairs, Groups, etc
1146 1144
1147 1145 #Parameters
1148 1146
1149 1147 data_param = None #Parameters obtained
1150 1148
1151 1149 data_pre = None #Data Pre Parametrization
1152 1150
1153 1151 data_SNR = None #Signal to Noise Ratio
1154 1152
1155 1153 # heightRange = None #Heights
1156 1154
1157 1155 abscissaList = None #Abscissa, can be velocities, lags or time
1158 1156
1159 1157 # noise = None #Noise Potency
1160 1158
1161 1159 utctimeInit = None #Initial UTC time
1162 1160
1163 1161 paramInterval = None #Time interval to calculate Parameters in seconds
1164 1162
1165 1163 useLocalTime = True
1166 1164
1167 1165 #Fitting
1168 1166
1169 1167 data_error = None #Error of the estimation
1170 1168
1171 1169 constants = None
1172 1170
1173 1171 library = None
1174 1172
1175 1173 #Output signal
1176 1174
1177 1175 outputInterval = None #Time interval to calculate output signal in seconds
1178 1176
1179 1177 data_output = None #Out signal
1180 1178
1181 1179 nAvg = None
1182 1180
1183 1181 noise_estimation = None
1184 1182
1185 1183 GauSPC = None #Fit gaussian SPC
1186 1184
1187 1185
1188 1186 def __init__(self):
1189 1187 '''
1190 1188 Constructor
1191 1189 '''
1192 1190 self.radarControllerHeaderObj = RadarControllerHeader()
1193 1191
1194 1192 self.systemHeaderObj = SystemHeader()
1195 1193
1196 1194 self.type = "Parameters"
1197 1195
1198 1196 def getTimeRange1(self, interval):
1199 1197
1200 1198 datatime = []
1201 1199
1202 1200 if self.useLocalTime:
1203 1201 time1 = self.utctimeInit - self.timeZone*60
1204 1202 else:
1205 1203 time1 = self.utctimeInit
1206 1204
1207 1205 datatime.append(time1)
1208 1206 datatime.append(time1 + interval)
1209 1207 datatime = numpy.array(datatime)
1210 1208
1211 1209 return datatime
1212 1210
1213 1211 def getTimeInterval(self):
1214 1212
1215 1213 if hasattr(self, 'timeInterval1'):
1216 1214 return self.timeInterval1
1217 1215 else:
1218 1216 return self.paramInterval
1219 1217
1220 1218 def setValue(self, value):
1221 1219
1222 1220 print "This property should not be initialized"
1223 1221
1224 1222 return
1225 1223
1226 1224 def getNoise(self):
1227 1225
1228 1226 return self.spc_noise
1229 1227
1230 1228 timeInterval = property(getTimeInterval)
1231 1229 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
Show file before | Show file after | Hide comments
@@ -1,2160 +1,2160
1 1 import os
2 2 import datetime
3 3 import numpy
4 4 import inspect
5 5 from figure import Figure, isRealtime, isTimeInHourRange
6 6 from plotting_codes import *
7 7
8 8
9 class FitGauPlot(Figure):
9 class SpcParamPlot(Figure):
10 10
11 11 isConfig = None
12 12 __nsubplots = None
13 13
14 14 WIDTHPROF = None
15 15 HEIGHTPROF = None
16 PREFIX = 'fitgau'
16 PREFIX = 'SpcParam'
17 17
18 18 def __init__(self, **kwargs):
19 19 Figure.__init__(self, **kwargs)
20 20 self.isConfig = False
21 21 self.__nsubplots = 1
22 22
23 23 self.WIDTH = 250
24 24 self.HEIGHT = 250
25 25 self.WIDTHPROF = 120
26 26 self.HEIGHTPROF = 0
27 27 self.counter_imagwr = 0
28 28
29 29 self.PLOT_CODE = SPEC_CODE
30 30
31 31 self.FTP_WEI = None
32 32 self.EXP_CODE = None
33 33 self.SUB_EXP_CODE = None
34 34 self.PLOT_POS = None
35 35
36 36 self.__xfilter_ena = False
37 37 self.__yfilter_ena = False
38 38
39 39 def getSubplots(self):
40 40
41 41 ncol = int(numpy.sqrt(self.nplots)+0.9)
42 42 nrow = int(self.nplots*1./ncol + 0.9)
43 43
44 44 return nrow, ncol
45 45
46 46 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
47 47
48 48 self.__showprofile = showprofile
49 49 self.nplots = nplots
50 50
51 51 ncolspan = 1
52 52 colspan = 1
53 53 if showprofile:
54 54 ncolspan = 3
55 55 colspan = 2
56 56 self.__nsubplots = 2
57 57
58 58 self.createFigure(id = id,
59 59 wintitle = wintitle,
60 60 widthplot = self.WIDTH + self.WIDTHPROF,
61 61 heightplot = self.HEIGHT + self.HEIGHTPROF,
62 62 show=show)
63 63
64 64 nrow, ncol = self.getSubplots()
65 65
66 66 counter = 0
67 67 for y in range(nrow):
68 68 for x in range(ncol):
69 69
70 70 if counter >= self.nplots:
71 71 break
72 72
73 73 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
74 74
75 75 if showprofile:
76 76 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
77 77
78 78 counter += 1
79 79
80 80 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
81 81 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
82 82 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
83 83 server=None, folder=None, username=None, password=None,
84 84 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
85 xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 0):
85 xaxis="frequency", colormap='jet', normFactor=None , Selector = 0):
86 86
87 87 """
88 88
89 89 Input:
90 90 dataOut :
91 91 id :
92 92 wintitle :
93 93 channelList :
94 94 showProfile :
95 95 xmin : None,
96 96 xmax : None,
97 97 ymin : None,
98 98 ymax : None,
99 99 zmin : None,
100 100 zmax : None
101 101 """
102 102 if realtime:
103 103 if not(isRealtime(utcdatatime = dataOut.utctime)):
104 104 print 'Skipping this plot function'
105 105 return
106 106
107 107 if channelList == None:
108 108 channelIndexList = dataOut.channelIndexList
109 109 else:
110 110 channelIndexList = []
111 111 for channel in channelList:
112 112 if channel not in dataOut.channelList:
113 113 raise ValueError, "Channel %d is not in dataOut.channelList" %channel
114 114 channelIndexList.append(dataOut.channelList.index(channel))
115 115
116 116 # if normFactor is None:
117 117 # factor = dataOut.normFactor
118 118 # else:
119 119 # factor = normFactor
120 120 if xaxis == "frequency":
121 x = dataOut.spc_range[0]
121 x = dataOut.spcparam_range[0]
122 122 xlabel = "Frequency (kHz)"
123 123
124 124 elif xaxis == "time":
125 x = dataOut.spc_range[1]
125 x = dataOut.spcparam_range[1]
126 126 xlabel = "Time (ms)"
127 127
128 128 else:
129 x = dataOut.spc_range[2]
129 x = dataOut.spcparam_range[2]
130 130 xlabel = "Velocity (m/s)"
131 131
132 132 ylabel = "Range (Km)"
133 133
134 134 y = dataOut.getHeiRange()
135 135
136 z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor
137 print 'GausSPC', z[0,32,10:40]
136 z = dataOut.SPCparam[Selector] #GauSelector] #dataOut.data_spc/factor
137 #print 'GausSPC', z[0,32,10:40]
138 138 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
139 139 zdB = 10*numpy.log10(z)
140 140
141 141 avg = numpy.average(z, axis=1)
142 142 avgdB = 10*numpy.log10(avg)
143 143
144 144 noise = dataOut.spc_noise
145 145 noisedB = 10*numpy.log10(noise)
146 146
147 147 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
148 148 title = wintitle + " Spectra"
149 149 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
150 150 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
151 151
152 152 if not self.isConfig:
153 153
154 154 nplots = len(channelIndexList)
155 155
156 156 self.setup(id=id,
157 157 nplots=nplots,
158 158 wintitle=wintitle,
159 159 showprofile=showprofile,
160 160 show=show)
161 161
162 162 if xmin == None: xmin = numpy.nanmin(x)
163 163 if xmax == None: xmax = numpy.nanmax(x)
164 164 if ymin == None: ymin = numpy.nanmin(y)
165 165 if ymax == None: ymax = numpy.nanmax(y)
166 166 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
167 167 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
168 168
169 169 self.FTP_WEI = ftp_wei
170 170 self.EXP_CODE = exp_code
171 171 self.SUB_EXP_CODE = sub_exp_code
172 172 self.PLOT_POS = plot_pos
173 173
174 174 self.isConfig = True
175 175
176 176 self.setWinTitle(title)
177 177
178 178 for i in range(self.nplots):
179 179 index = channelIndexList[i]
180 180 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
181 181 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
182 182 if len(dataOut.beam.codeList) != 0:
183 183 title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime)
184 184
185 185 axes = self.axesList[i*self.__nsubplots]
186 186 axes.pcolor(x, y, zdB[index,:,:],
187 187 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
188 188 xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap,
189 189 ticksize=9, cblabel='')
190 190
191 191 if self.__showprofile:
192 192 axes = self.axesList[i*self.__nsubplots +1]
193 193 axes.pline(avgdB[index,:], y,
194 194 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
195 195 xlabel='dB', ylabel='', title='',
196 196 ytick_visible=False,
197 197 grid='x')
198 198
199 199 noiseline = numpy.repeat(noisedB[index], len(y))
200 200 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
201 201
202 202 self.draw()
203 203
204 204 if figfile == None:
205 205 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
206 206 name = str_datetime
207 207 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
208 208 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
209 209 figfile = self.getFilename(name)
210 210
211 211 self.save(figpath=figpath,
212 212 figfile=figfile,
213 213 save=save,
214 214 ftp=ftp,
215 215 wr_period=wr_period,
216 216 thisDatetime=thisDatetime)
217 217
218 218
219 219
220 220 class MomentsPlot(Figure):
221 221
222 222 isConfig = None
223 223 __nsubplots = None
224 224
225 225 WIDTHPROF = None
226 226 HEIGHTPROF = None
227 227 PREFIX = 'prm'
228 228
229 229 def __init__(self, **kwargs):
230 230 Figure.__init__(self, **kwargs)
231 231 self.isConfig = False
232 232 self.__nsubplots = 1
233 233
234 234 self.WIDTH = 280
235 235 self.HEIGHT = 250
236 236 self.WIDTHPROF = 120
237 237 self.HEIGHTPROF = 0
238 238 self.counter_imagwr = 0
239 239
240 240 self.PLOT_CODE = MOMENTS_CODE
241 241
242 242 self.FTP_WEI = None
243 243 self.EXP_CODE = None
244 244 self.SUB_EXP_CODE = None
245 245 self.PLOT_POS = None
246 246
247 247 def getSubplots(self):
248 248
249 249 ncol = int(numpy.sqrt(self.nplots)+0.9)
250 250 nrow = int(self.nplots*1./ncol + 0.9)
251 251
252 252 return nrow, ncol
253 253
254 254 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
255 255
256 256 self.__showprofile = showprofile
257 257 self.nplots = nplots
258 258
259 259 ncolspan = 1
260 260 colspan = 1
261 261 if showprofile:
262 262 ncolspan = 3
263 263 colspan = 2
264 264 self.__nsubplots = 2
265 265
266 266 self.createFigure(id = id,
267 267 wintitle = wintitle,
268 268 widthplot = self.WIDTH + self.WIDTHPROF,
269 269 heightplot = self.HEIGHT + self.HEIGHTPROF,
270 270 show=show)
271 271
272 272 nrow, ncol = self.getSubplots()
273 273
274 274 counter = 0
275 275 for y in range(nrow):
276 276 for x in range(ncol):
277 277
278 278 if counter >= self.nplots:
279 279 break
280 280
281 281 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
282 282
283 283 if showprofile:
284 284 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
285 285
286 286 counter += 1
287 287
288 288 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
289 289 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
290 290 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
291 291 server=None, folder=None, username=None, password=None,
292 292 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
293 293
294 294 """
295 295
296 296 Input:
297 297 dataOut :
298 298 id :
299 299 wintitle :
300 300 channelList :
301 301 showProfile :
302 302 xmin : None,
303 303 xmax : None,
304 304 ymin : None,
305 305 ymax : None,
306 306 zmin : None,
307 307 zmax : None
308 308 """
309 309
310 310 if dataOut.flagNoData:
311 311 return None
312 312
313 313 if realtime:
314 314 if not(isRealtime(utcdatatime = dataOut.utctime)):
315 315 print 'Skipping this plot function'
316 316 return
317 317
318 318 if channelList == None:
319 319 channelIndexList = dataOut.channelIndexList
320 320 else:
321 321 channelIndexList = []
322 322 for channel in channelList:
323 323 if channel not in dataOut.channelList:
324 324 raise ValueError, "Channel %d is not in dataOut.channelList"
325 325 channelIndexList.append(dataOut.channelList.index(channel))
326 326
327 327 factor = dataOut.normFactor
328 328 x = dataOut.abscissaList
329 329 y = dataOut.heightList
330 330
331 331 z = dataOut.data_pre[channelIndexList,:,:]/factor
332 332 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
333 333 avg = numpy.average(z, axis=1)
334 334 noise = dataOut.noise/factor
335 335
336 336 zdB = 10*numpy.log10(z)
337 337 avgdB = 10*numpy.log10(avg)
338 338 noisedB = 10*numpy.log10(noise)
339 339
340 340 #thisDatetime = dataOut.datatime
341 341 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
342 342 title = wintitle + " Parameters"
343 343 xlabel = "Velocity (m/s)"
344 344 ylabel = "Range (Km)"
345 345
346 346 update_figfile = False
347 347
348 348 if not self.isConfig:
349 349
350 350 nplots = len(channelIndexList)
351 351
352 352 self.setup(id=id,
353 353 nplots=nplots,
354 354 wintitle=wintitle,
355 355 showprofile=showprofile,
356 356 show=show)
357 357
358 358 if xmin == None: xmin = numpy.nanmin(x)
359 359 if xmax == None: xmax = numpy.nanmax(x)
360 360 if ymin == None: ymin = numpy.nanmin(y)
361 361 if ymax == None: ymax = numpy.nanmax(y)
362 362 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
363 363 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
364 364
365 365 self.FTP_WEI = ftp_wei
366 366 self.EXP_CODE = exp_code
367 367 self.SUB_EXP_CODE = sub_exp_code
368 368 self.PLOT_POS = plot_pos
369 369
370 370 self.isConfig = True
371 371 update_figfile = True
372 372
373 373 self.setWinTitle(title)
374 374
375 375 for i in range(self.nplots):
376 376 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
377 377 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i], noisedB[i], str_datetime)
378 378 axes = self.axesList[i*self.__nsubplots]
379 379 axes.pcolor(x, y, zdB[i,:,:],
380 380 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
381 381 xlabel=xlabel, ylabel=ylabel, title=title,
382 382 ticksize=9, cblabel='')
383 383 #Mean Line
384 384 mean = dataOut.data_param[i, 1, :]
385 385 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
386 386
387 387 if self.__showprofile:
388 388 axes = self.axesList[i*self.__nsubplots +1]
389 389 axes.pline(avgdB[i], y,
390 390 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
391 391 xlabel='dB', ylabel='', title='',
392 392 ytick_visible=False,
393 393 grid='x')
394 394
395 395 noiseline = numpy.repeat(noisedB[i], len(y))
396 396 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
397 397
398 398 self.draw()
399 399
400 400 self.save(figpath=figpath,
401 401 figfile=figfile,
402 402 save=save,
403 403 ftp=ftp,
404 404 wr_period=wr_period,
405 405 thisDatetime=thisDatetime)
406 406
407 407
408 408
409 409 class SkyMapPlot(Figure):
410 410
411 411 __isConfig = None
412 412 __nsubplots = None
413 413
414 414 WIDTHPROF = None
415 415 HEIGHTPROF = None
416 416 PREFIX = 'mmap'
417 417
418 418 def __init__(self, **kwargs):
419 419 Figure.__init__(self, **kwargs)
420 420 self.isConfig = False
421 421 self.__nsubplots = 1
422 422
423 423 # self.WIDTH = 280
424 424 # self.HEIGHT = 250
425 425 self.WIDTH = 600
426 426 self.HEIGHT = 600
427 427 self.WIDTHPROF = 120
428 428 self.HEIGHTPROF = 0
429 429 self.counter_imagwr = 0
430 430
431 431 self.PLOT_CODE = MSKYMAP_CODE
432 432
433 433 self.FTP_WEI = None
434 434 self.EXP_CODE = None
435 435 self.SUB_EXP_CODE = None
436 436 self.PLOT_POS = None
437 437
438 438 def getSubplots(self):
439 439
440 440 ncol = int(numpy.sqrt(self.nplots)+0.9)
441 441 nrow = int(self.nplots*1./ncol + 0.9)
442 442
443 443 return nrow, ncol
444 444
445 445 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
446 446
447 447 self.__showprofile = showprofile
448 448 self.nplots = nplots
449 449
450 450 ncolspan = 1
451 451 colspan = 1
452 452
453 453 self.createFigure(id = id,
454 454 wintitle = wintitle,
455 455 widthplot = self.WIDTH, #+ self.WIDTHPROF,
456 456 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
457 457 show=show)
458 458
459 459 nrow, ncol = 1,1
460 460 counter = 0
461 461 x = 0
462 462 y = 0
463 463 self.addAxes(1, 1, 0, 0, 1, 1, True)
464 464
465 465 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
466 466 tmin=0, tmax=24, timerange=None,
467 467 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
468 468 server=None, folder=None, username=None, password=None,
469 469 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
470 470
471 471 """
472 472
473 473 Input:
474 474 dataOut :
475 475 id :
476 476 wintitle :
477 477 channelList :
478 478 showProfile :
479 479 xmin : None,
480 480 xmax : None,
481 481 ymin : None,
482 482 ymax : None,
483 483 zmin : None,
484 484 zmax : None
485 485 """
486 486
487 487 arrayParameters = dataOut.data_param
488 488 error = arrayParameters[:,-1]
489 489 indValid = numpy.where(error == 0)[0]
490 490 finalMeteor = arrayParameters[indValid,:]
491 491 finalAzimuth = finalMeteor[:,3]
492 492 finalZenith = finalMeteor[:,4]
493 493
494 494 x = finalAzimuth*numpy.pi/180
495 495 y = finalZenith
496 496 x1 = [dataOut.ltctime, dataOut.ltctime]
497 497
498 498 #thisDatetime = dataOut.datatime
499 499 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
500 500 title = wintitle + " Parameters"
501 501 xlabel = "Zonal Zenith Angle (deg) "
502 502 ylabel = "Meridional Zenith Angle (deg)"
503 503 update_figfile = False
504 504
505 505 if not self.isConfig:
506 506
507 507 nplots = 1
508 508
509 509 self.setup(id=id,
510 510 nplots=nplots,
511 511 wintitle=wintitle,
512 512 showprofile=showprofile,
513 513 show=show)
514 514
515 515 if self.xmin is None and self.xmax is None:
516 516 self.xmin, self.xmax = self.getTimeLim(x1, tmin, tmax, timerange)
517 517
518 518 if timerange != None:
519 519 self.timerange = timerange
520 520 else:
521 521 self.timerange = self.xmax - self.xmin
522 522
523 523 self.FTP_WEI = ftp_wei
524 524 self.EXP_CODE = exp_code
525 525 self.SUB_EXP_CODE = sub_exp_code
526 526 self.PLOT_POS = plot_pos
527 527 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
528 528 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
529 529 self.isConfig = True
530 530 update_figfile = True
531 531
532 532 self.setWinTitle(title)
533 533
534 534 i = 0
535 535 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
536 536
537 537 axes = self.axesList[i*self.__nsubplots]
538 538 nevents = axes.x_buffer.shape[0] + x.shape[0]
539 539 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
540 540 axes.polar(x, y,
541 541 title=title, xlabel=xlabel, ylabel=ylabel,
542 542 ticksize=9, cblabel='')
543 543
544 544 self.draw()
545 545
546 546 self.save(figpath=figpath,
547 547 figfile=figfile,
548 548 save=save,
549 549 ftp=ftp,
550 550 wr_period=wr_period,
551 551 thisDatetime=thisDatetime,
552 552 update_figfile=update_figfile)
553 553
554 554 if dataOut.ltctime >= self.xmax:
555 555 self.isConfigmagwr = wr_period
556 556 self.isConfig = False
557 557 update_figfile = True
558 558 axes.__firsttime = True
559 559 self.xmin += self.timerange
560 560 self.xmax += self.timerange
561 561
562 562
563 563
564 564
565 565 class WindProfilerPlot(Figure):
566 566
567 567 __isConfig = None
568 568 __nsubplots = None
569 569
570 570 WIDTHPROF = None
571 571 HEIGHTPROF = None
572 572 PREFIX = 'wind'
573 573
574 574 def __init__(self, **kwargs):
575 575 Figure.__init__(self, **kwargs)
576 576 self.timerange = None
577 577 self.isConfig = False
578 578 self.__nsubplots = 1
579 579
580 580 self.WIDTH = 800
581 581 self.HEIGHT = 300
582 582 self.WIDTHPROF = 120
583 583 self.HEIGHTPROF = 0
584 584 self.counter_imagwr = 0
585 585
586 586 self.PLOT_CODE = WIND_CODE
587 587
588 588 self.FTP_WEI = None
589 589 self.EXP_CODE = None
590 590 self.SUB_EXP_CODE = None
591 591 self.PLOT_POS = None
592 592 self.tmin = None
593 593 self.tmax = None
594 594
595 595 self.xmin = None
596 596 self.xmax = None
597 597
598 598 self.figfile = None
599 599
600 600 def getSubplots(self):
601 601
602 602 ncol = 1
603 603 nrow = self.nplots
604 604
605 605 return nrow, ncol
606 606
607 607 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
608 608
609 609 self.__showprofile = showprofile
610 610 self.nplots = nplots
611 611
612 612 ncolspan = 1
613 613 colspan = 1
614 614
615 615 self.createFigure(id = id,
616 616 wintitle = wintitle,
617 617 widthplot = self.WIDTH + self.WIDTHPROF,
618 618 heightplot = self.HEIGHT + self.HEIGHTPROF,
619 619 show=show)
620 620
621 621 nrow, ncol = self.getSubplots()
622 622
623 623 counter = 0
624 624 for y in range(nrow):
625 625 if counter >= self.nplots:
626 626 break
627 627
628 628 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
629 629 counter += 1
630 630
631 631 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='False',
632 632 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
633 633 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
634 634 timerange=None, SNRthresh = None,
635 635 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
636 636 server=None, folder=None, username=None, password=None,
637 637 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
638 638 """
639 639
640 640 Input:
641 641 dataOut :
642 642 id :
643 643 wintitle :
644 644 channelList :
645 645 showProfile :
646 646 xmin : None,
647 647 xmax : None,
648 648 ymin : None,
649 649 ymax : None,
650 650 zmin : None,
651 651 zmax : None
652 652 """
653 653
654 654 # if timerange is not None:
655 655 # self.timerange = timerange
656 656 #
657 657 # tmin = None
658 658 # tmax = None
659 659
660 660 x = dataOut.getTimeRange1(dataOut.paramInterval)
661 661 y = dataOut.heightList
662 662 z = dataOut.data_output.copy()
663 663 nplots = z.shape[0] #Number of wind dimensions estimated
664 664 nplotsw = nplots
665 665
666 666
667 667 #If there is a SNR function defined
668 668 if dataOut.data_SNR is not None:
669 669 nplots += 1
670 670 SNR = dataOut.data_SNR[0]
671 671 SNRavg = SNR#numpy.average(SNR, axis=0)
672 672
673 673 SNRdB = 10*numpy.log10(SNR)
674 674 SNRavgdB = 10*numpy.log10(SNRavg)
675 675
676 676 if SNRthresh == None:
677 677 SNRthresh = -5.0
678 678 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
679 679
680 680 for i in range(nplotsw):
681 681 z[i,ind] = numpy.nan
682 682
683 683 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
684 684 #thisDatetime = datetime.datetime.now()
685 685 title = wintitle + "Wind"
686 686 xlabel = ""
687 687 ylabel = "Height (km)"
688 688 update_figfile = False
689 689
690 690 if not self.isConfig:
691 691
692 692 self.setup(id=id,
693 693 nplots=nplots,
694 694 wintitle=wintitle,
695 695 showprofile=showprofile,
696 696 show=show)
697 697
698 698 if timerange is not None:
699 699 self.timerange = timerange
700 700
701 701 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
702 702
703 703 if ymin == None: ymin = numpy.nanmin(y)
704 704 if ymax == None: ymax = numpy.nanmax(y)
705 705
706 706 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
707 707 #if numpy.isnan(zmax): zmax = 50
708 708 if zmin == None: zmin = -zmax
709 709
710 710 if nplotsw == 3:
711 711 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
712 712 if zmin_ver == None: zmin_ver = -zmax_ver
713 713
714 714 if dataOut.data_SNR is not None:
715 715 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
716 716 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
717 717
718 718
719 719 self.FTP_WEI = ftp_wei
720 720 self.EXP_CODE = exp_code
721 721 self.SUB_EXP_CODE = sub_exp_code
722 722 self.PLOT_POS = plot_pos
723 723
724 724 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
725 725 self.isConfig = True
726 726 self.figfile = figfile
727 727 update_figfile = True
728 728
729 729 self.setWinTitle(title)
730 730
731 731 if ((self.xmax - x[1]) < (x[1]-x[0])):
732 732 x[1] = self.xmax
733 733
734 734 strWind = ['Zonal', 'Meridional', 'Vertical']
735 735 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
736 736 zmaxVector = [zmax, zmax, zmax_ver]
737 737 zminVector = [zmin, zmin, zmin_ver]
738 738 windFactor = [1,1,100]
739 739
740 740 for i in range(nplotsw):
741 741
742 742 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
743 743 axes = self.axesList[i*self.__nsubplots]
744 744
745 745 z1 = z[i,:].reshape((1,-1))*windFactor[i]
746 746
747 747 print 'x', x
748 748 print datetime.datetime.utcfromtimestamp(x[0])
749 749 print datetime.datetime.utcfromtimestamp(x[1])
750 750
751 751 #z1=numpy.ma.masked_where(z1==0.,z1)
752 752
753 753 axes.pcolorbuffer(x, y, z1,
754 754 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
755 755 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
756 756 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="seismic" )
757 757
758 758 if dataOut.data_SNR is not None:
759 759 i += 1
760 760 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
761 761 axes = self.axesList[i*self.__nsubplots]
762 762 SNRavgdB = SNRavgdB.reshape((1,-1))
763 763 axes.pcolorbuffer(x, y, SNRavgdB,
764 764 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
765 765 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
766 766 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
767 767
768 768 self.draw()
769 769
770 770 self.save(figpath=figpath,
771 771 figfile=figfile,
772 772 save=save,
773 773 ftp=ftp,
774 774 wr_period=wr_period,
775 775 thisDatetime=thisDatetime,
776 776 update_figfile=update_figfile)
777 777
778 778 if dataOut.ltctime + dataOut.paramInterval >= self.xmax:
779 779 self.counter_imagwr = wr_period
780 780 self.isConfig = False
781 781 update_figfile = True
782 782
783 783
784 784 class ParametersPlot(Figure):
785 785
786 786 __isConfig = None
787 787 __nsubplots = None
788 788
789 789 WIDTHPROF = None
790 790 HEIGHTPROF = None
791 791 PREFIX = 'param'
792 792
793 793 nplots = None
794 794 nchan = None
795 795
796 796 def __init__(self, **kwargs):
797 797 Figure.__init__(self, **kwargs)
798 798 self.timerange = None
799 799 self.isConfig = False
800 800 self.__nsubplots = 1
801 801
802 802 self.WIDTH = 800
803 803 self.HEIGHT = 180
804 804 self.WIDTHPROF = 120
805 805 self.HEIGHTPROF = 0
806 806 self.counter_imagwr = 0
807 807
808 808 self.PLOT_CODE = RTI_CODE
809 809
810 810 self.FTP_WEI = None
811 811 self.EXP_CODE = None
812 812 self.SUB_EXP_CODE = None
813 813 self.PLOT_POS = None
814 814 self.tmin = None
815 815 self.tmax = None
816 816
817 817 self.xmin = None
818 818 self.xmax = None
819 819
820 820 self.figfile = None
821 821
822 822 def getSubplots(self):
823 823
824 824 ncol = 1
825 825 nrow = self.nplots
826 826
827 827 return nrow, ncol
828 828
829 829 def setup(self, id, nplots, wintitle, show=True):
830 830
831 831 self.nplots = nplots
832 832
833 833 ncolspan = 1
834 834 colspan = 1
835 835
836 836 self.createFigure(id = id,
837 837 wintitle = wintitle,
838 838 widthplot = self.WIDTH + self.WIDTHPROF,
839 839 heightplot = self.HEIGHT + self.HEIGHTPROF,
840 840 show=show)
841 841
842 842 nrow, ncol = self.getSubplots()
843 843
844 844 counter = 0
845 845 for y in range(nrow):
846 846 for x in range(ncol):
847 847
848 848 if counter >= self.nplots:
849 849 break
850 850
851 851 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
852 852
853 853 counter += 1
854 854
855 855 def run(self, dataOut, id, wintitle="", channelList=None, paramIndex = 0, colormap="jet",
856 856 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, timerange=None,
857 857 showSNR=False, SNRthresh = -numpy.inf, SNRmin=None, SNRmax=None,
858 858 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
859 859 server=None, folder=None, username=None, password=None,
860 860 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, HEIGHT=None):
861 861 """
862 862
863 863 Input:
864 864 dataOut :
865 865 id :
866 866 wintitle :
867 867 channelList :
868 868 showProfile :
869 869 xmin : None,
870 870 xmax : None,
871 871 ymin : None,
872 872 ymax : None,
873 873 zmin : None,
874 874 zmax : None
875 875 """
876 876
877 877 if HEIGHT is not None:
878 878 self.HEIGHT = HEIGHT
879 879
880 880
881 881 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
882 882 return
883 883
884 884 if channelList == None:
885 885 channelIndexList = range(dataOut.data_param.shape[0])
886 886 else:
887 887 channelIndexList = []
888 888 for channel in channelList:
889 889 if channel not in dataOut.channelList:
890 890 raise ValueError, "Channel %d is not in dataOut.channelList"
891 891 channelIndexList.append(dataOut.channelList.index(channel))
892 892
893 893 x = dataOut.getTimeRange1(dataOut.paramInterval)
894 894 y = dataOut.getHeiRange()
895 895
896 896 if dataOut.data_param.ndim == 3:
897 897 z = dataOut.data_param[channelIndexList,paramIndex,:]
898 898 else:
899 899 z = dataOut.data_param[channelIndexList,:]
900 900
901 901 if showSNR:
902 902 #SNR data
903 903 SNRarray = dataOut.data_SNR[channelIndexList,:]
904 904 SNRdB = 10*numpy.log10(SNRarray)
905 905 ind = numpy.where(SNRdB < SNRthresh)
906 906 z[ind] = numpy.nan
907 907
908 908 thisDatetime = dataOut.datatime
909 909 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
910 910 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
911 911 xlabel = ""
912 912 ylabel = "Range (Km)"
913 913
914 914 update_figfile = False
915 915
916 916 if not self.isConfig:
917 917
918 918 nchan = len(channelIndexList)
919 919 self.nchan = nchan
920 920 self.plotFact = 1
921 921 nplots = nchan
922 922
923 923 if showSNR:
924 924 nplots = nchan*2
925 925 self.plotFact = 2
926 926 if SNRmin == None: SNRmin = numpy.nanmin(SNRdB)
927 927 if SNRmax == None: SNRmax = numpy.nanmax(SNRdB)
928 928
929 929 self.setup(id=id,
930 930 nplots=nplots,
931 931 wintitle=wintitle,
932 932 show=show)
933 933
934 934 if timerange != None:
935 935 self.timerange = timerange
936 936
937 937 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
938 938
939 939 if ymin == None: ymin = numpy.nanmin(y)
940 940 if ymax == None: ymax = numpy.nanmax(y)
941 941 if zmin == None: zmin = numpy.nanmin(z)
942 942 if zmax == None: zmax = numpy.nanmax(z)
943 943
944 944 self.FTP_WEI = ftp_wei
945 945 self.EXP_CODE = exp_code
946 946 self.SUB_EXP_CODE = sub_exp_code
947 947 self.PLOT_POS = plot_pos
948 948
949 949 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
950 950 self.isConfig = True
951 951 self.figfile = figfile
952 952 update_figfile = True
953 953
954 954 self.setWinTitle(title)
955 955
956 956 for i in range(self.nchan):
957 957 index = channelIndexList[i]
958 958 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
959 959 axes = self.axesList[i*self.plotFact]
960 960 z1 = z[i,:].reshape((1,-1))
961 961 axes.pcolorbuffer(x, y, z1,
962 962 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
963 963 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
964 964 ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
965 965
966 966 if showSNR:
967 967 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
968 968 axes = self.axesList[i*self.plotFact + 1]
969 969 SNRdB1 = SNRdB[i,:].reshape((1,-1))
970 970 axes.pcolorbuffer(x, y, SNRdB1,
971 971 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
972 972 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
973 973 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
974 974
975 975
976 976 self.draw()
977 977
978 978 if dataOut.ltctime >= self.xmax:
979 979 self.counter_imagwr = wr_period
980 980 self.isConfig = False
981 981 update_figfile = True
982 982
983 983 self.save(figpath=figpath,
984 984 figfile=figfile,
985 985 save=save,
986 986 ftp=ftp,
987 987 wr_period=wr_period,
988 988 thisDatetime=thisDatetime,
989 989 update_figfile=update_figfile)
990 990
991 991
992 992
993 993 class Parameters1Plot(Figure):
994 994
995 995 __isConfig = None
996 996 __nsubplots = None
997 997
998 998 WIDTHPROF = None
999 999 HEIGHTPROF = None
1000 1000 PREFIX = 'prm'
1001 1001
1002 1002 def __init__(self, **kwargs):
1003 1003 Figure.__init__(self, **kwargs)
1004 1004 self.timerange = 2*60*60
1005 1005 self.isConfig = False
1006 1006 self.__nsubplots = 1
1007 1007
1008 1008 self.WIDTH = 800
1009 1009 self.HEIGHT = 180
1010 1010 self.WIDTHPROF = 120
1011 1011 self.HEIGHTPROF = 0
1012 1012 self.counter_imagwr = 0
1013 1013
1014 1014 self.PLOT_CODE = PARMS_CODE
1015 1015
1016 1016 self.FTP_WEI = None
1017 1017 self.EXP_CODE = None
1018 1018 self.SUB_EXP_CODE = None
1019 1019 self.PLOT_POS = None
1020 1020 self.tmin = None
1021 1021 self.tmax = None
1022 1022
1023 1023 self.xmin = None
1024 1024 self.xmax = None
1025 1025
1026 1026 self.figfile = None
1027 1027
1028 1028 def getSubplots(self):
1029 1029
1030 1030 ncol = 1
1031 1031 nrow = self.nplots
1032 1032
1033 1033 return nrow, ncol
1034 1034
1035 1035 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1036 1036
1037 1037 self.__showprofile = showprofile
1038 1038 self.nplots = nplots
1039 1039
1040 1040 ncolspan = 1
1041 1041 colspan = 1
1042 1042
1043 1043 self.createFigure(id = id,
1044 1044 wintitle = wintitle,
1045 1045 widthplot = self.WIDTH + self.WIDTHPROF,
1046 1046 heightplot = self.HEIGHT + self.HEIGHTPROF,
1047 1047 show=show)
1048 1048
1049 1049 nrow, ncol = self.getSubplots()
1050 1050
1051 1051 counter = 0
1052 1052 for y in range(nrow):
1053 1053 for x in range(ncol):
1054 1054
1055 1055 if counter >= self.nplots:
1056 1056 break
1057 1057
1058 1058 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1059 1059
1060 1060 if showprofile:
1061 1061 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
1062 1062
1063 1063 counter += 1
1064 1064
1065 1065 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
1066 1066 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
1067 1067 parameterIndex = None, onlyPositive = False,
1068 1068 SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, onlySNR = False,
1069 1069 DOP = True,
1070 1070 zlabel = "", parameterName = "", parameterObject = "data_param",
1071 1071 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1072 1072 server=None, folder=None, username=None, password=None,
1073 1073 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1074 1074 #print inspect.getargspec(self.run).args
1075 1075 """
1076 1076
1077 1077 Input:
1078 1078 dataOut :
1079 1079 id :
1080 1080 wintitle :
1081 1081 channelList :
1082 1082 showProfile :
1083 1083 xmin : None,
1084 1084 xmax : None,
1085 1085 ymin : None,
1086 1086 ymax : None,
1087 1087 zmin : None,
1088 1088 zmax : None
1089 1089 """
1090 1090
1091 1091 data_param = getattr(dataOut, parameterObject)
1092 1092
1093 1093 if channelList == None:
1094 1094 channelIndexList = numpy.arange(data_param.shape[0])
1095 1095 else:
1096 1096 channelIndexList = numpy.array(channelList)
1097 1097
1098 1098 nchan = len(channelIndexList) #Number of channels being plotted
1099 1099
1100 1100 if nchan < 1:
1101 1101 return
1102 1102
1103 1103 nGraphsByChannel = 0
1104 1104
1105 1105 if SNR:
1106 1106 nGraphsByChannel += 1
1107 1107 if DOP:
1108 1108 nGraphsByChannel += 1
1109 1109
1110 1110 if nGraphsByChannel < 1:
1111 1111 return
1112 1112
1113 1113 nplots = nGraphsByChannel*nchan
1114 1114
1115 1115 if timerange is not None:
1116 1116 self.timerange = timerange
1117 1117
1118 1118 #tmin = None
1119 1119 #tmax = None
1120 1120 if parameterIndex == None:
1121 1121 parameterIndex = 1
1122 1122
1123 1123 x = dataOut.getTimeRange1(dataOut.paramInterval)
1124 1124 y = dataOut.heightList
1125 1125 z = data_param[channelIndexList,parameterIndex,:].copy()
1126 1126
1127 1127 zRange = dataOut.abscissaList
1128 1128 # nChannels = z.shape[0] #Number of wind dimensions estimated
1129 1129 # thisDatetime = dataOut.datatime
1130 1130
1131 1131 if dataOut.data_SNR is not None:
1132 1132 SNRarray = dataOut.data_SNR[channelIndexList,:]
1133 1133 SNRdB = 10*numpy.log10(SNRarray)
1134 1134 # SNRavgdB = 10*numpy.log10(SNRavg)
1135 1135 ind = numpy.where(SNRdB < 10**(SNRthresh/10))
1136 1136 z[ind] = numpy.nan
1137 1137
1138 1138 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1139 1139 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
1140 1140 xlabel = ""
1141 1141 ylabel = "Range (Km)"
1142 1142
1143 1143 if (SNR and not onlySNR): nplots = 2*nplots
1144 1144
1145 1145 if onlyPositive:
1146 1146 colormap = "jet"
1147 1147 zmin = 0
1148 1148 else: colormap = "RdBu_r"
1149 1149
1150 1150 if not self.isConfig:
1151 1151
1152 1152 self.setup(id=id,
1153 1153 nplots=nplots,
1154 1154 wintitle=wintitle,
1155 1155 showprofile=showprofile,
1156 1156 show=show)
1157 1157
1158 1158 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1159 1159
1160 1160 if ymin == None: ymin = numpy.nanmin(y)
1161 1161 if ymax == None: ymax = numpy.nanmax(y)
1162 1162 if zmin == None: zmin = numpy.nanmin(zRange)
1163 1163 if zmax == None: zmax = numpy.nanmax(zRange)
1164 1164
1165 1165 if SNR:
1166 1166 if SNRmin == None: SNRmin = numpy.nanmin(SNRdB)
1167 1167 if SNRmax == None: SNRmax = numpy.nanmax(SNRdB)
1168 1168
1169 1169 self.FTP_WEI = ftp_wei
1170 1170 self.EXP_CODE = exp_code
1171 1171 self.SUB_EXP_CODE = sub_exp_code
1172 1172 self.PLOT_POS = plot_pos
1173 1173
1174 1174 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1175 1175 self.isConfig = True
1176 1176 self.figfile = figfile
1177 1177
1178 1178 self.setWinTitle(title)
1179 1179
1180 1180 if ((self.xmax - x[1]) < (x[1]-x[0])):
1181 1181 x[1] = self.xmax
1182 1182
1183 1183 for i in range(nchan):
1184 1184
1185 1185 if (SNR and not onlySNR): j = 2*i
1186 1186 else: j = i
1187 1187
1188 1188 j = nGraphsByChannel*i
1189 1189
1190 1190 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
1191 1191 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
1192 1192
1193 1193 if not onlySNR:
1194 1194 axes = self.axesList[j*self.__nsubplots]
1195 1195 z1 = z[i,:].reshape((1,-1))
1196 1196 axes.pcolorbuffer(x, y, z1,
1197 1197 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1198 1198 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
1199 1199 ticksize=9, cblabel=zlabel, cbsize="1%")
1200 1200
1201 1201 if DOP:
1202 1202 title = "%s Channel %d: %s" %(parameterName, channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1203 1203
1204 1204 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
1205 1205 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
1206 1206 axes = self.axesList[j]
1207 1207 z1 = z[i,:].reshape((1,-1))
1208 1208 axes.pcolorbuffer(x, y, z1,
1209 1209 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1210 1210 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
1211 1211 ticksize=9, cblabel=zlabel, cbsize="1%")
1212 1212
1213 1213 if SNR:
1214 1214 title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1215 1215 axes = self.axesList[(j)*self.__nsubplots]
1216 1216 if not onlySNR:
1217 1217 axes = self.axesList[(j + 1)*self.__nsubplots]
1218 1218
1219 1219 axes = self.axesList[(j + nGraphsByChannel-1)]
1220 1220
1221 1221 z1 = SNRdB[i,:].reshape((1,-1))
1222 1222 axes.pcolorbuffer(x, y, z1,
1223 1223 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1224 1224 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap="jet",
1225 1225 ticksize=9, cblabel=zlabel, cbsize="1%")
1226 1226
1227 1227
1228 1228
1229 1229 self.draw()
1230 1230
1231 1231 if x[1] >= self.axesList[0].xmax:
1232 1232 self.counter_imagwr = wr_period
1233 1233 self.isConfig = False
1234 1234 self.figfile = None
1235 1235
1236 1236 self.save(figpath=figpath,
1237 1237 figfile=figfile,
1238 1238 save=save,
1239 1239 ftp=ftp,
1240 1240 wr_period=wr_period,
1241 1241 thisDatetime=thisDatetime,
1242 1242 update_figfile=False)
1243 1243
1244 1244 class SpectralFittingPlot(Figure):
1245 1245
1246 1246 __isConfig = None
1247 1247 __nsubplots = None
1248 1248
1249 1249 WIDTHPROF = None
1250 1250 HEIGHTPROF = None
1251 1251 PREFIX = 'prm'
1252 1252
1253 1253
1254 1254 N = None
1255 1255 ippSeconds = None
1256 1256
1257 1257 def __init__(self, **kwargs):
1258 1258 Figure.__init__(self, **kwargs)
1259 1259 self.isConfig = False
1260 1260 self.__nsubplots = 1
1261 1261
1262 1262 self.PLOT_CODE = SPECFIT_CODE
1263 1263
1264 1264 self.WIDTH = 450
1265 1265 self.HEIGHT = 250
1266 1266 self.WIDTHPROF = 0
1267 1267 self.HEIGHTPROF = 0
1268 1268
1269 1269 def getSubplots(self):
1270 1270
1271 1271 ncol = int(numpy.sqrt(self.nplots)+0.9)
1272 1272 nrow = int(self.nplots*1./ncol + 0.9)
1273 1273
1274 1274 return nrow, ncol
1275 1275
1276 1276 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
1277 1277
1278 1278 showprofile = False
1279 1279 self.__showprofile = showprofile
1280 1280 self.nplots = nplots
1281 1281
1282 1282 ncolspan = 5
1283 1283 colspan = 4
1284 1284 if showprofile:
1285 1285 ncolspan = 5
1286 1286 colspan = 4
1287 1287 self.__nsubplots = 2
1288 1288
1289 1289 self.createFigure(id = id,
1290 1290 wintitle = wintitle,
1291 1291 widthplot = self.WIDTH + self.WIDTHPROF,
1292 1292 heightplot = self.HEIGHT + self.HEIGHTPROF,
1293 1293 show=show)
1294 1294
1295 1295 nrow, ncol = self.getSubplots()
1296 1296
1297 1297 counter = 0
1298 1298 for y in range(nrow):
1299 1299 for x in range(ncol):
1300 1300
1301 1301 if counter >= self.nplots:
1302 1302 break
1303 1303
1304 1304 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1305 1305
1306 1306 if showprofile:
1307 1307 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
1308 1308
1309 1309 counter += 1
1310 1310
1311 1311 def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
1312 1312 xmin=None, xmax=None, ymin=None, ymax=None,
1313 1313 save=False, figpath='./', figfile=None, show=True):
1314 1314
1315 1315 """
1316 1316
1317 1317 Input:
1318 1318 dataOut :
1319 1319 id :
1320 1320 wintitle :
1321 1321 channelList :
1322 1322 showProfile :
1323 1323 xmin : None,
1324 1324 xmax : None,
1325 1325 zmin : None,
1326 1326 zmax : None
1327 1327 """
1328 1328
1329 1329 if cutHeight==None:
1330 1330 h=270
1331 1331 heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
1332 1332 cutHeight = dataOut.heightList[heightindex]
1333 1333
1334 1334 factor = dataOut.normFactor
1335 1335 x = dataOut.abscissaList[:-1]
1336 1336 #y = dataOut.getHeiRange()
1337 1337
1338 1338 z = dataOut.data_pre[:,:,heightindex]/factor
1339 1339 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1340 1340 avg = numpy.average(z, axis=1)
1341 1341 listChannels = z.shape[0]
1342 1342
1343 1343 #Reconstruct Function
1344 1344 if fit==True:
1345 1345 groupArray = dataOut.groupList
1346 1346 listChannels = groupArray.reshape((groupArray.size))
1347 1347 listChannels.sort()
1348 1348 spcFitLine = numpy.zeros(z.shape)
1349 1349 constants = dataOut.constants
1350 1350
1351 1351 nGroups = groupArray.shape[0]
1352 1352 nChannels = groupArray.shape[1]
1353 1353 nProfiles = z.shape[1]
1354 1354
1355 1355 for f in range(nGroups):
1356 1356 groupChann = groupArray[f,:]
1357 1357 p = dataOut.data_param[f,:,heightindex]
1358 1358 # p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
1359 1359 fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
1360 1360 fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
1361 1361 spcFitLine[groupChann,:] = fitLineAux
1362 1362 # spcFitLine = spcFitLine/factor
1363 1363
1364 1364 z = z[listChannels,:]
1365 1365 spcFitLine = spcFitLine[listChannels,:]
1366 1366 spcFitLinedB = 10*numpy.log10(spcFitLine)
1367 1367
1368 1368 zdB = 10*numpy.log10(z)
1369 1369 #thisDatetime = dataOut.datatime
1370 1370 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1371 1371 title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1372 1372 xlabel = "Velocity (m/s)"
1373 1373 ylabel = "Spectrum"
1374 1374
1375 1375 if not self.isConfig:
1376 1376
1377 1377 nplots = listChannels.size
1378 1378
1379 1379 self.setup(id=id,
1380 1380 nplots=nplots,
1381 1381 wintitle=wintitle,
1382 1382 showprofile=showprofile,
1383 1383 show=show)
1384 1384
1385 1385 if xmin == None: xmin = numpy.nanmin(x)
1386 1386 if xmax == None: xmax = numpy.nanmax(x)
1387 1387 if ymin == None: ymin = numpy.nanmin(zdB)
1388 1388 if ymax == None: ymax = numpy.nanmax(zdB)+2
1389 1389
1390 1390 self.isConfig = True
1391 1391
1392 1392 self.setWinTitle(title)
1393 1393 for i in range(self.nplots):
1394 1394 # title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
1395 1395 title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i])
1396 1396 axes = self.axesList[i*self.__nsubplots]
1397 1397 if fit == False:
1398 1398 axes.pline(x, zdB[i,:],
1399 1399 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1400 1400 xlabel=xlabel, ylabel=ylabel, title=title
1401 1401 )
1402 1402 if fit == True:
1403 1403 fitline=spcFitLinedB[i,:]
1404 1404 y=numpy.vstack([zdB[i,:],fitline] )
1405 1405 legendlabels=['Data','Fitting']
1406 1406 axes.pmultilineyaxis(x, y,
1407 1407 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1408 1408 xlabel=xlabel, ylabel=ylabel, title=title,
1409 1409 legendlabels=legendlabels, marker=None,
1410 1410 linestyle='solid', grid='both')
1411 1411
1412 1412 self.draw()
1413 1413
1414 1414 self.save(figpath=figpath,
1415 1415 figfile=figfile,
1416 1416 save=save,
1417 1417 ftp=ftp,
1418 1418 wr_period=wr_period,
1419 1419 thisDatetime=thisDatetime)
1420 1420
1421 1421
1422 1422 class EWDriftsPlot(Figure):
1423 1423
1424 1424 __isConfig = None
1425 1425 __nsubplots = None
1426 1426
1427 1427 WIDTHPROF = None
1428 1428 HEIGHTPROF = None
1429 1429 PREFIX = 'drift'
1430 1430
1431 1431 def __init__(self, **kwargs):
1432 1432 Figure.__init__(self, **kwargs)
1433 1433 self.timerange = 2*60*60
1434 1434 self.isConfig = False
1435 1435 self.__nsubplots = 1
1436 1436
1437 1437 self.WIDTH = 800
1438 1438 self.HEIGHT = 150
1439 1439 self.WIDTHPROF = 120
1440 1440 self.HEIGHTPROF = 0
1441 1441 self.counter_imagwr = 0
1442 1442
1443 1443 self.PLOT_CODE = EWDRIFT_CODE
1444 1444
1445 1445 self.FTP_WEI = None
1446 1446 self.EXP_CODE = None
1447 1447 self.SUB_EXP_CODE = None
1448 1448 self.PLOT_POS = None
1449 1449 self.tmin = None
1450 1450 self.tmax = None
1451 1451
1452 1452 self.xmin = None
1453 1453 self.xmax = None
1454 1454
1455 1455 self.figfile = None
1456 1456
1457 1457 def getSubplots(self):
1458 1458
1459 1459 ncol = 1
1460 1460 nrow = self.nplots
1461 1461
1462 1462 return nrow, ncol
1463 1463
1464 1464 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1465 1465
1466 1466 self.__showprofile = showprofile
1467 1467 self.nplots = nplots
1468 1468
1469 1469 ncolspan = 1
1470 1470 colspan = 1
1471 1471
1472 1472 self.createFigure(id = id,
1473 1473 wintitle = wintitle,
1474 1474 widthplot = self.WIDTH + self.WIDTHPROF,
1475 1475 heightplot = self.HEIGHT + self.HEIGHTPROF,
1476 1476 show=show)
1477 1477
1478 1478 nrow, ncol = self.getSubplots()
1479 1479
1480 1480 counter = 0
1481 1481 for y in range(nrow):
1482 1482 if counter >= self.nplots:
1483 1483 break
1484 1484
1485 1485 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
1486 1486 counter += 1
1487 1487
1488 1488 def run(self, dataOut, id, wintitle="", channelList=None,
1489 1489 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
1490 1490 zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
1491 1491 timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
1492 1492 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1493 1493 server=None, folder=None, username=None, password=None,
1494 1494 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1495 1495 """
1496 1496
1497 1497 Input:
1498 1498 dataOut :
1499 1499 id :
1500 1500 wintitle :
1501 1501 channelList :
1502 1502 showProfile :
1503 1503 xmin : None,
1504 1504 xmax : None,
1505 1505 ymin : None,
1506 1506 ymax : None,
1507 1507 zmin : None,
1508 1508 zmax : None
1509 1509 """
1510 1510
1511 1511 if timerange is not None:
1512 1512 self.timerange = timerange
1513 1513
1514 1514 tmin = None
1515 1515 tmax = None
1516 1516
1517 1517 x = dataOut.getTimeRange1(dataOut.outputInterval)
1518 1518 # y = dataOut.heightList
1519 1519 y = dataOut.heightList
1520 1520
1521 1521 z = dataOut.data_output
1522 1522 nplots = z.shape[0] #Number of wind dimensions estimated
1523 1523 nplotsw = nplots
1524 1524
1525 1525 #If there is a SNR function defined
1526 1526 if dataOut.data_SNR is not None:
1527 1527 nplots += 1
1528 1528 SNR = dataOut.data_SNR
1529 1529
1530 1530 if SNR_1:
1531 1531 SNR += 1
1532 1532
1533 1533 SNRavg = numpy.average(SNR, axis=0)
1534 1534
1535 1535 SNRdB = 10*numpy.log10(SNR)
1536 1536 SNRavgdB = 10*numpy.log10(SNRavg)
1537 1537
1538 1538 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
1539 1539
1540 1540 for i in range(nplotsw):
1541 1541 z[i,ind] = numpy.nan
1542 1542
1543 1543
1544 1544 showprofile = False
1545 1545 # thisDatetime = dataOut.datatime
1546 1546 thisDatetime = datetime.datetime.utcfromtimestamp(x[1])
1547 1547 title = wintitle + " EW Drifts"
1548 1548 xlabel = ""
1549 1549 ylabel = "Height (Km)"
1550 1550
1551 1551 if not self.isConfig:
1552 1552
1553 1553 self.setup(id=id,
1554 1554 nplots=nplots,
1555 1555 wintitle=wintitle,
1556 1556 showprofile=showprofile,
1557 1557 show=show)
1558 1558
1559 1559 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1560 1560
1561 1561 if ymin == None: ymin = numpy.nanmin(y)
1562 1562 if ymax == None: ymax = numpy.nanmax(y)
1563 1563
1564 1564 if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
1565 1565 if zminZonal == None: zminZonal = -zmaxZonal
1566 1566 if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
1567 1567 if zminVertical == None: zminVertical = -zmaxVertical
1568 1568
1569 1569 if dataOut.data_SNR is not None:
1570 1570 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
1571 1571 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
1572 1572
1573 1573 self.FTP_WEI = ftp_wei
1574 1574 self.EXP_CODE = exp_code
1575 1575 self.SUB_EXP_CODE = sub_exp_code
1576 1576 self.PLOT_POS = plot_pos
1577 1577
1578 1578 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1579 1579 self.isConfig = True
1580 1580
1581 1581
1582 1582 self.setWinTitle(title)
1583 1583
1584 1584 if ((self.xmax - x[1]) < (x[1]-x[0])):
1585 1585 x[1] = self.xmax
1586 1586
1587 1587 strWind = ['Zonal','Vertical']
1588 1588 strCb = 'Velocity (m/s)'
1589 1589 zmaxVector = [zmaxZonal, zmaxVertical]
1590 1590 zminVector = [zminZonal, zminVertical]
1591 1591
1592 1592 for i in range(nplotsw):
1593 1593
1594 1594 title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1595 1595 axes = self.axesList[i*self.__nsubplots]
1596 1596
1597 1597 z1 = z[i,:].reshape((1,-1))
1598 1598
1599 1599 axes.pcolorbuffer(x, y, z1,
1600 1600 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
1601 1601 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1602 1602 ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
1603 1603
1604 1604 if dataOut.data_SNR is not None:
1605 1605 i += 1
1606 1606 if SNR_1:
1607 1607 title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1608 1608 else:
1609 1609 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1610 1610 axes = self.axesList[i*self.__nsubplots]
1611 1611 SNRavgdB = SNRavgdB.reshape((1,-1))
1612 1612
1613 1613 axes.pcolorbuffer(x, y, SNRavgdB,
1614 1614 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1615 1615 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1616 1616 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
1617 1617
1618 1618 self.draw()
1619 1619
1620 1620 if x[1] >= self.axesList[0].xmax:
1621 1621 self.counter_imagwr = wr_period
1622 1622 self.isConfig = False
1623 1623 self.figfile = None
1624 1624
1625 1625
1626 1626
1627 1627
1628 1628 class PhasePlot(Figure):
1629 1629
1630 1630 __isConfig = None
1631 1631 __nsubplots = None
1632 1632
1633 1633 PREFIX = 'mphase'
1634 1634
1635 1635 def __init__(self, **kwargs):
1636 1636 Figure.__init__(self, **kwargs)
1637 1637 self.timerange = 24*60*60
1638 1638 self.isConfig = False
1639 1639 self.__nsubplots = 1
1640 1640 self.counter_imagwr = 0
1641 1641 self.WIDTH = 600
1642 1642 self.HEIGHT = 300
1643 1643 self.WIDTHPROF = 120
1644 1644 self.HEIGHTPROF = 0
1645 1645 self.xdata = None
1646 1646 self.ydata = None
1647 1647
1648 1648 self.PLOT_CODE = MPHASE_CODE
1649 1649
1650 1650 self.FTP_WEI = None
1651 1651 self.EXP_CODE = None
1652 1652 self.SUB_EXP_CODE = None
1653 1653 self.PLOT_POS = None
1654 1654
1655 1655
1656 1656 self.filename_phase = None
1657 1657
1658 1658 self.figfile = None
1659 1659
1660 1660 def getSubplots(self):
1661 1661
1662 1662 ncol = 1
1663 1663 nrow = 1
1664 1664
1665 1665 return nrow, ncol
1666 1666
1667 1667 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1668 1668
1669 1669 self.__showprofile = showprofile
1670 1670 self.nplots = nplots
1671 1671
1672 1672 ncolspan = 7
1673 1673 colspan = 6
1674 1674 self.__nsubplots = 2
1675 1675
1676 1676 self.createFigure(id = id,
1677 1677 wintitle = wintitle,
1678 1678 widthplot = self.WIDTH+self.WIDTHPROF,
1679 1679 heightplot = self.HEIGHT+self.HEIGHTPROF,
1680 1680 show=show)
1681 1681
1682 1682 nrow, ncol = self.getSubplots()
1683 1683
1684 1684 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1685 1685
1686 1686
1687 1687 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1688 1688 xmin=None, xmax=None, ymin=None, ymax=None,
1689 1689 timerange=None,
1690 1690 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1691 1691 server=None, folder=None, username=None, password=None,
1692 1692 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1693 1693
1694 1694
1695 1695 tmin = None
1696 1696 tmax = None
1697 1697 x = dataOut.getTimeRange1(dataOut.outputInterval)
1698 1698 y = dataOut.getHeiRange()
1699 1699
1700 1700
1701 1701 #thisDatetime = dataOut.datatime
1702 1702 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
1703 1703 title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1704 1704 xlabel = "Local Time"
1705 1705 ylabel = "Phase"
1706 1706
1707 1707
1708 1708 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1709 1709 phase_beacon = dataOut.data_output
1710 1710 update_figfile = False
1711 1711
1712 1712 if not self.isConfig:
1713 1713
1714 1714 self.nplots = phase_beacon.size
1715 1715
1716 1716 self.setup(id=id,
1717 1717 nplots=self.nplots,
1718 1718 wintitle=wintitle,
1719 1719 showprofile=showprofile,
1720 1720 show=show)
1721 1721
1722 1722 if timerange is not None:
1723 1723 self.timerange = timerange
1724 1724
1725 1725 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1726 1726
1727 1727 if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0
1728 1728 if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0
1729 1729
1730 1730 self.FTP_WEI = ftp_wei
1731 1731 self.EXP_CODE = exp_code
1732 1732 self.SUB_EXP_CODE = sub_exp_code
1733 1733 self.PLOT_POS = plot_pos
1734 1734
1735 1735 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1736 1736 self.isConfig = True
1737 1737 self.figfile = figfile
1738 1738 self.xdata = numpy.array([])
1739 1739 self.ydata = numpy.array([])
1740 1740
1741 1741 #open file beacon phase
1742 1742 path = '%s%03d' %(self.PREFIX, self.id)
1743 1743 beacon_file = os.path.join(path,'%s.txt'%self.name)
1744 1744 self.filename_phase = os.path.join(figpath,beacon_file)
1745 1745 update_figfile = True
1746 1746
1747 1747
1748 1748 #store data beacon phase
1749 1749 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1750 1750
1751 1751 self.setWinTitle(title)
1752 1752
1753 1753
1754 1754 title = "Phase Offset %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1755 1755
1756 1756 legendlabels = ["phase %d"%(chan) for chan in numpy.arange(self.nplots)]
1757 1757
1758 1758 axes = self.axesList[0]
1759 1759
1760 1760 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1761 1761
1762 1762 if len(self.ydata)==0:
1763 1763 self.ydata = phase_beacon.reshape(-1,1)
1764 1764 else:
1765 1765 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1766 1766
1767 1767
1768 1768 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1769 1769 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1770 1770 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1771 1771 XAxisAsTime=True, grid='both'
1772 1772 )
1773 1773
1774 1774 self.draw()
1775 1775
1776 1776 self.save(figpath=figpath,
1777 1777 figfile=figfile,
1778 1778 save=save,
1779 1779 ftp=ftp,
1780 1780 wr_period=wr_period,
1781 1781 thisDatetime=thisDatetime,
1782 1782 update_figfile=update_figfile)
1783 1783
1784 1784 if dataOut.ltctime + dataOut.outputInterval >= self.xmax:
1785 1785 self.counter_imagwr = wr_period
1786 1786 self.isConfig = False
1787 1787 update_figfile = True
1788 1788
1789 1789
1790 1790
1791 1791 class NSMeteorDetection1Plot(Figure):
1792 1792
1793 1793 isConfig = None
1794 1794 __nsubplots = None
1795 1795
1796 1796 WIDTHPROF = None
1797 1797 HEIGHTPROF = None
1798 1798 PREFIX = 'nsm'
1799 1799
1800 1800 zminList = None
1801 1801 zmaxList = None
1802 1802 cmapList = None
1803 1803 titleList = None
1804 1804 nPairs = None
1805 1805 nChannels = None
1806 1806 nParam = None
1807 1807
1808 1808 def __init__(self, **kwargs):
1809 1809 Figure.__init__(self, **kwargs)
1810 1810 self.isConfig = False
1811 1811 self.__nsubplots = 1
1812 1812
1813 1813 self.WIDTH = 750
1814 1814 self.HEIGHT = 250
1815 1815 self.WIDTHPROF = 120
1816 1816 self.HEIGHTPROF = 0
1817 1817 self.counter_imagwr = 0
1818 1818
1819 1819 self.PLOT_CODE = SPEC_CODE
1820 1820
1821 1821 self.FTP_WEI = None
1822 1822 self.EXP_CODE = None
1823 1823 self.SUB_EXP_CODE = None
1824 1824 self.PLOT_POS = None
1825 1825
1826 1826 self.__xfilter_ena = False
1827 1827 self.__yfilter_ena = False
1828 1828
1829 1829 def getSubplots(self):
1830 1830
1831 1831 ncol = 3
1832 1832 nrow = int(numpy.ceil(self.nplots/3.0))
1833 1833
1834 1834 return nrow, ncol
1835 1835
1836 1836 def setup(self, id, nplots, wintitle, show=True):
1837 1837
1838 1838 self.nplots = nplots
1839 1839
1840 1840 ncolspan = 1
1841 1841 colspan = 1
1842 1842
1843 1843 self.createFigure(id = id,
1844 1844 wintitle = wintitle,
1845 1845 widthplot = self.WIDTH + self.WIDTHPROF,
1846 1846 heightplot = self.HEIGHT + self.HEIGHTPROF,
1847 1847 show=show)
1848 1848
1849 1849 nrow, ncol = self.getSubplots()
1850 1850
1851 1851 counter = 0
1852 1852 for y in range(nrow):
1853 1853 for x in range(ncol):
1854 1854
1855 1855 if counter >= self.nplots:
1856 1856 break
1857 1857
1858 1858 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1859 1859
1860 1860 counter += 1
1861 1861
1862 1862 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
1863 1863 xmin=None, xmax=None, ymin=None, ymax=None, SNRmin=None, SNRmax=None,
1864 1864 vmin=None, vmax=None, wmin=None, wmax=None, mode = 'SA',
1865 1865 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1866 1866 server=None, folder=None, username=None, password=None,
1867 1867 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
1868 1868 xaxis="frequency"):
1869 1869
1870 1870 """
1871 1871
1872 1872 Input:
1873 1873 dataOut :
1874 1874 id :
1875 1875 wintitle :
1876 1876 channelList :
1877 1877 showProfile :
1878 1878 xmin : None,
1879 1879 xmax : None,
1880 1880 ymin : None,
1881 1881 ymax : None,
1882 1882 zmin : None,
1883 1883 zmax : None
1884 1884 """
1885 1885 #SEPARAR EN DOS PLOTS
1886 1886 nParam = dataOut.data_param.shape[1] - 3
1887 1887
1888 1888 utctime = dataOut.data_param[0,0]
1889 1889 tmet = dataOut.data_param[:,1].astype(int)
1890 1890 hmet = dataOut.data_param[:,2].astype(int)
1891 1891
1892 1892 x = dataOut.abscissaList
1893 1893 y = dataOut.heightList
1894 1894
1895 1895 z = numpy.zeros((nParam, y.size, x.size - 1))
1896 1896 z[:,:] = numpy.nan
1897 1897 z[:,hmet,tmet] = dataOut.data_param[:,3:].T
1898 1898 z[0,:,:] = 10*numpy.log10(z[0,:,:])
1899 1899
1900 1900 xlabel = "Time (s)"
1901 1901 ylabel = "Range (km)"
1902 1902
1903 1903 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
1904 1904
1905 1905 if not self.isConfig:
1906 1906
1907 1907 nplots = nParam
1908 1908
1909 1909 self.setup(id=id,
1910 1910 nplots=nplots,
1911 1911 wintitle=wintitle,
1912 1912 show=show)
1913 1913
1914 1914 if xmin is None: xmin = numpy.nanmin(x)
1915 1915 if xmax is None: xmax = numpy.nanmax(x)
1916 1916 if ymin is None: ymin = numpy.nanmin(y)
1917 1917 if ymax is None: ymax = numpy.nanmax(y)
1918 1918 if SNRmin is None: SNRmin = numpy.nanmin(z[0,:])
1919 1919 if SNRmax is None: SNRmax = numpy.nanmax(z[0,:])
1920 1920 if vmax is None: vmax = numpy.nanmax(numpy.abs(z[1,:]))
1921 1921 if vmin is None: vmin = -vmax
1922 1922 if wmin is None: wmin = 0
1923 1923 if wmax is None: wmax = 50
1924 1924
1925 1925 pairsList = dataOut.groupList
1926 1926 self.nPairs = len(dataOut.groupList)
1927 1927
1928 1928 zminList = [SNRmin, vmin, cmin] + [pmin]*self.nPairs
1929 1929 zmaxList = [SNRmax, vmax, cmax] + [pmax]*self.nPairs
1930 1930 titleList = ["SNR","Radial Velocity","Coherence"]
1931 1931 cmapList = ["jet","RdBu_r","jet"]
1932 1932
1933 1933 for i in range(self.nPairs):
1934 1934 strAux1 = "Phase Difference "+ str(pairsList[i][0]) + str(pairsList[i][1])
1935 1935 titleList = titleList + [strAux1]
1936 1936 cmapList = cmapList + ["RdBu_r"]
1937 1937
1938 1938 self.zminList = zminList
1939 1939 self.zmaxList = zmaxList
1940 1940 self.cmapList = cmapList
1941 1941 self.titleList = titleList
1942 1942
1943 1943 self.FTP_WEI = ftp_wei
1944 1944 self.EXP_CODE = exp_code
1945 1945 self.SUB_EXP_CODE = sub_exp_code
1946 1946 self.PLOT_POS = plot_pos
1947 1947
1948 1948 self.isConfig = True
1949 1949
1950 1950 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
1951 1951
1952 1952 for i in range(nParam):
1953 1953 title = self.titleList[i] + ": " +str_datetime
1954 1954 axes = self.axesList[i]
1955 1955 axes.pcolor(x, y, z[i,:].T,
1956 1956 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=self.zminList[i], zmax=self.zmaxList[i],
1957 1957 xlabel=xlabel, ylabel=ylabel, title=title, colormap=self.cmapList[i],ticksize=9, cblabel='')
1958 1958 self.draw()
1959 1959
1960 1960 if figfile == None:
1961 1961 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
1962 1962 name = str_datetime
1963 1963 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
1964 1964 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
1965 1965 figfile = self.getFilename(name)
1966 1966
1967 1967 self.save(figpath=figpath,
1968 1968 figfile=figfile,
1969 1969 save=save,
1970 1970 ftp=ftp,
1971 1971 wr_period=wr_period,
1972 1972 thisDatetime=thisDatetime)
1973 1973
1974 1974
1975 1975 class NSMeteorDetection2Plot(Figure):
1976 1976
1977 1977 isConfig = None
1978 1978 __nsubplots = None
1979 1979
1980 1980 WIDTHPROF = None
1981 1981 HEIGHTPROF = None
1982 1982 PREFIX = 'nsm'
1983 1983
1984 1984 zminList = None
1985 1985 zmaxList = None
1986 1986 cmapList = None
1987 1987 titleList = None
1988 1988 nPairs = None
1989 1989 nChannels = None
1990 1990 nParam = None
1991 1991
1992 1992 def __init__(self, **kwargs):
1993 1993 Figure.__init__(self, **kwargs)
1994 1994 self.isConfig = False
1995 1995 self.__nsubplots = 1
1996 1996
1997 1997 self.WIDTH = 750
1998 1998 self.HEIGHT = 250
1999 1999 self.WIDTHPROF = 120
2000 2000 self.HEIGHTPROF = 0
2001 2001 self.counter_imagwr = 0
2002 2002
2003 2003 self.PLOT_CODE = SPEC_CODE
2004 2004
2005 2005 self.FTP_WEI = None
2006 2006 self.EXP_CODE = None
2007 2007 self.SUB_EXP_CODE = None
2008 2008 self.PLOT_POS = None
2009 2009
2010 2010 self.__xfilter_ena = False
2011 2011 self.__yfilter_ena = False
2012 2012
2013 2013 def getSubplots(self):
2014 2014
2015 2015 ncol = 3
2016 2016 nrow = int(numpy.ceil(self.nplots/3.0))
2017 2017
2018 2018 return nrow, ncol
2019 2019
2020 2020 def setup(self, id, nplots, wintitle, show=True):
2021 2021
2022 2022 self.nplots = nplots
2023 2023
2024 2024 ncolspan = 1
2025 2025 colspan = 1
2026 2026
2027 2027 self.createFigure(id = id,
2028 2028 wintitle = wintitle,
2029 2029 widthplot = self.WIDTH + self.WIDTHPROF,
2030 2030 heightplot = self.HEIGHT + self.HEIGHTPROF,
2031 2031 show=show)
2032 2032
2033 2033 nrow, ncol = self.getSubplots()
2034 2034
2035 2035 counter = 0
2036 2036 for y in range(nrow):
2037 2037 for x in range(ncol):
2038 2038
2039 2039 if counter >= self.nplots:
2040 2040 break
2041 2041
2042 2042 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
2043 2043
2044 2044 counter += 1
2045 2045
2046 2046 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
2047 2047 xmin=None, xmax=None, ymin=None, ymax=None, SNRmin=None, SNRmax=None,
2048 2048 vmin=None, vmax=None, wmin=None, wmax=None, mode = 'SA',
2049 2049 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
2050 2050 server=None, folder=None, username=None, password=None,
2051 2051 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
2052 2052 xaxis="frequency"):
2053 2053
2054 2054 """
2055 2055
2056 2056 Input:
2057 2057 dataOut :
2058 2058 id :
2059 2059 wintitle :
2060 2060 channelList :
2061 2061 showProfile :
2062 2062 xmin : None,
2063 2063 xmax : None,
2064 2064 ymin : None,
2065 2065 ymax : None,
2066 2066 zmin : None,
2067 2067 zmax : None
2068 2068 """
2069 2069 #Rebuild matrix
2070 2070 utctime = dataOut.data_param[0,0]
2071 2071 cmet = dataOut.data_param[:,1].astype(int)
2072 2072 tmet = dataOut.data_param[:,2].astype(int)
2073 2073 hmet = dataOut.data_param[:,3].astype(int)
2074 2074
2075 2075 nParam = 3
2076 2076 nChan = len(dataOut.groupList)
2077 2077 x = dataOut.abscissaList
2078 2078 y = dataOut.heightList
2079 2079
2080 2080 z = numpy.full((nChan, nParam, y.size, x.size - 1),numpy.nan)
2081 2081 z[cmet,:,hmet,tmet] = dataOut.data_param[:,4:]
2082 2082 z[:,0,:,:] = 10*numpy.log10(z[:,0,:,:]) #logarithmic scale
2083 2083 z = numpy.reshape(z, (nChan*nParam, y.size, x.size-1))
2084 2084
2085 2085 xlabel = "Time (s)"
2086 2086 ylabel = "Range (km)"
2087 2087
2088 2088 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
2089 2089
2090 2090 if not self.isConfig:
2091 2091
2092 2092 nplots = nParam*nChan
2093 2093
2094 2094 self.setup(id=id,
2095 2095 nplots=nplots,
2096 2096 wintitle=wintitle,
2097 2097 show=show)
2098 2098
2099 2099 if xmin is None: xmin = numpy.nanmin(x)
2100 2100 if xmax is None: xmax = numpy.nanmax(x)
2101 2101 if ymin is None: ymin = numpy.nanmin(y)
2102 2102 if ymax is None: ymax = numpy.nanmax(y)
2103 2103 if SNRmin is None: SNRmin = numpy.nanmin(z[0,:])
2104 2104 if SNRmax is None: SNRmax = numpy.nanmax(z[0,:])
2105 2105 if vmax is None: vmax = numpy.nanmax(numpy.abs(z[1,:]))
2106 2106 if vmin is None: vmin = -vmax
2107 2107 if wmin is None: wmin = 0
2108 2108 if wmax is None: wmax = 50
2109 2109
2110 2110 self.nChannels = nChan
2111 2111
2112 2112 zminList = []
2113 2113 zmaxList = []
2114 2114 titleList = []
2115 2115 cmapList = []
2116 2116 for i in range(self.nChannels):
2117 2117 strAux1 = "SNR Channel "+ str(i)
2118 2118 strAux2 = "Radial Velocity Channel "+ str(i)
2119 2119 strAux3 = "Spectral Width Channel "+ str(i)
2120 2120
2121 2121 titleList = titleList + [strAux1,strAux2,strAux3]
2122 2122 cmapList = cmapList + ["jet","RdBu_r","jet"]
2123 2123 zminList = zminList + [SNRmin,vmin,wmin]
2124 2124 zmaxList = zmaxList + [SNRmax,vmax,wmax]
2125 2125
2126 2126 self.zminList = zminList
2127 2127 self.zmaxList = zmaxList
2128 2128 self.cmapList = cmapList
2129 2129 self.titleList = titleList
2130 2130
2131 2131 self.FTP_WEI = ftp_wei
2132 2132 self.EXP_CODE = exp_code
2133 2133 self.SUB_EXP_CODE = sub_exp_code
2134 2134 self.PLOT_POS = plot_pos
2135 2135
2136 2136 self.isConfig = True
2137 2137
2138 2138 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
2139 2139
2140 2140 for i in range(self.nplots):
2141 2141 title = self.titleList[i] + ": " +str_datetime
2142 2142 axes = self.axesList[i]
2143 2143 axes.pcolor(x, y, z[i,:].T,
2144 2144 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=self.zminList[i], zmax=self.zmaxList[i],
2145 2145 xlabel=xlabel, ylabel=ylabel, title=title, colormap=self.cmapList[i],ticksize=9, cblabel='')
2146 2146 self.draw()
2147 2147
2148 2148 if figfile == None:
2149 2149 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
2150 2150 name = str_datetime
2151 2151 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
2152 2152 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
2153 2153 figfile = self.getFilename(name)
2154 2154
2155 2155 self.save(figpath=figpath,
2156 2156 figfile=figfile,
2157 2157 save=save,
2158 2158 ftp=ftp,
2159 2159 wr_period=wr_period,
2160 2160 thisDatetime=thisDatetime)
Show file before | Show file after | Hide comments
@@ -1,1611 +1,1605
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9
10 10 import matplotlib.pyplot as plt
11 11
12 12 from figure import Figure, isRealtime, isTimeInHourRange
13 13 from plotting_codes import *
14 14 from matplotlib.pyplot import savefig
15 15
16 16 class SpectraPlot(Figure):
17 17
18 18 isConfig = None
19 19 __nsubplots = None
20 20
21 21 WIDTHPROF = None
22 22 HEIGHTPROF = None
23 23 PREFIX = 'spc'
24 24
25 25 def __init__(self, **kwargs):
26 26 Figure.__init__(self, **kwargs)
27 27 self.isConfig = False
28 28 self.__nsubplots = 1
29 29
30 30 self.WIDTH = 300
31 31 self.HEIGHT = 300
32 32 self.WIDTHPROF = 120
33 33 self.HEIGHTPROF = 0
34 34 self.counter_imagwr = 0
35 35
36 36 self.PLOT_CODE = SPEC_CODE
37 37
38 38 self.FTP_WEI = None
39 39 self.EXP_CODE = None
40 40 self.SUB_EXP_CODE = None
41 41 self.PLOT_POS = None
42 42
43 43 self.__xfilter_ena = False
44 44 self.__yfilter_ena = False
45 45
46 46 self.indice=1
47 47
48 48 def getSubplots(self):
49 49
50 50 ncol = int(numpy.sqrt(self.nplots)+0.9)
51 51 nrow = int(self.nplots*1./ncol + 0.9)
52 52
53 53 return nrow, ncol
54 54
55 55 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
56 56
57 57 self.__showprofile = showprofile
58 58 self.nplots = nplots
59 59
60 60 ncolspan = 1
61 61 colspan = 1
62 62 if showprofile:
63 63 ncolspan = 3
64 64 colspan = 2
65 65 self.__nsubplots = 2
66 66
67 67 self.createFigure(id = id,
68 68 wintitle = wintitle,
69 69 widthplot = self.WIDTH + self.WIDTHPROF,
70 70 heightplot = self.HEIGHT + self.HEIGHTPROF,
71 71 show=show)
72 72
73 73 nrow, ncol = self.getSubplots()
74 74
75 75 counter = 0
76 76 for y in range(nrow):
77 77 for x in range(ncol):
78 78
79 79 if counter >= self.nplots:
80 80 break
81 81
82 82 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
83 83
84 84 if showprofile:
85 85 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
86 86
87 87 counter += 1
88 88
89 89 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
90 90 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
91 91 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
92 92 server=None, folder=None, username=None, password=None,
93 93 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
94 94 xaxis="frequency", colormap='jet', normFactor=None):
95 95
96 96 """
97 97
98 98 Input:
99 99 dataOut :
100 100 id :
101 101 wintitle :
102 102 channelList :
103 103 showProfile :
104 104 xmin : None,
105 105 xmax : None,
106 106 ymin : None,
107 107 ymax : None,
108 108 zmin : None,
109 109 zmax : None
110 110 """
111 111 if realtime:
112 112 if not(isRealtime(utcdatatime = dataOut.utctime)):
113 113 print 'Skipping this plot function'
114 114 return
115 115
116 116 if channelList == None:
117 117 channelIndexList = dataOut.channelIndexList
118 118 else:
119 119 channelIndexList = []
120 120 for channel in channelList:
121 121 if channel not in dataOut.channelList:
122 122 raise ValueError, "Channel %d is not in dataOut.channelList" %channel
123 123 channelIndexList.append(dataOut.channelList.index(channel))
124 124
125 125 if normFactor is None:
126 126 factor = dataOut.normFactor
127 127 else:
128 128 factor = normFactor
129 129 if xaxis == "frequency":
130 130 x = dataOut.getFreqRange(1)/1000.
131 print 'FRECUENCIA MAXIMA', numpy.amax(x)
132 asfasfasdfaf
133 print '#######################################################'
134 print 'xlen', len(x)
135 print x
136 print '#######################################################'
137 131 xlabel = "Frequency (kHz)"
138 132
139 133 elif xaxis == "time":
140 134 x = dataOut.getAcfRange(1)
141 135 xlabel = "Time (ms)"
142 136
143 137 else:
144 138 x = dataOut.getVelRange(1)
145 139 xlabel = "Velocity (m/s)"
146 140
147 141 ylabel = "Range (Km)"
148 142
149 143 y = dataOut.getHeiRange()
150 144
151 145 z = dataOut.data_spc/factor
152 146 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
153 147 zdB = 10*numpy.log10(z)
154 148
155 149 avg = numpy.average(z, axis=1)
156 150 avgdB = 10*numpy.log10(avg)
157 151
158 152 noise = dataOut.getNoise()/factor
159 153 noisedB = 10*numpy.log10(noise)
160 154
161 155 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
162 156 title = wintitle + " Spectra"
163 157
164 158
165 159
166 160 print 'len de X',len(x), numpy.shape(x), 'len de spc line',len(dataOut.data_spc[1,:,15]), numpy.shape(dataOut.data_spc)
167 161 print 'Altura:', y[0], y[1], y[13], y[14], y[10]
168 162 #a=z[1,:,15]
169 163
170 164 # fig = plt.figure(10+self.indice)
171 165 # plt.plot( x[0:128], zdB[0,:,10] )
172 166 # plt.axis([-12, 12, 15, 50])
173 167 # plt.title(" %s" %( '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))) )
174 168 # plt.ylabel('Intensidad [dB]')
175 169 # plt.xlabel('Velocidad [m/s]')
176 170 # fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
177 171 #
178 172 # plt.show()
179 173 #
180 174 # self.indice=self.indice+1
181 175
182 176
183 177
184 178 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
185 179 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
186 180
187 181 if not self.isConfig:
188 182
189 183 nplots = len(channelIndexList)
190 184
191 185 self.setup(id=id,
192 186 nplots=nplots,
193 187 wintitle=wintitle,
194 188 showprofile=showprofile,
195 189 show=show)
196 190
197 191 if xmin == None: xmin = numpy.nanmin(x)
198 192 if xmax == None: xmax = numpy.nanmax(x)
199 193 if ymin == None: ymin = numpy.nanmin(y)
200 194 if ymax == None: ymax = numpy.nanmax(y)
201 195 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
202 196 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
203 197
204 198 self.FTP_WEI = ftp_wei
205 199 self.EXP_CODE = exp_code
206 200 self.SUB_EXP_CODE = sub_exp_code
207 201 self.PLOT_POS = plot_pos
208 202
209 203 self.isConfig = True
210 204
211 205 self.setWinTitle(title)
212 206
213 207 for i in range(self.nplots):
214 208 index = channelIndexList[i]
215 209 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
216 210 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
217 211 if len(dataOut.beam.codeList) != 0:
218 212 title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime)
219 213
220 214 axes = self.axesList[i*self.__nsubplots]
221 215 axes.pcolor(x, y, zdB[index,:,:],
222 216 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
223 217 xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap,
224 218 ticksize=9, cblabel='')
225 219
226 220 if self.__showprofile:
227 221 axes = self.axesList[i*self.__nsubplots +1]
228 222 axes.pline(avgdB[index,:], y,
229 223 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
230 224 xlabel='dB', ylabel='', title='',
231 225 ytick_visible=False,
232 226 grid='x')
233 227
234 228 noiseline = numpy.repeat(noisedB[index], len(y))
235 229 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
236 230
237 231 self.draw()
238 232
239 233 if figfile == None:
240 234 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
241 235 name = str_datetime
242 236 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
243 237 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
244 238 figfile = self.getFilename(name)
245 239
246 240 self.save(figpath=figpath,
247 241 figfile=figfile,
248 242 save=save,
249 243 ftp=ftp,
250 244 wr_period=wr_period,
251 245 thisDatetime=thisDatetime)
252 246
253 247
254 248 class CrossSpectraPlot(Figure):
255 249
256 250 isConfig = None
257 251 __nsubplots = None
258 252
259 253 WIDTH = None
260 254 HEIGHT = None
261 255 WIDTHPROF = None
262 256 HEIGHTPROF = None
263 257 PREFIX = 'cspc'
264 258
265 259 def __init__(self, **kwargs):
266 260 Figure.__init__(self, **kwargs)
267 261 self.isConfig = False
268 262 self.__nsubplots = 4
269 263 self.counter_imagwr = 0
270 264 self.WIDTH = 250
271 265 self.HEIGHT = 250
272 266 self.WIDTHPROF = 0
273 267 self.HEIGHTPROF = 0
274 268
275 269 self.PLOT_CODE = CROSS_CODE
276 270 self.FTP_WEI = None
277 271 self.EXP_CODE = None
278 272 self.SUB_EXP_CODE = None
279 273 self.PLOT_POS = None
280 274
281 275 self.indice=0
282 276
283 277 def getSubplots(self):
284 278
285 279 ncol = 4
286 280 nrow = self.nplots
287 281
288 282 return nrow, ncol
289 283
290 284 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
291 285
292 286 self.__showprofile = showprofile
293 287 self.nplots = nplots
294 288
295 289 ncolspan = 1
296 290 colspan = 1
297 291
298 292 self.createFigure(id = id,
299 293 wintitle = wintitle,
300 294 widthplot = self.WIDTH + self.WIDTHPROF,
301 295 heightplot = self.HEIGHT + self.HEIGHTPROF,
302 296 show=True)
303 297
304 298 nrow, ncol = self.getSubplots()
305 299
306 300 counter = 0
307 301 for y in range(nrow):
308 302 for x in range(ncol):
309 303 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
310 304
311 305 counter += 1
312 306
313 307 def run(self, dataOut, id, wintitle="", pairsList=None,
314 308 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
315 309 coh_min=None, coh_max=None, phase_min=None, phase_max=None,
316 310 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
317 311 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
318 312 server=None, folder=None, username=None, password=None,
319 313 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None,
320 314 xaxis='frequency'):
321 315
322 316 """
323 317
324 318 Input:
325 319 dataOut :
326 320 id :
327 321 wintitle :
328 322 channelList :
329 323 showProfile :
330 324 xmin : None,
331 325 xmax : None,
332 326 ymin : None,
333 327 ymax : None,
334 328 zmin : None,
335 329 zmax : None
336 330 """
337 331
338 332 if pairsList == None:
339 333 pairsIndexList = dataOut.pairsIndexList
340 334 else:
341 335 pairsIndexList = []
342 336 for pair in pairsList:
343 337 if pair not in dataOut.pairsList:
344 338 raise ValueError, "Pair %s is not in dataOut.pairsList" %str(pair)
345 339 pairsIndexList.append(dataOut.pairsList.index(pair))
346 340
347 341 if not pairsIndexList:
348 342 return
349 343
350 344 if len(pairsIndexList) > 4:
351 345 pairsIndexList = pairsIndexList[0:4]
352 346
353 347 if normFactor is None:
354 348 factor = dataOut.normFactor
355 349 else:
356 350 factor = normFactor
357 351 x = dataOut.getVelRange(1)
358 352 y = dataOut.getHeiRange()
359 353 z = dataOut.data_spc[:,:,:]/factor
360 354 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
361 355
362 356 noise = dataOut.noise/factor
363 357
364 358 zdB = 10*numpy.log10(z)
365 359 noisedB = 10*numpy.log10(noise)
366 360
367 361 if coh_min == None:
368 362 coh_min = 0.0
369 363 if coh_max == None:
370 364 coh_max = 1.0
371 365
372 366 if phase_min == None:
373 367 phase_min = -180
374 368 if phase_max == None:
375 369 phase_max = 180
376 370
377 371 #thisDatetime = dataOut.datatime
378 372 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
379 373 title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
380 374 # xlabel = "Velocity (m/s)"
381 375 ylabel = "Range (Km)"
382 376
383 377 if xaxis == "frequency":
384 378 x = dataOut.getFreqRange(1)/1000.
385 379 xlabel = "Frequency (kHz)"
386 380
387 381 elif xaxis == "time":
388 382 x = dataOut.getAcfRange(1)
389 383 xlabel = "Time (ms)"
390 384
391 385 else:
392 386 x = dataOut.getVelRange(1)
393 387 xlabel = "Velocity (m/s)"
394 388
395 389 if not self.isConfig:
396 390
397 391 nplots = len(pairsIndexList)
398 392
399 393 self.setup(id=id,
400 394 nplots=nplots,
401 395 wintitle=wintitle,
402 396 showprofile=False,
403 397 show=show)
404 398
405 399 avg = numpy.abs(numpy.average(z, axis=1))
406 400 avgdB = 10*numpy.log10(avg)
407 401
408 402 if xmin == None: xmin = numpy.nanmin(x)
409 403 if xmax == None: xmax = numpy.nanmax(x)
410 404 if ymin == None: ymin = numpy.nanmin(y)
411 405 if ymax == None: ymax = numpy.nanmax(y)
412 406 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
413 407 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
414 408
415 409 self.FTP_WEI = ftp_wei
416 410 self.EXP_CODE = exp_code
417 411 self.SUB_EXP_CODE = sub_exp_code
418 412 self.PLOT_POS = plot_pos
419 413
420 414 self.isConfig = True
421 415
422 416 self.setWinTitle(title)
423 417
424 418
425 419 for i in range(self.nplots):
426 420 pair = dataOut.pairsList[pairsIndexList[i]]
427 421
428 422 chan_index0 = dataOut.channelList.index(pair[0])
429 423 chan_index1 = dataOut.channelList.index(pair[1])
430 424
431 425 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
432 426 title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[chan_index0], str_datetime)
433 427 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index0,:,:]/factor)
434 428 axes0 = self.axesList[i*self.__nsubplots]
435 429 axes0.pcolor(x, y, zdB,
436 430 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
437 431 xlabel=xlabel, ylabel=ylabel, title=title,
438 432 ticksize=9, colormap=power_cmap, cblabel='')
439 433
440 434 title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[chan_index1], str_datetime)
441 435 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index1,:,:]/factor)
442 436 axes0 = self.axesList[i*self.__nsubplots+1]
443 437 axes0.pcolor(x, y, zdB,
444 438 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
445 439 xlabel=xlabel, ylabel=ylabel, title=title,
446 440 ticksize=9, colormap=power_cmap, cblabel='')
447 441
448 442 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:] / numpy.sqrt( dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:] )
449 443 coherence = numpy.abs(coherenceComplex)
450 444 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
451 445 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
452 446
453 447
454 448
455 449
456 450 # #print 'FASE', numpy.shape(phase), y[25]
457 451 # fig = plt.figure(10+self.indice)
458 452 # #plt.plot( x[0:256],coherence[:,25] )
459 453 # cohAv = numpy.average(coherence,1)
460 454 #
461 455 # plt.plot( x[0:256],cohAv )
462 456 # #plt.axis([-12, 12, 15, 50])
463 457 # plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
464 458 # plt.ylabel('Desfase [grados]')
465 459 # plt.xlabel('Velocidad [m/s]')
466 460 # fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
467 461 #
468 462 # plt.show()
469 463 # self.indice=self.indice+1
470 464
471 465
472 466 # print 'FASE', numpy.shape(phase), y[25]
473 467 # fig = plt.figure(10+self.indice)
474 468 # plt.plot( x[0:256],phase[:,25] )
475 469 # #plt.axis([-12, 12, 15, 50])
476 470 # plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
477 471 # plt.ylabel('Desfase [grados]')
478 472 # plt.xlabel('Velocidad [m/s]')
479 473 # fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
480 474 #
481 475 # plt.show()
482 476 # self.indice=self.indice+1
483 477
484 478
485 479
486 480
487 481 title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1])
488 482 axes0 = self.axesList[i*self.__nsubplots+2]
489 483 axes0.pcolor(x, y, coherence,
490 484 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=coh_min, zmax=coh_max,
491 485 xlabel=xlabel, ylabel=ylabel, title=title,
492 486 ticksize=9, colormap=coherence_cmap, cblabel='')
493 487
494 488 title = "Phase Ch%d * Ch%d" %(pair[0], pair[1])
495 489 axes0 = self.axesList[i*self.__nsubplots+3]
496 490 axes0.pcolor(x, y, phase,
497 491 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
498 492 xlabel=xlabel, ylabel=ylabel, title=title,
499 493 ticksize=9, colormap=phase_cmap, cblabel='')
500 494
501 495
502 496
503 497 self.draw()
504 498
505 499 self.save(figpath=figpath,
506 500 figfile=figfile,
507 501 save=save,
508 502 ftp=ftp,
509 503 wr_period=wr_period,
510 504 thisDatetime=thisDatetime)
511 505
512 506
513 507 class RTIPlot(Figure):
514 508
515 509 __isConfig = None
516 510 __nsubplots = None
517 511
518 512 WIDTHPROF = None
519 513 HEIGHTPROF = None
520 514 PREFIX = 'rti'
521 515
522 516 def __init__(self, **kwargs):
523 517
524 518 Figure.__init__(self, **kwargs)
525 519 self.timerange = None
526 520 self.isConfig = False
527 521 self.__nsubplots = 1
528 522
529 523 self.WIDTH = 800
530 524 self.HEIGHT = 250
531 525 self.WIDTHPROF = 120
532 526 self.HEIGHTPROF = 0
533 527 self.counter_imagwr = 0
534 528
535 529 self.PLOT_CODE = RTI_CODE
536 530
537 531 self.FTP_WEI = None
538 532 self.EXP_CODE = None
539 533 self.SUB_EXP_CODE = None
540 534 self.PLOT_POS = None
541 535 self.tmin = None
542 536 self.tmax = None
543 537
544 538 self.xmin = None
545 539 self.xmax = None
546 540
547 541 self.figfile = None
548 542
549 543 def getSubplots(self):
550 544
551 545 ncol = 1
552 546 nrow = self.nplots
553 547
554 548 return nrow, ncol
555 549
556 550 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
557 551
558 552 self.__showprofile = showprofile
559 553 self.nplots = nplots
560 554
561 555 ncolspan = 1
562 556 colspan = 1
563 557 if showprofile:
564 558 ncolspan = 7
565 559 colspan = 6
566 560 self.__nsubplots = 2
567 561
568 562 self.createFigure(id = id,
569 563 wintitle = wintitle,
570 564 widthplot = self.WIDTH + self.WIDTHPROF,
571 565 heightplot = self.HEIGHT + self.HEIGHTPROF,
572 566 show=show)
573 567
574 568 nrow, ncol = self.getSubplots()
575 569
576 570 counter = 0
577 571 for y in range(nrow):
578 572 for x in range(ncol):
579 573
580 574 if counter >= self.nplots:
581 575 break
582 576
583 577 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
584 578
585 579 if showprofile:
586 580 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
587 581
588 582 counter += 1
589 583
590 584 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
591 585 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
592 586 timerange=None, colormap='jet',
593 587 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
594 588 server=None, folder=None, username=None, password=None,
595 589 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, HEIGHT=None):
596 590
597 591 """
598 592
599 593 Input:
600 594 dataOut :
601 595 id :
602 596 wintitle :
603 597 channelList :
604 598 showProfile :
605 599 xmin : None,
606 600 xmax : None,
607 601 ymin : None,
608 602 ymax : None,
609 603 zmin : None,
610 604 zmax : None
611 605 """
612 606
613 607 #colormap = kwargs.get('colormap', 'jet')
614 608 if HEIGHT is not None:
615 609 self.HEIGHT = HEIGHT
616 610
617 611 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
618 612 return
619 613
620 614 if channelList == None:
621 615 channelIndexList = dataOut.channelIndexList
622 616 else:
623 617 channelIndexList = []
624 618 for channel in channelList:
625 619 if channel not in dataOut.channelList:
626 620 raise ValueError, "Channel %d is not in dataOut.channelList"
627 621 channelIndexList.append(dataOut.channelList.index(channel))
628 622
629 623 if normFactor is None:
630 624 factor = dataOut.normFactor
631 625 else:
632 626 factor = normFactor
633 627
634 628 # factor = dataOut.normFactor
635 629 x = dataOut.getTimeRange()
636 630 y = dataOut.getHeiRange()
637 631
638 632 z = dataOut.data_spc/factor
639 633 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
640 634 avg = numpy.average(z, axis=1)
641 635 avgdB = 10.*numpy.log10(avg)
642 636 # avgdB = dataOut.getPower()
643 637
644 638
645 639 thisDatetime = dataOut.datatime
646 640 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
647 641 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
648 642 xlabel = ""
649 643 ylabel = "Range (Km)"
650 644
651 645 update_figfile = False
652 646
653 647 if dataOut.ltctime >= self.xmax:
654 648 self.counter_imagwr = wr_period
655 649 self.isConfig = False
656 650 update_figfile = True
657 651
658 652 if not self.isConfig:
659 653
660 654 nplots = len(channelIndexList)
661 655
662 656 self.setup(id=id,
663 657 nplots=nplots,
664 658 wintitle=wintitle,
665 659 showprofile=showprofile,
666 660 show=show)
667 661
668 662 if timerange != None:
669 663 self.timerange = timerange
670 664
671 665 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
672 666
673 667 noise = dataOut.noise/factor
674 668 noisedB = 10*numpy.log10(noise)
675 669
676 670 if ymin == None: ymin = numpy.nanmin(y)
677 671 if ymax == None: ymax = numpy.nanmax(y)
678 672 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
679 673 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
680 674
681 675 self.FTP_WEI = ftp_wei
682 676 self.EXP_CODE = exp_code
683 677 self.SUB_EXP_CODE = sub_exp_code
684 678 self.PLOT_POS = plot_pos
685 679
686 680 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
687 681 self.isConfig = True
688 682 self.figfile = figfile
689 683 update_figfile = True
690 684
691 685 self.setWinTitle(title)
692 686
693 687 for i in range(self.nplots):
694 688 index = channelIndexList[i]
695 689 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
696 690 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
697 691 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
698 692 axes = self.axesList[i*self.__nsubplots]
699 693 zdB = avgdB[index].reshape((1,-1))
700 694 axes.pcolorbuffer(x, y, zdB,
701 695 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
702 696 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
703 697 ticksize=9, cblabel='', cbsize="1%", colormap=colormap)
704 698
705 699 if self.__showprofile:
706 700 axes = self.axesList[i*self.__nsubplots +1]
707 701 axes.pline(avgdB[index], y,
708 702 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
709 703 xlabel='dB', ylabel='', title='',
710 704 ytick_visible=False,
711 705 grid='x')
712 706
713 707 self.draw()
714 708
715 709 self.save(figpath=figpath,
716 710 figfile=figfile,
717 711 save=save,
718 712 ftp=ftp,
719 713 wr_period=wr_period,
720 714 thisDatetime=thisDatetime,
721 715 update_figfile=update_figfile)
722 716
723 717 class CoherenceMap(Figure):
724 718 isConfig = None
725 719 __nsubplots = None
726 720
727 721 WIDTHPROF = None
728 722 HEIGHTPROF = None
729 723 PREFIX = 'cmap'
730 724
731 725 def __init__(self, **kwargs):
732 726 Figure.__init__(self, **kwargs)
733 727 self.timerange = 2*60*60
734 728 self.isConfig = False
735 729 self.__nsubplots = 1
736 730
737 731 self.WIDTH = 800
738 732 self.HEIGHT = 180
739 733 self.WIDTHPROF = 120
740 734 self.HEIGHTPROF = 0
741 735 self.counter_imagwr = 0
742 736
743 737 self.PLOT_CODE = COH_CODE
744 738
745 739 self.FTP_WEI = None
746 740 self.EXP_CODE = None
747 741 self.SUB_EXP_CODE = None
748 742 self.PLOT_POS = None
749 743 self.counter_imagwr = 0
750 744
751 745 self.xmin = None
752 746 self.xmax = None
753 747
754 748 def getSubplots(self):
755 749 ncol = 1
756 750 nrow = self.nplots*2
757 751
758 752 return nrow, ncol
759 753
760 754 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
761 755 self.__showprofile = showprofile
762 756 self.nplots = nplots
763 757
764 758 ncolspan = 1
765 759 colspan = 1
766 760 if showprofile:
767 761 ncolspan = 7
768 762 colspan = 6
769 763 self.__nsubplots = 2
770 764
771 765 self.createFigure(id = id,
772 766 wintitle = wintitle,
773 767 widthplot = self.WIDTH + self.WIDTHPROF,
774 768 heightplot = self.HEIGHT + self.HEIGHTPROF,
775 769 show=True)
776 770
777 771 nrow, ncol = self.getSubplots()
778 772
779 773 for y in range(nrow):
780 774 for x in range(ncol):
781 775
782 776 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
783 777
784 778 if showprofile:
785 779 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
786 780
787 781 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
788 782 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
789 783 timerange=None, phase_min=None, phase_max=None,
790 784 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
791 785 coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
792 786 server=None, folder=None, username=None, password=None,
793 787 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
794 788
795 789 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
796 790 return
797 791
798 792 if pairsList == None:
799 793 pairsIndexList = dataOut.pairsIndexList
800 794 else:
801 795 pairsIndexList = []
802 796 for pair in pairsList:
803 797 if pair not in dataOut.pairsList:
804 798 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
805 799 pairsIndexList.append(dataOut.pairsList.index(pair))
806 800
807 801 if pairsIndexList == []:
808 802 return
809 803
810 804 if len(pairsIndexList) > 4:
811 805 pairsIndexList = pairsIndexList[0:4]
812 806
813 807 if phase_min == None:
814 808 phase_min = -180
815 809 if phase_max == None:
816 810 phase_max = 180
817 811
818 812 x = dataOut.getTimeRange()
819 813 y = dataOut.getHeiRange()
820 814
821 815 thisDatetime = dataOut.datatime
822 816
823 817 title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
824 818 xlabel = ""
825 819 ylabel = "Range (Km)"
826 820 update_figfile = False
827 821
828 822 if not self.isConfig:
829 823 nplots = len(pairsIndexList)
830 824 self.setup(id=id,
831 825 nplots=nplots,
832 826 wintitle=wintitle,
833 827 showprofile=showprofile,
834 828 show=show)
835 829
836 830 if timerange != None:
837 831 self.timerange = timerange
838 832
839 833 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
840 834
841 835 if ymin == None: ymin = numpy.nanmin(y)
842 836 if ymax == None: ymax = numpy.nanmax(y)
843 837 if zmin == None: zmin = 0.
844 838 if zmax == None: zmax = 1.
845 839
846 840 self.FTP_WEI = ftp_wei
847 841 self.EXP_CODE = exp_code
848 842 self.SUB_EXP_CODE = sub_exp_code
849 843 self.PLOT_POS = plot_pos
850 844
851 845 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
852 846
853 847 self.isConfig = True
854 848 update_figfile = True
855 849
856 850 self.setWinTitle(title)
857 851
858 852 for i in range(self.nplots):
859 853
860 854 pair = dataOut.pairsList[pairsIndexList[i]]
861 855
862 856 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
863 857 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
864 858 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
865 859
866 860
867 861 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
868 862 coherence = numpy.abs(avgcoherenceComplex)
869 863
870 864 z = coherence.reshape((1,-1))
871 865
872 866 counter = 0
873 867
874 868 title = "Coherence Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
875 869 axes = self.axesList[i*self.__nsubplots*2]
876 870 axes.pcolorbuffer(x, y, z,
877 871 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
878 872 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
879 873 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
880 874
881 875 if self.__showprofile:
882 876 counter += 1
883 877 axes = self.axesList[i*self.__nsubplots*2 + counter]
884 878 axes.pline(coherence, y,
885 879 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
886 880 xlabel='', ylabel='', title='', ticksize=7,
887 881 ytick_visible=False, nxticks=5,
888 882 grid='x')
889 883
890 884 counter += 1
891 885
892 886 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
893 887
894 888 z = phase.reshape((1,-1))
895 889
896 890 title = "Phase Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
897 891 axes = self.axesList[i*self.__nsubplots*2 + counter]
898 892 axes.pcolorbuffer(x, y, z,
899 893 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
900 894 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
901 895 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
902 896
903 897 if self.__showprofile:
904 898 counter += 1
905 899 axes = self.axesList[i*self.__nsubplots*2 + counter]
906 900 axes.pline(phase, y,
907 901 xmin=phase_min, xmax=phase_max, ymin=ymin, ymax=ymax,
908 902 xlabel='', ylabel='', title='', ticksize=7,
909 903 ytick_visible=False, nxticks=4,
910 904 grid='x')
911 905
912 906 self.draw()
913 907
914 908 if dataOut.ltctime >= self.xmax:
915 909 self.counter_imagwr = wr_period
916 910 self.isConfig = False
917 911 update_figfile = True
918 912
919 913 self.save(figpath=figpath,
920 914 figfile=figfile,
921 915 save=save,
922 916 ftp=ftp,
923 917 wr_period=wr_period,
924 918 thisDatetime=thisDatetime,
925 919 update_figfile=update_figfile)
926 920
927 921 class PowerProfilePlot(Figure):
928 922
929 923 isConfig = None
930 924 __nsubplots = None
931 925
932 926 WIDTHPROF = None
933 927 HEIGHTPROF = None
934 928 PREFIX = 'spcprofile'
935 929
936 930 def __init__(self, **kwargs):
937 931 Figure.__init__(self, **kwargs)
938 932 self.isConfig = False
939 933 self.__nsubplots = 1
940 934
941 935 self.PLOT_CODE = POWER_CODE
942 936
943 937 self.WIDTH = 300
944 938 self.HEIGHT = 500
945 939 self.counter_imagwr = 0
946 940
947 941 def getSubplots(self):
948 942 ncol = 1
949 943 nrow = 1
950 944
951 945 return nrow, ncol
952 946
953 947 def setup(self, id, nplots, wintitle, show):
954 948
955 949 self.nplots = nplots
956 950
957 951 ncolspan = 1
958 952 colspan = 1
959 953
960 954 self.createFigure(id = id,
961 955 wintitle = wintitle,
962 956 widthplot = self.WIDTH,
963 957 heightplot = self.HEIGHT,
964 958 show=show)
965 959
966 960 nrow, ncol = self.getSubplots()
967 961
968 962 counter = 0
969 963 for y in range(nrow):
970 964 for x in range(ncol):
971 965 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
972 966
973 967 def run(self, dataOut, id, wintitle="", channelList=None,
974 968 xmin=None, xmax=None, ymin=None, ymax=None,
975 969 save=False, figpath='./', figfile=None, show=True,
976 970 ftp=False, wr_period=1, server=None,
977 971 folder=None, username=None, password=None):
978 972
979 973
980 974 if channelList == None:
981 975 channelIndexList = dataOut.channelIndexList
982 976 channelList = dataOut.channelList
983 977 else:
984 978 channelIndexList = []
985 979 for channel in channelList:
986 980 if channel not in dataOut.channelList:
987 981 raise ValueError, "Channel %d is not in dataOut.channelList"
988 982 channelIndexList.append(dataOut.channelList.index(channel))
989 983
990 984 factor = dataOut.normFactor
991 985
992 986 y = dataOut.getHeiRange()
993 987
994 988 #for voltage
995 989 if dataOut.type == 'Voltage':
996 990 x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
997 991 x = x.real
998 992 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
999 993
1000 994 #for spectra
1001 995 if dataOut.type == 'Spectra':
1002 996 x = dataOut.data_spc[channelIndexList,:,:]/factor
1003 997 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
1004 998 x = numpy.average(x, axis=1)
1005 999
1006 1000
1007 1001 xdB = 10*numpy.log10(x)
1008 1002
1009 1003 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1010 1004 title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
1011 1005 xlabel = "dB"
1012 1006 ylabel = "Range (Km)"
1013 1007
1014 1008 if not self.isConfig:
1015 1009
1016 1010 nplots = 1
1017 1011
1018 1012 self.setup(id=id,
1019 1013 nplots=nplots,
1020 1014 wintitle=wintitle,
1021 1015 show=show)
1022 1016
1023 1017 if ymin == None: ymin = numpy.nanmin(y)
1024 1018 if ymax == None: ymax = numpy.nanmax(y)
1025 1019 if xmin == None: xmin = numpy.nanmin(xdB)*0.9
1026 1020 if xmax == None: xmax = numpy.nanmax(xdB)*1.1
1027 1021
1028 1022 self.isConfig = True
1029 1023
1030 1024 self.setWinTitle(title)
1031 1025
1032 1026 title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1033 1027 axes = self.axesList[0]
1034 1028
1035 1029 legendlabels = ["channel %d"%x for x in channelList]
1036 1030 axes.pmultiline(xdB, y,
1037 1031 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1038 1032 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1039 1033 ytick_visible=True, nxticks=5,
1040 1034 grid='x')
1041 1035
1042 1036 self.draw()
1043 1037
1044 1038 self.save(figpath=figpath,
1045 1039 figfile=figfile,
1046 1040 save=save,
1047 1041 ftp=ftp,
1048 1042 wr_period=wr_period,
1049 1043 thisDatetime=thisDatetime)
1050 1044
1051 1045 class SpectraCutPlot(Figure):
1052 1046
1053 1047 isConfig = None
1054 1048 __nsubplots = None
1055 1049
1056 1050 WIDTHPROF = None
1057 1051 HEIGHTPROF = None
1058 1052 PREFIX = 'spc_cut'
1059 1053
1060 1054 def __init__(self, **kwargs):
1061 1055 Figure.__init__(self, **kwargs)
1062 1056 self.isConfig = False
1063 1057 self.__nsubplots = 1
1064 1058
1065 1059 self.PLOT_CODE = POWER_CODE
1066 1060
1067 1061 self.WIDTH = 700
1068 1062 self.HEIGHT = 500
1069 1063 self.counter_imagwr = 0
1070 1064
1071 1065 def getSubplots(self):
1072 1066 ncol = 1
1073 1067 nrow = 1
1074 1068
1075 1069 return nrow, ncol
1076 1070
1077 1071 def setup(self, id, nplots, wintitle, show):
1078 1072
1079 1073 self.nplots = nplots
1080 1074
1081 1075 ncolspan = 1
1082 1076 colspan = 1
1083 1077
1084 1078 self.createFigure(id = id,
1085 1079 wintitle = wintitle,
1086 1080 widthplot = self.WIDTH,
1087 1081 heightplot = self.HEIGHT,
1088 1082 show=show)
1089 1083
1090 1084 nrow, ncol = self.getSubplots()
1091 1085
1092 1086 counter = 0
1093 1087 for y in range(nrow):
1094 1088 for x in range(ncol):
1095 1089 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1096 1090
1097 1091 def run(self, dataOut, id, wintitle="", channelList=None,
1098 1092 xmin=None, xmax=None, ymin=None, ymax=None,
1099 1093 save=False, figpath='./', figfile=None, show=True,
1100 1094 ftp=False, wr_period=1, server=None,
1101 1095 folder=None, username=None, password=None,
1102 1096 xaxis="frequency"):
1103 1097
1104 1098
1105 1099 if channelList == None:
1106 1100 channelIndexList = dataOut.channelIndexList
1107 1101 channelList = dataOut.channelList
1108 1102 else:
1109 1103 channelIndexList = []
1110 1104 for channel in channelList:
1111 1105 if channel not in dataOut.channelList:
1112 1106 raise ValueError, "Channel %d is not in dataOut.channelList"
1113 1107 channelIndexList.append(dataOut.channelList.index(channel))
1114 1108
1115 1109 factor = dataOut.normFactor
1116 1110
1117 1111 y = dataOut.getHeiRange()
1118 1112
1119 1113 z = dataOut.data_spc/factor
1120 1114 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1121 1115
1122 1116 hei_index = numpy.arange(25)*3 + 20
1123 1117
1124 1118 if xaxis == "frequency":
1125 1119 x = dataOut.getFreqRange()/1000.
1126 1120 zdB = 10*numpy.log10(z[0,:,hei_index])
1127 1121 xlabel = "Frequency (kHz)"
1128 1122 ylabel = "Power (dB)"
1129 1123
1130 1124 elif xaxis == "time":
1131 1125 x = dataOut.getAcfRange()
1132 1126 zdB = z[0,:,hei_index]
1133 1127 xlabel = "Time (ms)"
1134 1128 ylabel = "ACF"
1135 1129
1136 1130 else:
1137 1131 x = dataOut.getVelRange()
1138 1132 zdB = 10*numpy.log10(z[0,:,hei_index])
1139 1133 xlabel = "Velocity (m/s)"
1140 1134 ylabel = "Power (dB)"
1141 1135
1142 1136 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1143 1137 title = wintitle + " Range Cuts %s" %(thisDatetime.strftime("%d-%b-%Y"))
1144 1138
1145 1139 if not self.isConfig:
1146 1140
1147 1141 nplots = 1
1148 1142
1149 1143 self.setup(id=id,
1150 1144 nplots=nplots,
1151 1145 wintitle=wintitle,
1152 1146 show=show)
1153 1147
1154 1148 if xmin == None: xmin = numpy.nanmin(x)*0.9
1155 1149 if xmax == None: xmax = numpy.nanmax(x)*1.1
1156 1150 if ymin == None: ymin = numpy.nanmin(zdB)
1157 1151 if ymax == None: ymax = numpy.nanmax(zdB)
1158 1152
1159 1153 self.isConfig = True
1160 1154
1161 1155 self.setWinTitle(title)
1162 1156
1163 1157 title = "Spectra Cuts: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1164 1158 axes = self.axesList[0]
1165 1159
1166 1160 legendlabels = ["Range = %dKm" %y[i] for i in hei_index]
1167 1161
1168 1162 axes.pmultilineyaxis( x, zdB,
1169 1163 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1170 1164 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1171 1165 ytick_visible=True, nxticks=5,
1172 1166 grid='x')
1173 1167
1174 1168 self.draw()
1175 1169
1176 1170 self.save(figpath=figpath,
1177 1171 figfile=figfile,
1178 1172 save=save,
1179 1173 ftp=ftp,
1180 1174 wr_period=wr_period,
1181 1175 thisDatetime=thisDatetime)
1182 1176
1183 1177 class Noise(Figure):
1184 1178
1185 1179 isConfig = None
1186 1180 __nsubplots = None
1187 1181
1188 1182 PREFIX = 'noise'
1189 1183
1190 1184
1191 1185 def __init__(self, **kwargs):
1192 1186 Figure.__init__(self, **kwargs)
1193 1187 self.timerange = 24*60*60
1194 1188 self.isConfig = False
1195 1189 self.__nsubplots = 1
1196 1190 self.counter_imagwr = 0
1197 1191 self.WIDTH = 800
1198 1192 self.HEIGHT = 400
1199 1193 self.WIDTHPROF = 120
1200 1194 self.HEIGHTPROF = 0
1201 1195 self.xdata = None
1202 1196 self.ydata = None
1203 1197
1204 1198 self.PLOT_CODE = NOISE_CODE
1205 1199
1206 1200 self.FTP_WEI = None
1207 1201 self.EXP_CODE = None
1208 1202 self.SUB_EXP_CODE = None
1209 1203 self.PLOT_POS = None
1210 1204 self.figfile = None
1211 1205
1212 1206 self.xmin = None
1213 1207 self.xmax = None
1214 1208
1215 1209 def getSubplots(self):
1216 1210
1217 1211 ncol = 1
1218 1212 nrow = 1
1219 1213
1220 1214 return nrow, ncol
1221 1215
1222 1216 def openfile(self, filename):
1223 1217 dirname = os.path.dirname(filename)
1224 1218
1225 1219 if not os.path.exists(dirname):
1226 1220 os.mkdir(dirname)
1227 1221
1228 1222 f = open(filename,'w+')
1229 1223 f.write('\n\n')
1230 1224 f.write('JICAMARCA RADIO OBSERVATORY - Noise \n')
1231 1225 f.write('DD MM YYYY HH MM SS Channel0 Channel1 Channel2 Channel3\n\n' )
1232 1226 f.close()
1233 1227
1234 1228 def save_data(self, filename_phase, data, data_datetime):
1235 1229
1236 1230 f=open(filename_phase,'a')
1237 1231
1238 1232 timetuple_data = data_datetime.timetuple()
1239 1233 day = str(timetuple_data.tm_mday)
1240 1234 month = str(timetuple_data.tm_mon)
1241 1235 year = str(timetuple_data.tm_year)
1242 1236 hour = str(timetuple_data.tm_hour)
1243 1237 minute = str(timetuple_data.tm_min)
1244 1238 second = str(timetuple_data.tm_sec)
1245 1239
1246 1240 data_msg = ''
1247 1241 for i in range(len(data)):
1248 1242 data_msg += str(data[i]) + ' '
1249 1243
1250 1244 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' ' + data_msg + '\n')
1251 1245 f.close()
1252 1246
1253 1247
1254 1248 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1255 1249
1256 1250 self.__showprofile = showprofile
1257 1251 self.nplots = nplots
1258 1252
1259 1253 ncolspan = 7
1260 1254 colspan = 6
1261 1255 self.__nsubplots = 2
1262 1256
1263 1257 self.createFigure(id = id,
1264 1258 wintitle = wintitle,
1265 1259 widthplot = self.WIDTH+self.WIDTHPROF,
1266 1260 heightplot = self.HEIGHT+self.HEIGHTPROF,
1267 1261 show=show)
1268 1262
1269 1263 nrow, ncol = self.getSubplots()
1270 1264
1271 1265 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1272 1266
1273 1267
1274 1268 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1275 1269 xmin=None, xmax=None, ymin=None, ymax=None,
1276 1270 timerange=None,
1277 1271 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1278 1272 server=None, folder=None, username=None, password=None,
1279 1273 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1280 1274
1281 1275 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1282 1276 return
1283 1277
1284 1278 if channelList == None:
1285 1279 channelIndexList = dataOut.channelIndexList
1286 1280 channelList = dataOut.channelList
1287 1281 else:
1288 1282 channelIndexList = []
1289 1283 for channel in channelList:
1290 1284 if channel not in dataOut.channelList:
1291 1285 raise ValueError, "Channel %d is not in dataOut.channelList"
1292 1286 channelIndexList.append(dataOut.channelList.index(channel))
1293 1287
1294 1288 x = dataOut.getTimeRange()
1295 1289 #y = dataOut.getHeiRange()
1296 1290 factor = dataOut.normFactor
1297 1291 noise = dataOut.noise[channelIndexList]/factor
1298 1292 noisedB = 10*numpy.log10(noise)
1299 1293
1300 1294 thisDatetime = dataOut.datatime
1301 1295
1302 1296 title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1303 1297 xlabel = ""
1304 1298 ylabel = "Intensity (dB)"
1305 1299 update_figfile = False
1306 1300
1307 1301 if not self.isConfig:
1308 1302
1309 1303 nplots = 1
1310 1304
1311 1305 self.setup(id=id,
1312 1306 nplots=nplots,
1313 1307 wintitle=wintitle,
1314 1308 showprofile=showprofile,
1315 1309 show=show)
1316 1310
1317 1311 if timerange != None:
1318 1312 self.timerange = timerange
1319 1313
1320 1314 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1321 1315
1322 1316 if ymin == None: ymin = numpy.floor(numpy.nanmin(noisedB)) - 10.0
1323 1317 if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
1324 1318
1325 1319 self.FTP_WEI = ftp_wei
1326 1320 self.EXP_CODE = exp_code
1327 1321 self.SUB_EXP_CODE = sub_exp_code
1328 1322 self.PLOT_POS = plot_pos
1329 1323
1330 1324
1331 1325 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1332 1326 self.isConfig = True
1333 1327 self.figfile = figfile
1334 1328 self.xdata = numpy.array([])
1335 1329 self.ydata = numpy.array([])
1336 1330
1337 1331 update_figfile = True
1338 1332
1339 1333 #open file beacon phase
1340 1334 path = '%s%03d' %(self.PREFIX, self.id)
1341 1335 noise_file = os.path.join(path,'%s.txt'%self.name)
1342 1336 self.filename_noise = os.path.join(figpath,noise_file)
1343 1337
1344 1338 self.setWinTitle(title)
1345 1339
1346 1340 title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1347 1341
1348 1342 legendlabels = ["channel %d"%(idchannel) for idchannel in channelList]
1349 1343 axes = self.axesList[0]
1350 1344
1351 1345 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1352 1346
1353 1347 if len(self.ydata)==0:
1354 1348 self.ydata = noisedB.reshape(-1,1)
1355 1349 else:
1356 1350 self.ydata = numpy.hstack((self.ydata, noisedB.reshape(-1,1)))
1357 1351
1358 1352
1359 1353 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1360 1354 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1361 1355 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1362 1356 XAxisAsTime=True, grid='both'
1363 1357 )
1364 1358
1365 1359 self.draw()
1366 1360
1367 1361 if dataOut.ltctime >= self.xmax:
1368 1362 self.counter_imagwr = wr_period
1369 1363 self.isConfig = False
1370 1364 update_figfile = True
1371 1365
1372 1366 self.save(figpath=figpath,
1373 1367 figfile=figfile,
1374 1368 save=save,
1375 1369 ftp=ftp,
1376 1370 wr_period=wr_period,
1377 1371 thisDatetime=thisDatetime,
1378 1372 update_figfile=update_figfile)
1379 1373
1380 1374 #store data beacon phase
1381 1375 if save:
1382 1376 self.save_data(self.filename_noise, noisedB, thisDatetime)
1383 1377
1384 1378 class BeaconPhase(Figure):
1385 1379
1386 1380 __isConfig = None
1387 1381 __nsubplots = None
1388 1382
1389 1383 PREFIX = 'beacon_phase'
1390 1384
1391 1385 def __init__(self, **kwargs):
1392 1386 Figure.__init__(self, **kwargs)
1393 1387 self.timerange = 24*60*60
1394 1388 self.isConfig = False
1395 1389 self.__nsubplots = 1
1396 1390 self.counter_imagwr = 0
1397 1391 self.WIDTH = 800
1398 1392 self.HEIGHT = 400
1399 1393 self.WIDTHPROF = 120
1400 1394 self.HEIGHTPROF = 0
1401 1395 self.xdata = None
1402 1396 self.ydata = None
1403 1397
1404 1398 self.PLOT_CODE = BEACON_CODE
1405 1399
1406 1400 self.FTP_WEI = None
1407 1401 self.EXP_CODE = None
1408 1402 self.SUB_EXP_CODE = None
1409 1403 self.PLOT_POS = None
1410 1404
1411 1405 self.filename_phase = None
1412 1406
1413 1407 self.figfile = None
1414 1408
1415 1409 self.xmin = None
1416 1410 self.xmax = None
1417 1411
1418 1412 def getSubplots(self):
1419 1413
1420 1414 ncol = 1
1421 1415 nrow = 1
1422 1416
1423 1417 return nrow, ncol
1424 1418
1425 1419 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1426 1420
1427 1421 self.__showprofile = showprofile
1428 1422 self.nplots = nplots
1429 1423
1430 1424 ncolspan = 7
1431 1425 colspan = 6
1432 1426 self.__nsubplots = 2
1433 1427
1434 1428 self.createFigure(id = id,
1435 1429 wintitle = wintitle,
1436 1430 widthplot = self.WIDTH+self.WIDTHPROF,
1437 1431 heightplot = self.HEIGHT+self.HEIGHTPROF,
1438 1432 show=show)
1439 1433
1440 1434 nrow, ncol = self.getSubplots()
1441 1435
1442 1436 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1443 1437
1444 1438 def save_phase(self, filename_phase):
1445 1439 f = open(filename_phase,'w+')
1446 1440 f.write('\n\n')
1447 1441 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1448 1442 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1449 1443 f.close()
1450 1444
1451 1445 def save_data(self, filename_phase, data, data_datetime):
1452 1446 f=open(filename_phase,'a')
1453 1447 timetuple_data = data_datetime.timetuple()
1454 1448 day = str(timetuple_data.tm_mday)
1455 1449 month = str(timetuple_data.tm_mon)
1456 1450 year = str(timetuple_data.tm_year)
1457 1451 hour = str(timetuple_data.tm_hour)
1458 1452 minute = str(timetuple_data.tm_min)
1459 1453 second = str(timetuple_data.tm_sec)
1460 1454 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1461 1455 f.close()
1462 1456
1463 1457
1464 1458 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1465 1459 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1466 1460 timerange=None,
1467 1461 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1468 1462 server=None, folder=None, username=None, password=None,
1469 1463 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1470 1464
1471 1465 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1472 1466 return
1473 1467
1474 1468 if pairsList == None:
1475 1469 pairsIndexList = dataOut.pairsIndexList[:10]
1476 1470 else:
1477 1471 pairsIndexList = []
1478 1472 for pair in pairsList:
1479 1473 if pair not in dataOut.pairsList:
1480 1474 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
1481 1475 pairsIndexList.append(dataOut.pairsList.index(pair))
1482 1476
1483 1477 if pairsIndexList == []:
1484 1478 return
1485 1479
1486 1480 # if len(pairsIndexList) > 4:
1487 1481 # pairsIndexList = pairsIndexList[0:4]
1488 1482
1489 1483 hmin_index = None
1490 1484 hmax_index = None
1491 1485
1492 1486 if hmin != None and hmax != None:
1493 1487 indexes = numpy.arange(dataOut.nHeights)
1494 1488 hmin_list = indexes[dataOut.heightList >= hmin]
1495 1489 hmax_list = indexes[dataOut.heightList <= hmax]
1496 1490
1497 1491 if hmin_list.any():
1498 1492 hmin_index = hmin_list[0]
1499 1493
1500 1494 if hmax_list.any():
1501 1495 hmax_index = hmax_list[-1]+1
1502 1496
1503 1497 x = dataOut.getTimeRange()
1504 1498 #y = dataOut.getHeiRange()
1505 1499
1506 1500
1507 1501 thisDatetime = dataOut.datatime
1508 1502
1509 1503 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1510 1504 xlabel = "Local Time"
1511 1505 ylabel = "Phase (degrees)"
1512 1506
1513 1507 update_figfile = False
1514 1508
1515 1509 nplots = len(pairsIndexList)
1516 1510 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1517 1511 phase_beacon = numpy.zeros(len(pairsIndexList))
1518 1512 for i in range(nplots):
1519 1513 pair = dataOut.pairsList[pairsIndexList[i]]
1520 1514 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1521 1515 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1522 1516 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1523 1517 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1524 1518 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1525 1519
1526 1520 #print "Phase %d%d" %(pair[0], pair[1])
1527 1521 #print phase[dataOut.beacon_heiIndexList]
1528 1522
1529 1523 if dataOut.beacon_heiIndexList:
1530 1524 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1531 1525 else:
1532 1526 phase_beacon[i] = numpy.average(phase)
1533 1527
1534 1528 if not self.isConfig:
1535 1529
1536 1530 nplots = len(pairsIndexList)
1537 1531
1538 1532 self.setup(id=id,
1539 1533 nplots=nplots,
1540 1534 wintitle=wintitle,
1541 1535 showprofile=showprofile,
1542 1536 show=show)
1543 1537
1544 1538 if timerange != None:
1545 1539 self.timerange = timerange
1546 1540
1547 1541 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1548 1542
1549 1543 if ymin == None: ymin = 0
1550 1544 if ymax == None: ymax = 360
1551 1545
1552 1546 self.FTP_WEI = ftp_wei
1553 1547 self.EXP_CODE = exp_code
1554 1548 self.SUB_EXP_CODE = sub_exp_code
1555 1549 self.PLOT_POS = plot_pos
1556 1550
1557 1551 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1558 1552 self.isConfig = True
1559 1553 self.figfile = figfile
1560 1554 self.xdata = numpy.array([])
1561 1555 self.ydata = numpy.array([])
1562 1556
1563 1557 update_figfile = True
1564 1558
1565 1559 #open file beacon phase
1566 1560 path = '%s%03d' %(self.PREFIX, self.id)
1567 1561 beacon_file = os.path.join(path,'%s.txt'%self.name)
1568 1562 self.filename_phase = os.path.join(figpath,beacon_file)
1569 1563 #self.save_phase(self.filename_phase)
1570 1564
1571 1565
1572 1566 #store data beacon phase
1573 1567 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1574 1568
1575 1569 self.setWinTitle(title)
1576 1570
1577 1571
1578 1572 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1579 1573
1580 1574 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1581 1575
1582 1576 axes = self.axesList[0]
1583 1577
1584 1578 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1585 1579
1586 1580 if len(self.ydata)==0:
1587 1581 self.ydata = phase_beacon.reshape(-1,1)
1588 1582 else:
1589 1583 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1590 1584
1591 1585
1592 1586 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1593 1587 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1594 1588 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1595 1589 XAxisAsTime=True, grid='both'
1596 1590 )
1597 1591
1598 1592 self.draw()
1599 1593
1600 1594 if dataOut.ltctime >= self.xmax:
1601 1595 self.counter_imagwr = wr_period
1602 1596 self.isConfig = False
1603 1597 update_figfile = True
1604 1598
1605 1599 self.save(figpath=figpath,
1606 1600 figfile=figfile,
1607 1601 save=save,
1608 1602 ftp=ftp,
1609 1603 wr_period=wr_period,
1610 1604 thisDatetime=thisDatetime,
1611 1605 update_figfile=update_figfile)
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,1068 +1,1067
1 1 import numpy
2 2
3 3 from jroproc_base import ProcessingUnit, Operation
4 4 from schainpy.model.data.jrodata import Spectra
5 5 from schainpy.model.data.jrodata import hildebrand_sekhon
6 6
7 7 import matplotlib.pyplot as plt
8 8
9 9 class SpectraProc(ProcessingUnit):
10 10
11 11 def __init__(self, **kwargs):
12 12
13 13 ProcessingUnit.__init__(self, **kwargs)
14 14
15 15 self.buffer = None
16 16 self.firstdatatime = None
17 17 self.profIndex = 0
18 18 self.dataOut = Spectra()
19 19 self.id_min = None
20 20 self.id_max = None
21 21
22 22 def __updateSpecFromVoltage(self):
23 23
24 24 self.dataOut.timeZone = self.dataIn.timeZone
25 25 self.dataOut.dstFlag = self.dataIn.dstFlag
26 26 self.dataOut.errorCount = self.dataIn.errorCount
27 27 self.dataOut.useLocalTime = self.dataIn.useLocalTime
28 28
29 29 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
30 30 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
31 31 self.dataOut.channelList = self.dataIn.channelList
32 32 self.dataOut.heightList = self.dataIn.heightList
33 33 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
34 34
35 35 self.dataOut.nBaud = self.dataIn.nBaud
36 36 self.dataOut.nCode = self.dataIn.nCode
37 37 self.dataOut.code = self.dataIn.code
38 38 self.dataOut.nProfiles = self.dataOut.nFFTPoints
39 39
40 40 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
41 41 self.dataOut.utctime = self.firstdatatime
42 42 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
43 43 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
44 44 self.dataOut.flagShiftFFT = False
45 45
46 46 self.dataOut.nCohInt = self.dataIn.nCohInt
47 47 self.dataOut.nIncohInt = 1
48 48
49 49 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
50 50
51 51 self.dataOut.frequency = self.dataIn.frequency
52 52 self.dataOut.realtime = self.dataIn.realtime
53 53
54 54 self.dataOut.azimuth = self.dataIn.azimuth
55 55 self.dataOut.zenith = self.dataIn.zenith
56 56
57 57 self.dataOut.beam.codeList = self.dataIn.beam.codeList
58 58 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
59 59 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
60 60
61 61 def __getFft(self):
62 62 """
63 63 Convierte valores de Voltaje a Spectra
64 64
65 65 Affected:
66 66 self.dataOut.data_spc
67 67 self.dataOut.data_cspc
68 68 self.dataOut.data_dc
69 69 self.dataOut.heightList
70 70 self.profIndex
71 71 self.buffer
72 72 self.dataOut.flagNoData
73 73 """
74 74 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
75 75
76 76 fft_volt = fft_volt.astype(numpy.dtype('complex'))
77 77 dc = fft_volt[:,0,:]
78 78
79 79 #calculo de self-spectra
80 80 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
81 81 spc = fft_volt * numpy.conjugate(fft_volt)
82 82 spc = spc.real
83 83
84 84 blocksize = 0
85 85 blocksize += dc.size
86 86 blocksize += spc.size
87 87
88 88 cspc = None
89 89 pairIndex = 0
90 90 if self.dataOut.pairsList != None:
91 91 #calculo de cross-spectra
92 92 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
93 93 for pair in self.dataOut.pairsList:
94 94 if pair[0] not in self.dataOut.channelList:
95 95 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
96 96 if pair[1] not in self.dataOut.channelList:
97 97 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
98 98
99 99 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
100 100 pairIndex += 1
101 101 blocksize += cspc.size
102 102
103 103 self.dataOut.data_spc = spc
104 104 self.dataOut.data_cspc = cspc
105 105 self.dataOut.data_dc = dc
106 106 self.dataOut.blockSize = blocksize
107 107 self.dataOut.flagShiftFFT = True
108 108
109 109 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None):
110 110
111 111 self.dataOut.flagNoData = True
112 112
113 113 if self.dataIn.type == "Spectra":
114 114 self.dataOut.copy(self.dataIn)
115 115 # self.__selectPairs(pairsList)
116 116 return True
117 117
118 118 if self.dataIn.type == "Voltage":
119 119
120 120 if nFFTPoints == None:
121 121 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
122 122
123 123 if nProfiles == None:
124 124 nProfiles = nFFTPoints
125 125
126 126 if ippFactor == None:
127 127 ippFactor = 1
128 128
129 129 self.dataOut.ippFactor = ippFactor
130 130
131 131 self.dataOut.nFFTPoints = nFFTPoints
132 132 self.dataOut.pairsList = pairsList
133 133
134 134 if self.buffer is None:
135 135 self.buffer = numpy.zeros( (self.dataIn.nChannels,
136 136 nProfiles,
137 137 self.dataIn.nHeights),
138 138 dtype='complex')
139 139
140 140 if self.dataIn.flagDataAsBlock:
141 141 #data dimension: [nChannels, nProfiles, nSamples]
142 142
143 143 nVoltProfiles = self.dataIn.data.shape[1]
144 144 # nVoltProfiles = self.dataIn.nProfiles
145 145
146 146 if nVoltProfiles == nProfiles:
147 147 self.buffer = self.dataIn.data.copy()
148 148 self.profIndex = nVoltProfiles
149 149
150 150 elif nVoltProfiles < nProfiles:
151 151
152 152 if self.profIndex == 0:
153 153 self.id_min = 0
154 154 self.id_max = nVoltProfiles
155 155
156 156 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
157 157 self.profIndex += nVoltProfiles
158 158 self.id_min += nVoltProfiles
159 159 self.id_max += nVoltProfiles
160 160 else:
161 161 raise ValueError, "The type object %s has %d profiles, it should just has %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
162 162 self.dataOut.flagNoData = True
163 163 return 0
164 164 else:
165 print 'DATA shape', self.dataIn.data.shape
166 sadsdf
165
167 166 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
168 167 self.profIndex += 1
169 168
170 169 if self.firstdatatime == None:
171 170 self.firstdatatime = self.dataIn.utctime
172 171
173 172 if self.profIndex == nProfiles:
174 173 self.__updateSpecFromVoltage()
175 174 self.__getFft()
176 175
177 176 self.dataOut.flagNoData = False
178 177 self.firstdatatime = None
179 178 self.profIndex = 0
180 179
181 180 return True
182 181
183 182 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
184 183
185 184 def __selectPairs(self, pairsList):
186 185
187 186 if channelList == None:
188 187 return
189 188
190 189 pairsIndexListSelected = []
191 190
192 191 for thisPair in pairsList:
193 192
194 193 if thisPair not in self.dataOut.pairsList:
195 194 continue
196 195
197 196 pairIndex = self.dataOut.pairsList.index(thisPair)
198 197
199 198 pairsIndexListSelected.append(pairIndex)
200 199
201 200 if not pairsIndexListSelected:
202 201 self.dataOut.data_cspc = None
203 202 self.dataOut.pairsList = []
204 203 return
205 204
206 205 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
207 206 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
208 207
209 208 return
210 209
211 210 def __selectPairsByChannel(self, channelList=None):
212 211
213 212 if channelList == None:
214 213 return
215 214
216 215 pairsIndexListSelected = []
217 216 for pairIndex in self.dataOut.pairsIndexList:
218 217 #First pair
219 218 if self.dataOut.pairsList[pairIndex][0] not in channelList:
220 219 continue
221 220 #Second pair
222 221 if self.dataOut.pairsList[pairIndex][1] not in channelList:
223 222 continue
224 223
225 224 pairsIndexListSelected.append(pairIndex)
226 225
227 226 if not pairsIndexListSelected:
228 227 self.dataOut.data_cspc = None
229 228 self.dataOut.pairsList = []
230 229 return
231 230
232 231 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
233 232 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
234 233
235 234 return
236 235
237 236 def selectChannels(self, channelList):
238 237
239 238 channelIndexList = []
240 239
241 240 for channel in channelList:
242 241 if channel not in self.dataOut.channelList:
243 242 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
244 243
245 244 index = self.dataOut.channelList.index(channel)
246 245 channelIndexList.append(index)
247 246
248 247 self.selectChannelsByIndex(channelIndexList)
249 248
250 249 def selectChannelsByIndex(self, channelIndexList):
251 250 """
252 251 Selecciona un bloque de datos en base a canales segun el channelIndexList
253 252
254 253 Input:
255 254 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
256 255
257 256 Affected:
258 257 self.dataOut.data_spc
259 258 self.dataOut.channelIndexList
260 259 self.dataOut.nChannels
261 260
262 261 Return:
263 262 None
264 263 """
265 264
266 265 for channelIndex in channelIndexList:
267 266 if channelIndex not in self.dataOut.channelIndexList:
268 267 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
269 268
270 269 # nChannels = len(channelIndexList)
271 270
272 271 data_spc = self.dataOut.data_spc[channelIndexList,:]
273 272 data_dc = self.dataOut.data_dc[channelIndexList,:]
274 273
275 274 self.dataOut.data_spc = data_spc
276 275 self.dataOut.data_dc = data_dc
277 276
278 277 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
279 278 # self.dataOut.nChannels = nChannels
280 279
281 280 self.__selectPairsByChannel(self.dataOut.channelList)
282 281
283 282 return 1
284 283
285 284
286 285 def selectFFTs(self, minFFT, maxFFT ):
287 286 """
288 287 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
289 288 minFFT<= FFT <= maxFFT
290 289 """
291 290
292 291 if (minFFT > maxFFT):
293 292 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT)
294 293
295 294 if (minFFT < self.dataOut.getFreqRange()[0]):
296 295 minFFT = self.dataOut.getFreqRange()[0]
297 296
298 297 if (maxFFT > self.dataOut.getFreqRange()[-1]):
299 298 maxFFT = self.dataOut.getFreqRange()[-1]
300 299
301 300 minIndex = 0
302 301 maxIndex = 0
303 302 FFTs = self.dataOut.getFreqRange()
304 303
305 304 inda = numpy.where(FFTs >= minFFT)
306 305 indb = numpy.where(FFTs <= maxFFT)
307 306
308 307 try:
309 308 minIndex = inda[0][0]
310 309 except:
311 310 minIndex = 0
312 311
313 312 try:
314 313 maxIndex = indb[0][-1]
315 314 except:
316 315 maxIndex = len(FFTs)
317 316
318 317 self.selectFFTsByIndex(minIndex, maxIndex)
319 318
320 319 return 1
321 320
322 321
323 322
324 323 def selectHeights(self, minHei, maxHei):
325 324 """
326 325 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
327 326 minHei <= height <= maxHei
328 327
329 328 Input:
330 329 minHei : valor minimo de altura a considerar
331 330 maxHei : valor maximo de altura a considerar
332 331
333 332 Affected:
334 333 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
335 334
336 335 Return:
337 336 1 si el metodo se ejecuto con exito caso contrario devuelve 0
338 337 """
339 338
340 339
341 340 if (minHei > maxHei):
342 341 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
343 342
344 343 if (minHei < self.dataOut.heightList[0]):
345 344 minHei = self.dataOut.heightList[0]
346 345
347 346 if (maxHei > self.dataOut.heightList[-1]):
348 347 maxHei = self.dataOut.heightList[-1]
349 348
350 349 minIndex = 0
351 350 maxIndex = 0
352 351 heights = self.dataOut.heightList
353 352
354 353 inda = numpy.where(heights >= minHei)
355 354 indb = numpy.where(heights <= maxHei)
356 355
357 356 try:
358 357 minIndex = inda[0][0]
359 358 except:
360 359 minIndex = 0
361 360
362 361 try:
363 362 maxIndex = indb[0][-1]
364 363 except:
365 364 maxIndex = len(heights)
366 365
367 366 self.selectHeightsByIndex(minIndex, maxIndex)
368 367
369 368
370 369 return 1
371 370
372 371 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
373 372 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
374 373
375 374 if hei_ref != None:
376 375 newheis = numpy.where(self.dataOut.heightList>hei_ref)
377 376
378 377 minIndex = min(newheis[0])
379 378 maxIndex = max(newheis[0])
380 379 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
381 380 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
382 381
383 382 # determina indices
384 383 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
385 384 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
386 385 beacon_dB = numpy.sort(avg_dB)[-nheis:]
387 386 beacon_heiIndexList = []
388 387 for val in avg_dB.tolist():
389 388 if val >= beacon_dB[0]:
390 389 beacon_heiIndexList.append(avg_dB.tolist().index(val))
391 390
392 391 #data_spc = data_spc[:,:,beacon_heiIndexList]
393 392 data_cspc = None
394 393 if self.dataOut.data_cspc is not None:
395 394 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
396 395 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
397 396
398 397 data_dc = None
399 398 if self.dataOut.data_dc is not None:
400 399 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
401 400 #data_dc = data_dc[:,beacon_heiIndexList]
402 401
403 402 self.dataOut.data_spc = data_spc
404 403 self.dataOut.data_cspc = data_cspc
405 404 self.dataOut.data_dc = data_dc
406 405 self.dataOut.heightList = heightList
407 406 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
408 407
409 408 return 1
410 409
411 410 def selectFFTsByIndex(self, minIndex, maxIndex):
412 411 """
413 412
414 413 """
415 414
416 415 if (minIndex < 0) or (minIndex > maxIndex):
417 416 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
418 417
419 418 if (maxIndex >= self.dataOut.nProfiles):
420 419 maxIndex = self.dataOut.nProfiles-1
421 420
422 421 #Spectra
423 422 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
424 423
425 424 data_cspc = None
426 425 if self.dataOut.data_cspc is not None:
427 426 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
428 427
429 428 data_dc = None
430 429 if self.dataOut.data_dc is not None:
431 430 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
432 431
433 432 self.dataOut.data_spc = data_spc
434 433 self.dataOut.data_cspc = data_cspc
435 434 self.dataOut.data_dc = data_dc
436 435
437 436 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
438 437 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
439 438 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
440 439
441 440 #self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
442 441
443 442 return 1
444 443
445 444
446 445
447 446 def selectHeightsByIndex(self, minIndex, maxIndex):
448 447 """
449 448 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
450 449 minIndex <= index <= maxIndex
451 450
452 451 Input:
453 452 minIndex : valor de indice minimo de altura a considerar
454 453 maxIndex : valor de indice maximo de altura a considerar
455 454
456 455 Affected:
457 456 self.dataOut.data_spc
458 457 self.dataOut.data_cspc
459 458 self.dataOut.data_dc
460 459 self.dataOut.heightList
461 460
462 461 Return:
463 462 1 si el metodo se ejecuto con exito caso contrario devuelve 0
464 463 """
465 464
466 465 if (minIndex < 0) or (minIndex > maxIndex):
467 466 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
468 467
469 468 if (maxIndex >= self.dataOut.nHeights):
470 469 maxIndex = self.dataOut.nHeights-1
471 470
472 471 #Spectra
473 472 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
474 473
475 474 data_cspc = None
476 475 if self.dataOut.data_cspc is not None:
477 476 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
478 477
479 478 data_dc = None
480 479 if self.dataOut.data_dc is not None:
481 480 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
482 481
483 482 self.dataOut.data_spc = data_spc
484 483 self.dataOut.data_cspc = data_cspc
485 484 self.dataOut.data_dc = data_dc
486 485
487 486 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
488 487
489 488 return 1
490 489
491 490
492 491 def removeDC(self, mode = 2):
493 492 jspectra = self.dataOut.data_spc
494 493 jcspectra = self.dataOut.data_cspc
495 494
496 495
497 496 num_chan = jspectra.shape[0]
498 497 num_hei = jspectra.shape[2]
499 498
500 499 if jcspectra is not None:
501 500 jcspectraExist = True
502 501 num_pairs = jcspectra.shape[0]
503 502 else: jcspectraExist = False
504 503
505 504 freq_dc = jspectra.shape[1]/2
506 505 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
507 506
508 507 if ind_vel[0]<0:
509 508 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
510 509
511 510 if mode == 1:
512 511 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
513 512
514 513 if jcspectraExist:
515 514 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
516 515
517 516 if mode == 2:
518 517
519 518 vel = numpy.array([-2,-1,1,2])
520 519 xx = numpy.zeros([4,4])
521 520
522 521 for fil in range(4):
523 522 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
524 523
525 524 xx_inv = numpy.linalg.inv(xx)
526 525 xx_aux = xx_inv[0,:]
527 526
528 527 for ich in range(num_chan):
529 528 yy = jspectra[ich,ind_vel,:]
530 529 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
531 530
532 531 junkid = jspectra[ich,freq_dc,:]<=0
533 532 cjunkid = sum(junkid)
534 533
535 534 if cjunkid.any():
536 535 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
537 536
538 537 if jcspectraExist:
539 538 for ip in range(num_pairs):
540 539 yy = jcspectra[ip,ind_vel,:]
541 540 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
542 541
543 542
544 543 self.dataOut.data_spc = jspectra
545 544 self.dataOut.data_cspc = jcspectra
546 545
547 546 return 1
548 547
549 548 def removeInterference2(self):
550 549
551 550 cspc = self.dataOut.data_cspc
552 551 spc = self.dataOut.data_spc
553 552 print numpy.shape(spc)
554 553 Heights = numpy.arange(cspc.shape[2])
555 554 realCspc = numpy.abs(cspc)
556 555
557 556 for i in range(cspc.shape[0]):
558 557 LinePower= numpy.sum(realCspc[i], axis=0)
559 558 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
560 559 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
561 560 #print numpy.shape(realCspc)
562 561 #print '',SelectedHeights, '', numpy.shape(realCspc[i,:,SelectedHeights])
563 562 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
564 563 print SelectedHeights
565 564 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
566 565 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
567 566
568 567
569 568 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
570 569 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
571 570 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
572 571 cspc[i,InterferenceRange,:] = numpy.NaN
573 572
574 573 print '########################################################################################'
575 574 print 'Len interference sum',len(InterferenceSum)
576 575 print 'InterferenceThresholdMin', InterferenceThresholdMin, 'InterferenceThresholdMax', InterferenceThresholdMax
577 576 print 'InterferenceRange',InterferenceRange
578 577 print '########################################################################################'
579 578
580 579 ''' Ploteo '''
581 580
582 581 #for i in range(3):
583 582 #print 'FASE', numpy.shape(phase), y[25]
584 583 #print numpy.shape(coherence)
585 584 #fig = plt.figure(10+ int(numpy.random.rand()*100))
586 585 #plt.plot( x[0:256],coherence[:,25] )
587 586 #cohAv = numpy.average(coherence[i],1)
588 587 #Pendiente = FrecRange * PhaseSlope[i]
589 588 #plt.plot( InterferenceSum)
590 589 #plt.plot( numpy.sort(InterferenceSum))
591 590 #plt.plot( LinePower )
592 591 #plt.plot( xFrec,phase[i])
593 592
594 593 #CSPCmean = numpy.mean(numpy.abs(CSPCSamples),0)
595 594 #plt.plot(xFrec, FitGauss01)
596 595 #plt.plot(xFrec, CSPCmean)
597 596 #plt.plot(xFrec, numpy.abs(CSPCSamples[0]))
598 597 #plt.plot(xFrec, FitGauss)
599 598 #plt.plot(xFrec, yMean)
600 599 #plt.plot(xFrec, numpy.abs(coherence[0]))
601 600
602 601 #plt.axis([-12, 12, 15, 50])
603 602 #plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
604 603
605 604
606 605 #fig.savefig('/home/erick/Documents/Pics/nom{}.png'.format(int(numpy.random.rand()*100)))
607 606
608 607 #plt.show()
609 608 #self.indice=self.indice+1
610 609 #raise
611 610
612 611
613 612 self.dataOut.data_cspc = cspc
614 613
615 614 # for i in range(spc.shape[0]):
616 615 # LinePower= numpy.sum(spc[i], axis=0)
617 616 # Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
618 617 # SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
619 618 # #print numpy.shape(realCspc)
620 619 # #print '',SelectedHeights, '', numpy.shape(realCspc[i,:,SelectedHeights])
621 620 # InterferenceSum = numpy.sum( spc[i,:,SelectedHeights], axis=0 )
622 621 # InterferenceThreshold = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
623 622 # InterferenceRange = numpy.where( InterferenceSum > InterferenceThreshold )
624 623 # if len(InterferenceRange)<int(spc.shape[1]*0.03):
625 624 # spc[i,InterferenceRange,:] = numpy.NaN
626 625
627 626 #self.dataOut.data_spc = spc
628 627
629 628 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
630 629
631 630 jspectra = self.dataOut.data_spc
632 631 jcspectra = self.dataOut.data_cspc
633 632 jnoise = self.dataOut.getNoise()
634 633 num_incoh = self.dataOut.nIncohInt
635 634
636 635 num_channel = jspectra.shape[0]
637 636 num_prof = jspectra.shape[1]
638 637 num_hei = jspectra.shape[2]
639 638
640 639 #hei_interf
641 640 if hei_interf is None:
642 641 count_hei = num_hei/2 #Como es entero no importa
643 642 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
644 643 hei_interf = numpy.asarray(hei_interf)[0]
645 644 #nhei_interf
646 645 if (nhei_interf == None):
647 646 nhei_interf = 5
648 647 if (nhei_interf < 1):
649 648 nhei_interf = 1
650 649 if (nhei_interf > count_hei):
651 650 nhei_interf = count_hei
652 651 if (offhei_interf == None):
653 652 offhei_interf = 0
654 653
655 654 ind_hei = range(num_hei)
656 655 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
657 656 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
658 657 mask_prof = numpy.asarray(range(num_prof))
659 658 num_mask_prof = mask_prof.size
660 659 comp_mask_prof = [0, num_prof/2]
661 660
662 661
663 662 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
664 663 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
665 664 jnoise = numpy.nan
666 665 noise_exist = jnoise[0] < numpy.Inf
667 666
668 667 #Subrutina de Remocion de la Interferencia
669 668 for ich in range(num_channel):
670 669 #Se ordena los espectros segun su potencia (menor a mayor)
671 670 power = jspectra[ich,mask_prof,:]
672 671 power = power[:,hei_interf]
673 672 power = power.sum(axis = 0)
674 673 psort = power.ravel().argsort()
675 674
676 675 #Se estima la interferencia promedio en los Espectros de Potencia empleando
677 676 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
678 677
679 678 if noise_exist:
680 679 # tmp_noise = jnoise[ich] / num_prof
681 680 tmp_noise = jnoise[ich]
682 681 junkspc_interf = junkspc_interf - tmp_noise
683 682 #junkspc_interf[:,comp_mask_prof] = 0
684 683
685 684 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
686 685 jspc_interf = jspc_interf.transpose()
687 686 #Calculando el espectro de interferencia promedio
688 687 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
689 688 noiseid = noiseid[0]
690 689 cnoiseid = noiseid.size
691 690 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
692 691 interfid = interfid[0]
693 692 cinterfid = interfid.size
694 693
695 694 if (cnoiseid > 0): jspc_interf[noiseid] = 0
696 695
697 696 #Expandiendo los perfiles a limpiar
698 697 if (cinterfid > 0):
699 698 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
700 699 new_interfid = numpy.asarray(new_interfid)
701 700 new_interfid = {x for x in new_interfid}
702 701 new_interfid = numpy.array(list(new_interfid))
703 702 new_cinterfid = new_interfid.size
704 703 else: new_cinterfid = 0
705 704
706 705 for ip in range(new_cinterfid):
707 706 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
708 707 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
709 708
710 709
711 710 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
712 711
713 712 #Removiendo la interferencia del punto de mayor interferencia
714 713 ListAux = jspc_interf[mask_prof].tolist()
715 714 maxid = ListAux.index(max(ListAux))
716 715
717 716
718 717 if cinterfid > 0:
719 718 for ip in range(cinterfid*(interf == 2) - 1):
720 719 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
721 720 cind = len(ind)
722 721
723 722 if (cind > 0):
724 723 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
725 724
726 725 ind = numpy.array([-2,-1,1,2])
727 726 xx = numpy.zeros([4,4])
728 727
729 728 for id1 in range(4):
730 729 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
731 730
732 731 xx_inv = numpy.linalg.inv(xx)
733 732 xx = xx_inv[:,0]
734 733 ind = (ind + maxid + num_mask_prof)%num_mask_prof
735 734 yy = jspectra[ich,mask_prof[ind],:]
736 735 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
737 736
738 737
739 738 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
740 739 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
741 740
742 741 #Remocion de Interferencia en el Cross Spectra
743 742 if jcspectra is None: return jspectra, jcspectra
744 743 num_pairs = jcspectra.size/(num_prof*num_hei)
745 744 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
746 745
747 746 for ip in range(num_pairs):
748 747
749 748 #-------------------------------------------
750 749
751 750 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
752 751 cspower = cspower[:,hei_interf]
753 752 cspower = cspower.sum(axis = 0)
754 753
755 754 cspsort = cspower.ravel().argsort()
756 755 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
757 756 junkcspc_interf = junkcspc_interf.transpose()
758 757 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
759 758
760 759 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
761 760
762 761 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
763 762 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
764 763 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
765 764
766 765 for iprof in range(num_prof):
767 766 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
768 767 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
769 768
770 769 #Removiendo la Interferencia
771 770 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
772 771
773 772 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
774 773 maxid = ListAux.index(max(ListAux))
775 774
776 775 ind = numpy.array([-2,-1,1,2])
777 776 xx = numpy.zeros([4,4])
778 777
779 778 for id1 in range(4):
780 779 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
781 780
782 781 xx_inv = numpy.linalg.inv(xx)
783 782 xx = xx_inv[:,0]
784 783
785 784 ind = (ind + maxid + num_mask_prof)%num_mask_prof
786 785 yy = jcspectra[ip,mask_prof[ind],:]
787 786 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
788 787
789 788 #Guardar Resultados
790 789 self.dataOut.data_spc = jspectra
791 790 self.dataOut.data_cspc = jcspectra
792 791
793 792 return 1
794 793
795 794 def setRadarFrequency(self, frequency=None):
796 795
797 796 if frequency != None:
798 797 self.dataOut.frequency = frequency
799 798
800 799 return 1
801 800
802 801 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
803 802 #validacion de rango
804 803 if minHei == None:
805 804 minHei = self.dataOut.heightList[0]
806 805
807 806 if maxHei == None:
808 807 maxHei = self.dataOut.heightList[-1]
809 808
810 809 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
811 810 print 'minHei: %.2f is out of the heights range'%(minHei)
812 811 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
813 812 minHei = self.dataOut.heightList[0]
814 813
815 814 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
816 815 print 'maxHei: %.2f is out of the heights range'%(maxHei)
817 816 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
818 817 maxHei = self.dataOut.heightList[-1]
819 818
820 819 # validacion de velocidades
821 820 velrange = self.dataOut.getVelRange(1)
822 821
823 822 if minVel == None:
824 823 minVel = velrange[0]
825 824
826 825 if maxVel == None:
827 826 maxVel = velrange[-1]
828 827
829 828 if (minVel < velrange[0]) or (minVel > maxVel):
830 829 print 'minVel: %.2f is out of the velocity range'%(minVel)
831 830 print 'minVel is setting to %.2f'%(velrange[0])
832 831 minVel = velrange[0]
833 832
834 833 if (maxVel > velrange[-1]) or (maxVel < minVel):
835 834 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
836 835 print 'maxVel is setting to %.2f'%(velrange[-1])
837 836 maxVel = velrange[-1]
838 837
839 838 # seleccion de indices para rango
840 839 minIndex = 0
841 840 maxIndex = 0
842 841 heights = self.dataOut.heightList
843 842
844 843 inda = numpy.where(heights >= minHei)
845 844 indb = numpy.where(heights <= maxHei)
846 845
847 846 try:
848 847 minIndex = inda[0][0]
849 848 except:
850 849 minIndex = 0
851 850
852 851 try:
853 852 maxIndex = indb[0][-1]
854 853 except:
855 854 maxIndex = len(heights)
856 855
857 856 if (minIndex < 0) or (minIndex > maxIndex):
858 857 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
859 858
860 859 if (maxIndex >= self.dataOut.nHeights):
861 860 maxIndex = self.dataOut.nHeights-1
862 861
863 862 # seleccion de indices para velocidades
864 863 indminvel = numpy.where(velrange >= minVel)
865 864 indmaxvel = numpy.where(velrange <= maxVel)
866 865 try:
867 866 minIndexVel = indminvel[0][0]
868 867 except:
869 868 minIndexVel = 0
870 869
871 870 try:
872 871 maxIndexVel = indmaxvel[0][-1]
873 872 except:
874 873 maxIndexVel = len(velrange)
875 874
876 875 #seleccion del espectro
877 876 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
878 877 #estimacion de ruido
879 878 noise = numpy.zeros(self.dataOut.nChannels)
880 879
881 880 for channel in range(self.dataOut.nChannels):
882 881 daux = data_spc[channel,:,:]
883 882 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
884 883
885 884 self.dataOut.noise_estimation = noise.copy()
886 885
887 886 return 1
888 887
889 888 class IncohInt(Operation):
890 889
891 890
892 891 __profIndex = 0
893 892 __withOverapping = False
894 893
895 894 __byTime = False
896 895 __initime = None
897 896 __lastdatatime = None
898 897 __integrationtime = None
899 898
900 899 __buffer_spc = None
901 900 __buffer_cspc = None
902 901 __buffer_dc = None
903 902
904 903 __dataReady = False
905 904
906 905 __timeInterval = None
907 906
908 907 n = None
909 908
910 909
911 910
912 911 def __init__(self, **kwargs):
913 912
914 913 Operation.__init__(self, **kwargs)
915 914 # self.isConfig = False
916 915
917 916 def setup(self, n=None, timeInterval=None, overlapping=False):
918 917 """
919 918 Set the parameters of the integration class.
920 919
921 920 Inputs:
922 921
923 922 n : Number of coherent integrations
924 923 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
925 924 overlapping :
926 925
927 926 """
928 927
929 928 self.__initime = None
930 929 self.__lastdatatime = 0
931 930
932 931 self.__buffer_spc = 0
933 932 self.__buffer_cspc = 0
934 933 self.__buffer_dc = 0
935 934
936 935 self.__profIndex = 0
937 936 self.__dataReady = False
938 937 self.__byTime = False
939 938
940 939 if n is None and timeInterval is None:
941 940 raise ValueError, "n or timeInterval should be specified ..."
942 941
943 942 if n is not None:
944 943 self.n = int(n)
945 944 else:
946 945 self.__integrationtime = int(timeInterval) #if (type(timeInterval)!=integer) -> change this line
947 946 self.n = None
948 947 self.__byTime = True
949 948
950 949 def putData(self, data_spc, data_cspc, data_dc):
951 950
952 951 """
953 952 Add a profile to the __buffer_spc and increase in one the __profileIndex
954 953
955 954 """
956 955
957 956 self.__buffer_spc += data_spc
958 957
959 958 if data_cspc is None:
960 959 self.__buffer_cspc = None
961 960 else:
962 961 self.__buffer_cspc += data_cspc
963 962
964 963 if data_dc is None:
965 964 self.__buffer_dc = None
966 965 else:
967 966 self.__buffer_dc += data_dc
968 967
969 968 self.__profIndex += 1
970 969
971 970 return
972 971
973 972 def pushData(self):
974 973 """
975 974 Return the sum of the last profiles and the profiles used in the sum.
976 975
977 976 Affected:
978 977
979 978 self.__profileIndex
980 979
981 980 """
982 981
983 982 data_spc = self.__buffer_spc
984 983 data_cspc = self.__buffer_cspc
985 984 data_dc = self.__buffer_dc
986 985 n = self.__profIndex
987 986
988 987 self.__buffer_spc = 0
989 988 self.__buffer_cspc = 0
990 989 self.__buffer_dc = 0
991 990 self.__profIndex = 0
992 991
993 992 return data_spc, data_cspc, data_dc, n
994 993
995 994 def byProfiles(self, *args):
996 995
997 996 self.__dataReady = False
998 997 avgdata_spc = None
999 998 avgdata_cspc = None
1000 999 avgdata_dc = None
1001 1000
1002 1001 self.putData(*args)
1003 1002
1004 1003 if self.__profIndex == self.n:
1005 1004
1006 1005 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1007 1006 self.n = n
1008 1007 self.__dataReady = True
1009 1008
1010 1009 return avgdata_spc, avgdata_cspc, avgdata_dc
1011 1010
1012 1011 def byTime(self, datatime, *args):
1013 1012
1014 1013 self.__dataReady = False
1015 1014 avgdata_spc = None
1016 1015 avgdata_cspc = None
1017 1016 avgdata_dc = None
1018 1017
1019 1018 self.putData(*args)
1020 1019
1021 1020 if (datatime - self.__initime) >= self.__integrationtime:
1022 1021 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1023 1022 self.n = n
1024 1023 self.__dataReady = True
1025 1024
1026 1025 return avgdata_spc, avgdata_cspc, avgdata_dc
1027 1026
1028 1027 def integrate(self, datatime, *args):
1029 1028
1030 1029 if self.__profIndex == 0:
1031 1030 self.__initime = datatime
1032 1031
1033 1032 if self.__byTime:
1034 1033 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
1035 1034 else:
1036 1035 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1037 1036
1038 1037 if not self.__dataReady:
1039 1038 return None, None, None, None
1040 1039
1041 1040 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1042 1041
1043 1042 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1044 1043
1045 1044 if n==1:
1046 1045 return
1047 1046
1048 1047 dataOut.flagNoData = True
1049 1048
1050 1049 if not self.isConfig:
1051 1050 self.setup(n, timeInterval, overlapping)
1052 1051 self.isConfig = True
1053 1052
1054 1053 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1055 1054 dataOut.data_spc,
1056 1055 dataOut.data_cspc,
1057 1056 dataOut.data_dc)
1058 1057
1059 1058 if self.__dataReady:
1060 1059
1061 1060 dataOut.data_spc = avgdata_spc
1062 1061 dataOut.data_cspc = avgdata_cspc
1063 1062 dataOut.data_dc = avgdata_dc
1064 1063
1065 1064 dataOut.nIncohInt *= self.n
1066 1065 dataOut.utctime = avgdatatime
1067 1066 dataOut.flagNoData = False
1068 1067
Show file before | Show file after | Hide comments
@@ -1,1 +1,1
1 <Project description="Segundo Test" id="191" name="test01"><ReadUnit datatype="SpectraReader" id="1911" inputId="0" name="SpectraReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="SpectraReader" /><Parameter format="str" id="191112" name="path" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdatatest/test1024/d2017234" /><Parameter format="date" id="191113" name="startDate" value="2017/08/22" /><Parameter format="date" id="191114" name="endDate" value="2017/08/22" /><Parameter format="time" id="191115" name="startTime" value="02:00:00" /><Parameter format="time" id="191116" name="endTime" value="06:00:00" /><Parameter format="int" id="191118" name="online" value="0" /><Parameter format="int" id="191119" name="walk" value="0" /></Operation><Operation id="19112" name="printInfo" priority="2" type="self" /><Operation id="19113" name="printNumberOfBlock" priority="3" type="self" /></ReadUnit><ProcUnit datatype="Parameters" id="1913" inputId="1912" name="ParametersProc"><Operation id="19131" name="run" priority="1" type="self" /><Operation id="19132" name="SpectralMoments" priority="2" type="other" /><Operation id="19133" name="FullSpectralAnalysis" priority="3" type="other"><Parameter format="float" id="191331" name="SNRlimit" value="-16" /><Parameter format="float" id="191332" name="E01" value="1.500" /><Parameter format="float" id="191333" name="E02" value="1.500" /><Parameter format="float" id="191334" name="E12" value="0" /><Parameter format="float" id="191335" name="N01" value="0.875" /><Parameter format="float" id="191336" name="N02" value="-0.875" /><Parameter format="float" id="191337" name="N12" value="-1.750" /></Operation><Operation id="19134" name="ParamWriter" priority="4" type="other"><Parameter format="str" id="191341" name="path" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdatatest/test1024" /><Parameter format="int" id="191342" name="blocksPerFile" value="100" /><Parameter format="list" id="191343" name="metadataList" value="heightList,timeZone,paramInterval" /><Parameter format="list" id="191344" name="dataList" value="data_output,data_SNR,utctime,utctimeInit" /></Operation></ProcUnit><ProcUnit datatype="SpectraProc" id="1912" inputId="1911" name="SpectraProc"><Operation id="19121" name="run" priority="1" type="self"><Parameter format="pairslist" id="191211" name="pairsList" value="(0,1),(0,2),(1,2)" /></Operation><Operation id="19122" name="setRadarFrequency" priority="2" type="self"><Parameter format="float" id="191221" name="frequency" value="445.09e6" /></Operation><Operation id="19123" name="IncohInt" priority="3" type="external"><Parameter format="float" id="191231" name="n" value="5" /></Operation><Operation id="19124" name="SpectraPlot" priority="4" type="external"><Parameter format="int" id="191241" name="id" value="11" /><Parameter format="str" id="191242" name="wintitle" value="SpectraPlot" /><Parameter format="str" id="191243" name="xaxis" value="frequency" /><Parameter format="float" id="191244" name="ymin" value="1" /><Parameter format="int" id="191245" name="ymax" value="5" /><Parameter format="int" id="191246" name="zmin" value="15" /><Parameter format="int" id="191247" name="zmax" value="50" /><Parameter format="int" id="191248" name="save" value="2" /><Parameter format="int" id="191249" name="save" value="5" /><Parameter format="str" id="191250" name="figpath" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/Images/FirstResults1024" /></Operation><Operation id="19125" name="CrossSpectraPlot" priority="5" type="other"><Parameter format="int" id="191251" name="id" value="2005" /><Parameter format="str" id="191252" name="wintitle" value="CrossSpectraPlot_ShortPulse" /><Parameter format="str" id="191253" name="xaxis" value="Velocity" /><Parameter format="float" id="191254" name="xmin" value="-10" /><Parameter format="float" id="191255" name="xmax" value="10" /><Parameter format="int" id="191256" name="zmin" value="15" /><Parameter format="int" id="191257" name="zmax" value="50" /><Parameter format="str" id="191258" name="phase_cmap" value="bwr" /><Parameter format="float" id="191259" name="ymin" value="1" /><Parameter format="float" id="191260" name="ymax" value="5" /></Operation></ProcUnit></Project> No newline at end of file
1 <Project description="Segundo Test" id="191" name="test01"><ReadUnit datatype="SpectraReader" id="1911" inputId="0" name="SpectraReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="SpectraReader" /><Parameter format="str" id="191112" name="path" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /><Parameter format="date" id="191113" name="startDate" value="2018/02/01" /><Parameter format="date" id="191114" name="endDate" value="2018/02/01" /><Parameter format="time" id="191115" name="startTime" value="17:00:00" /><Parameter format="time" id="191116" name="endTime" value="20:00:00" /><Parameter format="int" id="191118" name="online" value="0" /><Parameter format="int" id="191119" name="walk" value="1" /></Operation><Operation id="19112" name="printInfo" priority="2" type="self" /><Operation id="19113" name="printNumberOfBlock" priority="3" type="self" /></ReadUnit><ProcUnit datatype="Parameters" id="1913" inputId="1912" name="ParametersProc"><Operation id="19131" name="run" priority="1" type="self" /><Operation id="19132" name="SpectralFilters" priority="2" type="other"><Parameter format="float" id="191321" name="Wind_Velocity_Limit" value="2.5" /><Parameter format="float" id="191322" name="Rain_Velocity_Limit" value="1.5" /></Operation><Operation id="19133" name="FullSpectralAnalysis" priority="3" type="other"><Parameter format="float" id="191331" name="SNRlimit" value="-16" /><Parameter format="float" id="191332" name="E01" value="1.500" /><Parameter format="float" id="191333" name="E02" value="1.500" /><Parameter format="float" id="191334" name="E12" value="0" /><Parameter format="float" id="191335" name="N01" value="0.875" /><Parameter format="float" id="191336" name="N02" value="-0.875" /><Parameter format="float" id="191337" name="N12" value="-1.750" /></Operation><Operation id="19134" name="WindProfilerPlot" priority="4" type="other"><Parameter format="int" id="191341" name="id" value="4" /><Parameter format="str" id="191342" name="wintitle" value="Wind Profiler" /><Parameter format="float" id="191343" name="xmin" value="17" /><Parameter format="float" id="191344" name="xmax" value="20" /><Parameter format="float" id="191345" name="ymin" value="0" /><Parameter format="int" id="191346" name="ymax" value="8" /><Parameter format="float" id="191347" name="zmin" value="-20" /><Parameter format="float" id="191348" name="zmax" value="20" /><Parameter format="float" id="191349" name="SNRmin" value="-20" /><Parameter format="float" id="191350" name="SNRmax" value="20" /><Parameter format="float" id="191351" name="zmin_ver" value="-200" /><Parameter format="float" id="191352" name="zmax_ver" value="200" /><Parameter format="float" id="191353" name="SNRthresh" value="-20" /><Parameter format="int" id="191354" name="save" value="1" /><Parameter format="str" id="191355" name="figpath" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /></Operation><Operation id="19135" name="PrecipitationProc" priority="5" type="other" /><Operation id="19136" name="ParametersPlot" priority="6" type="other"><Parameter format="int" id="191361" name="id" value="10" /><Parameter format="str" id="191362" name="wintitle" value="First_gg" /><Parameter format="int" id="191363" name="zmin" value="-20" /><Parameter format="int" id="191364" name="zmax" value="60" /><Parameter format="int" id="191365" name="ymin" value="0" /><Parameter format="int" id="191366" name="ymax" value="8" /><Parameter format="int" id="191367" name="xmin" value="17" /><Parameter format="int" id="191368" name="xmax" value="20" /><Parameter format="int" id="191369" name="save" value="1" /><Parameter format="str" id="191370" name="figpath" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /></Operation><Operation id="19137" name="ParamWriter" priority="7" type="other"><Parameter format="str" id="191371" name="path" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdatatest/test1024" /><Parameter format="int" id="191372" name="blocksPerFile" value="100" /><Parameter format="list" id="191373" name="metadataList" value="heightList,timeZone,paramInterval" /><Parameter format="list" id="191374" name="dataList" value="data_output,data_SNR,utctime,utctimeInit" /></Operation></ProcUnit><ProcUnit datatype="SpectraProc" id="1912" inputId="1911" name="SpectraProc"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="setRadarFrequency" priority="2" type="self"><Parameter format="float" id="191221" name="frequency" value="445.09e6" /></Operation></ProcUnit></Project> No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,238 +1,239
1 1 import os, sys
2 2
3 3 from schainpy.controller import Project
4 4
5 5 if __name__ == '__main__':
6 6
7 7 desc = "Segundo Test"
8 8 filename = "schain.xml"
9 9
10 pathW='/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdatatest/test1024'
11 figpath = '/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/Images/test1024'
10 pathW='/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra'
11 figpath = '/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra'
12 12
13 13 controllerObj = Project()
14 14
15 15 controllerObj.setup(id='191', name='test01', description=desc)
16 16
17 17 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
18 18 path='/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/',
19 19 #path='/home/erick/Documents/Data/Claire_Data/raw',
20 startDate='2017/08/22',
21 endDate='2017/08/22',
22 startTime='01:00:00',
23 endTime='06:00:00',
20 startDate='2018/02/01',
21 endDate='2018/02/01',
22 startTime='12:00:00',
23 endTime='20:00:00',
24 24 online=0,
25 25 walk=1)
26 26
27 27 opObj00 = readUnitConfObj.addOperation(name='printInfo')
28 28 #
29 29 # procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc',
30 30 # inputId=readUnitConfObj.getId())
31 31 #
32 32 # opObj10 = procUnitConfObj0.addOperation(name='selectHeights')
33 33 # opObj10.addParameter(name='minHei', value='0', format='float')
34 34 # opObj10.addParameter(name='maxHei', value='8', format='float')
35 35 #
36 36 # opObj10 = procUnitConfObj0.addOperation(name='filterByHeights')
37 37 # opObj10.addParameter(name='window', value='2', format='float')
38 38 #
39 39 # opObj10 = procUnitConfObj0.addOperation(name='Decoder', optype='external')
40 40 # opObj10.addParameter(name='code', value='1,-1', format='intlist')
41 41 # opObj10.addParameter(name='nCode', value='2', format='float')
42 42 # opObj10.addParameter(name='nBaud', value='1', format='float')
43 43 #
44 44 #
45 45 # opObj10 = procUnitConfObj0.addOperation(name='CohInt', optype='external')
46 46 # opObj10.addParameter(name='n', value='1296', format='float')
47 47
48 48 opObj00 = readUnitConfObj.addOperation(name='printNumberOfBlock')
49 49
50 50 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc',
51 51 inputId=readUnitConfObj.getId())
52 52
53 53
54 54 opObj10 = procUnitConfObj0.addOperation(name='setRadarFrequency')
55 55 opObj10.addParameter(name='frequency', value='445.09e6', format='float')
56 56
57 #opObj10 = procUnitConfObj0.addOperation(name='CohInt', optype='external')
58 #opObj10.addParameter(name='n', value='1', format='float')
57 opObj10 = procUnitConfObj0.addOperation(name='CohInt', optype='external')
58 opObj10.addParameter(name='n', value='2', format='float')
59 59
60 60 #opObj10 = procUnitConfObj0.addOperation(name='selectHeights')
61 61 #opObj10.addParameter(name='minHei', value='1', format='float')
62 62 #opObj10.addParameter(name='maxHei', value='15', format='float')
63 63
64 64 #opObj10 = procUnitConfObj0.addOperation(name='selectFFTs')
65 65 #opObj10.addParameter(name='minHei', value='', format='float')
66 66 #opObj10.addParameter(name='maxHei', value='', format='float')
67 67
68 68 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc',
69 69 inputId=procUnitConfObj0.getId())
70 70
71 71 # Creating a processing object with its parameters
72 72 # schainpy.model.proc.jroproc_spectra.SpectraProc.run()
73 73 # If you need to add more parameters can use the "addParameter method"
74 procUnitConfObj1.addParameter(name='nFFTPoints', value='1024', format='int')
74 procUnitConfObj1.addParameter(name='nFFTPoints', value='256', format='int')
75 75
76 76
77 77 opObj10 = procUnitConfObj1.addOperation(name='removeDC')
78 78 #opObj10 = procUnitConfObj1.addOperation(name='removeInterference')
79 #opObj10 = procUnitConfObj1.addOperation(name='IncohInt', optype='external')
80 #opObj10.addParameter(name='n', value='30', format='float')
79 opObj10 = procUnitConfObj1.addOperation(name='IncohInt', optype='external')
80 opObj10.addParameter(name='n', value='10', format='float')
81 81
82 82
83 83
84 84 #opObj10 = procUnitConfObj1.addOperation(name='selectFFTs')
85 85 #opObj10.addParameter(name='minFFT', value='-15', format='float')
86 86 #opObj10.addParameter(name='maxFFT', value='15', format='float')
87 87
88 88
89 89
90 90 opObj10 = procUnitConfObj1.addOperation(name='SpectraWriter', optype='other')
91 91 opObj10.addParameter(name='blocksPerFile', value='64', format = 'int')
92 92 opObj10.addParameter(name='path', value=pathW)
93 93 # Using internal methods
94 94 # schainpy.model.proc.jroproc_spectra.SpectraProc.selectChannels()
95 95 # opObj10 = procUnitConfObj1.addOperation(name='selectChannels')
96 96 # opObj10.addParameter(name='channelList', value='0,1', format='intlist')
97 97
98 98 # Using internal methods
99 99 # schainpy.model.proc.jroproc_spectra.SpectraProc.selectHeights()
100 100 # opObj10 = procUnitConfObj1.addOperation(name='selectHeights')
101 101 # opObj10.addParameter(name='minHei', value='90', format='float')
102 102 # opObj10.addParameter(name='maxHei', value='180', format='float')
103 103
104 104 # Using external methods (new modules)
105 105 # #schainpy.model.proc.jroproc_spectra.IncohInt.setup()
106 106 # opObj12 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
107 107 # opObj12.addParameter(name='n', value='1', format='int')
108 108
109 109 # Using external methods (new modules)
110 110 # schainpy.model.graphics.jroplot_spectra.SpectraPlot.setup()
111 111 opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='external')
112 112 opObj11.addParameter(name='id', value='11', format='int')
113 113 opObj11.addParameter(name='wintitle', value='SpectraPlot', format='str')
114 114 opObj11.addParameter(name='xaxis', value='velocity', format='str')
115 115 # opObj11.addParameter(name='xmin', value='-10', format='int')
116 116 # opObj11.addParameter(name='xmax', value='10', format='int')
117 117
118 118 # opObj11.addParameter(name='ymin', value='1', format='float')
119 119 # opObj11.addParameter(name='ymax', value='3', format='int')
120 120 #opObj11.addParameter(name='zmin', value='10', format='int')
121 121 #opObj11.addParameter(name='zmax', value='35', format='int')
122 122 # opObj11.addParameter(name='save', value='2', format='int')
123 123 # opObj11.addParameter(name='save', value='5', format='int')
124 124 # opObj11.addParameter(name='figpath', value=figpath, format='str')
125 125
126 126
127 127 opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='other')
128 128 procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList')
129 129 opObj11.addParameter(name='id', value='2005', format='int')
130 130 #opObj11.addParameter(name='wintitle', value='CrossSpectraPlot_ShortPulse', format='str')
131 131 #opObj11.addParameter(name='exp_code', value='13', format='int')
132 132 opObj11.addParameter(name='xaxis', value='Velocity', format='str')
133 133 #opObj11.addParameter(name='xmin', value='-6', format='float')
134 134 #opObj11.addParameter(name='xmax', value='6', format='float')
135 135 opObj11.addParameter(name='zmin', value='15', format='float')
136 136 opObj11.addParameter(name='zmax', value='50', format='float')
137 137 opObj11.addParameter(name='ymin', value='0', format='float')
138 138 opObj11.addParameter(name='ymax', value='7', format='float')
139 139 #opObj11.addParameter(name='phase_min', value='-4', format='int')
140 140 #opObj11.addParameter(name='phase_max', value='4', format='int')
141 141 #
142 142
143 143 # Using external methods (new modules)
144 144 # schainpy.model.graphics.jroplot_spectra.RTIPlot.setup()
145 145 opObj11 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
146 146 opObj11.addParameter(name='id', value='30', format='int')
147 147 opObj11.addParameter(name='wintitle', value='RTI', format='str')
148 148 opObj11.addParameter(name='zmin', value='15', format='int')
149 149 opObj11.addParameter(name='zmax', value='40', format='int')
150 150 opObj11.addParameter(name='ymin', value='1', format='int')
151 151 opObj11.addParameter(name='ymax', value='7', format='int')
152 152 opObj11.addParameter(name='showprofile', value='1', format='int')
153 153 # opObj11.addParameter(name='timerange', value=str(5*60*60*60), format='int')
154 opObj11.addParameter(name='xmin', value='1', format='float')
155 opObj11.addParameter(name='xmax', value='6', format='float')
154 #opObj11.addParameter(name='xmin', value='1', format='float')
155 #opObj11.addParameter(name='xmax', value='6', format='float')
156 156 opObj11.addParameter(name='save', value='1', format='int')
157 opObj11.addParameter(name='figpath', value=figpath, format='str')
157 158
158 159
159 160 # '''#########################################################################################'''
160 161 #
161 162 #
162 163 # procUnitConfObj2 = controllerObj.addProcUnit(datatype='Parameters', inputId=procUnitConfObj1.getId())
163 164 # opObj11 = procUnitConfObj2.addOperation(name='SpectralMoments', optype='other')
164 165 #
165 166 # '''
166 167 # # Discriminacion de ecos
167 168 # opObj11 = procUnitConfObj2.addOperation(name='GaussianFit', optype='other')
168 169 # opObj11.addParameter(name='SNRlimit', value='0', format='int')
169 170 # '''
170 171 #
171 172 # '''
172 173 # # Estimacion de Precipitacion
173 174 # opObj11 = procUnitConfObj2.addOperation(name='PrecipitationProc', optype='other')
174 175 # '''
175 176 #
176 177 # opObj22 = procUnitConfObj2.addOperation(name='FullSpectralAnalysis', optype='other')
177 178 #
178 179 # opObj22.addParameter(name='SNRlimit', value='-10', format='float')
179 180 # opObj22.addParameter(name='E01', value='1.500', format='float')
180 181 # opObj22.addParameter(name='E02', value='1.500', format='float')
181 182 # opObj22.addParameter(name='E12', value='0', format='float')
182 183 # opObj22.addParameter(name='N01', value='0.875', format='float')
183 184 # opObj22.addParameter(name='N02', value='-0.875', format='float')
184 185 # opObj22.addParameter(name='N12', value='-1.750', format='float')
185 186 #
186 187 #
187 188 # opObj22 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
188 189 # opObj22.addParameter(name='id', value='4', format='int')
189 190 # opObj22.addParameter(name='wintitle', value='Wind Profiler', format='str')
190 191 # opObj22.addParameter(name='save', value='1', format='bool')
191 192 # opObj22.addParameter(name='xmin', value='0', format='float')
192 193 # opObj22.addParameter(name='xmax', value='6', format='float')
193 194 # opObj22.addParameter(name='ymin', value='1', format='float')
194 195 # opObj22.addParameter(name='ymax', value='3.5', format='float')
195 196 # opObj22.addParameter(name='zmin', value='-1', format='float')
196 197 # opObj22.addParameter(name='zmax', value='1', format='float')
197 198 # opObj22.addParameter(name='SNRmin', value='-15', format='float')
198 199 # opObj22.addParameter(name='SNRmax', value='20', format='float')
199 200 # opObj22.addParameter(name='zmin_ver', value='-200', format='float')
200 201 # opObj22.addParameter(name='zmax_ver', value='200', format='float')
201 202 # opObj22.addParameter(name='save', value='1', format='int')
202 203 # opObj22.addParameter(name='figpath', value=figpath, format='str')
203 204 #
204 205 #
205 206 #
206 207 # #opObj11.addParameter(name='zmin', value='75', format='int')
207 208 #
208 209 # #opObj12 = procUnitConfObj2.addOperation(name='ParametersPlot', optype='other')
209 210 # #opObj12.addParameter(name='id',value='4',format='int')
210 211 # #opObj12.addParameter(name='wintitle',value='First_gg',format='str')
211 212 # '''
212 213 # #Ploteo de Discriminacion de Gaussianas
213 214 #
214 215 # opObj11 = procUnitConfObj2.addOperation(name='FitGauPlot', optype='other')
215 216 # opObj11.addParameter(name='id', value='21', format='int')
216 217 # opObj11.addParameter(name='wintitle', value='Rainfall Gaussian', format='str')
217 218 # opObj11.addParameter(name='xaxis', value='velocity', format='str')
218 219 # opObj11.addParameter(name='showprofile', value='1', format='int')
219 220 # opObj11.addParameter(name='zmin', value='75', format='int')
220 221 # opObj11.addParameter(name='zmax', value='100', format='int')
221 222 # opObj11.addParameter(name='GauSelector', value='1', format='int')
222 223 # #opObj11.addParameter(name='save', value='1', format='int')
223 224 # #opObj11.addParameter(name='figpath', value='/home/erick/Documents/Data/d2015106')
224 225 #
225 226 # opObj11 = procUnitConfObj2.addOperation(name='FitGauPlot', optype='other')
226 227 # opObj11.addParameter(name='id', value='22', format='int')
227 228 # opObj11.addParameter(name='wintitle', value='Wind Gaussian', format='str')
228 229 # opObj11.addParameter(name='xaxis', value='velocity', format='str')
229 230 # opObj11.addParameter(name='showprofile', value='1', format='int')
230 231 # opObj11.addParameter(name='zmin', value='75', format='int')
231 232 # opObj11.addParameter(name='zmax', value='100', format='int')
232 233 # opObj11.addParameter(name='GauSelector', value='0', format='int')
233 234 # '''
234 235 #
235 236 #
236 237
237 238
238 239 controllerObj.start()
General Comments 0
You need to be logged in to leave comments. Login now