##// END OF EJS Templates
schain-jc
RSS
BLTR ok
Juan C. Espinoza -
r1018:f15cf5631d3c
parent child
Show More
Side by Side Unified
Context file:
r1018:f15cf5631d3c Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1229 +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 640 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
641 641 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
642 642
643 643 return velrange
644 644
645 645 def getNPairs(self):
646 646
647 647 return len(self.pairsList)
648 648
649 649 def getPairsIndexList(self):
650 650
651 651 return range(self.nPairs)
652 652
653 653 def getNormFactor(self):
654 654
655 655 pwcode = 1
656 656
657 657 if self.flagDecodeData:
658 658 pwcode = numpy.sum(self.code[0]**2)
659 659 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
660 660 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
661 661
662 662 return normFactor
663 663
664 664 def getFlagCspc(self):
665 665
666 666 if self.data_cspc is None:
667 667 return True
668 668
669 669 return False
670 670
671 671 def getFlagDc(self):
672 672
673 673 if self.data_dc is None:
674 674 return True
675 675
676 676 return False
677 677
678 678 def getTimeInterval(self):
679 679
680 680 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
681 681
682 682 return timeInterval
683 683
684 684 def getPower(self):
685 685
686 686 factor = self.normFactor
687 687 z = self.data_spc/factor
688 688 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
689 689 avg = numpy.average(z, axis=1)
690 690
691 691 return 10*numpy.log10(avg)
692 692
693 693 def getCoherence(self, pairsList=None, phase=False):
694 694
695 695 z = []
696 696 if pairsList is None:
697 697 pairsIndexList = self.pairsIndexList
698 698 else:
699 699 pairsIndexList = []
700 700 for pair in pairsList:
701 701 if pair not in self.pairsList:
702 702 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
703 703 pairsIndexList.append(self.pairsList.index(pair))
704 704 for i in range(len(pairsIndexList)):
705 705 pair = self.pairsList[pairsIndexList[i]]
706 706 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
707 707 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
708 708 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
709 709 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
710 710 if phase:
711 711 data = numpy.arctan2(avgcoherenceComplex.imag,
712 712 avgcoherenceComplex.real)*180/numpy.pi
713 713 else:
714 714 data = numpy.abs(avgcoherenceComplex)
715 715
716 716 z.append(data)
717 717
718 718 return numpy.array(z)
719 719
720 720 def setValue(self, value):
721 721
722 722 print "This property should not be initialized"
723 723
724 724 return
725 725
726 726 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
727 727 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
728 728 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
729 729 flag_cspc = property(getFlagCspc, setValue)
730 730 flag_dc = property(getFlagDc, setValue)
731 731 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
732 732 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
733 733
734 734 class SpectraHeis(Spectra):
735 735
736 736 data_spc = None
737 737
738 738 data_cspc = None
739 739
740 740 data_dc = None
741 741
742 742 nFFTPoints = None
743 743
744 744 # nPairs = None
745 745
746 746 pairsList = None
747 747
748 748 nCohInt = None
749 749
750 750 nIncohInt = None
751 751
752 752 def __init__(self):
753 753
754 754 self.radarControllerHeaderObj = RadarControllerHeader()
755 755
756 756 self.systemHeaderObj = SystemHeader()
757 757
758 758 self.type = "SpectraHeis"
759 759
760 760 # self.dtype = None
761 761
762 762 # self.nChannels = 0
763 763
764 764 # self.nHeights = 0
765 765
766 766 self.nProfiles = None
767 767
768 768 self.heightList = None
769 769
770 770 self.channelList = None
771 771
772 772 # self.channelIndexList = None
773 773
774 774 self.flagNoData = True
775 775
776 776 self.flagDiscontinuousBlock = False
777 777
778 778 # self.nPairs = 0
779 779
780 780 self.utctime = None
781 781
782 782 self.blocksize = None
783 783
784 784 self.profileIndex = 0
785 785
786 786 self.nCohInt = 1
787 787
788 788 self.nIncohInt = 1
789 789
790 790 def getNormFactor(self):
791 791 pwcode = 1
792 792 if self.flagDecodeData:
793 793 pwcode = numpy.sum(self.code[0]**2)
794 794
795 795 normFactor = self.nIncohInt*self.nCohInt*pwcode
796 796
797 797 return normFactor
798 798
799 799 def getTimeInterval(self):
800 800
801 801 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
802 802
803 803 return timeInterval
804 804
805 805 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
806 806 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
807 807
808 808 class Fits(JROData):
809 809
810 810 heightList = None
811 811
812 812 channelList = None
813 813
814 814 flagNoData = True
815 815
816 816 flagDiscontinuousBlock = False
817 817
818 818 useLocalTime = False
819 819
820 820 utctime = None
821 821
822 822 timeZone = None
823 823
824 824 # ippSeconds = None
825 825
826 826 # timeInterval = None
827 827
828 828 nCohInt = None
829 829
830 830 nIncohInt = None
831 831
832 832 noise = None
833 833
834 834 windowOfFilter = 1
835 835
836 836 #Speed of ligth
837 837 C = 3e8
838 838
839 839 frequency = 49.92e6
840 840
841 841 realtime = False
842 842
843 843
844 844 def __init__(self):
845 845
846 846 self.type = "Fits"
847 847
848 848 self.nProfiles = None
849 849
850 850 self.heightList = None
851 851
852 852 self.channelList = None
853 853
854 854 # self.channelIndexList = None
855 855
856 856 self.flagNoData = True
857 857
858 858 self.utctime = None
859 859
860 860 self.nCohInt = 1
861 861
862 862 self.nIncohInt = 1
863 863
864 864 self.useLocalTime = True
865 865
866 866 self.profileIndex = 0
867 867
868 868 # self.utctime = None
869 869 # self.timeZone = None
870 870 # self.ltctime = None
871 871 # self.timeInterval = None
872 872 # self.header = None
873 873 # self.data_header = None
874 874 # self.data = None
875 875 # self.datatime = None
876 876 # self.flagNoData = False
877 877 # self.expName = ''
878 878 # self.nChannels = None
879 879 # self.nSamples = None
880 880 # self.dataBlocksPerFile = None
881 881 # self.comments = ''
882 882 #
883 883
884 884
885 885 def getltctime(self):
886 886
887 887 if self.useLocalTime:
888 888 return self.utctime - self.timeZone*60
889 889
890 890 return self.utctime
891 891
892 892 def getDatatime(self):
893 893
894 894 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
895 895 return datatime
896 896
897 897 def getTimeRange(self):
898 898
899 899 datatime = []
900 900
901 901 datatime.append(self.ltctime)
902 902 datatime.append(self.ltctime + self.timeInterval)
903 903
904 904 datatime = numpy.array(datatime)
905 905
906 906 return datatime
907 907
908 908 def getHeiRange(self):
909 909
910 910 heis = self.heightList
911 911
912 912 return heis
913 913
914 914 def getNHeights(self):
915 915
916 916 return len(self.heightList)
917 917
918 918 def getNChannels(self):
919 919
920 920 return len(self.channelList)
921 921
922 922 def getChannelIndexList(self):
923 923
924 924 return range(self.nChannels)
925 925
926 926 def getNoise(self, type = 1):
927 927
928 928 #noise = numpy.zeros(self.nChannels)
929 929
930 930 if type == 1:
931 931 noise = self.getNoisebyHildebrand()
932 932
933 933 if type == 2:
934 934 noise = self.getNoisebySort()
935 935
936 936 if type == 3:
937 937 noise = self.getNoisebyWindow()
938 938
939 939 return noise
940 940
941 941 def getTimeInterval(self):
942 942
943 943 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
944 944
945 945 return timeInterval
946 946
947 947 datatime = property(getDatatime, "I'm the 'datatime' property")
948 948 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
949 949 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
950 950 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
951 951 noise = property(getNoise, "I'm the 'nHeights' property.")
952 952
953 953 ltctime = property(getltctime, "I'm the 'ltctime' property")
954 954 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
955 955
956 956
957 957 class Correlation(JROData):
958 958
959 959 noise = None
960 960
961 961 SNR = None
962 962
963 963 #--------------------------------------------------
964 964
965 965 mode = None
966 966
967 967 split = False
968 968
969 969 data_cf = None
970 970
971 971 lags = None
972 972
973 973 lagRange = None
974 974
975 975 pairsList = None
976 976
977 977 normFactor = None
978 978
979 979 #--------------------------------------------------
980 980
981 981 # calculateVelocity = None
982 982
983 983 nLags = None
984 984
985 985 nPairs = None
986 986
987 987 nAvg = None
988 988
989 989
990 990 def __init__(self):
991 991 '''
992 992 Constructor
993 993 '''
994 994 self.radarControllerHeaderObj = RadarControllerHeader()
995 995
996 996 self.systemHeaderObj = SystemHeader()
997 997
998 998 self.type = "Correlation"
999 999
1000 1000 self.data = None
1001 1001
1002 1002 self.dtype = None
1003 1003
1004 1004 self.nProfiles = None
1005 1005
1006 1006 self.heightList = None
1007 1007
1008 1008 self.channelList = None
1009 1009
1010 1010 self.flagNoData = True
1011 1011
1012 1012 self.flagDiscontinuousBlock = False
1013 1013
1014 1014 self.utctime = None
1015 1015
1016 1016 self.timeZone = None
1017 1017
1018 1018 self.dstFlag = None
1019 1019
1020 1020 self.errorCount = None
1021 1021
1022 1022 self.blocksize = None
1023 1023
1024 1024 self.flagDecodeData = False #asumo q la data no esta decodificada
1025 1025
1026 1026 self.flagDeflipData = False #asumo q la data no esta sin flip
1027 1027
1028 1028 self.pairsList = None
1029 1029
1030 1030 self.nPoints = None
1031 1031
1032 1032 def getPairsList(self):
1033 1033
1034 1034 return self.pairsList
1035 1035
1036 1036 def getNoise(self, mode = 2):
1037 1037
1038 1038 indR = numpy.where(self.lagR == 0)[0][0]
1039 1039 indT = numpy.where(self.lagT == 0)[0][0]
1040 1040
1041 1041 jspectra0 = self.data_corr[:,:,indR,:]
1042 1042 jspectra = copy.copy(jspectra0)
1043 1043
1044 1044 num_chan = jspectra.shape[0]
1045 1045 num_hei = jspectra.shape[2]
1046 1046
1047 1047 freq_dc = jspectra.shape[1]/2
1048 1048 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1049 1049
1050 1050 if ind_vel[0]<0:
1051 1051 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1052 1052
1053 1053 if mode == 1:
1054 1054 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1055 1055
1056 1056 if mode == 2:
1057 1057
1058 1058 vel = numpy.array([-2,-1,1,2])
1059 1059 xx = numpy.zeros([4,4])
1060 1060
1061 1061 for fil in range(4):
1062 1062 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1063 1063
1064 1064 xx_inv = numpy.linalg.inv(xx)
1065 1065 xx_aux = xx_inv[0,:]
1066 1066
1067 1067 for ich in range(num_chan):
1068 1068 yy = jspectra[ich,ind_vel,:]
1069 1069 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1070 1070
1071 1071 junkid = jspectra[ich,freq_dc,:]<=0
1072 1072 cjunkid = sum(junkid)
1073 1073
1074 1074 if cjunkid.any():
1075 1075 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1076 1076
1077 1077 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1078 1078
1079 1079 return noise
1080 1080
1081 1081 def getTimeInterval(self):
1082 1082
1083 1083 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1084 1084
1085 1085 return timeInterval
1086 1086
1087 1087 def splitFunctions(self):
1088 1088
1089 1089 pairsList = self.pairsList
1090 1090 ccf_pairs = []
1091 1091 acf_pairs = []
1092 1092 ccf_ind = []
1093 1093 acf_ind = []
1094 1094 for l in range(len(pairsList)):
1095 1095 chan0 = pairsList[l][0]
1096 1096 chan1 = pairsList[l][1]
1097 1097
1098 1098 #Obteniendo pares de Autocorrelacion
1099 1099 if chan0 == chan1:
1100 1100 acf_pairs.append(chan0)
1101 1101 acf_ind.append(l)
1102 1102 else:
1103 1103 ccf_pairs.append(pairsList[l])
1104 1104 ccf_ind.append(l)
1105 1105
1106 1106 data_acf = self.data_cf[acf_ind]
1107 1107 data_ccf = self.data_cf[ccf_ind]
1108 1108
1109 1109 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1110 1110
1111 1111 def getNormFactor(self):
1112 1112 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1113 1113 acf_pairs = numpy.array(acf_pairs)
1114 1114 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1115 1115
1116 1116 for p in range(self.nPairs):
1117 1117 pair = self.pairsList[p]
1118 1118
1119 1119 ch0 = pair[0]
1120 1120 ch1 = pair[1]
1121 1121
1122 1122 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1123 1123 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1124 1124 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1125 1125
1126 1126 return normFactor
1127 1127
1128 1128 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1129 1129 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1130 1130
1131 1131 class Parameters(Spectra):
1132 1132
1133 1133 experimentInfo = None #Information about the experiment
1134 1134
1135 1135 #Information from previous data
1136 1136
1137 1137 inputUnit = None #Type of data to be processed
1138 1138
1139 1139 operation = None #Type of operation to parametrize
1140 1140
1141 1141 #normFactor = None #Normalization Factor
1142 1142
1143 1143 groupList = None #List of Pairs, Groups, etc
1144 1144
1145 1145 #Parameters
1146 1146
1147 1147 data_param = None #Parameters obtained
1148 1148
1149 1149 data_pre = None #Data Pre Parametrization
1150 1150
1151 1151 data_SNR = None #Signal to Noise Ratio
1152 1152
1153 1153 # heightRange = None #Heights
1154 1154
1155 1155 abscissaList = None #Abscissa, can be velocities, lags or time
1156 1156
1157 1157 # noise = None #Noise Potency
1158 1158
1159 1159 utctimeInit = None #Initial UTC time
1160 1160
1161 1161 paramInterval = None #Time interval to calculate Parameters in seconds
1162 1162
1163 1163 useLocalTime = True
1164 1164
1165 1165 #Fitting
1166 1166
1167 1167 data_error = None #Error of the estimation
1168 1168
1169 1169 constants = None
1170 1170
1171 1171 library = None
1172 1172
1173 1173 #Output signal
1174 1174
1175 1175 outputInterval = None #Time interval to calculate output signal in seconds
1176 1176
1177 1177 data_output = None #Out signal
1178 1178
1179 1179 nAvg = None
1180 1180
1181 1181 noise_estimation = None
1182 1182
1183 1183 GauSPC = None #Fit gaussian SPC
1184 1184
1185 1185
1186 1186 def __init__(self):
1187 1187 '''
1188 1188 Constructor
1189 1189 '''
1190 1190 self.radarControllerHeaderObj = RadarControllerHeader()
1191 1191
1192 1192 self.systemHeaderObj = SystemHeader()
1193 1193
1194 1194 self.type = "Parameters"
1195 1195
1196 1196 def getTimeRange1(self, interval):
1197 1197
1198 1198 datatime = []
1199 1199
1200 1200 if self.useLocalTime:
1201 1201 time1 = self.utctimeInit - self.timeZone*60
1202 1202 else:
1203 1203 time1 = self.utctimeInit
1204 print 'interval',interval
1204
1205 1205 datatime.append(time1)
1206 1206 datatime.append(time1 + interval)
1207 1207 datatime = numpy.array(datatime)
1208 1208
1209 1209 return datatime
1210 1210
1211 1211 def getTimeInterval(self):
1212 1212
1213 1213 if hasattr(self, 'timeInterval1'):
1214 1214 return self.timeInterval1
1215 1215 else:
1216 1216 return self.paramInterval
1217 1217
1218 1218 def setValue(self, value):
1219 1219
1220 1220 print "This property should not be initialized"
1221 1221
1222 1222 return
1223 1223
1224 1224 def getNoise(self):
1225 1225
1226 1226 return self.spc_noise
1227 1227
1228 1228 timeInterval = property(getTimeInterval)
1229 1229 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file, binary diff hidden
Show file before | Show file after | Hide comments
@@ -1,2159 +1,2154
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 9 class FitGauPlot(Figure):
10 10
11 11 isConfig = None
12 12 __nsubplots = None
13 13
14 14 WIDTHPROF = None
15 15 HEIGHTPROF = None
16 16 PREFIX = 'fitgau'
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 85 xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 1):
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 121 x = dataOut.spc_range[0]
122 122 xlabel = "Frequency (kHz)"
123 123
124 124 elif xaxis == "time":
125 125 x = dataOut.spc_range[1]
126 126 xlabel = "Time (ms)"
127 127
128 128 else:
129 129 x = dataOut.spc_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 136 z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor
137 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 y = dataOut.heightList
662 z = dataOut.data_output.copy()
663 print ' '
664 print 'Xvel',z[0]
665 print ' '
666 print 'Yvel',z[1]
667 print ' '
661 y = dataOut.heightList
662 z = dataOut.data_output.copy()
668 663 nplots = z.shape[0] #Number of wind dimensions estimated
669 664 nplotsw = nplots
670 665
671 666
672 667 #If there is a SNR function defined
673 668 if dataOut.data_SNR is not None:
674 669 nplots += 1
675 670 SNR = dataOut.data_SNR
676 671 SNRavg = numpy.average(SNR, axis=0)
677 672
678 673 SNRdB = 10*numpy.log10(SNR)
679 674 SNRavgdB = 10*numpy.log10(SNRavg)
680 675
681 676 if SNRthresh == None: SNRthresh = -5.0
682 677 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
683 678
684 679 for i in range(nplotsw):
685 680 z[i,ind] = numpy.nan
686 681
687 682 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
688 683 #thisDatetime = datetime.datetime.now()
689 684 title = wintitle + "Wind"
690 685 xlabel = ""
691 686 ylabel = "Height (km)"
692 687 update_figfile = False
693 688
694 689 if not self.isConfig:
695 690
696 691 self.setup(id=id,
697 692 nplots=nplots,
698 693 wintitle=wintitle,
699 694 showprofile=showprofile,
700 695 show=show)
701 696
702 697 if timerange is not None:
703 698 self.timerange = timerange
704 699
705 700 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
706 701
707 702 if ymin == None: ymin = numpy.nanmin(y)
708 703 if ymax == None: ymax = numpy.nanmax(y)
709 704
710 705 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
711 706 #if numpy.isnan(zmax): zmax = 50
712 707 if zmin == None: zmin = -zmax
713 708
714 709 if nplotsw == 3:
715 710 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
716 711 if zmin_ver == None: zmin_ver = -zmax_ver
717 712
718 713 if dataOut.data_SNR is not None:
719 714 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
720 715 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
721 716
722 717
723 718 self.FTP_WEI = ftp_wei
724 719 self.EXP_CODE = exp_code
725 720 self.SUB_EXP_CODE = sub_exp_code
726 721 self.PLOT_POS = plot_pos
727 722
728 723 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
729 724 self.isConfig = True
730 725 self.figfile = figfile
731 726 update_figfile = True
732 727
733 728 self.setWinTitle(title)
734 729
735 730 if ((self.xmax - x[1]) < (x[1]-x[0])):
736 731 x[1] = self.xmax
737 732
738 733 strWind = ['Zonal', 'Meridional', 'Vertical']
739 734 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
740 735 zmaxVector = [zmax, zmax, zmax_ver]
741 736 zminVector = [zmin, zmin, zmin_ver]
742 737 windFactor = [1,1,100]
743 738
744 739 for i in range(nplotsw):
745 740
746 741 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
747 742 axes = self.axesList[i*self.__nsubplots]
748 743
749 744 z1 = z[i,:].reshape((1,-1))*windFactor[i]
750 745 #z1=numpy.ma.masked_where(z1==0.,z1)
751 746
752 747 axes.pcolorbuffer(x, y, z1,
753 748 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
754 749 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
755 750 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="seismic" )
756 751
757 752 if dataOut.data_SNR is not None:
758 753 i += 1
759 754 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
760 755 axes = self.axesList[i*self.__nsubplots]
761 756 SNRavgdB = SNRavgdB.reshape((1,-1))
762 757 axes.pcolorbuffer(x, y, SNRavgdB,
763 758 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
764 759 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
765 760 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
766 761
767 762 self.draw()
768 763
769 764 self.save(figpath=figpath,
770 765 figfile=figfile,
771 766 save=save,
772 767 ftp=ftp,
773 768 wr_period=wr_period,
774 769 thisDatetime=thisDatetime,
775 770 update_figfile=update_figfile)
776 771
777 772 if dataOut.ltctime + dataOut.paramInterval >= self.xmax:
778 773 self.counter_imagwr = wr_period
779 774 self.isConfig = False
780 775 update_figfile = True
781 776
782 777
783 778 class ParametersPlot(Figure):
784 779
785 780 __isConfig = None
786 781 __nsubplots = None
787 782
788 783 WIDTHPROF = None
789 784 HEIGHTPROF = None
790 785 PREFIX = 'param'
791 786
792 787 nplots = None
793 788 nchan = None
794 789
795 790 def __init__(self, **kwargs):
796 791 Figure.__init__(self, **kwargs)
797 792 self.timerange = None
798 793 self.isConfig = False
799 794 self.__nsubplots = 1
800 795
801 796 self.WIDTH = 800
802 797 self.HEIGHT = 180
803 798 self.WIDTHPROF = 120
804 799 self.HEIGHTPROF = 0
805 800 self.counter_imagwr = 0
806 801
807 802 self.PLOT_CODE = RTI_CODE
808 803
809 804 self.FTP_WEI = None
810 805 self.EXP_CODE = None
811 806 self.SUB_EXP_CODE = None
812 807 self.PLOT_POS = None
813 808 self.tmin = None
814 809 self.tmax = None
815 810
816 811 self.xmin = None
817 812 self.xmax = None
818 813
819 814 self.figfile = None
820 815
821 816 def getSubplots(self):
822 817
823 818 ncol = 1
824 819 nrow = self.nplots
825 820
826 821 return nrow, ncol
827 822
828 823 def setup(self, id, nplots, wintitle, show=True):
829 824
830 825 self.nplots = nplots
831 826
832 827 ncolspan = 1
833 828 colspan = 1
834 829
835 830 self.createFigure(id = id,
836 831 wintitle = wintitle,
837 832 widthplot = self.WIDTH + self.WIDTHPROF,
838 833 heightplot = self.HEIGHT + self.HEIGHTPROF,
839 834 show=show)
840 835
841 836 nrow, ncol = self.getSubplots()
842 837
843 838 counter = 0
844 839 for y in range(nrow):
845 840 for x in range(ncol):
846 841
847 842 if counter >= self.nplots:
848 843 break
849 844
850 845 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
851 846
852 847 counter += 1
853 848
854 849 def run(self, dataOut, id, wintitle="", channelList=None, paramIndex = 0, colormap="jet",
855 850 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, timerange=None,
856 851 showSNR=False, SNRthresh = -numpy.inf, SNRmin=None, SNRmax=None,
857 852 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
858 853 server=None, folder=None, username=None, password=None,
859 854 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, HEIGHT=None):
860 855 """
861 856
862 857 Input:
863 858 dataOut :
864 859 id :
865 860 wintitle :
866 861 channelList :
867 862 showProfile :
868 863 xmin : None,
869 864 xmax : None,
870 865 ymin : None,
871 866 ymax : None,
872 867 zmin : None,
873 868 zmax : None
874 869 """
875 870
876 871 if HEIGHT is not None:
877 872 self.HEIGHT = HEIGHT
878 873
879 874
880 875 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
881 876 return
882 877
883 878 if channelList == None:
884 879 channelIndexList = range(dataOut.data_param.shape[0])
885 880 else:
886 881 channelIndexList = []
887 882 for channel in channelList:
888 883 if channel not in dataOut.channelList:
889 884 raise ValueError, "Channel %d is not in dataOut.channelList"
890 885 channelIndexList.append(dataOut.channelList.index(channel))
891 886
892 887 x = dataOut.getTimeRange1(dataOut.paramInterval)
893 888 y = dataOut.getHeiRange()
894 889
895 890 if dataOut.data_param.ndim == 3:
896 891 z = dataOut.data_param[channelIndexList,paramIndex,:]
897 892 else:
898 893 z = dataOut.data_param[channelIndexList,:]
899 894
900 895 if showSNR:
901 896 #SNR data
902 897 SNRarray = dataOut.data_SNR[channelIndexList,:]
903 898 SNRdB = 10*numpy.log10(SNRarray)
904 899 ind = numpy.where(SNRdB < SNRthresh)
905 900 z[ind] = numpy.nan
906 901
907 902 thisDatetime = dataOut.datatime
908 903 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
909 904 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
910 905 xlabel = ""
911 906 ylabel = "Range (Km)"
912 907
913 908 update_figfile = False
914 909
915 910 if not self.isConfig:
916 911
917 912 nchan = len(channelIndexList)
918 913 self.nchan = nchan
919 914 self.plotFact = 1
920 915 nplots = nchan
921 916
922 917 if showSNR:
923 918 nplots = nchan*2
924 919 self.plotFact = 2
925 920 if SNRmin == None: SNRmin = numpy.nanmin(SNRdB)
926 921 if SNRmax == None: SNRmax = numpy.nanmax(SNRdB)
927 922
928 923 self.setup(id=id,
929 924 nplots=nplots,
930 925 wintitle=wintitle,
931 926 show=show)
932 927
933 928 if timerange != None:
934 929 self.timerange = timerange
935 930
936 931 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
937 932
938 933 if ymin == None: ymin = numpy.nanmin(y)
939 934 if ymax == None: ymax = numpy.nanmax(y)
940 935 if zmin == None: zmin = numpy.nanmin(z)
941 936 if zmax == None: zmax = numpy.nanmax(z)
942 937
943 938 self.FTP_WEI = ftp_wei
944 939 self.EXP_CODE = exp_code
945 940 self.SUB_EXP_CODE = sub_exp_code
946 941 self.PLOT_POS = plot_pos
947 942
948 943 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
949 944 self.isConfig = True
950 945 self.figfile = figfile
951 946 update_figfile = True
952 947
953 948 self.setWinTitle(title)
954 949
955 950 for i in range(self.nchan):
956 951 index = channelIndexList[i]
957 952 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
958 953 axes = self.axesList[i*self.plotFact]
959 954 z1 = z[i,:].reshape((1,-1))
960 955 axes.pcolorbuffer(x, y, z1,
961 956 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
962 957 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
963 958 ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
964 959
965 960 if showSNR:
966 961 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
967 962 axes = self.axesList[i*self.plotFact + 1]
968 963 SNRdB1 = SNRdB[i,:].reshape((1,-1))
969 964 axes.pcolorbuffer(x, y, SNRdB1,
970 965 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
971 966 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
972 967 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
973 968
974 969
975 970 self.draw()
976 971
977 972 if dataOut.ltctime >= self.xmax:
978 973 self.counter_imagwr = wr_period
979 974 self.isConfig = False
980 975 update_figfile = True
981 976
982 977 self.save(figpath=figpath,
983 978 figfile=figfile,
984 979 save=save,
985 980 ftp=ftp,
986 981 wr_period=wr_period,
987 982 thisDatetime=thisDatetime,
988 983 update_figfile=update_figfile)
989 984
990 985
991 986
992 987 class Parameters1Plot(Figure):
993 988
994 989 __isConfig = None
995 990 __nsubplots = None
996 991
997 992 WIDTHPROF = None
998 993 HEIGHTPROF = None
999 994 PREFIX = 'prm'
1000 995
1001 996 def __init__(self, **kwargs):
1002 997 Figure.__init__(self, **kwargs)
1003 998 self.timerange = 2*60*60
1004 999 self.isConfig = False
1005 1000 self.__nsubplots = 1
1006 1001
1007 1002 self.WIDTH = 800
1008 1003 self.HEIGHT = 180
1009 1004 self.WIDTHPROF = 120
1010 1005 self.HEIGHTPROF = 0
1011 1006 self.counter_imagwr = 0
1012 1007
1013 1008 self.PLOT_CODE = PARMS_CODE
1014 1009
1015 1010 self.FTP_WEI = None
1016 1011 self.EXP_CODE = None
1017 1012 self.SUB_EXP_CODE = None
1018 1013 self.PLOT_POS = None
1019 1014 self.tmin = None
1020 1015 self.tmax = None
1021 1016
1022 1017 self.xmin = None
1023 1018 self.xmax = None
1024 1019
1025 1020 self.figfile = None
1026 1021
1027 1022 def getSubplots(self):
1028 1023
1029 1024 ncol = 1
1030 1025 nrow = self.nplots
1031 1026
1032 1027 return nrow, ncol
1033 1028
1034 1029 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1035 1030
1036 1031 self.__showprofile = showprofile
1037 1032 self.nplots = nplots
1038 1033
1039 1034 ncolspan = 1
1040 1035 colspan = 1
1041 1036
1042 1037 self.createFigure(id = id,
1043 1038 wintitle = wintitle,
1044 1039 widthplot = self.WIDTH + self.WIDTHPROF,
1045 1040 heightplot = self.HEIGHT + self.HEIGHTPROF,
1046 1041 show=show)
1047 1042
1048 1043 nrow, ncol = self.getSubplots()
1049 1044
1050 1045 counter = 0
1051 1046 for y in range(nrow):
1052 1047 for x in range(ncol):
1053 1048
1054 1049 if counter >= self.nplots:
1055 1050 break
1056 1051
1057 1052 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1058 1053
1059 1054 if showprofile:
1060 1055 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
1061 1056
1062 1057 counter += 1
1063 1058
1064 1059 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
1065 1060 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
1066 1061 parameterIndex = None, onlyPositive = False,
1067 1062 SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, onlySNR = False,
1068 1063 DOP = True,
1069 1064 zlabel = "", parameterName = "", parameterObject = "data_param",
1070 1065 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1071 1066 server=None, folder=None, username=None, password=None,
1072 1067 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1073 1068 #print inspect.getargspec(self.run).args
1074 1069 """
1075 1070
1076 1071 Input:
1077 1072 dataOut :
1078 1073 id :
1079 1074 wintitle :
1080 1075 channelList :
1081 1076 showProfile :
1082 1077 xmin : None,
1083 1078 xmax : None,
1084 1079 ymin : None,
1085 1080 ymax : None,
1086 1081 zmin : None,
1087 1082 zmax : None
1088 1083 """
1089 1084
1090 1085 data_param = getattr(dataOut, parameterObject)
1091 1086
1092 1087 if channelList == None:
1093 1088 channelIndexList = numpy.arange(data_param.shape[0])
1094 1089 else:
1095 1090 channelIndexList = numpy.array(channelList)
1096 1091
1097 1092 nchan = len(channelIndexList) #Number of channels being plotted
1098 1093
1099 1094 if nchan < 1:
1100 1095 return
1101 1096
1102 1097 nGraphsByChannel = 0
1103 1098
1104 1099 if SNR:
1105 1100 nGraphsByChannel += 1
1106 1101 if DOP:
1107 1102 nGraphsByChannel += 1
1108 1103
1109 1104 if nGraphsByChannel < 1:
1110 1105 return
1111 1106
1112 1107 nplots = nGraphsByChannel*nchan
1113 1108
1114 1109 if timerange is not None:
1115 1110 self.timerange = timerange
1116 1111
1117 1112 #tmin = None
1118 1113 #tmax = None
1119 1114 if parameterIndex == None:
1120 1115 parameterIndex = 1
1121 1116
1122 1117 x = dataOut.getTimeRange1(dataOut.paramInterval)
1123 1118 y = dataOut.heightList
1124 1119 z = data_param[channelIndexList,parameterIndex,:].copy()
1125 1120
1126 1121 zRange = dataOut.abscissaList
1127 1122 # nChannels = z.shape[0] #Number of wind dimensions estimated
1128 1123 # thisDatetime = dataOut.datatime
1129 1124
1130 1125 if dataOut.data_SNR is not None:
1131 1126 SNRarray = dataOut.data_SNR[channelIndexList,:]
1132 1127 SNRdB = 10*numpy.log10(SNRarray)
1133 1128 # SNRavgdB = 10*numpy.log10(SNRavg)
1134 1129 ind = numpy.where(SNRdB < 10**(SNRthresh/10))
1135 1130 z[ind] = numpy.nan
1136 1131
1137 1132 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1138 1133 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
1139 1134 xlabel = ""
1140 1135 ylabel = "Range (Km)"
1141 1136
1142 1137 if (SNR and not onlySNR): nplots = 2*nplots
1143 1138
1144 1139 if onlyPositive:
1145 1140 colormap = "jet"
1146 1141 zmin = 0
1147 1142 else: colormap = "RdBu_r"
1148 1143
1149 1144 if not self.isConfig:
1150 1145
1151 1146 self.setup(id=id,
1152 1147 nplots=nplots,
1153 1148 wintitle=wintitle,
1154 1149 showprofile=showprofile,
1155 1150 show=show)
1156 1151
1157 1152 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1158 1153
1159 1154 if ymin == None: ymin = numpy.nanmin(y)
1160 1155 if ymax == None: ymax = numpy.nanmax(y)
1161 1156 if zmin == None: zmin = numpy.nanmin(zRange)
1162 1157 if zmax == None: zmax = numpy.nanmax(zRange)
1163 1158
1164 1159 if SNR:
1165 1160 if SNRmin == None: SNRmin = numpy.nanmin(SNRdB)
1166 1161 if SNRmax == None: SNRmax = numpy.nanmax(SNRdB)
1167 1162
1168 1163 self.FTP_WEI = ftp_wei
1169 1164 self.EXP_CODE = exp_code
1170 1165 self.SUB_EXP_CODE = sub_exp_code
1171 1166 self.PLOT_POS = plot_pos
1172 1167
1173 1168 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1174 1169 self.isConfig = True
1175 1170 self.figfile = figfile
1176 1171
1177 1172 self.setWinTitle(title)
1178 1173
1179 1174 if ((self.xmax - x[1]) < (x[1]-x[0])):
1180 1175 x[1] = self.xmax
1181 1176
1182 1177 for i in range(nchan):
1183 1178
1184 1179 if (SNR and not onlySNR): j = 2*i
1185 1180 else: j = i
1186 1181
1187 1182 j = nGraphsByChannel*i
1188 1183
1189 1184 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
1190 1185 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
1191 1186
1192 1187 if not onlySNR:
1193 1188 axes = self.axesList[j*self.__nsubplots]
1194 1189 z1 = z[i,:].reshape((1,-1))
1195 1190 axes.pcolorbuffer(x, y, z1,
1196 1191 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1197 1192 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
1198 1193 ticksize=9, cblabel=zlabel, cbsize="1%")
1199 1194
1200 1195 if DOP:
1201 1196 title = "%s Channel %d: %s" %(parameterName, channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1202 1197
1203 1198 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
1204 1199 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
1205 1200 axes = self.axesList[j]
1206 1201 z1 = z[i,:].reshape((1,-1))
1207 1202 axes.pcolorbuffer(x, y, z1,
1208 1203 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1209 1204 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
1210 1205 ticksize=9, cblabel=zlabel, cbsize="1%")
1211 1206
1212 1207 if SNR:
1213 1208 title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1214 1209 axes = self.axesList[(j)*self.__nsubplots]
1215 1210 if not onlySNR:
1216 1211 axes = self.axesList[(j + 1)*self.__nsubplots]
1217 1212
1218 1213 axes = self.axesList[(j + nGraphsByChannel-1)]
1219 1214
1220 1215 z1 = SNRdB[i,:].reshape((1,-1))
1221 1216 axes.pcolorbuffer(x, y, z1,
1222 1217 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1223 1218 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap="jet",
1224 1219 ticksize=9, cblabel=zlabel, cbsize="1%")
1225 1220
1226 1221
1227 1222
1228 1223 self.draw()
1229 1224
1230 1225 if x[1] >= self.axesList[0].xmax:
1231 1226 self.counter_imagwr = wr_period
1232 1227 self.isConfig = False
1233 1228 self.figfile = None
1234 1229
1235 1230 self.save(figpath=figpath,
1236 1231 figfile=figfile,
1237 1232 save=save,
1238 1233 ftp=ftp,
1239 1234 wr_period=wr_period,
1240 1235 thisDatetime=thisDatetime,
1241 1236 update_figfile=False)
1242 1237
1243 1238 class SpectralFittingPlot(Figure):
1244 1239
1245 1240 __isConfig = None
1246 1241 __nsubplots = None
1247 1242
1248 1243 WIDTHPROF = None
1249 1244 HEIGHTPROF = None
1250 1245 PREFIX = 'prm'
1251 1246
1252 1247
1253 1248 N = None
1254 1249 ippSeconds = None
1255 1250
1256 1251 def __init__(self, **kwargs):
1257 1252 Figure.__init__(self, **kwargs)
1258 1253 self.isConfig = False
1259 1254 self.__nsubplots = 1
1260 1255
1261 1256 self.PLOT_CODE = SPECFIT_CODE
1262 1257
1263 1258 self.WIDTH = 450
1264 1259 self.HEIGHT = 250
1265 1260 self.WIDTHPROF = 0
1266 1261 self.HEIGHTPROF = 0
1267 1262
1268 1263 def getSubplots(self):
1269 1264
1270 1265 ncol = int(numpy.sqrt(self.nplots)+0.9)
1271 1266 nrow = int(self.nplots*1./ncol + 0.9)
1272 1267
1273 1268 return nrow, ncol
1274 1269
1275 1270 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
1276 1271
1277 1272 showprofile = False
1278 1273 self.__showprofile = showprofile
1279 1274 self.nplots = nplots
1280 1275
1281 1276 ncolspan = 5
1282 1277 colspan = 4
1283 1278 if showprofile:
1284 1279 ncolspan = 5
1285 1280 colspan = 4
1286 1281 self.__nsubplots = 2
1287 1282
1288 1283 self.createFigure(id = id,
1289 1284 wintitle = wintitle,
1290 1285 widthplot = self.WIDTH + self.WIDTHPROF,
1291 1286 heightplot = self.HEIGHT + self.HEIGHTPROF,
1292 1287 show=show)
1293 1288
1294 1289 nrow, ncol = self.getSubplots()
1295 1290
1296 1291 counter = 0
1297 1292 for y in range(nrow):
1298 1293 for x in range(ncol):
1299 1294
1300 1295 if counter >= self.nplots:
1301 1296 break
1302 1297
1303 1298 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1304 1299
1305 1300 if showprofile:
1306 1301 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
1307 1302
1308 1303 counter += 1
1309 1304
1310 1305 def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
1311 1306 xmin=None, xmax=None, ymin=None, ymax=None,
1312 1307 save=False, figpath='./', figfile=None, show=True):
1313 1308
1314 1309 """
1315 1310
1316 1311 Input:
1317 1312 dataOut :
1318 1313 id :
1319 1314 wintitle :
1320 1315 channelList :
1321 1316 showProfile :
1322 1317 xmin : None,
1323 1318 xmax : None,
1324 1319 zmin : None,
1325 1320 zmax : None
1326 1321 """
1327 1322
1328 1323 if cutHeight==None:
1329 1324 h=270
1330 1325 heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
1331 1326 cutHeight = dataOut.heightList[heightindex]
1332 1327
1333 1328 factor = dataOut.normFactor
1334 1329 x = dataOut.abscissaList[:-1]
1335 1330 #y = dataOut.getHeiRange()
1336 1331
1337 1332 z = dataOut.data_pre[:,:,heightindex]/factor
1338 1333 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1339 1334 avg = numpy.average(z, axis=1)
1340 1335 listChannels = z.shape[0]
1341 1336
1342 1337 #Reconstruct Function
1343 1338 if fit==True:
1344 1339 groupArray = dataOut.groupList
1345 1340 listChannels = groupArray.reshape((groupArray.size))
1346 1341 listChannels.sort()
1347 1342 spcFitLine = numpy.zeros(z.shape)
1348 1343 constants = dataOut.constants
1349 1344
1350 1345 nGroups = groupArray.shape[0]
1351 1346 nChannels = groupArray.shape[1]
1352 1347 nProfiles = z.shape[1]
1353 1348
1354 1349 for f in range(nGroups):
1355 1350 groupChann = groupArray[f,:]
1356 1351 p = dataOut.data_param[f,:,heightindex]
1357 1352 # p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
1358 1353 fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
1359 1354 fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
1360 1355 spcFitLine[groupChann,:] = fitLineAux
1361 1356 # spcFitLine = spcFitLine/factor
1362 1357
1363 1358 z = z[listChannels,:]
1364 1359 spcFitLine = spcFitLine[listChannels,:]
1365 1360 spcFitLinedB = 10*numpy.log10(spcFitLine)
1366 1361
1367 1362 zdB = 10*numpy.log10(z)
1368 1363 #thisDatetime = dataOut.datatime
1369 1364 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1370 1365 title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1371 1366 xlabel = "Velocity (m/s)"
1372 1367 ylabel = "Spectrum"
1373 1368
1374 1369 if not self.isConfig:
1375 1370
1376 1371 nplots = listChannels.size
1377 1372
1378 1373 self.setup(id=id,
1379 1374 nplots=nplots,
1380 1375 wintitle=wintitle,
1381 1376 showprofile=showprofile,
1382 1377 show=show)
1383 1378
1384 1379 if xmin == None: xmin = numpy.nanmin(x)
1385 1380 if xmax == None: xmax = numpy.nanmax(x)
1386 1381 if ymin == None: ymin = numpy.nanmin(zdB)
1387 1382 if ymax == None: ymax = numpy.nanmax(zdB)+2
1388 1383
1389 1384 self.isConfig = True
1390 1385
1391 1386 self.setWinTitle(title)
1392 1387 for i in range(self.nplots):
1393 1388 # title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
1394 1389 title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i])
1395 1390 axes = self.axesList[i*self.__nsubplots]
1396 1391 if fit == False:
1397 1392 axes.pline(x, zdB[i,:],
1398 1393 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1399 1394 xlabel=xlabel, ylabel=ylabel, title=title
1400 1395 )
1401 1396 if fit == True:
1402 1397 fitline=spcFitLinedB[i,:]
1403 1398 y=numpy.vstack([zdB[i,:],fitline] )
1404 1399 legendlabels=['Data','Fitting']
1405 1400 axes.pmultilineyaxis(x, y,
1406 1401 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1407 1402 xlabel=xlabel, ylabel=ylabel, title=title,
1408 1403 legendlabels=legendlabels, marker=None,
1409 1404 linestyle='solid', grid='both')
1410 1405
1411 1406 self.draw()
1412 1407
1413 1408 self.save(figpath=figpath,
1414 1409 figfile=figfile,
1415 1410 save=save,
1416 1411 ftp=ftp,
1417 1412 wr_period=wr_period,
1418 1413 thisDatetime=thisDatetime)
1419 1414
1420 1415
1421 1416 class EWDriftsPlot(Figure):
1422 1417
1423 1418 __isConfig = None
1424 1419 __nsubplots = None
1425 1420
1426 1421 WIDTHPROF = None
1427 1422 HEIGHTPROF = None
1428 1423 PREFIX = 'drift'
1429 1424
1430 1425 def __init__(self, **kwargs):
1431 1426 Figure.__init__(self, **kwargs)
1432 1427 self.timerange = 2*60*60
1433 1428 self.isConfig = False
1434 1429 self.__nsubplots = 1
1435 1430
1436 1431 self.WIDTH = 800
1437 1432 self.HEIGHT = 150
1438 1433 self.WIDTHPROF = 120
1439 1434 self.HEIGHTPROF = 0
1440 1435 self.counter_imagwr = 0
1441 1436
1442 1437 self.PLOT_CODE = EWDRIFT_CODE
1443 1438
1444 1439 self.FTP_WEI = None
1445 1440 self.EXP_CODE = None
1446 1441 self.SUB_EXP_CODE = None
1447 1442 self.PLOT_POS = None
1448 1443 self.tmin = None
1449 1444 self.tmax = None
1450 1445
1451 1446 self.xmin = None
1452 1447 self.xmax = None
1453 1448
1454 1449 self.figfile = None
1455 1450
1456 1451 def getSubplots(self):
1457 1452
1458 1453 ncol = 1
1459 1454 nrow = self.nplots
1460 1455
1461 1456 return nrow, ncol
1462 1457
1463 1458 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1464 1459
1465 1460 self.__showprofile = showprofile
1466 1461 self.nplots = nplots
1467 1462
1468 1463 ncolspan = 1
1469 1464 colspan = 1
1470 1465
1471 1466 self.createFigure(id = id,
1472 1467 wintitle = wintitle,
1473 1468 widthplot = self.WIDTH + self.WIDTHPROF,
1474 1469 heightplot = self.HEIGHT + self.HEIGHTPROF,
1475 1470 show=show)
1476 1471
1477 1472 nrow, ncol = self.getSubplots()
1478 1473
1479 1474 counter = 0
1480 1475 for y in range(nrow):
1481 1476 if counter >= self.nplots:
1482 1477 break
1483 1478
1484 1479 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
1485 1480 counter += 1
1486 1481
1487 1482 def run(self, dataOut, id, wintitle="", channelList=None,
1488 1483 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
1489 1484 zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
1490 1485 timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
1491 1486 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1492 1487 server=None, folder=None, username=None, password=None,
1493 1488 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1494 1489 """
1495 1490
1496 1491 Input:
1497 1492 dataOut :
1498 1493 id :
1499 1494 wintitle :
1500 1495 channelList :
1501 1496 showProfile :
1502 1497 xmin : None,
1503 1498 xmax : None,
1504 1499 ymin : None,
1505 1500 ymax : None,
1506 1501 zmin : None,
1507 1502 zmax : None
1508 1503 """
1509 1504
1510 1505 if timerange is not None:
1511 1506 self.timerange = timerange
1512 1507
1513 1508 tmin = None
1514 1509 tmax = None
1515 1510
1516 1511 x = dataOut.getTimeRange1(dataOut.outputInterval)
1517 1512 # y = dataOut.heightList
1518 1513 y = dataOut.heightList
1519 1514
1520 1515 z = dataOut.data_output
1521 1516 nplots = z.shape[0] #Number of wind dimensions estimated
1522 1517 nplotsw = nplots
1523 1518
1524 1519 #If there is a SNR function defined
1525 1520 if dataOut.data_SNR is not None:
1526 1521 nplots += 1
1527 1522 SNR = dataOut.data_SNR
1528 1523
1529 1524 if SNR_1:
1530 1525 SNR += 1
1531 1526
1532 1527 SNRavg = numpy.average(SNR, axis=0)
1533 1528
1534 1529 SNRdB = 10*numpy.log10(SNR)
1535 1530 SNRavgdB = 10*numpy.log10(SNRavg)
1536 1531
1537 1532 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
1538 1533
1539 1534 for i in range(nplotsw):
1540 1535 z[i,ind] = numpy.nan
1541 1536
1542 1537
1543 1538 showprofile = False
1544 1539 # thisDatetime = dataOut.datatime
1545 1540 thisDatetime = datetime.datetime.utcfromtimestamp(x[1])
1546 1541 title = wintitle + " EW Drifts"
1547 1542 xlabel = ""
1548 1543 ylabel = "Height (Km)"
1549 1544
1550 1545 if not self.isConfig:
1551 1546
1552 1547 self.setup(id=id,
1553 1548 nplots=nplots,
1554 1549 wintitle=wintitle,
1555 1550 showprofile=showprofile,
1556 1551 show=show)
1557 1552
1558 1553 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1559 1554
1560 1555 if ymin == None: ymin = numpy.nanmin(y)
1561 1556 if ymax == None: ymax = numpy.nanmax(y)
1562 1557
1563 1558 if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
1564 1559 if zminZonal == None: zminZonal = -zmaxZonal
1565 1560 if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
1566 1561 if zminVertical == None: zminVertical = -zmaxVertical
1567 1562
1568 1563 if dataOut.data_SNR is not None:
1569 1564 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
1570 1565 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
1571 1566
1572 1567 self.FTP_WEI = ftp_wei
1573 1568 self.EXP_CODE = exp_code
1574 1569 self.SUB_EXP_CODE = sub_exp_code
1575 1570 self.PLOT_POS = plot_pos
1576 1571
1577 1572 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1578 1573 self.isConfig = True
1579 1574
1580 1575
1581 1576 self.setWinTitle(title)
1582 1577
1583 1578 if ((self.xmax - x[1]) < (x[1]-x[0])):
1584 1579 x[1] = self.xmax
1585 1580
1586 1581 strWind = ['Zonal','Vertical']
1587 1582 strCb = 'Velocity (m/s)'
1588 1583 zmaxVector = [zmaxZonal, zmaxVertical]
1589 1584 zminVector = [zminZonal, zminVertical]
1590 1585
1591 1586 for i in range(nplotsw):
1592 1587
1593 1588 title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1594 1589 axes = self.axesList[i*self.__nsubplots]
1595 1590
1596 1591 z1 = z[i,:].reshape((1,-1))
1597 1592
1598 1593 axes.pcolorbuffer(x, y, z1,
1599 1594 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
1600 1595 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1601 1596 ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
1602 1597
1603 1598 if dataOut.data_SNR is not None:
1604 1599 i += 1
1605 1600 if SNR_1:
1606 1601 title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1607 1602 else:
1608 1603 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1609 1604 axes = self.axesList[i*self.__nsubplots]
1610 1605 SNRavgdB = SNRavgdB.reshape((1,-1))
1611 1606
1612 1607 axes.pcolorbuffer(x, y, SNRavgdB,
1613 1608 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1614 1609 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1615 1610 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
1616 1611
1617 1612 self.draw()
1618 1613
1619 1614 if x[1] >= self.axesList[0].xmax:
1620 1615 self.counter_imagwr = wr_period
1621 1616 self.isConfig = False
1622 1617 self.figfile = None
1623 1618
1624 1619
1625 1620
1626 1621
1627 1622 class PhasePlot(Figure):
1628 1623
1629 1624 __isConfig = None
1630 1625 __nsubplots = None
1631 1626
1632 1627 PREFIX = 'mphase'
1633 1628
1634 1629 def __init__(self, **kwargs):
1635 1630 Figure.__init__(self, **kwargs)
1636 1631 self.timerange = 24*60*60
1637 1632 self.isConfig = False
1638 1633 self.__nsubplots = 1
1639 1634 self.counter_imagwr = 0
1640 1635 self.WIDTH = 600
1641 1636 self.HEIGHT = 300
1642 1637 self.WIDTHPROF = 120
1643 1638 self.HEIGHTPROF = 0
1644 1639 self.xdata = None
1645 1640 self.ydata = None
1646 1641
1647 1642 self.PLOT_CODE = MPHASE_CODE
1648 1643
1649 1644 self.FTP_WEI = None
1650 1645 self.EXP_CODE = None
1651 1646 self.SUB_EXP_CODE = None
1652 1647 self.PLOT_POS = None
1653 1648
1654 1649
1655 1650 self.filename_phase = None
1656 1651
1657 1652 self.figfile = None
1658 1653
1659 1654 def getSubplots(self):
1660 1655
1661 1656 ncol = 1
1662 1657 nrow = 1
1663 1658
1664 1659 return nrow, ncol
1665 1660
1666 1661 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1667 1662
1668 1663 self.__showprofile = showprofile
1669 1664 self.nplots = nplots
1670 1665
1671 1666 ncolspan = 7
1672 1667 colspan = 6
1673 1668 self.__nsubplots = 2
1674 1669
1675 1670 self.createFigure(id = id,
1676 1671 wintitle = wintitle,
1677 1672 widthplot = self.WIDTH+self.WIDTHPROF,
1678 1673 heightplot = self.HEIGHT+self.HEIGHTPROF,
1679 1674 show=show)
1680 1675
1681 1676 nrow, ncol = self.getSubplots()
1682 1677
1683 1678 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1684 1679
1685 1680
1686 1681 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1687 1682 xmin=None, xmax=None, ymin=None, ymax=None,
1688 1683 timerange=None,
1689 1684 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1690 1685 server=None, folder=None, username=None, password=None,
1691 1686 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1692 1687
1693 1688
1694 1689 tmin = None
1695 1690 tmax = None
1696 1691 x = dataOut.getTimeRange1(dataOut.outputInterval)
1697 1692 y = dataOut.getHeiRange()
1698 1693
1699 1694
1700 1695 #thisDatetime = dataOut.datatime
1701 1696 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
1702 1697 title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1703 1698 xlabel = "Local Time"
1704 1699 ylabel = "Phase"
1705 1700
1706 1701
1707 1702 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1708 1703 phase_beacon = dataOut.data_output
1709 1704 update_figfile = False
1710 1705
1711 1706 if not self.isConfig:
1712 1707
1713 1708 self.nplots = phase_beacon.size
1714 1709
1715 1710 self.setup(id=id,
1716 1711 nplots=self.nplots,
1717 1712 wintitle=wintitle,
1718 1713 showprofile=showprofile,
1719 1714 show=show)
1720 1715
1721 1716 if timerange is not None:
1722 1717 self.timerange = timerange
1723 1718
1724 1719 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1725 1720
1726 1721 if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0
1727 1722 if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0
1728 1723
1729 1724 self.FTP_WEI = ftp_wei
1730 1725 self.EXP_CODE = exp_code
1731 1726 self.SUB_EXP_CODE = sub_exp_code
1732 1727 self.PLOT_POS = plot_pos
1733 1728
1734 1729 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1735 1730 self.isConfig = True
1736 1731 self.figfile = figfile
1737 1732 self.xdata = numpy.array([])
1738 1733 self.ydata = numpy.array([])
1739 1734
1740 1735 #open file beacon phase
1741 1736 path = '%s%03d' %(self.PREFIX, self.id)
1742 1737 beacon_file = os.path.join(path,'%s.txt'%self.name)
1743 1738 self.filename_phase = os.path.join(figpath,beacon_file)
1744 1739 update_figfile = True
1745 1740
1746 1741
1747 1742 #store data beacon phase
1748 1743 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1749 1744
1750 1745 self.setWinTitle(title)
1751 1746
1752 1747
1753 1748 title = "Phase Offset %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1754 1749
1755 1750 legendlabels = ["phase %d"%(chan) for chan in numpy.arange(self.nplots)]
1756 1751
1757 1752 axes = self.axesList[0]
1758 1753
1759 1754 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1760 1755
1761 1756 if len(self.ydata)==0:
1762 1757 self.ydata = phase_beacon.reshape(-1,1)
1763 1758 else:
1764 1759 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1765 1760
1766 1761
1767 1762 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1768 1763 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1769 1764 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1770 1765 XAxisAsTime=True, grid='both'
1771 1766 )
1772 1767
1773 1768 self.draw()
1774 1769
1775 1770 self.save(figpath=figpath,
1776 1771 figfile=figfile,
1777 1772 save=save,
1778 1773 ftp=ftp,
1779 1774 wr_period=wr_period,
1780 1775 thisDatetime=thisDatetime,
1781 1776 update_figfile=update_figfile)
1782 1777
1783 1778 if dataOut.ltctime + dataOut.outputInterval >= self.xmax:
1784 1779 self.counter_imagwr = wr_period
1785 1780 self.isConfig = False
1786 1781 update_figfile = True
1787 1782
1788 1783
1789 1784
1790 1785 class NSMeteorDetection1Plot(Figure):
1791 1786
1792 1787 isConfig = None
1793 1788 __nsubplots = None
1794 1789
1795 1790 WIDTHPROF = None
1796 1791 HEIGHTPROF = None
1797 1792 PREFIX = 'nsm'
1798 1793
1799 1794 zminList = None
1800 1795 zmaxList = None
1801 1796 cmapList = None
1802 1797 titleList = None
1803 1798 nPairs = None
1804 1799 nChannels = None
1805 1800 nParam = None
1806 1801
1807 1802 def __init__(self, **kwargs):
1808 1803 Figure.__init__(self, **kwargs)
1809 1804 self.isConfig = False
1810 1805 self.__nsubplots = 1
1811 1806
1812 1807 self.WIDTH = 750
1813 1808 self.HEIGHT = 250
1814 1809 self.WIDTHPROF = 120
1815 1810 self.HEIGHTPROF = 0
1816 1811 self.counter_imagwr = 0
1817 1812
1818 1813 self.PLOT_CODE = SPEC_CODE
1819 1814
1820 1815 self.FTP_WEI = None
1821 1816 self.EXP_CODE = None
1822 1817 self.SUB_EXP_CODE = None
1823 1818 self.PLOT_POS = None
1824 1819
1825 1820 self.__xfilter_ena = False
1826 1821 self.__yfilter_ena = False
1827 1822
1828 1823 def getSubplots(self):
1829 1824
1830 1825 ncol = 3
1831 1826 nrow = int(numpy.ceil(self.nplots/3.0))
1832 1827
1833 1828 return nrow, ncol
1834 1829
1835 1830 def setup(self, id, nplots, wintitle, show=True):
1836 1831
1837 1832 self.nplots = nplots
1838 1833
1839 1834 ncolspan = 1
1840 1835 colspan = 1
1841 1836
1842 1837 self.createFigure(id = id,
1843 1838 wintitle = wintitle,
1844 1839 widthplot = self.WIDTH + self.WIDTHPROF,
1845 1840 heightplot = self.HEIGHT + self.HEIGHTPROF,
1846 1841 show=show)
1847 1842
1848 1843 nrow, ncol = self.getSubplots()
1849 1844
1850 1845 counter = 0
1851 1846 for y in range(nrow):
1852 1847 for x in range(ncol):
1853 1848
1854 1849 if counter >= self.nplots:
1855 1850 break
1856 1851
1857 1852 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1858 1853
1859 1854 counter += 1
1860 1855
1861 1856 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
1862 1857 xmin=None, xmax=None, ymin=None, ymax=None, SNRmin=None, SNRmax=None,
1863 1858 vmin=None, vmax=None, wmin=None, wmax=None, mode = 'SA',
1864 1859 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1865 1860 server=None, folder=None, username=None, password=None,
1866 1861 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
1867 1862 xaxis="frequency"):
1868 1863
1869 1864 """
1870 1865
1871 1866 Input:
1872 1867 dataOut :
1873 1868 id :
1874 1869 wintitle :
1875 1870 channelList :
1876 1871 showProfile :
1877 1872 xmin : None,
1878 1873 xmax : None,
1879 1874 ymin : None,
1880 1875 ymax : None,
1881 1876 zmin : None,
1882 1877 zmax : None
1883 1878 """
1884 1879 #SEPARAR EN DOS PLOTS
1885 1880 nParam = dataOut.data_param.shape[1] - 3
1886 1881
1887 1882 utctime = dataOut.data_param[0,0]
1888 1883 tmet = dataOut.data_param[:,1].astype(int)
1889 1884 hmet = dataOut.data_param[:,2].astype(int)
1890 1885
1891 1886 x = dataOut.abscissaList
1892 1887 y = dataOut.heightList
1893 1888
1894 1889 z = numpy.zeros((nParam, y.size, x.size - 1))
1895 1890 z[:,:] = numpy.nan
1896 1891 z[:,hmet,tmet] = dataOut.data_param[:,3:].T
1897 1892 z[0,:,:] = 10*numpy.log10(z[0,:,:])
1898 1893
1899 1894 xlabel = "Time (s)"
1900 1895 ylabel = "Range (km)"
1901 1896
1902 1897 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
1903 1898
1904 1899 if not self.isConfig:
1905 1900
1906 1901 nplots = nParam
1907 1902
1908 1903 self.setup(id=id,
1909 1904 nplots=nplots,
1910 1905 wintitle=wintitle,
1911 1906 show=show)
1912 1907
1913 1908 if xmin is None: xmin = numpy.nanmin(x)
1914 1909 if xmax is None: xmax = numpy.nanmax(x)
1915 1910 if ymin is None: ymin = numpy.nanmin(y)
1916 1911 if ymax is None: ymax = numpy.nanmax(y)
1917 1912 if SNRmin is None: SNRmin = numpy.nanmin(z[0,:])
1918 1913 if SNRmax is None: SNRmax = numpy.nanmax(z[0,:])
1919 1914 if vmax is None: vmax = numpy.nanmax(numpy.abs(z[1,:]))
1920 1915 if vmin is None: vmin = -vmax
1921 1916 if wmin is None: wmin = 0
1922 1917 if wmax is None: wmax = 50
1923 1918
1924 1919 pairsList = dataOut.groupList
1925 1920 self.nPairs = len(dataOut.groupList)
1926 1921
1927 1922 zminList = [SNRmin, vmin, cmin] + [pmin]*self.nPairs
1928 1923 zmaxList = [SNRmax, vmax, cmax] + [pmax]*self.nPairs
1929 1924 titleList = ["SNR","Radial Velocity","Coherence"]
1930 1925 cmapList = ["jet","RdBu_r","jet"]
1931 1926
1932 1927 for i in range(self.nPairs):
1933 1928 strAux1 = "Phase Difference "+ str(pairsList[i][0]) + str(pairsList[i][1])
1934 1929 titleList = titleList + [strAux1]
1935 1930 cmapList = cmapList + ["RdBu_r"]
1936 1931
1937 1932 self.zminList = zminList
1938 1933 self.zmaxList = zmaxList
1939 1934 self.cmapList = cmapList
1940 1935 self.titleList = titleList
1941 1936
1942 1937 self.FTP_WEI = ftp_wei
1943 1938 self.EXP_CODE = exp_code
1944 1939 self.SUB_EXP_CODE = sub_exp_code
1945 1940 self.PLOT_POS = plot_pos
1946 1941
1947 1942 self.isConfig = True
1948 1943
1949 1944 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
1950 1945
1951 1946 for i in range(nParam):
1952 1947 title = self.titleList[i] + ": " +str_datetime
1953 1948 axes = self.axesList[i]
1954 1949 axes.pcolor(x, y, z[i,:].T,
1955 1950 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=self.zminList[i], zmax=self.zmaxList[i],
1956 1951 xlabel=xlabel, ylabel=ylabel, title=title, colormap=self.cmapList[i],ticksize=9, cblabel='')
1957 1952 self.draw()
1958 1953
1959 1954 if figfile == None:
1960 1955 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
1961 1956 name = str_datetime
1962 1957 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
1963 1958 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
1964 1959 figfile = self.getFilename(name)
1965 1960
1966 1961 self.save(figpath=figpath,
1967 1962 figfile=figfile,
1968 1963 save=save,
1969 1964 ftp=ftp,
1970 1965 wr_period=wr_period,
1971 1966 thisDatetime=thisDatetime)
1972 1967
1973 1968
1974 1969 class NSMeteorDetection2Plot(Figure):
1975 1970
1976 1971 isConfig = None
1977 1972 __nsubplots = None
1978 1973
1979 1974 WIDTHPROF = None
1980 1975 HEIGHTPROF = None
1981 1976 PREFIX = 'nsm'
1982 1977
1983 1978 zminList = None
1984 1979 zmaxList = None
1985 1980 cmapList = None
1986 1981 titleList = None
1987 1982 nPairs = None
1988 1983 nChannels = None
1989 1984 nParam = None
1990 1985
1991 1986 def __init__(self, **kwargs):
1992 1987 Figure.__init__(self, **kwargs)
1993 1988 self.isConfig = False
1994 1989 self.__nsubplots = 1
1995 1990
1996 1991 self.WIDTH = 750
1997 1992 self.HEIGHT = 250
1998 1993 self.WIDTHPROF = 120
1999 1994 self.HEIGHTPROF = 0
2000 1995 self.counter_imagwr = 0
2001 1996
2002 1997 self.PLOT_CODE = SPEC_CODE
2003 1998
2004 1999 self.FTP_WEI = None
2005 2000 self.EXP_CODE = None
2006 2001 self.SUB_EXP_CODE = None
2007 2002 self.PLOT_POS = None
2008 2003
2009 2004 self.__xfilter_ena = False
2010 2005 self.__yfilter_ena = False
2011 2006
2012 2007 def getSubplots(self):
2013 2008
2014 2009 ncol = 3
2015 2010 nrow = int(numpy.ceil(self.nplots/3.0))
2016 2011
2017 2012 return nrow, ncol
2018 2013
2019 2014 def setup(self, id, nplots, wintitle, show=True):
2020 2015
2021 2016 self.nplots = nplots
2022 2017
2023 2018 ncolspan = 1
2024 2019 colspan = 1
2025 2020
2026 2021 self.createFigure(id = id,
2027 2022 wintitle = wintitle,
2028 2023 widthplot = self.WIDTH + self.WIDTHPROF,
2029 2024 heightplot = self.HEIGHT + self.HEIGHTPROF,
2030 2025 show=show)
2031 2026
2032 2027 nrow, ncol = self.getSubplots()
2033 2028
2034 2029 counter = 0
2035 2030 for y in range(nrow):
2036 2031 for x in range(ncol):
2037 2032
2038 2033 if counter >= self.nplots:
2039 2034 break
2040 2035
2041 2036 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
2042 2037
2043 2038 counter += 1
2044 2039
2045 2040 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
2046 2041 xmin=None, xmax=None, ymin=None, ymax=None, SNRmin=None, SNRmax=None,
2047 2042 vmin=None, vmax=None, wmin=None, wmax=None, mode = 'SA',
2048 2043 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
2049 2044 server=None, folder=None, username=None, password=None,
2050 2045 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
2051 2046 xaxis="frequency"):
2052 2047
2053 2048 """
2054 2049
2055 2050 Input:
2056 2051 dataOut :
2057 2052 id :
2058 2053 wintitle :
2059 2054 channelList :
2060 2055 showProfile :
2061 2056 xmin : None,
2062 2057 xmax : None,
2063 2058 ymin : None,
2064 2059 ymax : None,
2065 2060 zmin : None,
2066 2061 zmax : None
2067 2062 """
2068 2063 #Rebuild matrix
2069 2064 utctime = dataOut.data_param[0,0]
2070 2065 cmet = dataOut.data_param[:,1].astype(int)
2071 2066 tmet = dataOut.data_param[:,2].astype(int)
2072 2067 hmet = dataOut.data_param[:,3].astype(int)
2073 2068
2074 2069 nParam = 3
2075 2070 nChan = len(dataOut.groupList)
2076 2071 x = dataOut.abscissaList
2077 2072 y = dataOut.heightList
2078 2073
2079 2074 z = numpy.full((nChan, nParam, y.size, x.size - 1),numpy.nan)
2080 2075 z[cmet,:,hmet,tmet] = dataOut.data_param[:,4:]
2081 2076 z[:,0,:,:] = 10*numpy.log10(z[:,0,:,:]) #logarithmic scale
2082 2077 z = numpy.reshape(z, (nChan*nParam, y.size, x.size-1))
2083 2078
2084 2079 xlabel = "Time (s)"
2085 2080 ylabel = "Range (km)"
2086 2081
2087 2082 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
2088 2083
2089 2084 if not self.isConfig:
2090 2085
2091 2086 nplots = nParam*nChan
2092 2087
2093 2088 self.setup(id=id,
2094 2089 nplots=nplots,
2095 2090 wintitle=wintitle,
2096 2091 show=show)
2097 2092
2098 2093 if xmin is None: xmin = numpy.nanmin(x)
2099 2094 if xmax is None: xmax = numpy.nanmax(x)
2100 2095 if ymin is None: ymin = numpy.nanmin(y)
2101 2096 if ymax is None: ymax = numpy.nanmax(y)
2102 2097 if SNRmin is None: SNRmin = numpy.nanmin(z[0,:])
2103 2098 if SNRmax is None: SNRmax = numpy.nanmax(z[0,:])
2104 2099 if vmax is None: vmax = numpy.nanmax(numpy.abs(z[1,:]))
2105 2100 if vmin is None: vmin = -vmax
2106 2101 if wmin is None: wmin = 0
2107 2102 if wmax is None: wmax = 50
2108 2103
2109 2104 self.nChannels = nChan
2110 2105
2111 2106 zminList = []
2112 2107 zmaxList = []
2113 2108 titleList = []
2114 2109 cmapList = []
2115 2110 for i in range(self.nChannels):
2116 2111 strAux1 = "SNR Channel "+ str(i)
2117 2112 strAux2 = "Radial Velocity Channel "+ str(i)
2118 2113 strAux3 = "Spectral Width Channel "+ str(i)
2119 2114
2120 2115 titleList = titleList + [strAux1,strAux2,strAux3]
2121 2116 cmapList = cmapList + ["jet","RdBu_r","jet"]
2122 2117 zminList = zminList + [SNRmin,vmin,wmin]
2123 2118 zmaxList = zmaxList + [SNRmax,vmax,wmax]
2124 2119
2125 2120 self.zminList = zminList
2126 2121 self.zmaxList = zmaxList
2127 2122 self.cmapList = cmapList
2128 2123 self.titleList = titleList
2129 2124
2130 2125 self.FTP_WEI = ftp_wei
2131 2126 self.EXP_CODE = exp_code
2132 2127 self.SUB_EXP_CODE = sub_exp_code
2133 2128 self.PLOT_POS = plot_pos
2134 2129
2135 2130 self.isConfig = True
2136 2131
2137 2132 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
2138 2133
2139 2134 for i in range(self.nplots):
2140 2135 title = self.titleList[i] + ": " +str_datetime
2141 2136 axes = self.axesList[i]
2142 2137 axes.pcolor(x, y, z[i,:].T,
2143 2138 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=self.zminList[i], zmax=self.zmaxList[i],
2144 2139 xlabel=xlabel, ylabel=ylabel, title=title, colormap=self.cmapList[i],ticksize=9, cblabel='')
2145 2140 self.draw()
2146 2141
2147 2142 if figfile == None:
2148 2143 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
2149 2144 name = str_datetime
2150 2145 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
2151 2146 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
2152 2147 figfile = self.getFilename(name)
2153 2148
2154 2149 self.save(figpath=figpath,
2155 2150 figfile=figfile,
2156 2151 save=save,
2157 2152 ftp=ftp,
2158 2153 wr_period=wr_period,
2159 2154 thisDatetime=thisDatetime)
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file, binary diff hidden
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file, binary diff hidden
Show file before | Show file after | Hide comments
@@ -1,364 +1,364
1 1 '''
2 2 Created on Nov 9, 2016
3 3
4 4 @author: roj- LouVD
5 5 '''
6 6
7 7
8 8 import os
9 9 import sys
10 10 import time
11 11 import glob
12 12 import datetime
13
13 14 import numpy
14 15
15 16 from schainpy.model.proc.jroproc_base import ProcessingUnit
16 17 from schainpy.model.data.jrodata import Parameters
17 18 from schainpy.model.io.jroIO_base import JRODataReader, isNumber
18 19
19 20 FILE_HEADER_STRUCTURE = numpy.dtype([
20 21 ('FMN', '<u4'),
21 22 ('nrec', '<u4'),
22 23 ('fr_offset', '<u4'),
23 24 ('id', '<u4'),
24 25 ('site', 'u1', (32,))
25 26 ])
26 27
27 28 REC_HEADER_STRUCTURE = numpy.dtype([
28 29 ('rmn', '<u4'),
29 30 ('rcounter', '<u4'),
30 31 ('nr_offset', '<u4'),
31 32 ('tr_offset', '<u4'),
32 33 ('time', '<u4'),
33 34 ('time_msec', '<u4'),
34 35 ('tag', 'u1', (32,)),
35 36 ('comments', 'u1', (32,)),
36 37 ('lat', '<f4'),
37 38 ('lon', '<f4'),
38 39 ('gps_status', '<u4'),
39 40 ('freq', '<u4'),
40 41 ('freq0', '<u4'),
41 42 ('nchan', '<u4'),
42 43 ('delta_r', '<u4'),
43 44 ('nranges', '<u4'),
44 45 ('r0', '<u4'),
45 46 ('prf', '<u4'),
46 47 ('ncoh', '<u4'),
47 48 ('npoints', '<u4'),
48 49 ('polarization', '<i4'),
49 50 ('rx_filter', '<u4'),
50 51 ('nmodes', '<u4'),
51 52 ('dmode_index', '<u4'),
52 53 ('dmode_rngcorr', '<u4'),
53 54 ('nrxs', '<u4'),
54 55 ('acf_length', '<u4'),
55 56 ('acf_lags', '<u4'),
56 57 ('sea_to_atmos', '<f4'),
57 58 ('sea_notch', '<u4'),
58 59 ('lh_sea', '<u4'),
59 60 ('hh_sea', '<u4'),
60 61 ('nbins_sea', '<u4'),
61 62 ('min_snr', '<f4'),
62 63 ('min_cc', '<f4'),
63 64 ('max_time_diff', '<f4')
64 65 ])
65 66
66 67 DATA_STRUCTURE = numpy.dtype([
67 68 ('range', '<u4'),
68 69 ('status', '<u4'),
69 70 ('zonal', '<f4'),
70 71 ('meridional', '<f4'),
71 72 ('vertical', '<f4'),
72 73 ('zonal_a', '<f4'),
73 74 ('meridional_a', '<f4'),
74 75 ('corrected_fading', '<f4'), # seconds
75 76 ('uncorrected_fading', '<f4'), # seconds
76 77 ('time_diff', '<f4'),
77 78 ('major_axis', '<f4'),
78 79 ('axial_ratio', '<f4'),
79 80 ('orientation', '<f4'),
80 81 ('sea_power', '<u4'),
81 82 ('sea_algorithm', '<u4')
82 83 ])
83 84
84 85 class BLTRParamReader(JRODataReader, ProcessingUnit):
85 86 '''
86 87 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR from *.sswma files
87 88 '''
88 89
89 90 ext = '.sswma'
90 91
91 92 def __init__(self, **kwargs):
92 93
93 94 ProcessingUnit.__init__(self , **kwargs)
94 95
95 96 self.dataOut = Parameters()
96 97 self.counter_records = 0
97 98 self.flagNoMoreFiles = 0
98 99 self.isConfig = False
99 100 self.filename = None
100 101
101 102 def setup(self,
102 103 path=None,
103 104 startDate=None,
104 105 endDate=None,
105 106 ext=None,
106 107 startTime=datetime.time(0, 0, 0),
107 108 endTime=datetime.time(23, 59, 59),
108 109 timezone=0,
109 110 status_value=0,
110 111 **kwargs):
111 112
112 113 self.path = path
113 114 self.startTime = startTime
114 115 self.endTime = endTime
115 116 self.status_value = status_value
116 117
117 118 if self.path is None:
118 119 raise ValueError, "The path is not valid"
119 120
120 121 if ext is None:
121 122 ext = self.ext
122 123
123 124 self.search_files(self.path, startDate, endDate, ext)
124 125 self.timezone = timezone
125 126 self.fileIndex = 0
126 127
127 128 if not self.fileList:
128 129 raise Warning, "There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' "%(path)
129 130
130 131 self.setNextFile()
131 132
132 133 def search_files(self, path, startDate, endDate, ext):
133 134 '''
134 135 Searching for BLTR rawdata file in path
135 136 Creating a list of file to proces included in [startDate,endDate]
136 137
137 138 Input:
138 139 path - Path to find BLTR rawdata files
139 140 startDate - Select file from this date
140 141 enDate - Select file until this date
141 142 ext - Extension of the file to read
142 143
143 144 '''
144 145
145 146 print 'Searching file in %s ' % (path)
146 147 foldercounter = 0
147 148 fileList0 = glob.glob1(path, "*%s" % ext)
148 149 fileList0.sort()
149 150
150 151 self.fileList = []
151 152 self.dateFileList = []
152 153
153 154 for thisFile in fileList0:
154 155 year = thisFile[-14:-10]
155 156 if not isNumber(year):
156 157 continue
157 158
158 159 month = thisFile[-10:-8]
159 160 if not isNumber(month):
160 161 continue
161 162
162 163 day = thisFile[-8:-6]
163 164 if not isNumber(day):
164 165 continue
165 166
166 167 year, month, day = int(year), int(month), int(day)
167 168 dateFile = datetime.date(year, month, day)
168 169
169 170 if (startDate > dateFile) or (endDate < dateFile):
170 171 continue
171 172
172 173 self.fileList.append(thisFile)
173 174 self.dateFileList.append(dateFile)
174 175
175 176 return
176 177
177 178 def setNextFile(self):
178 179
179 180 file_id = self.fileIndex
180 181
181 182 if file_id == len(self.fileList):
182 183 print '\nNo more files in the folder'
183 184 print 'Total number of file(s) read : {}'.format(self.fileIndex + 1)
184 185 self.flagNoMoreFiles = 1
185 186 return 0
186 187
187 188 print '\n[Setting file] (%s) ...' % self.fileList[file_id]
188 189 filename = os.path.join(self.path, self.fileList[file_id])
189 190
190 191 dirname, name = os.path.split(filename)
191 192 self.siteFile = name.split('.')[0] # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya
192 193 if self.filename is not None:
193 194 self.fp.close()
194 195 self.filename = filename
195 196 self.fp = open(self.filename, 'rb')
196 197 self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
197 198 self.nrecords = self.header_file['nrec'][0]
198 199 self.sizeOfFile = os.path.getsize(self.filename)
199 200 self.counter_records = 0
200 201 self.flagIsNewFile = 0
201 202 self.fileIndex += 1
202 203
203 204 return 1
204 205
205 206 def readNextBlock(self):
206 207
207 208 while True:
208 209 if self.counter_records == self.nrecords:
209 210 self.flagIsNewFile = 1
210 211 if not self.setNextFile():
211 212 return 0
212 213
213 214 self.readBlock()
214 215
215 216 if (self.datatime.time() < self.startTime) or (self.datatime.time() > self.endTime):
216 217 print "[Reading] Record No. %d/%d -> %s [Skipping]" %(
217 218 self.counter_records,
218 219 self.nrecords,
219 220 self.datatime.ctime())
220 221 continue
221 222 break
222 223
223 224 print "[Reading] Record No. %d/%d -> %s" %(
224 225 self.counter_records,
225 226 self.nrecords,
226 227 self.datatime.ctime())
227 228
228 229 return 1
229 230
230 231 def readBlock(self):
231 232
232 233 pointer = self.fp.tell()
233 234 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
234 235 self.nchannels = header_rec['nchan'][0]/2
235 236 self.kchan = header_rec['nrxs'][0]
236 237 self.nmodes = header_rec['nmodes'][0]
237 238 self.nranges = header_rec['nranges'][0]
238 239 self.fp.seek(pointer)
239 240 self.height = numpy.empty((self.nmodes, self.nranges))
240 241 self.snr = numpy.empty((self.nmodes, self.nchannels, self.nranges))
241 242 self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
242 243
243 244 for mode in range(self.nmodes):
244 self.readHeader()
245 self.readHeader()
245 246 data = self.readData()
246 247 self.height[mode] = (data[0] - self.correction) / 1000.
247 248 self.buffer[mode] = data[1]
248 249 self.snr[mode] = data[2]
249 250
250 251 self.counter_records = self.counter_records + self.nmodes
251 252
252 253 return
253 254
254 255 def readHeader(self):
255 256 '''
256 257 RecordHeader of BLTR rawdata file
257 258 '''
258 259
259 260 header_structure = numpy.dtype(
260 261 REC_HEADER_STRUCTURE.descr + [
261 262 ('antenna_coord', 'f4', (2, self.nchannels)),
262 263 ('rx_gains', 'u4', (self.nchannels,)),
263 264 ('rx_analysis', 'u4', (self.nchannels,))
264 265 ]
265 266 )
266 267
267 268 self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
268 269 self.lat = self.header_rec['lat'][0]
269 270 self.lon = self.header_rec['lon'][0]
270 271 self.delta = self.header_rec['delta_r'][0]
271 272 self.correction = self.header_rec['dmode_rngcorr'][0]
272 273 self.imode = self.header_rec['dmode_index'][0]
273 274 self.antenna = self.header_rec['antenna_coord']
274 275 self.rx_gains = self.header_rec['rx_gains']
275 self.time1 = self.header_rec['time'][0]
276 self.time = self.header_rec['time'][0]
276 277 tseconds = self.header_rec['time'][0]
277 278 local_t1 = time.localtime(tseconds)
278 279 self.year = local_t1.tm_year
279 280 self.month = local_t1.tm_mon
280 281 self.day = local_t1.tm_mday
281 282 self.t = datetime.datetime(self.year, self.month, self.day)
282 self.datatime = datetime.datetime.utcfromtimestamp(self.time1)
283 self.datatime = datetime.datetime.utcfromtimestamp(self.time)
283 284
284 285 def readData(self):
285 286 '''
286 287 Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
287 288
288 289 Input:
289 290 status_value - Array data is set to NAN for values that are not equal to status_value
290 291
291 292 '''
292 293
293 294 data_structure = numpy.dtype(
294 295 DATA_STRUCTURE.descr + [
295 296 ('rx_saturation', 'u4', (self.nchannels,)),
296 297 ('chan_offset', 'u4', (2 * self.nchannels,)),
297 298 ('rx_amp', 'u4', (self.nchannels,)),
298 299 ('rx_snr', 'f4', (self.nchannels,)),
299 300 ('cross_snr', 'f4', (self.kchan,)),
300 301 ('sea_power_relative', 'f4', (self.kchan,))]
301 302 )
302 303
303 304 data = numpy.fromfile(self.fp, data_structure, self.nranges)
304 305
305 306 height = data['range']
306 307 winds = numpy.array((data['zonal'], data['meridional'], data['vertical']))
307 308 snr = data['rx_snr'].T
308 309
309 winds[numpy.where(winds == -9999.)] = numpy.nan
310 winds[numpy.where(winds == -9999.)] = numpy.nan
310 311 winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
311 312 snr[numpy.where(snr == -9999.)] = numpy.nan
312 313 snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
313 314 snr = numpy.power(10, snr / 10)
314 315
315 316 return height, winds, snr
316 317
317 318 def set_output(self):
318 319 '''
319 320 Storing data from databuffer to dataOut object
320 321 '''
321
322 self.dataOut.time1 = self.time1
322
323 323 self.dataOut.data_SNR = self.snr
324 self.dataOut.height= self.height
324 self.dataOut.height = self.height
325 325 self.dataOut.data_output = self.buffer
326 self.dataOut.utctimeInit = self.time1
326 self.dataOut.utctimeInit = self.time
327 327 self.dataOut.utctime = self.dataOut.utctimeInit
328 328 self.dataOut.counter_records = self.counter_records
329 329 self.dataOut.nrecords = self.nrecords
330 330 self.dataOut.useLocalTime = False
331 331 self.dataOut.paramInterval = 157
332 332 self.dataOut.timezone = self.timezone
333 333 self.dataOut.site = self.siteFile
334 334 self.dataOut.nrecords = self.nrecords
335 335 self.dataOut.sizeOfFile = self.sizeOfFile
336 336 self.dataOut.lat = self.lat
337 self.dataOut.lon = self.lon
337 self.dataOut.lon = self.lon
338 338 self.dataOut.channelList = range(self.nchannels)
339 339 self.dataOut.kchan = self.kchan
340 340 # self.dataOut.nHeights = self.nranges
341 341 self.dataOut.delta = self.delta
342 342 self.dataOut.correction = self.correction
343 343 self.dataOut.nmodes = self.nmodes
344 344 self.dataOut.imode = self.imode
345 345 self.dataOut.antenna = self.antenna
346 346 self.dataOut.rx_gains = self.rx_gains
347 347 self.dataOut.flagNoData = False
348 348
349 349 def getData(self):
350 350 '''
351 351 Storing data from databuffer to dataOut object
352 352 '''
353 353 if self.flagNoMoreFiles:
354 354 self.dataOut.flagNoData = True
355 355 print 'No file left to process'
356 356 return 0
357 357
358 if not(self.readNextBlock()):
358 if not self.readNextBlock():
359 359 self.dataOut.flagNoData = True
360 360 return 0
361 361
362 362 self.set_output()
363 363
364 364 return 1
Show file before | Show file after | Hide comments
@@ -1,16 +1,15
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: Processor.py 1 2012-11-12 18:56:07Z murco $
5 5 '''
6 6
7 7 from jroproc_voltage import *
8 8 from jroproc_spectra import *
9 9 from jroproc_heispectra import *
10 10 from jroproc_amisr import *
11 11 from jroproc_correlation import *
12 12 from jroproc_parameters import *
13 13 from jroproc_spectra_lags import *
14 14 from jroproc_spectra_acf import *
15 from jroproc_bltr import *
16
15 from bltrproc_parameters import *
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file, binary diff hidden
Show file before | Show file after | Hide comments
@@ -1,564 +1,393
1 1 '''
2 2 Created on Oct 24, 2016
3 3
4 4 @author: roj- LouVD
5 5 '''
6 6
7 7 import numpy
8 8 import copy
9 9 import datetime
10 10 import time
11 11 from time import gmtime
12 12
13 from jroproc_base import ProcessingUnit
14 from schainpy.model.data.jrodata import Parameters
15 13 from numpy import transpose
16 14
17 from matplotlib import cm
18 import matplotlib.pyplot as plt
19 from matplotlib.mlab import griddata
15 from jroproc_base import ProcessingUnit, Operation
16 from schainpy.model.data.jrodata import Parameters
20 17
21 18
22 19 class BLTRParametersProc(ProcessingUnit):
23 20 '''
24 21 Processing unit for BLTR parameters data (winds)
25 22
26 23 Inputs:
27 24 self.dataOut.nmodes - Number of operation modes
28 25 self.dataOut.nchannels - Number of channels
29 26 self.dataOut.nranges - Number of ranges
30 27
31 28 self.dataOut.data_SNR - SNR array
32 29 self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
33 30 self.dataOut.height - Height array (km)
34 31 self.dataOut.time - Time array (seconds)
35 32
36 33 self.dataOut.fileIndex -Index of the file currently read
37 34 self.dataOut.lat - Latitude coordinate of BLTR location
38 35
39 36 self.dataOut.doy - Experiment doy (number of the day in the current year)
40 37 self.dataOut.month - Experiment month
41 38 self.dataOut.day - Experiment day
42 39 self.dataOut.year - Experiment year
43 40 '''
44 41
45 42 def __init__(self, **kwargs):
46 43 '''
47 44 Inputs: None
48 45 '''
49 46 ProcessingUnit.__init__(self, **kwargs)
50 47 self.dataOut = Parameters()
51 48
52 def run (self, mode):
49 def run(self, mode, snr_threshold=None):
53 50 '''
51
52 Inputs:
53 mode = High resolution (0) or Low resolution (1) data
54 snr_threshold = snr filter value
54 55 '''
55 if self.dataIn.type == "Parameters":
56 if self.dataIn.type == 'Parameters':
56 57 self.dataOut.copy(self.dataIn)
57
58
58 59 self.dataOut.data_output = self.dataOut.data_output[mode]
59 self.dataOut.heightList = self.dataOut.height[mode]
60 self.dataOut.heightList = self.dataOut.height[0]
61 self.dataOut.data_SNR = self.dataOut.data_SNR[mode]
60 62
61 def TimeSelect(self):
62 '''
63 Selecting the time array according to the day of the experiment with a duration of 24 hours
64 '''
65
66 k1 = datetime.datetime(self.dataOut.year, self.dataOut.month, self.dataOut.day) - datetime.timedelta(hours=5)
67 k2 = datetime.datetime(self.dataOut.year, self.dataOut.month, self.dataOut.day) + datetime.timedelta(hours=25) - datetime.timedelta(hours=5)
68 limit_sec1 = time.mktime(k1.timetuple())
69 limit_sec2 = time.mktime(k2.timetuple())
70 valid_data = 0
71
72 doy = self.dataOut.doy
73 t1 = numpy.where(self.dataOut.time[0, :] >= limit_sec1)
74 t2 = numpy.where(self.dataOut.time[0, :] < limit_sec2)
75 time_select = []
76 for val_sec in t1[0]:
77 if val_sec in t2[0]:
78 time_select.append(val_sec)
79
80 time_select = numpy.array(time_select, dtype='int')
81 valid_data = valid_data + len(time_select)
82
83
84 if len(time_select) > 0:
85 self.f_timesec = self.dataOut.time[:, time_select]
86 snr = self.dataOut.data_SNR[time_select, :, :, :]
87 zon = self.dataOut.data_output[0][time_select, :, :]
88 mer = self.dataOut.data_output[1][time_select, :, :]
89 ver = self.dataOut.data_output[2][time_select, :, :]
90
91 if valid_data > 0:
92 self.timesec1 = self.f_timesec[0, :]
93 self.f_height = self.dataOut.height
94 self.f_zon = zon
95 self.f_mer = mer
96 self.f_ver = ver
97 self.f_snr = snr
98 self.f_timedate = []
99 self.f_time = []
100
101 for valuet in self.timesec1:
102 time_t = time.gmtime(valuet)
103 year = time_t.tm_year
104 month = time_t.tm_mon
105 day = time_t.tm_mday
106 hour = time_t.tm_hour
107 minute = time_t.tm_min
108 second = time_t.tm_sec
109 f_timedate_0 = datetime.datetime(year, month, day, hour, minute, second)
110 self.f_timedate.append(f_timedate_0)
111
112 return self.f_timedate, self.f_timesec, self.f_height, self.f_zon, self.f_mer, self.f_ver, self.f_snr
113
114 else:
115 self.f_timesec = None
116 self.f_timedate = None
117 self.f_height = None
118 self.f_zon = None
119 self.f_mer = None
120 self.f_ver = None
121 self.f_snr = None
122 print 'Invalid time'
123
124 return self.f_timedate, self.f_height, self.f_zon, self.f_mer, self.f_ver, self.f_snr
125
126 def SnrFilter(self, snr_val,modetofilter):
63 if snr_threshold is not None:
64 SNRavg = numpy.average(self.dataOut.data_SNR, axis=0)
65 SNRavgdB = 10*numpy.log10(SNRavg)
66 for i in range(3):
67 self.dataOut.data_output[i][SNRavgdB <= snr_threshold] = numpy.nan
68
69 # TODO
70 class OutliersFilter(Operation):
71
72 def __init__(self, **kwargs):
127 73 '''
128 Inputs: snr_val - Threshold value
129
130 74 '''
131 if modetofilter!=2 and modetofilter!=1 :
132 raise ValueError,'Mode to filter should be "1" or "2". {} is not valid, check "Modetofilter" value.'.format(modetofilter)
133 m = modetofilter-1
134
135 print ' SNR filter [mode {}]: SNR <= {}: data_output = NA'.format(modetofilter,snr_val)
136 for k in range(self.dataOut.nchannels):
137 for r in range(self.dataOut.nranges):
138 if self.dataOut.data_SNR[r,k,m] <= snr_val:
139 self.dataOut.data_output[2][r,m] = numpy.nan
140 self.dataOut.data_output[1][r,m] = numpy.nan
141 self.dataOut.data_output[0][r,m] = numpy.nan
142
75 Operation.__init__(self, **kwargs)
143 76
144
145 def OutliersFilter(self,modetofilter,svalue,svalue2,method,factor,filter,npoints):
77 def run(self, svalue2, method, factor, filter, npoints=9):
146 78 '''
147 79 Inputs:
148 80 svalue - string to select array velocity
149 81 svalue2 - string to choose axis filtering
150 82 method - 0 for SMOOTH or 1 for MEDIAN
151 factor - number used to set threshold
83 factor - number used to set threshold
152 84 filter - 1 for data filtering using the standard deviation criteria else 0
153 85 npoints - number of points for mask filter
154
155 '''
156 if modetofilter!=2 and modetofilter!=1 :
157 raise ValueError,'Mode to filter should be "1" or "2". {} is not valid, check "Modetofilter" value.'.format(modetofilter)
158
159 m = modetofilter-1
160
161 print ' Outliers Filter [mode {}]: {} {} / threshold = {}'.format(modetofilter,svalue,svalue,factor)
162
163 npoints = 9
164 novalid = 0.1
165 if svalue == 'zonal':
166 value = self.dataOut.data_output[0]
167
168 elif svalue == 'meridional':
169 value = self.dataOut.data_output[1]
170
171 elif svalue == 'vertical':
172 value = self.dataOut.data_output[2]
173
174 else:
175 print 'value is not defined'
176 return
86 '''
87
88 print ' Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor)
89
177 90
178 if svalue2 == 'inTime':
179 yaxis = self.dataOut.height
180 xaxis = numpy.array([[self.dataOut.time1],[self.dataOut.time1]])
181
182 elif svalue2 == 'inHeight':
183 yaxis = numpy.array([[self.dataOut.time1],[self.dataOut.time1]])
184 xaxis = self.dataOut.height
185
186 else:
187 print 'svalue2 is required, either inHeight or inTime'
188 return
91 yaxis = self.dataOut.heightList
92 xaxis = numpy.array([[self.dataOut.utctime]])
189 93
190 output_array = value
94 # Zonal
95 value_temp = self.dataOut.data_output[0]
191 96
192 value_temp = value[:,m]
193 error = numpy.zeros(len(self.dataOut.time[m,:]))
194 if svalue2 == 'inHeight':
195 value_temp = numpy.transpose(value_temp)
196 error = numpy.zeros(len(self.dataOut.height))
97 # Zonal
98 value_temp = self.dataOut.data_output[1]
197 99
198 htemp = yaxis[m,:]
100 # Vertical
101 value_temp = numpy.transpose(self.dataOut.data_output[2])
102
103 htemp = yaxis
199 104 std = value_temp
200 105 for h in range(len(htemp)):
201 if filter: #standard deviation filtering
202 std[h] = numpy.std(value_temp[h],ddof = npoints)
203 value_temp[numpy.where(std[h] > 5),h] = numpy.nan
204 error[numpy.where(std[h] > 5)] = error[numpy.where(std[h] > 5)] + 1
205
206
207 106 nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
208 minvalid = novalid*len(xaxis[m,:])
209 if minvalid <= npoints:
210 minvalid = npoints
107 minvalid = npoints
211 108
212 109 #only if valid values greater than the minimum required (10%)
213 110 if nvalues_valid > minvalid:
214 111
215 112 if method == 0:
216 113 #SMOOTH
217 114 w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
218 115
219 116
220 117 if method == 1:
221 118 #MEDIAN
222 119 w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
223 120
224 121 dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
225 122
226 123 threshold = dw*factor
227 124 value_temp[numpy.where(w > threshold),h] = numpy.nan
228 125 value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
229 126
230 127
231 128 #At the end
232 129 if svalue2 == 'inHeight':
233 130 value_temp = numpy.transpose(value_temp)
234 131 output_array[:,m] = value_temp
235 132
236 133 if svalue == 'zonal':
237 134 self.dataOut.data_output[0] = output_array
238 135
239 136 elif svalue == 'meridional':
240 137 self.dataOut.data_output[1] = output_array
241 138
242 139 elif svalue == 'vertical':
243 140 self.dataOut.data_output[2] = output_array
244 141
245 142 return self.dataOut.data_output
246 143
247 144
248 145 def Median(self,input,width):
249 146 '''
250 147 Inputs:
251 148 input - Velocity array
252 149 width - Number of points for mask filter
253 150
254 151 '''
255 152
256 153 if numpy.mod(width,2) == 1:
257 154 pc = int((width - 1) / 2)
258 155 cont = 0
259 156 output = []
260 157
261 158 for i in range(len(input)):
262 159 if i >= pc and i < len(input) - pc:
263 160 new2 = input[i-pc:i+pc+1]
264 161 temp = numpy.where(numpy.isfinite(new2))
265 162 new = new2[temp]
266 163 value = numpy.median(new)
267 164 output.append(value)
268 165
269 166 output = numpy.array(output)
270 167 output = numpy.hstack((input[0:pc],output))
271 168 output = numpy.hstack((output,input[-pc:len(input)]))
272 169
273 170 return output
274 171
275 172 def Smooth(self,input,width,edge_truncate = None):
276 173 '''
277 174 Inputs:
278 175 input - Velocity array
279 176 width - Number of points for mask filter
280 177 edge_truncate - 1 for truncate the convolution product else
281 178
282 179 '''
283 180
284 181 if numpy.mod(width,2) == 0:
285 182 real_width = width + 1
286 183 nzeros = width / 2
287 184 else:
288 185 real_width = width
289 186 nzeros = (width - 1) / 2
290 187
291 188 half_width = int(real_width)/2
292 189 length = len(input)
293 190
294 191 gate = numpy.ones(real_width,dtype='float')
295 192 norm_of_gate = numpy.sum(gate)
296 193
297 194 nan_process = 0
298 195 nan_id = numpy.where(numpy.isnan(input))
299 196 if len(nan_id[0]) > 0:
300 197 nan_process = 1
301 198 pb = numpy.zeros(len(input))
302 199 pb[nan_id] = 1.
303 200 input[nan_id] = 0.
304 201
305 202 if edge_truncate == True:
306 203 output = numpy.convolve(input/norm_of_gate,gate,mode='same')
307 204 elif edge_truncate == False or edge_truncate == None:
308 205 output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
309 206 output = numpy.hstack((input[0:half_width],output))
310 207 output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
311 208
312 209 if nan_process:
313 210 pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
314 211 pb = numpy.hstack((numpy.zeros(half_width),pb))
315 212 pb = numpy.hstack((pb,numpy.zeros(half_width)))
316 213 output[numpy.where(pb > 0.9999)] = numpy.nan
317 214 input[nan_id] = numpy.nan
318 215 return output
319 216
320 217 def Average(self,aver=0,nhaver=1):
321 218 '''
322 219 Inputs:
323 220 aver - Indicates the time period over which is averaged or consensus data
324 221 nhaver - Indicates the decimation factor in heights
325 222
326 223 '''
327 224 nhpoints = 48
328 225
329 226 lat_piura = -5.17
330 227 lat_huancayo = -12.04
331 228 lat_porcuya = -5.8
332 229
333 230 if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
334 231 hcm = 3.
335 232 if self.dataOut.year == 2003 :
336 233 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
337 234 nhpoints = 12
338 235
339 236 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
340 237 hcm = 3.
341 238 if self.dataOut.year == 2003 :
342 239 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
343 240 nhpoints = 12
344 241
345 242
346 243 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
347 244 hcm = 5.#2
348 245
349 246 pdata = 0.2
350 247 taver = [1,2,3,4,6,8,12,24]
351 248 t0 = 0
352 249 tf = 24
353 250 ntime =(tf-t0)/taver[aver]
354 251 ti = numpy.arange(ntime)
355 252 tf = numpy.arange(ntime) + taver[aver]
356 253
357 254
358 255 old_height = self.dataOut.heightList
359 256
360 257 if nhaver > 1:
361 258 num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
362 259 deltha = 0.05*nhaver
363 260 minhvalid = pdata*nhaver
364 261 for im in range(self.dataOut.nmodes):
365 262 new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
366 263
367 264
368 265 data_fHeigths_List = []
369 266 data_fZonal_List = []
370 267 data_fMeridional_List = []
371 268 data_fVertical_List = []
372 269 startDTList = []
373 270
374 271
375 272 for i in range(ntime):
376 273 height = old_height
377 274
378 275 start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
379 276 stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
380 277
381 278
382 279 limit_sec1 = time.mktime(start.timetuple())
383 280 limit_sec2 = time.mktime(stop.timetuple())
384 281
385 282 t1 = numpy.where(self.f_timesec >= limit_sec1)
386 283 t2 = numpy.where(self.f_timesec < limit_sec2)
387 284 time_select = []
388 285 for val_sec in t1[0]:
389 286 if val_sec in t2[0]:
390 287 time_select.append(val_sec)
391 288
392 289
393 290 time_select = numpy.array(time_select,dtype = 'int')
394 291 minvalid = numpy.ceil(pdata*nhpoints)
395 292
396 293 zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
397 294 mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
398 295 ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
399 296
400 297 if nhaver > 1:
401 298 new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
402 299 new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
403 300 new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
404 301
405 302 if len(time_select) > minvalid:
406 303 time_average = self.f_timesec[time_select]
407 304
408 305 for im in range(self.dataOut.nmodes):
409 306
410 307 for ih in range(self.dataOut.nranges):
411 308 if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
412 309 zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
413 310
414 311 if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
415 312 mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
416 313
417 314 if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
418 315 ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
419 316
420 317 if nhaver > 1:
421 318 for ih in range(num_hei):
422 319 hvalid = numpy.arange(nhaver) + nhaver*ih
423 320
424 321 if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
425 322 new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) / numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
426 323
427 324 if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
428 325 new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) / numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
429 326
430 327 if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
431 328 new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) / numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
432 329 if nhaver > 1:
433 330 zon_aver = new_zon_aver
434 331 mer_aver = new_mer_aver
435 332 ver_aver = new_ver_aver
436 333 height = new_height
437 334
438 335
439 336 tstart = time_average[0]
440 337 tend = time_average[-1]
441 338 startTime = time.gmtime(tstart)
442 339
443 340 year = startTime.tm_year
444 341 month = startTime.tm_mon
445 342 day = startTime.tm_mday
446 343 hour = startTime.tm_hour
447 344 minute = startTime.tm_min
448 345 second = startTime.tm_sec
449 346
450 347 startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
451 348
452 349
453 350 o_height = numpy.array([])
454 351 o_zon_aver = numpy.array([])
455 352 o_mer_aver = numpy.array([])
456 353 o_ver_aver = numpy.array([])
457 354 if self.dataOut.nmodes > 1:
458 355 for im in range(self.dataOut.nmodes):
459 356
460 357 if im == 0:
461 358 h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
462 359 else:
463 360 h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
464 361
465 362
466 363 ht = h_select[0]
467 364
468 365 o_height = numpy.hstack((o_height,height[im,ht]))
469 366 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
470 367 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
471 368 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
472 369
473 370 data_fHeigths_List.append(o_height)
474 371 data_fZonal_List.append(o_zon_aver)
475 372 data_fMeridional_List.append(o_mer_aver)
476 373 data_fVertical_List.append(o_ver_aver)
477 374
478 375
479 376 else:
480 377 h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
481 378 ht = h_select[0]
482 379 o_height = numpy.hstack((o_height,height[im,ht]))
483 380 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
484 381 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
485 382 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
486 383
487 384 data_fHeigths_List.append(o_height)
488 385 data_fZonal_List.append(o_zon_aver)
489 386 data_fMeridional_List.append(o_mer_aver)
490 387 data_fVertical_List.append(o_ver_aver)
491 388
492 389
493 390 return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List
494
495
496 def prePlot(self,modeselect=None):
497 391
498 '''
499 Inputs:
500
501 self.dataOut.data_output - Zonal, Meridional and Vertical velocity array
502 self.dataOut.height - height array
503 self.dataOut.time - Time array (seconds)
504 self.dataOut.data_SNR - SNR array
505
506 '''
507
508 m = modeselect -1
509
510 print ' [Plotting mode {}]'.format(modeselect)
511 if not (m ==1 or m==0):
512 raise IndexError("'Mode' must be egual to : 1 or 2")
513 #
514 if self.flagfirstmode==0:
515 #copy of the data
516 self.data_output_copy = self.dataOut.data_output.copy()
517 self.data_height_copy = self.dataOut.height.copy()
518 self.data_time_copy = self.dataOut.time.copy()
519 self.data_SNR_copy = self.dataOut.data_SNR.copy()
520 self.flagfirstmode = 1
521
522 else:
523 self.dataOut.data_output = self.data_output_copy
524 self.dataOut.height = self.data_height_copy
525 self.dataOut.time = self.data_time_copy
526 self.dataOut.data_SNR = self.data_SNR_copy
527 self.flagfirstmode = 0
528
529
530 #select data for mode m
531 #self.dataOut.data_output = self.dataOut.data_output[:,:,m]
532 self.dataOut.heightList = self.dataOut.height[0,:]
533
534 data_SNR = self.dataOut.data_SNR[:,:,m]
535 self.dataOut.data_SNR= transpose(data_SNR)
536
537 if m==1 and self.dataOut.counter_records%2==0:
538 print '*********'
539 print 'MODO 2'
540 #print 'Zonal', self.dataOut.data_output[0]
541 #print 'Meridional', self.dataOut.data_output[1]
542 #print 'Vertical', self.dataOut.data_output[2]
543
544 print '*********'
545
546 Vx=self.dataOut.data_output[0,:,m]
547 Vy=self.dataOut.data_output[1,:,m]
548
549 Vmag=numpy.sqrt(Vx**2+Vy**2)
550 Vang=numpy.arctan2(Vy,Vx)
551 #print 'Vmag', Vmag
552 #print 'Vang', Vang
553
554 self.dataOut.data_output[0,:,m]=Vmag
555 self.dataOut.data_output[1,:,m]=Vang
556
557 prin= self.dataOut.data_output[0,:,m][~numpy.isnan(self.dataOut.data_output[0,:,m])]
558 print ' '
559 print 'VmagAverage',numpy.mean(prin)
560 print ' '
561 self.dataOut.data_output = self.dataOut.data_output[:,:,m]
562
563 392
564 393 No newline at end of file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file, binary diff hidden
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed, binary diff hidden
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed, binary diff hidden
General Comments 0
You need to be logged in to leave comments. Login now