##// END OF EJS Templates
schain-jc
RSS
-Functional HDF5 Format Writing and Reading Unit...
Julio Valdez -
r804:4c41460abf8c
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,1169 +1,1169
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
13 13 def getNumpyDtype(dataTypeCode):
14 14
15 15 if dataTypeCode == 0:
16 16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 17 elif dataTypeCode == 1:
18 18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 19 elif dataTypeCode == 2:
20 20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 21 elif dataTypeCode == 3:
22 22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 23 elif dataTypeCode == 4:
24 24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 25 elif dataTypeCode == 5:
26 26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 27 else:
28 28 raise ValueError, 'dataTypeCode was not defined'
29 29
30 30 return numpyDtype
31 31
32 32 def getDataTypeCode(numpyDtype):
33 33
34 34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 35 datatype = 0
36 36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 37 datatype = 1
38 38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 39 datatype = 2
40 40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 41 datatype = 3
42 42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 43 datatype = 4
44 44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 45 datatype = 5
46 46 else:
47 47 datatype = None
48 48
49 49 return datatype
50 50
51 51 def hildebrand_sekhon(data, navg):
52 52 """
53 53 This method is for the objective determination of the noise level in Doppler spectra. This
54 54 implementation technique is based on the fact that the standard deviation of the spectral
55 55 densities is equal to the mean spectral density for white Gaussian noise
56 56
57 57 Inputs:
58 58 Data : heights
59 59 navg : numbers of averages
60 60
61 61 Return:
62 62 -1 : any error
63 63 anoise : noise's level
64 64 """
65 65
66 66 sortdata = numpy.sort(data,axis=None)
67 67 lenOfData = len(sortdata)
68 68 nums_min = lenOfData*0.2
69 69
70 70 if nums_min <= 5:
71 71 nums_min = 5
72 72
73 73 sump = 0.
74 74
75 75 sumq = 0.
76 76
77 77 j = 0
78 78
79 79 cont = 1
80 80
81 81 while((cont==1)and(j<lenOfData)):
82 82
83 83 sump += sortdata[j]
84 84
85 85 sumq += sortdata[j]**2
86 86
87 87 if j > nums_min:
88 88 rtest = float(j)/(j-1) + 1.0/navg
89 89 if ((sumq*j) > (rtest*sump**2)):
90 90 j = j - 1
91 91 sump = sump - sortdata[j]
92 92 sumq = sumq - sortdata[j]**2
93 93 cont = 0
94 94
95 95 j += 1
96 96
97 97 lnoise = sump /j
98 98 # stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
99 99 return lnoise
100 100
101 101 class Beam:
102 102 def __init__(self):
103 103 self.codeList = []
104 104 self.azimuthList = []
105 105 self.zenithList = []
106 106
107 107 class GenericData(object):
108 108
109 109 flagNoData = True
110 110
111 111 def __init__(self):
112 112
113 113 raise NotImplementedError
114 114
115 115 def copy(self, inputObj=None):
116 116
117 117 if inputObj == None:
118 118 return copy.deepcopy(self)
119 119
120 120 for key in inputObj.__dict__.keys():
121 121
122 122 attribute = inputObj.__dict__[key]
123 123
124 124 #If this attribute is a tuple or list
125 125 if type(inputObj.__dict__[key]) in (tuple, list):
126 126 self.__dict__[key] = attribute[:]
127 127 continue
128 128
129 129 #If this attribute is another object or instance
130 130 if hasattr(attribute, '__dict__'):
131 131 self.__dict__[key] = attribute.copy()
132 132 continue
133 133
134 134 self.__dict__[key] = inputObj.__dict__[key]
135 135
136 136 def deepcopy(self):
137 137
138 138 return copy.deepcopy(self)
139 139
140 140 def isEmpty(self):
141 141
142 142 return self.flagNoData
143 143
144 144 class JROData(GenericData):
145 145
146 146 # m_BasicHeader = BasicHeader()
147 147 # m_ProcessingHeader = ProcessingHeader()
148 148
149 149 systemHeaderObj = SystemHeader()
150 150
151 151 radarControllerHeaderObj = RadarControllerHeader()
152 152
153 153 # data = None
154 154
155 155 type = None
156 156
157 157 datatype = None #dtype but in string
158 158
159 159 # dtype = None
160 160
161 161 # nChannels = None
162 162
163 163 # nHeights = None
164 164
165 165 nProfiles = None
166 166
167 167 heightList = None
168 168
169 169 channelList = None
170 170
171 171 flagDiscontinuousBlock = False
172 172
173 173 useLocalTime = False
174 174
175 175 utctime = None
176 176
177 177 timeZone = None
178 178
179 179 dstFlag = None
180 180
181 181 errorCount = None
182 182
183 183 blocksize = None
184 184
185 185 # nCode = None
186 186 #
187 187 # nBaud = None
188 188 #
189 189 # code = None
190 190
191 191 flagDecodeData = False #asumo q la data no esta decodificada
192 192
193 193 flagDeflipData = False #asumo q la data no esta sin flip
194 194
195 195 flagShiftFFT = False
196 196
197 197 # ippSeconds = None
198 198
199 199 # timeInterval = None
200 200
201 201 nCohInt = None
202 202
203 203 # noise = None
204 204
205 205 windowOfFilter = 1
206 206
207 207 #Speed of ligth
208 208 C = 3e8
209 209
210 210 frequency = 49.92e6
211 211
212 212 realtime = False
213 213
214 214 beacon_heiIndexList = None
215 215
216 216 last_block = None
217 217
218 218 blocknow = None
219 219
220 220 azimuth = None
221 221
222 222 zenith = None
223 223
224 224 beam = Beam()
225 225
226 226 profileIndex = None
227 227
228 228 def __init__(self):
229 229
230 230 raise NotImplementedError
231 231
232 232 def getNoise(self):
233 233
234 234 raise NotImplementedError
235 235
236 236 def getNChannels(self):
237 237
238 238 return len(self.channelList)
239 239
240 240 def getChannelIndexList(self):
241 241
242 242 return range(self.nChannels)
243 243
244 244 def getNHeights(self):
245 245
246 246 return len(self.heightList)
247 247
248 248 def getHeiRange(self, extrapoints=0):
249 249
250 250 heis = self.heightList
251 251 # deltah = self.heightList[1] - self.heightList[0]
252 252 #
253 253 # heis.append(self.heightList[-1])
254 254
255 255 return heis
256 256
257 257 def getDeltaH(self):
258 258
259 259 delta = self.heightList[1] - self.heightList[0]
260 260
261 261 return delta
262 262
263 263 def getltctime(self):
264 264
265 265 if self.useLocalTime:
266 266 return self.utctime - self.timeZone*60
267 267
268 268 return self.utctime
269 269
270 270 def getDatatime(self):
271 271
272 272 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
273 273 return datatimeValue
274 274
275 275 def getTimeRange(self):
276 276
277 277 datatime = []
278 278
279 279 datatime.append(self.ltctime)
280 280 datatime.append(self.ltctime + self.timeInterval+60)
281 281
282 282 datatime = numpy.array(datatime)
283 283
284 284 return datatime
285 285
286 286 def getFmaxTimeResponse(self):
287 287
288 288 period = (10**-6)*self.getDeltaH()/(0.15)
289 289
290 290 PRF = 1./(period * self.nCohInt)
291 291
292 292 fmax = PRF
293 293
294 294 return fmax
295 295
296 296 def getFmax(self):
297 297
298 298 PRF = 1./(self.ippSeconds * self.nCohInt)
299 299
300 300 fmax = PRF
301 301
302 302 return fmax
303 303
304 304 def getVmax(self):
305 305
306 306 _lambda = self.C/self.frequency
307 307
308 308 vmax = self.getFmax() * _lambda
309 309
310 310 return vmax
311 311
312 312 def get_ippSeconds(self):
313 313 '''
314 314 '''
315 315 return self.radarControllerHeaderObj.ippSeconds
316 316
317 317 def set_ippSeconds(self, ippSeconds):
318 318 '''
319 319 '''
320 320
321 321 self.radarControllerHeaderObj.ippSeconds = ippSeconds
322 322
323 323 return
324 324
325 325 def get_dtype(self):
326 326 '''
327 327 '''
328 328 return getNumpyDtype(self.datatype)
329 329
330 330 def set_dtype(self, numpyDtype):
331 331 '''
332 332 '''
333 333
334 334 self.datatype = getDataTypeCode(numpyDtype)
335 335
336 336 def get_code(self):
337 337 '''
338 338 '''
339 339 return self.radarControllerHeaderObj.code
340 340
341 341 def set_code(self, code):
342 342 '''
343 343 '''
344 344 self.radarControllerHeaderObj.code = code
345 345
346 346 return
347 347
348 348 def get_ncode(self):
349 349 '''
350 350 '''
351 351 return self.radarControllerHeaderObj.nCode
352 352
353 353 def set_ncode(self, nCode):
354 354 '''
355 355 '''
356 356 self.radarControllerHeaderObj.nCode = nCode
357 357
358 358 return
359 359
360 360 def get_nbaud(self):
361 361 '''
362 362 '''
363 363 return self.radarControllerHeaderObj.nBaud
364 364
365 365 def set_nbaud(self, nBaud):
366 366 '''
367 367 '''
368 368 self.radarControllerHeaderObj.nBaud = nBaud
369 369
370 370 return
371 371
372 372 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
373 373 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
374 374 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
375 375 #noise = property(getNoise, "I'm the 'nHeights' property.")
376 376 datatime = property(getDatatime, "I'm the 'datatime' property")
377 377 ltctime = property(getltctime, "I'm the 'ltctime' property")
378 378 ippSeconds = property(get_ippSeconds, set_ippSeconds)
379 379 dtype = property(get_dtype, set_dtype)
380 380 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
381 381 code = property(get_code, set_code)
382 382 nCode = property(get_ncode, set_ncode)
383 383 nBaud = property(get_nbaud, set_nbaud)
384 384
385 385 class Voltage(JROData):
386 386
387 387 #data es un numpy array de 2 dmensiones (canales, alturas)
388 388 data = None
389 389
390 390 def __init__(self):
391 391 '''
392 392 Constructor
393 393 '''
394 394
395 395 self.useLocalTime = True
396 396
397 397 self.radarControllerHeaderObj = RadarControllerHeader()
398 398
399 399 self.systemHeaderObj = SystemHeader()
400 400
401 401 self.type = "Voltage"
402 402
403 403 self.data = None
404 404
405 405 # self.dtype = None
406 406
407 407 # self.nChannels = 0
408 408
409 409 # self.nHeights = 0
410 410
411 411 self.nProfiles = None
412 412
413 413 self.heightList = None
414 414
415 415 self.channelList = None
416 416
417 417 # self.channelIndexList = None
418 418
419 419 self.flagNoData = True
420 420
421 421 self.flagDiscontinuousBlock = False
422 422
423 423 self.utctime = None
424 424
425 425 self.timeZone = None
426 426
427 427 self.dstFlag = None
428 428
429 429 self.errorCount = None
430 430
431 431 self.nCohInt = None
432 432
433 433 self.blocksize = None
434 434
435 435 self.flagDecodeData = False #asumo q la data no esta decodificada
436 436
437 437 self.flagDeflipData = False #asumo q la data no esta sin flip
438 438
439 439 self.flagShiftFFT = False
440 440
441 441 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
442 442
443 443 self.profileIndex = 0
444 444
445 445 def getNoisebyHildebrand(self, channel = None):
446 446 """
447 447 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
448 448
449 449 Return:
450 450 noiselevel
451 451 """
452 452
453 453 if channel != None:
454 454 data = self.data[channel]
455 455 nChannels = 1
456 456 else:
457 457 data = self.data
458 458 nChannels = self.nChannels
459 459
460 460 noise = numpy.zeros(nChannels)
461 461 power = data * numpy.conjugate(data)
462 462
463 463 for thisChannel in range(nChannels):
464 464 if nChannels == 1:
465 465 daux = power[:].real
466 466 else:
467 467 daux = power[thisChannel,:].real
468 468 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
469 469
470 470 return noise
471 471
472 472 def getNoise(self, type = 1, channel = None):
473 473
474 474 if type == 1:
475 475 noise = self.getNoisebyHildebrand(channel)
476 476
477 477 return noise
478 478
479 479 def getPower(self, channel = None):
480 480
481 481 if channel != None:
482 482 data = self.data[channel]
483 483 else:
484 484 data = self.data
485 485
486 486 power = data * numpy.conjugate(data)
487 487
488 488 return 10*numpy.log10(power.real)
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 es un numpy array de 2 dmensiones (canales, perfiles, alturas)
502 502 data_spc = None
503 503
504 504 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
505 505 data_cspc = None
506 506
507 507 #data es un numpy array de 2 dmensiones (canales, alturas)
508 508 data_dc = None
509 509
510 510 nFFTPoints = None
511 511
512 512 # nPairs = None
513 513
514 514 pairsList = None
515 515
516 516 nIncohInt = None
517 517
518 518 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
519 519
520 520 nCohInt = None #se requiere para determinar el valor de timeInterval
521 521
522 522 ippFactor = None
523 523
524 524 profileIndex = 0
525 525
526 526 plotting = "spectra"
527 527
528 528 def __init__(self):
529 529 '''
530 530 Constructor
531 531 '''
532 532
533 533 self.useLocalTime = True
534 534
535 535 self.radarControllerHeaderObj = RadarControllerHeader()
536 536
537 537 self.systemHeaderObj = SystemHeader()
538 538
539 539 self.type = "Spectra"
540 540
541 541 # self.data = None
542 542
543 543 # self.dtype = None
544 544
545 545 # self.nChannels = 0
546 546
547 547 # self.nHeights = 0
548 548
549 549 self.nProfiles = None
550 550
551 551 self.heightList = None
552 552
553 553 self.channelList = None
554 554
555 555 # self.channelIndexList = None
556 556
557 557 self.pairsList = None
558 558
559 559 self.flagNoData = True
560 560
561 561 self.flagDiscontinuousBlock = False
562 562
563 563 self.utctime = None
564 564
565 565 self.nCohInt = None
566 566
567 567 self.nIncohInt = None
568 568
569 569 self.blocksize = None
570 570
571 571 self.nFFTPoints = None
572 572
573 573 self.wavelength = None
574 574
575 575 self.flagDecodeData = False #asumo q la data no esta decodificada
576 576
577 577 self.flagDeflipData = False #asumo q la data no esta sin flip
578 578
579 579 self.flagShiftFFT = False
580 580
581 581 self.ippFactor = 1
582 582
583 583 #self.noise = None
584 584
585 585 self.beacon_heiIndexList = []
586 586
587 587 self.noise_estimation = None
588 588
589 589
590 590 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
591 591 """
592 592 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
593 593
594 594 Return:
595 595 noiselevel
596 596 """
597 597
598 598 noise = numpy.zeros(self.nChannels)
599 599
600 600 for channel in range(self.nChannels):
601 601 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
602 602 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
603 603
604 604 return noise
605 605
606 606 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
607 607
608 608 if self.noise_estimation is not None:
609 609 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
610 610 else:
611 611 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
612 612 return noise
613 613
614 614 def getFreqRangeTimeResponse(self, extrapoints=0):
615 615
616 616 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
617 617 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
618 618
619 619 return freqrange
620 620
621 621 def getAcfRange(self, extrapoints=0):
622 622
623 623 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
624 624 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
625 625
626 626 return freqrange
627 627
628 628 def getFreqRange(self, extrapoints=0):
629 629
630 630 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
631 631 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
632 632
633 633 return freqrange
634 634
635 635 def getVelRange(self, extrapoints=0):
636 636
637 637 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
638 638 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
639 639
640 640 return velrange
641 641
642 642 def getNPairs(self):
643 643
644 644 return len(self.pairsList)
645 645
646 646 def getPairsIndexList(self):
647 647
648 648 return range(self.nPairs)
649 649
650 650 def getNormFactor(self):
651 651
652 652 pwcode = 1
653 653
654 654 if self.flagDecodeData:
655 655 pwcode = numpy.sum(self.code[0]**2)
656 656 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
657 657 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
658 658
659 659 return normFactor
660 660
661 661 def getFlagCspc(self):
662 662
663 663 if self.data_cspc is None:
664 664 return True
665 665
666 666 return False
667 667
668 668 def getFlagDc(self):
669 669
670 670 if self.data_dc is None:
671 671 return True
672 672
673 673 return False
674 674
675 675 def getTimeInterval(self):
676 676
677 677 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
678 678
679 679 return timeInterval
680 680
681 681 def setValue(self, value):
682 682
683 683 print "This property should not be initialized"
684 684
685 685 return
686 686
687 687 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
688 688 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
689 689 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
690 690 flag_cspc = property(getFlagCspc, setValue)
691 691 flag_dc = property(getFlagDc, setValue)
692 692 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
693 693 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
694 694
695 695 class SpectraHeis(Spectra):
696 696
697 697 data_spc = None
698 698
699 699 data_cspc = None
700 700
701 701 data_dc = None
702 702
703 703 nFFTPoints = None
704 704
705 705 # nPairs = None
706 706
707 707 pairsList = None
708 708
709 709 nCohInt = None
710 710
711 711 nIncohInt = None
712 712
713 713 def __init__(self):
714 714
715 715 self.radarControllerHeaderObj = RadarControllerHeader()
716 716
717 717 self.systemHeaderObj = SystemHeader()
718 718
719 719 self.type = "SpectraHeis"
720 720
721 721 # self.dtype = None
722 722
723 723 # self.nChannels = 0
724 724
725 725 # self.nHeights = 0
726 726
727 727 self.nProfiles = None
728 728
729 729 self.heightList = None
730 730
731 731 self.channelList = None
732 732
733 733 # self.channelIndexList = None
734 734
735 735 self.flagNoData = True
736 736
737 737 self.flagDiscontinuousBlock = False
738 738
739 739 # self.nPairs = 0
740 740
741 741 self.utctime = None
742 742
743 743 self.blocksize = None
744 744
745 745 self.profileIndex = 0
746 746
747 747 self.nCohInt = 1
748 748
749 749 self.nIncohInt = 1
750 750
751 751 def getNormFactor(self):
752 752 pwcode = 1
753 753 if self.flagDecodeData:
754 754 pwcode = numpy.sum(self.code[0]**2)
755 755
756 756 normFactor = self.nIncohInt*self.nCohInt*pwcode
757 757
758 758 return normFactor
759 759
760 760 def getTimeInterval(self):
761 761
762 762 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
763 763
764 764 return timeInterval
765 765
766 766 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
767 767 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
768 768
769 769 class Fits(JROData):
770 770
771 771 heightList = None
772 772
773 773 channelList = None
774 774
775 775 flagNoData = True
776 776
777 777 flagDiscontinuousBlock = False
778 778
779 779 useLocalTime = False
780 780
781 781 utctime = None
782 782
783 783 timeZone = None
784 784
785 785 # ippSeconds = None
786 786
787 787 # timeInterval = None
788 788
789 789 nCohInt = None
790 790
791 791 nIncohInt = None
792 792
793 793 noise = None
794 794
795 795 windowOfFilter = 1
796 796
797 797 #Speed of ligth
798 798 C = 3e8
799 799
800 800 frequency = 49.92e6
801 801
802 802 realtime = False
803 803
804 804
805 805 def __init__(self):
806 806
807 807 self.type = "Fits"
808 808
809 809 self.nProfiles = None
810 810
811 811 self.heightList = None
812 812
813 813 self.channelList = None
814 814
815 815 # self.channelIndexList = None
816 816
817 817 self.flagNoData = True
818 818
819 819 self.utctime = None
820 820
821 821 self.nCohInt = 1
822 822
823 823 self.nIncohInt = 1
824 824
825 825 self.useLocalTime = True
826 826
827 827 self.profileIndex = 0
828 828
829 829 # self.utctime = None
830 830 # self.timeZone = None
831 831 # self.ltctime = None
832 832 # self.timeInterval = None
833 833 # self.header = None
834 834 # self.data_header = None
835 835 # self.data = None
836 836 # self.datatime = None
837 837 # self.flagNoData = False
838 838 # self.expName = ''
839 839 # self.nChannels = None
840 840 # self.nSamples = None
841 841 # self.dataBlocksPerFile = None
842 842 # self.comments = ''
843 843 #
844 844
845 845
846 846 def getltctime(self):
847 847
848 848 if self.useLocalTime:
849 849 return self.utctime - self.timeZone*60
850 850
851 851 return self.utctime
852 852
853 853 def getDatatime(self):
854 854
855 855 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
856 856 return datatime
857 857
858 858 def getTimeRange(self):
859 859
860 860 datatime = []
861 861
862 862 datatime.append(self.ltctime)
863 863 datatime.append(self.ltctime + self.timeInterval)
864 864
865 865 datatime = numpy.array(datatime)
866 866
867 867 return datatime
868 868
869 869 def getHeiRange(self):
870 870
871 871 heis = self.heightList
872 872
873 873 return heis
874 874
875 875 def getNHeights(self):
876 876
877 877 return len(self.heightList)
878 878
879 879 def getNChannels(self):
880 880
881 881 return len(self.channelList)
882 882
883 883 def getChannelIndexList(self):
884 884
885 885 return range(self.nChannels)
886 886
887 887 def getNoise(self, type = 1):
888 888
889 889 #noise = numpy.zeros(self.nChannels)
890 890
891 891 if type == 1:
892 892 noise = self.getNoisebyHildebrand()
893 893
894 894 if type == 2:
895 895 noise = self.getNoisebySort()
896 896
897 897 if type == 3:
898 898 noise = self.getNoisebyWindow()
899 899
900 900 return noise
901 901
902 902 def getTimeInterval(self):
903 903
904 904 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
905 905
906 906 return timeInterval
907 907
908 908 datatime = property(getDatatime, "I'm the 'datatime' property")
909 909 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
910 910 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
911 911 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
912 912 noise = property(getNoise, "I'm the 'nHeights' property.")
913 913
914 914 ltctime = property(getltctime, "I'm the 'ltctime' property")
915 915 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
916 916
917 917 class Correlation(JROData):
918 918
919 919 noise = None
920 920
921 921 SNR = None
922 922
923 923 pairsAutoCorr = None #Pairs of Autocorrelation
924 924
925 925 #--------------------------------------------------
926 926
927 927 data_corr = None
928 928
929 929 data_volt = None
930 930
931 931 lagT = None # each element value is a profileIndex
932 932
933 933 lagR = None # each element value is in km
934 934
935 935 pairsList = None
936 936
937 937 calculateVelocity = None
938 938
939 939 nPoints = None
940 940
941 941 nAvg = None
942 942
943 943 bufferSize = None
944 944
945 945 def __init__(self):
946 946 '''
947 947 Constructor
948 948 '''
949 949 self.radarControllerHeaderObj = RadarControllerHeader()
950 950
951 951 self.systemHeaderObj = SystemHeader()
952 952
953 953 self.type = "Correlation"
954 954
955 955 self.data = None
956 956
957 957 self.dtype = None
958 958
959 959 self.nProfiles = None
960 960
961 961 self.heightList = None
962 962
963 963 self.channelList = None
964 964
965 965 self.flagNoData = True
966 966
967 967 self.flagDiscontinuousBlock = False
968 968
969 969 self.utctime = None
970 970
971 971 self.timeZone = None
972 972
973 973 self.dstFlag = None
974 974
975 975 self.errorCount = None
976 976
977 977 self.blocksize = None
978 978
979 979 self.flagDecodeData = False #asumo q la data no esta decodificada
980 980
981 981 self.flagDeflipData = False #asumo q la data no esta sin flip
982 982
983 983 self.pairsList = None
984 984
985 985 self.nPoints = None
986 986
987 987 def getLagTRange(self, extrapoints=0):
988 988
989 989 lagTRange = self.lagT
990 990 diff = lagTRange[1] - lagTRange[0]
991 991 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
992 992 lagTRange = numpy.hstack((lagTRange, extra))
993 993
994 994 return lagTRange
995 995
996 996 def getLagRRange(self, extrapoints=0):
997 997
998 998 return self.lagR
999 999
1000 1000 def getPairsList(self):
1001 1001
1002 1002 return self.pairsList
1003 1003
1004 1004 def getCalculateVelocity(self):
1005 1005
1006 1006 return self.calculateVelocity
1007 1007
1008 1008 def getNPoints(self):
1009 1009
1010 1010 return self.nPoints
1011 1011
1012 1012 def getNAvg(self):
1013 1013
1014 1014 return self.nAvg
1015 1015
1016 1016 def getBufferSize(self):
1017 1017
1018 1018 return self.bufferSize
1019 1019
1020 1020 def getPairsAutoCorr(self):
1021 1021 pairsList = self.pairsList
1022 1022 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
1023 1023
1024 1024 for l in range(len(pairsList)):
1025 1025 firstChannel = pairsList[l][0]
1026 1026 secondChannel = pairsList[l][1]
1027 1027
1028 1028 #Obteniendo pares de Autocorrelacion
1029 1029 if firstChannel == secondChannel:
1030 1030 pairsAutoCorr[firstChannel] = int(l)
1031 1031
1032 1032 pairsAutoCorr = pairsAutoCorr.astype(int)
1033 1033
1034 1034 return pairsAutoCorr
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.nPoints
1084 1084
1085 1085 return timeInterval
1086 1086
1087 1087 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1088 1088 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
1089 1089 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
1090 1090 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
1091 1091 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
1092 1092 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
1093 1093
1094 1094
1095 1095 class Parameters(JROData):
1096 1096
1097 1097 #Information from previous data
1098 1098
1099 1099 inputUnit = None #Type of data to be processed
1100 1100
1101 1101 operation = None #Type of operation to parametrize
1102 1102
1103 1103 normFactor = None #Normalization Factor
1104 1104
1105 1105 groupList = None #List of Pairs, Groups, etc
1106 1106
1107 1107 #Parameters
1108 1108
1109 1109 data_param = None #Parameters obtained
1110 1110
1111 1111 data_pre = None #Data Pre Parametrization
1112 1112
1113 1113 data_SNR = None #Signal to Noise Ratio
1114 1114
1115 1115 # heightRange = None #Heights
1116 1116
1117 1117 abscissaList = None #Abscissa, can be velocities, lags or time
1118 1118
1119 1119 noise = None #Noise Potency
1120 1120
1121 1121 utctimeInit = None #Initial UTC time
1122 1122
1123 1123 paramInterval = None #Time interval to calculate Parameters in seconds
1124 1124
1125 1125 useLocalTime = True
1126 1126
1127 1127 #Fitting
1128 1128
1129 1129 data_error = None #Error of the estimation
1130 1130
1131 1131 constants = None
1132 1132
1133 1133 library = None
1134 1134
1135 1135 #Output signal
1136 1136
1137 1137 outputInterval = None #Time interval to calculate output signal in seconds
1138 1138
1139 1139 data_output = None #Out signal
1140 1140
1141 1141
1142 1142
1143 1143 def __init__(self):
1144 1144 '''
1145 1145 Constructor
1146 1146 '''
1147 1147 self.radarControllerHeaderObj = RadarControllerHeader()
1148 1148
1149 1149 self.systemHeaderObj = SystemHeader()
1150 1150
1151 1151 self.type = "Parameters"
1152 1152
1153 def getTimeRange1(self):
1153 def getTimeRange1(self, interval):
1154 1154
1155 1155 datatime = []
1156 1156
1157 1157 if self.useLocalTime:
1158 1158 time1 = self.utctimeInit - self.timeZone*60
1159 1159 else:
1160 1160 time1 = self.utctimeInit
1161 1161
1162 1162 # datatime.append(self.utctimeInit)
1163 1163 # datatime.append(self.utctimeInit + self.outputInterval)
1164 1164 datatime.append(time1)
1165 datatime.append(time1 + self.outputInterval)
1165 datatime.append(time1 + interval)
1166 1166
1167 1167 datatime = numpy.array(datatime)
1168 1168
1169 1169 return datatime
Show file before | Show file after | Hide comments
@@ -1,1373 +1,1364
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from figure import Figure, isRealtime
6 6 from plotting_codes import *
7 7
8 8 class MomentsPlot(Figure):
9 9
10 10 isConfig = None
11 11 __nsubplots = None
12 12
13 13 WIDTHPROF = None
14 14 HEIGHTPROF = None
15 15 PREFIX = 'prm'
16 16
17 17 def __init__(self):
18 18
19 19 self.isConfig = False
20 20 self.__nsubplots = 1
21 21
22 22 self.WIDTH = 280
23 23 self.HEIGHT = 250
24 24 self.WIDTHPROF = 120
25 25 self.HEIGHTPROF = 0
26 26 self.counter_imagwr = 0
27 27
28 28 self.PLOT_CODE = MOMENTS_CODE
29 29
30 30 self.FTP_WEI = None
31 31 self.EXP_CODE = None
32 32 self.SUB_EXP_CODE = None
33 33 self.PLOT_POS = None
34 34
35 35 def getSubplots(self):
36 36
37 37 ncol = int(numpy.sqrt(self.nplots)+0.9)
38 38 nrow = int(self.nplots*1./ncol + 0.9)
39 39
40 40 return nrow, ncol
41 41
42 42 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
43 43
44 44 self.__showprofile = showprofile
45 45 self.nplots = nplots
46 46
47 47 ncolspan = 1
48 48 colspan = 1
49 49 if showprofile:
50 50 ncolspan = 3
51 51 colspan = 2
52 52 self.__nsubplots = 2
53 53
54 54 self.createFigure(id = id,
55 55 wintitle = wintitle,
56 56 widthplot = self.WIDTH + self.WIDTHPROF,
57 57 heightplot = self.HEIGHT + self.HEIGHTPROF,
58 58 show=show)
59 59
60 60 nrow, ncol = self.getSubplots()
61 61
62 62 counter = 0
63 63 for y in range(nrow):
64 64 for x in range(ncol):
65 65
66 66 if counter >= self.nplots:
67 67 break
68 68
69 69 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
70 70
71 71 if showprofile:
72 72 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
73 73
74 74 counter += 1
75 75
76 76 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
77 77 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
78 78 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
79 79 server=None, folder=None, username=None, password=None,
80 80 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
81 81
82 82 """
83 83
84 84 Input:
85 85 dataOut :
86 86 id :
87 87 wintitle :
88 88 channelList :
89 89 showProfile :
90 90 xmin : None,
91 91 xmax : None,
92 92 ymin : None,
93 93 ymax : None,
94 94 zmin : None,
95 95 zmax : None
96 96 """
97 97
98 98 if dataOut.flagNoData:
99 99 return None
100 100
101 101 if realtime:
102 102 if not(isRealtime(utcdatatime = dataOut.utctime)):
103 103 print 'Skipping this plot function'
104 104 return
105 105
106 106 if channelList == None:
107 107 channelIndexList = dataOut.channelIndexList
108 108 else:
109 109 channelIndexList = []
110 110 for channel in channelList:
111 111 if channel not in dataOut.channelList:
112 112 raise ValueError, "Channel %d is not in dataOut.channelList"
113 113 channelIndexList.append(dataOut.channelList.index(channel))
114 114
115 115 factor = dataOut.normFactor
116 116 x = dataOut.abscissaList
117 117 y = dataOut.heightList
118 118
119 119 z = dataOut.data_pre[channelIndexList,:,:]/factor
120 120 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
121 121 avg = numpy.average(z, axis=1)
122 122 noise = dataOut.noise/factor
123 123
124 124 zdB = 10*numpy.log10(z)
125 125 avgdB = 10*numpy.log10(avg)
126 126 noisedB = 10*numpy.log10(noise)
127 127
128 128 #thisDatetime = dataOut.datatime
129 129 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
130 130 title = wintitle + " Parameters"
131 131 xlabel = "Velocity (m/s)"
132 132 ylabel = "Range (Km)"
133 133
134 134 update_figfile = False
135 135
136 136 if not self.isConfig:
137 137
138 138 nplots = len(channelIndexList)
139 139
140 140 self.setup(id=id,
141 141 nplots=nplots,
142 142 wintitle=wintitle,
143 143 showprofile=showprofile,
144 144 show=show)
145 145
146 146 if xmin == None: xmin = numpy.nanmin(x)
147 147 if xmax == None: xmax = numpy.nanmax(x)
148 148 if ymin == None: ymin = numpy.nanmin(y)
149 149 if ymax == None: ymax = numpy.nanmax(y)
150 150 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
151 151 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
152 152
153 153 self.FTP_WEI = ftp_wei
154 154 self.EXP_CODE = exp_code
155 155 self.SUB_EXP_CODE = sub_exp_code
156 156 self.PLOT_POS = plot_pos
157 157
158 158 self.isConfig = True
159 159 update_figfile = True
160 160
161 161 self.setWinTitle(title)
162 162
163 163 for i in range(self.nplots):
164 164 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
165 165 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i], noisedB[i], str_datetime)
166 166 axes = self.axesList[i*self.__nsubplots]
167 167 axes.pcolor(x, y, zdB[i,:,:],
168 168 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
169 169 xlabel=xlabel, ylabel=ylabel, title=title,
170 170 ticksize=9, cblabel='')
171 171 #Mean Line
172 172 mean = dataOut.data_param[i, 1, :]
173 173 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
174 174
175 175 if self.__showprofile:
176 176 axes = self.axesList[i*self.__nsubplots +1]
177 177 axes.pline(avgdB[i], y,
178 178 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
179 179 xlabel='dB', ylabel='', title='',
180 180 ytick_visible=False,
181 181 grid='x')
182 182
183 183 noiseline = numpy.repeat(noisedB[i], len(y))
184 184 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
185 185
186 186 self.draw()
187 187
188 188 self.save(figpath=figpath,
189 189 figfile=figfile,
190 190 save=save,
191 191 ftp=ftp,
192 192 wr_period=wr_period,
193 193 thisDatetime=thisDatetime)
194 194
195 195
196 196
197 197 class SkyMapPlot(Figure):
198 198
199 199 __isConfig = None
200 200 __nsubplots = None
201 201
202 202 WIDTHPROF = None
203 203 HEIGHTPROF = None
204 204 PREFIX = 'mmap'
205 205
206 206 def __init__(self):
207 207
208 208 self.isConfig = False
209 209 self.__nsubplots = 1
210 210
211 211 # self.WIDTH = 280
212 212 # self.HEIGHT = 250
213 213 self.WIDTH = 600
214 214 self.HEIGHT = 600
215 215 self.WIDTHPROF = 120
216 216 self.HEIGHTPROF = 0
217 217 self.counter_imagwr = 0
218 218
219 219 self.PLOT_CODE = MSKYMAP_CODE
220 220
221 221 self.FTP_WEI = None
222 222 self.EXP_CODE = None
223 223 self.SUB_EXP_CODE = None
224 224 self.PLOT_POS = None
225 225
226 226 def getSubplots(self):
227 227
228 228 ncol = int(numpy.sqrt(self.nplots)+0.9)
229 229 nrow = int(self.nplots*1./ncol + 0.9)
230 230
231 231 return nrow, ncol
232 232
233 233 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
234 234
235 235 self.__showprofile = showprofile
236 236 self.nplots = nplots
237 237
238 238 ncolspan = 1
239 239 colspan = 1
240 240
241 241 self.createFigure(id = id,
242 242 wintitle = wintitle,
243 243 widthplot = self.WIDTH, #+ self.WIDTHPROF,
244 244 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
245 245 show=show)
246 246
247 247 nrow, ncol = 1,1
248 248 counter = 0
249 249 x = 0
250 250 y = 0
251 251 self.addAxes(1, 1, 0, 0, 1, 1, True)
252 252
253 253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
254 tmin=None, tmax=None, timerange=None,
254 tmin=0, tmax=24, timerange=None,
255 255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
256 256 server=None, folder=None, username=None, password=None,
257 257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
258 258
259 259 """
260 260
261 261 Input:
262 262 dataOut :
263 263 id :
264 264 wintitle :
265 265 channelList :
266 266 showProfile :
267 267 xmin : None,
268 268 xmax : None,
269 269 ymin : None,
270 270 ymax : None,
271 271 zmin : None,
272 272 zmax : None
273 273 """
274 274
275 arrayParameters = dataOut.data_param[0,:]
275 arrayParameters = dataOut.data_param
276 276 error = arrayParameters[:,-1]
277 277 indValid = numpy.where(error == 0)[0]
278 278 finalMeteor = arrayParameters[indValid,:]
279 finalAzimuth = finalMeteor[:,4]
280 finalZenith = finalMeteor[:,5]
279 finalAzimuth = finalMeteor[:,3]
280 finalZenith = finalMeteor[:,4]
281 281
282 282 x = finalAzimuth*numpy.pi/180
283 283 y = finalZenith
284 x1 = dataOut.getTimeRange()
284 x1 = [dataOut.ltctime, dataOut.ltctime]
285 285
286 286 #thisDatetime = dataOut.datatime
287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
288 288 title = wintitle + " Parameters"
289 289 xlabel = "Zonal Zenith Angle (deg) "
290 290 ylabel = "Meridional Zenith Angle (deg)"
291 291 update_figfile = False
292 292
293 293 if not self.isConfig:
294 294
295 295 nplots = 1
296 296
297 297 self.setup(id=id,
298 298 nplots=nplots,
299 299 wintitle=wintitle,
300 300 showprofile=showprofile,
301 301 show=show)
302 302
303 303 if self.xmin is None and self.xmax is None:
304 304 self.xmin, self.xmax = self.getTimeLim(x1, tmin, tmax, timerange)
305 305
306 306 if timerange != None:
307 307 self.timerange = timerange
308 308 else:
309 309 self.timerange = self.xmax - self.xmin
310 310
311 311 self.FTP_WEI = ftp_wei
312 312 self.EXP_CODE = exp_code
313 313 self.SUB_EXP_CODE = sub_exp_code
314 314 self.PLOT_POS = plot_pos
315 315 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
316 316 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
317 317 self.isConfig = True
318 318 update_figfile = True
319 319
320 320 self.setWinTitle(title)
321 321
322 322 i = 0
323 323 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
324 324
325 325 axes = self.axesList[i*self.__nsubplots]
326 326 nevents = axes.x_buffer.shape[0] + x.shape[0]
327 327 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
328 328 axes.polar(x, y,
329 329 title=title, xlabel=xlabel, ylabel=ylabel,
330 330 ticksize=9, cblabel='')
331 331
332 332 self.draw()
333 333
334 334 self.save(figpath=figpath,
335 335 figfile=figfile,
336 336 save=save,
337 337 ftp=ftp,
338 338 wr_period=wr_period,
339 339 thisDatetime=thisDatetime,
340 340 update_figfile=update_figfile)
341 341
342 342 if dataOut.ltctime >= self.xmax:
343 343 self.isConfigmagwr = wr_period
344 344 self.isConfig = False
345 345 update_figfile = True
346 346 axes.__firsttime = True
347 347 self.xmin += self.timerange
348 348 self.xmax += self.timerange
349 349
350 350
351 351
352 352
353 353 class WindProfilerPlot(Figure):
354 354
355 355 __isConfig = None
356 356 __nsubplots = None
357 357
358 358 WIDTHPROF = None
359 359 HEIGHTPROF = None
360 360 PREFIX = 'wind'
361 361
362 362 def __init__(self):
363 363
364 364 self.timerange = None
365 365 self.isConfig = False
366 366 self.__nsubplots = 1
367 367
368 368 self.WIDTH = 800
369 369 self.HEIGHT = 150
370 370 self.WIDTHPROF = 120
371 371 self.HEIGHTPROF = 0
372 372 self.counter_imagwr = 0
373 373
374 374 self.PLOT_CODE = WIND_CODE
375 375
376 376 self.FTP_WEI = None
377 377 self.EXP_CODE = None
378 378 self.SUB_EXP_CODE = None
379 379 self.PLOT_POS = None
380 380 self.tmin = None
381 381 self.tmax = None
382 382
383 383 self.xmin = None
384 384 self.xmax = None
385 385
386 386 self.figfile = None
387 387
388 388 def getSubplots(self):
389 389
390 390 ncol = 1
391 391 nrow = self.nplots
392 392
393 393 return nrow, ncol
394 394
395 395 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
396 396
397 397 self.__showprofile = showprofile
398 398 self.nplots = nplots
399 399
400 400 ncolspan = 1
401 401 colspan = 1
402 402
403 403 self.createFigure(id = id,
404 404 wintitle = wintitle,
405 405 widthplot = self.WIDTH + self.WIDTHPROF,
406 406 heightplot = self.HEIGHT + self.HEIGHTPROF,
407 407 show=show)
408 408
409 409 nrow, ncol = self.getSubplots()
410 410
411 411 counter = 0
412 412 for y in range(nrow):
413 413 if counter >= self.nplots:
414 414 break
415 415
416 416 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
417 417 counter += 1
418 418
419 419 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='False',
420 420 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
421 421 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
422 422 timerange=None, SNRthresh = None,
423 423 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
424 424 server=None, folder=None, username=None, password=None,
425 425 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
426 426 """
427 427
428 428 Input:
429 429 dataOut :
430 430 id :
431 431 wintitle :
432 432 channelList :
433 433 showProfile :
434 434 xmin : None,
435 435 xmax : None,
436 436 ymin : None,
437 437 ymax : None,
438 438 zmin : None,
439 439 zmax : None
440 440 """
441 441
442 if channelList == None:
443 channelIndexList = dataOut.channelIndexList
444 else:
445 channelIndexList = []
446 for channel in channelList:
447 if channel not in dataOut.channelList:
448 raise ValueError, "Channel %d is not in dataOut.channelList"
449 channelIndexList.append(dataOut.channelList.index(channel))
450
451 442 # if timerange is not None:
452 443 # self.timerange = timerange
453 444 #
454 445 # tmin = None
455 446 # tmax = None
456 447
457 x = dataOut.getTimeRange1()
448 x = dataOut.getTimeRange1(dataOut.outputInterval)
458 449 # y = dataOut.heightList
459 450 y = dataOut.heightList
460 451
461 452 z = dataOut.data_output.copy()
462 453 nplots = z.shape[0] #Number of wind dimensions estimated
463 454 nplotsw = nplots
464 455
465 456 #If there is a SNR function defined
466 457 if dataOut.data_SNR is not None:
467 458 nplots += 1
468 459 SNR = dataOut.data_SNR
469 460 SNRavg = numpy.average(SNR, axis=0)
470 461
471 462 SNRdB = 10*numpy.log10(SNR)
472 463 SNRavgdB = 10*numpy.log10(SNRavg)
473 464
474 465 if SNRthresh == None: SNRthresh = -5.0
475 466 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
476 467
477 468 for i in range(nplotsw):
478 469 z[i,ind] = numpy.nan
479 470
480 471
481 472 # showprofile = False
482 473 # thisDatetime = dataOut.datatime
483 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
474 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
484 475 title = wintitle + "Wind"
485 476 xlabel = ""
486 477 ylabel = "Range (Km)"
487 478 update_figfile = False
488 479
489 480 if not self.isConfig:
490 481
491 482 self.setup(id=id,
492 483 nplots=nplots,
493 484 wintitle=wintitle,
494 485 showprofile=showprofile,
495 486 show=show)
496 487
497 488 if timerange is not None:
498 489 self.timerange = timerange
499 490
500 491 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
501 492
502 493 if ymin == None: ymin = numpy.nanmin(y)
503 494 if ymax == None: ymax = numpy.nanmax(y)
504 495
505 496 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
506 497 #if numpy.isnan(zmax): zmax = 50
507 498 if zmin == None: zmin = -zmax
508 499
509 500 if nplotsw == 3:
510 501 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
511 502 if zmin_ver == None: zmin_ver = -zmax_ver
512 503
513 504 if dataOut.data_SNR is not None:
514 505 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
515 506 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
516 507
517 508
518 509 self.FTP_WEI = ftp_wei
519 510 self.EXP_CODE = exp_code
520 511 self.SUB_EXP_CODE = sub_exp_code
521 512 self.PLOT_POS = plot_pos
522 513
523 514 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
524 515 self.isConfig = True
525 516 self.figfile = figfile
526 517 update_figfile = True
527 518
528 519 self.setWinTitle(title)
529 520
530 521 if ((self.xmax - x[1]) < (x[1]-x[0])):
531 522 x[1] = self.xmax
532 523
533 524 strWind = ['Zonal', 'Meridional', 'Vertical']
534 525 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
535 526 zmaxVector = [zmax, zmax, zmax_ver]
536 527 zminVector = [zmin, zmin, zmin_ver]
537 528 windFactor = [1,1,100]
538 529
539 530 for i in range(nplotsw):
540 531
541 532 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
542 533 axes = self.axesList[i*self.__nsubplots]
543 534
544 535 z1 = z[i,:].reshape((1,-1))*windFactor[i]
545 536
546 537 axes.pcolorbuffer(x, y, z1,
547 538 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
548 539 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
549 540 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
550 541
551 542 if dataOut.data_SNR is not None:
552 543 i += 1
553 544 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
554 545 axes = self.axesList[i*self.__nsubplots]
555 546
556 547 SNRavgdB = SNRavgdB.reshape((1,-1))
557 548
558 549 axes.pcolorbuffer(x, y, SNRavgdB,
559 550 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
560 551 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
561 552 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
562 553
563 554 self.draw()
564 555
565 556 if dataOut.ltctime >= self.xmax:
566 557 self.counter_imagwr = wr_period
567 558 self.isConfig = False
568 559 update_figfile = True
569 560
570 561 self.save(figpath=figpath,
571 562 figfile=figfile,
572 563 save=save,
573 564 ftp=ftp,
574 565 wr_period=wr_period,
575 566 thisDatetime=thisDatetime,
576 567 update_figfile=update_figfile)
577 568
578 569
579 570
580 571 class ParametersPlot(Figure):
581 572
582 573 __isConfig = None
583 574 __nsubplots = None
584 575
585 576 WIDTHPROF = None
586 577 HEIGHTPROF = None
587 578 PREFIX = 'prm'
588 579
589 580 def __init__(self):
590 581
591 582 self.timerange = 2*60*60
592 583 self.isConfig = False
593 584 self.__nsubplots = 1
594 585
595 586 self.WIDTH = 800
596 587 self.HEIGHT = 150
597 588 self.WIDTHPROF = 120
598 589 self.HEIGHTPROF = 0
599 590 self.counter_imagwr = 0
600 591
601 592 self.PLOT_CODE = PARMS_CODE
602 593
603 594 self.FTP_WEI = None
604 595 self.EXP_CODE = None
605 596 self.SUB_EXP_CODE = None
606 597 self.PLOT_POS = None
607 598 self.tmin = None
608 599 self.tmax = None
609 600
610 601 self.xmin = None
611 602 self.xmax = None
612 603
613 604 self.figfile = None
614 605
615 606 def getSubplots(self):
616 607
617 608 ncol = 1
618 609 nrow = self.nplots
619 610
620 611 return nrow, ncol
621 612
622 613 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
623 614
624 615 self.__showprofile = showprofile
625 616 self.nplots = nplots
626 617
627 618 ncolspan = 1
628 619 colspan = 1
629 620
630 621 self.createFigure(id = id,
631 622 wintitle = wintitle,
632 623 widthplot = self.WIDTH + self.WIDTHPROF,
633 624 heightplot = self.HEIGHT + self.HEIGHTPROF,
634 625 show=show)
635 626
636 627 nrow, ncol = self.getSubplots()
637 628
638 629 counter = 0
639 630 for y in range(nrow):
640 631 for x in range(ncol):
641 632
642 633 if counter >= self.nplots:
643 634 break
644 635
645 636 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
646 637
647 638 if showprofile:
648 639 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
649 640
650 641 counter += 1
651 642
652 643 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
653 644 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
654 645 parameterIndex = None, onlyPositive = False,
655 646 SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, onlySNR = False,
656 647 DOP = True,
657 648 zlabel = "", parameterName = "", parameterObject = "data_param",
658 649 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
659 650 server=None, folder=None, username=None, password=None,
660 651 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
661 652
662 653 """
663 654
664 655 Input:
665 656 dataOut :
666 657 id :
667 658 wintitle :
668 659 channelList :
669 660 showProfile :
670 661 xmin : None,
671 662 xmax : None,
672 663 ymin : None,
673 664 ymax : None,
674 665 zmin : None,
675 666 zmax : None
676 667 """
677 668
678 669 data_param = getattr(dataOut, parameterObject)
679 670
680 671 if channelList == None:
681 672 channelIndexList = numpy.arange(data_param.shape[0])
682 673 else:
683 674 channelIndexList = numpy.array(channelList)
684 675
685 676 nchan = len(channelIndexList) #Number of channels being plotted
686 677
687 678 if nchan < 1:
688 679 return
689 680
690 681 nGraphsByChannel = 0
691 682
692 683 if SNR:
693 684 nGraphsByChannel += 1
694 685 if DOP:
695 686 nGraphsByChannel += 1
696 687
697 688 if nGraphsByChannel < 1:
698 689 return
699 690
700 691 nplots = nGraphsByChannel*nchan
701 692
702 693 if timerange is not None:
703 694 self.timerange = timerange
704 695
705 696 #tmin = None
706 697 #tmax = None
707 698 if parameterIndex == None:
708 699 parameterIndex = 1
709 700
710 x = dataOut.getTimeRange1()
701 x = dataOut.getTimeRange1(dataOut.paramInterval)
711 702 y = dataOut.heightList
712 703 z = data_param[channelIndexList,parameterIndex,:].copy()
713 704
714 705 zRange = dataOut.abscissaList
715 706 # nChannels = z.shape[0] #Number of wind dimensions estimated
716 707 # thisDatetime = dataOut.datatime
717 708
718 709 if dataOut.data_SNR is not None:
719 710 SNRarray = dataOut.data_SNR[channelIndexList,:]
720 711 SNRdB = 10*numpy.log10(SNRarray)
721 712 # SNRavgdB = 10*numpy.log10(SNRavg)
722 713 ind = numpy.where(SNRdB < 10**(SNRthresh/10))
723 714 z[ind] = numpy.nan
724 715
725 716 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
726 717 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
727 718 xlabel = ""
728 719 ylabel = "Range (Km)"
729 720
730 721 if (SNR and not onlySNR): nplots = 2*nplots
731 722
732 723 if onlyPositive:
733 724 colormap = "jet"
734 725 zmin = 0
735 726 else: colormap = "RdBu_r"
736 727
737 728 if not self.isConfig:
738 729
739 730 self.setup(id=id,
740 731 nplots=nplots,
741 732 wintitle=wintitle,
742 733 showprofile=showprofile,
743 734 show=show)
744 735
745 736 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
746 737
747 738 if ymin == None: ymin = numpy.nanmin(y)
748 739 if ymax == None: ymax = numpy.nanmax(y)
749 740 if zmin == None: zmin = numpy.nanmin(zRange)
750 741 if zmax == None: zmax = numpy.nanmax(zRange)
751 742
752 743 if SNR:
753 744 if SNRmin == None: SNRmin = numpy.nanmin(SNRdB)
754 745 if SNRmax == None: SNRmax = numpy.nanmax(SNRdB)
755 746
756 747 self.FTP_WEI = ftp_wei
757 748 self.EXP_CODE = exp_code
758 749 self.SUB_EXP_CODE = sub_exp_code
759 750 self.PLOT_POS = plot_pos
760 751
761 752 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
762 753 self.isConfig = True
763 754 self.figfile = figfile
764 755
765 756 self.setWinTitle(title)
766 757
767 758 if ((self.xmax - x[1]) < (x[1]-x[0])):
768 759 x[1] = self.xmax
769 760
770 761 for i in range(nchan):
771 762
772 763 if (SNR and not onlySNR): j = 2*i
773 764 else: j = i
774 765
775 766 j = nGraphsByChannel*i
776 767
777 768 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
778 769 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
779 770
780 771 if not onlySNR:
781 772 axes = self.axesList[j*self.__nsubplots]
782 773 z1 = z[i,:].reshape((1,-1))
783 774 axes.pcolorbuffer(x, y, z1,
784 775 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
785 776 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
786 777 ticksize=9, cblabel=zlabel, cbsize="1%")
787 778
788 779 if DOP:
789 780 title = "%s Channel %d: %s" %(parameterName, channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
790 781
791 782 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
792 783 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
793 784 axes = self.axesList[j]
794 785 z1 = z[i,:].reshape((1,-1))
795 786 axes.pcolorbuffer(x, y, z1,
796 787 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
797 788 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
798 789 ticksize=9, cblabel=zlabel, cbsize="1%")
799 790
800 791 if SNR:
801 792 title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
802 793 axes = self.axesList[(j)*self.__nsubplots]
803 794 if not onlySNR:
804 795 axes = self.axesList[(j + 1)*self.__nsubplots]
805 796
806 797 axes = self.axesList[(j + nGraphsByChannel-1)]
807 798
808 799 z1 = SNRdB[i,:].reshape((1,-1))
809 800 axes.pcolorbuffer(x, y, z1,
810 801 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
811 802 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap="jet",
812 803 ticksize=9, cblabel=zlabel, cbsize="1%")
813 804
814 805
815 806
816 807 self.draw()
817 808
818 809 if x[1] >= self.axesList[0].xmax:
819 810 self.counter_imagwr = wr_period
820 811 self.isConfig = False
821 812 self.figfile = None
822 813
823 814 self.save(figpath=figpath,
824 815 figfile=figfile,
825 816 save=save,
826 817 ftp=ftp,
827 818 wr_period=wr_period,
828 819 thisDatetime=thisDatetime,
829 820 update_figfile=False)
830 821
831 822 class SpectralFittingPlot(Figure):
832 823
833 824 __isConfig = None
834 825 __nsubplots = None
835 826
836 827 WIDTHPROF = None
837 828 HEIGHTPROF = None
838 829 PREFIX = 'prm'
839 830
840 831
841 832 N = None
842 833 ippSeconds = None
843 834
844 835 def __init__(self):
845 836 self.isConfig = False
846 837 self.__nsubplots = 1
847 838
848 839 self.PLOT_CODE = SPECFIT_CODE
849 840
850 841 self.WIDTH = 450
851 842 self.HEIGHT = 250
852 843 self.WIDTHPROF = 0
853 844 self.HEIGHTPROF = 0
854 845
855 846 def getSubplots(self):
856 847
857 848 ncol = int(numpy.sqrt(self.nplots)+0.9)
858 849 nrow = int(self.nplots*1./ncol + 0.9)
859 850
860 851 return nrow, ncol
861 852
862 853 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
863 854
864 855 showprofile = False
865 856 self.__showprofile = showprofile
866 857 self.nplots = nplots
867 858
868 859 ncolspan = 5
869 860 colspan = 4
870 861 if showprofile:
871 862 ncolspan = 5
872 863 colspan = 4
873 864 self.__nsubplots = 2
874 865
875 866 self.createFigure(id = id,
876 867 wintitle = wintitle,
877 868 widthplot = self.WIDTH + self.WIDTHPROF,
878 869 heightplot = self.HEIGHT + self.HEIGHTPROF,
879 870 show=show)
880 871
881 872 nrow, ncol = self.getSubplots()
882 873
883 874 counter = 0
884 875 for y in range(nrow):
885 876 for x in range(ncol):
886 877
887 878 if counter >= self.nplots:
888 879 break
889 880
890 881 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
891 882
892 883 if showprofile:
893 884 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
894 885
895 886 counter += 1
896 887
897 888 def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
898 889 xmin=None, xmax=None, ymin=None, ymax=None,
899 890 save=False, figpath='./', figfile=None, show=True):
900 891
901 892 """
902 893
903 894 Input:
904 895 dataOut :
905 896 id :
906 897 wintitle :
907 898 channelList :
908 899 showProfile :
909 900 xmin : None,
910 901 xmax : None,
911 902 zmin : None,
912 903 zmax : None
913 904 """
914 905
915 906 if cutHeight==None:
916 907 h=270
917 908 heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
918 909 cutHeight = dataOut.heightList[heightindex]
919 910
920 911 factor = dataOut.normFactor
921 912 x = dataOut.abscissaList[:-1]
922 913 #y = dataOut.getHeiRange()
923 914
924 915 z = dataOut.data_pre[:,:,heightindex]/factor
925 916 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
926 917 avg = numpy.average(z, axis=1)
927 918 listChannels = z.shape[0]
928 919
929 920 #Reconstruct Function
930 921 if fit==True:
931 922 groupArray = dataOut.groupList
932 923 listChannels = groupArray.reshape((groupArray.size))
933 924 listChannels.sort()
934 925 spcFitLine = numpy.zeros(z.shape)
935 926 constants = dataOut.constants
936 927
937 928 nGroups = groupArray.shape[0]
938 929 nChannels = groupArray.shape[1]
939 930 nProfiles = z.shape[1]
940 931
941 932 for f in range(nGroups):
942 933 groupChann = groupArray[f,:]
943 934 p = dataOut.data_param[f,:,heightindex]
944 935 # p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
945 936 fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
946 937 fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
947 938 spcFitLine[groupChann,:] = fitLineAux
948 939 # spcFitLine = spcFitLine/factor
949 940
950 941 z = z[listChannels,:]
951 942 spcFitLine = spcFitLine[listChannels,:]
952 943 spcFitLinedB = 10*numpy.log10(spcFitLine)
953 944
954 945 zdB = 10*numpy.log10(z)
955 946 #thisDatetime = dataOut.datatime
956 947 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
957 948 title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
958 949 xlabel = "Velocity (m/s)"
959 950 ylabel = "Spectrum"
960 951
961 952 if not self.isConfig:
962 953
963 954 nplots = listChannels.size
964 955
965 956 self.setup(id=id,
966 957 nplots=nplots,
967 958 wintitle=wintitle,
968 959 showprofile=showprofile,
969 960 show=show)
970 961
971 962 if xmin == None: xmin = numpy.nanmin(x)
972 963 if xmax == None: xmax = numpy.nanmax(x)
973 964 if ymin == None: ymin = numpy.nanmin(zdB)
974 965 if ymax == None: ymax = numpy.nanmax(zdB)+2
975 966
976 967 self.isConfig = True
977 968
978 969 self.setWinTitle(title)
979 970 for i in range(self.nplots):
980 971 # title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
981 972 title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i])
982 973 axes = self.axesList[i*self.__nsubplots]
983 974 if fit == False:
984 975 axes.pline(x, zdB[i,:],
985 976 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
986 977 xlabel=xlabel, ylabel=ylabel, title=title
987 978 )
988 979 if fit == True:
989 980 fitline=spcFitLinedB[i,:]
990 981 y=numpy.vstack([zdB[i,:],fitline] )
991 982 legendlabels=['Data','Fitting']
992 983 axes.pmultilineyaxis(x, y,
993 984 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
994 985 xlabel=xlabel, ylabel=ylabel, title=title,
995 986 legendlabels=legendlabels, marker=None,
996 987 linestyle='solid', grid='both')
997 988
998 989 self.draw()
999 990
1000 991 self.save(figpath=figpath,
1001 992 figfile=figfile,
1002 993 save=save,
1003 994 ftp=ftp,
1004 995 wr_period=wr_period,
1005 996 thisDatetime=thisDatetime)
1006 997
1007 998
1008 999 class EWDriftsPlot(Figure):
1009 1000
1010 1001 __isConfig = None
1011 1002 __nsubplots = None
1012 1003
1013 1004 WIDTHPROF = None
1014 1005 HEIGHTPROF = None
1015 1006 PREFIX = 'drift'
1016 1007
1017 1008 def __init__(self):
1018 1009
1019 1010 self.timerange = 2*60*60
1020 1011 self.isConfig = False
1021 1012 self.__nsubplots = 1
1022 1013
1023 1014 self.WIDTH = 800
1024 1015 self.HEIGHT = 150
1025 1016 self.WIDTHPROF = 120
1026 1017 self.HEIGHTPROF = 0
1027 1018 self.counter_imagwr = 0
1028 1019
1029 1020 self.PLOT_CODE = EWDRIFT_CODE
1030 1021
1031 1022 self.FTP_WEI = None
1032 1023 self.EXP_CODE = None
1033 1024 self.SUB_EXP_CODE = None
1034 1025 self.PLOT_POS = None
1035 1026 self.tmin = None
1036 1027 self.tmax = None
1037 1028
1038 1029 self.xmin = None
1039 1030 self.xmax = None
1040 1031
1041 1032 self.figfile = None
1042 1033
1043 1034 def getSubplots(self):
1044 1035
1045 1036 ncol = 1
1046 1037 nrow = self.nplots
1047 1038
1048 1039 return nrow, ncol
1049 1040
1050 1041 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1051 1042
1052 1043 self.__showprofile = showprofile
1053 1044 self.nplots = nplots
1054 1045
1055 1046 ncolspan = 1
1056 1047 colspan = 1
1057 1048
1058 1049 self.createFigure(id = id,
1059 1050 wintitle = wintitle,
1060 1051 widthplot = self.WIDTH + self.WIDTHPROF,
1061 1052 heightplot = self.HEIGHT + self.HEIGHTPROF,
1062 1053 show=show)
1063 1054
1064 1055 nrow, ncol = self.getSubplots()
1065 1056
1066 1057 counter = 0
1067 1058 for y in range(nrow):
1068 1059 if counter >= self.nplots:
1069 1060 break
1070 1061
1071 1062 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
1072 1063 counter += 1
1073 1064
1074 1065 def run(self, dataOut, id, wintitle="", channelList=None,
1075 1066 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
1076 1067 zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
1077 1068 timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
1078 1069 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1079 1070 server=None, folder=None, username=None, password=None,
1080 1071 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1081 1072 """
1082 1073
1083 1074 Input:
1084 1075 dataOut :
1085 1076 id :
1086 1077 wintitle :
1087 1078 channelList :
1088 1079 showProfile :
1089 1080 xmin : None,
1090 1081 xmax : None,
1091 1082 ymin : None,
1092 1083 ymax : None,
1093 1084 zmin : None,
1094 1085 zmax : None
1095 1086 """
1096 1087
1097 1088 if timerange is not None:
1098 1089 self.timerange = timerange
1099 1090
1100 1091 tmin = None
1101 1092 tmax = None
1102 1093
1103 x = dataOut.getTimeRange1()
1094 x = dataOut.getTimeRange1(dataOut.outputInterval)
1104 1095 # y = dataOut.heightList
1105 1096 y = dataOut.heightList
1106 1097
1107 1098 z = dataOut.data_output
1108 1099 nplots = z.shape[0] #Number of wind dimensions estimated
1109 1100 nplotsw = nplots
1110 1101
1111 1102 #If there is a SNR function defined
1112 1103 if dataOut.data_SNR is not None:
1113 1104 nplots += 1
1114 1105 SNR = dataOut.data_SNR
1115 1106
1116 1107 if SNR_1:
1117 1108 SNR += 1
1118 1109
1119 1110 SNRavg = numpy.average(SNR, axis=0)
1120 1111
1121 1112 SNRdB = 10*numpy.log10(SNR)
1122 1113 SNRavgdB = 10*numpy.log10(SNRavg)
1123 1114
1124 1115 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
1125 1116
1126 1117 for i in range(nplotsw):
1127 1118 z[i,ind] = numpy.nan
1128 1119
1129 1120
1130 1121 showprofile = False
1131 1122 # thisDatetime = dataOut.datatime
1132 1123 thisDatetime = datetime.datetime.utcfromtimestamp(x[1])
1133 1124 title = wintitle + " EW Drifts"
1134 1125 xlabel = ""
1135 1126 ylabel = "Height (Km)"
1136 1127
1137 1128 if not self.isConfig:
1138 1129
1139 1130 self.setup(id=id,
1140 1131 nplots=nplots,
1141 1132 wintitle=wintitle,
1142 1133 showprofile=showprofile,
1143 1134 show=show)
1144 1135
1145 1136 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1146 1137
1147 1138 if ymin == None: ymin = numpy.nanmin(y)
1148 1139 if ymax == None: ymax = numpy.nanmax(y)
1149 1140
1150 1141 if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
1151 1142 if zminZonal == None: zminZonal = -zmaxZonal
1152 1143 if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
1153 1144 if zminVertical == None: zminVertical = -zmaxVertical
1154 1145
1155 1146 if dataOut.data_SNR is not None:
1156 1147 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
1157 1148 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
1158 1149
1159 1150 self.FTP_WEI = ftp_wei
1160 1151 self.EXP_CODE = exp_code
1161 1152 self.SUB_EXP_CODE = sub_exp_code
1162 1153 self.PLOT_POS = plot_pos
1163 1154
1164 1155 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1165 1156 self.isConfig = True
1166 1157
1167 1158
1168 1159 self.setWinTitle(title)
1169 1160
1170 1161 if ((self.xmax - x[1]) < (x[1]-x[0])):
1171 1162 x[1] = self.xmax
1172 1163
1173 1164 strWind = ['Zonal','Vertical']
1174 1165 strCb = 'Velocity (m/s)'
1175 1166 zmaxVector = [zmaxZonal, zmaxVertical]
1176 1167 zminVector = [zminZonal, zminVertical]
1177 1168
1178 1169 for i in range(nplotsw):
1179 1170
1180 1171 title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1181 1172 axes = self.axesList[i*self.__nsubplots]
1182 1173
1183 1174 z1 = z[i,:].reshape((1,-1))
1184 1175
1185 1176 axes.pcolorbuffer(x, y, z1,
1186 1177 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
1187 1178 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1188 1179 ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
1189 1180
1190 1181 if dataOut.data_SNR is not None:
1191 1182 i += 1
1192 1183 if SNR_1:
1193 1184 title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1194 1185 else:
1195 1186 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1196 1187 axes = self.axesList[i*self.__nsubplots]
1197 1188 SNRavgdB = SNRavgdB.reshape((1,-1))
1198 1189
1199 1190 axes.pcolorbuffer(x, y, SNRavgdB,
1200 1191 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1201 1192 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1202 1193 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
1203 1194
1204 1195 self.draw()
1205 1196
1206 1197 if x[1] >= self.axesList[0].xmax:
1207 1198 self.counter_imagwr = wr_period
1208 1199 self.isConfig = False
1209 1200 self.figfile = None
1210 1201
1211 1202
1212 1203
1213 1204
1214 1205 class PhasePlot(Figure):
1215 1206
1216 1207 __isConfig = None
1217 1208 __nsubplots = None
1218 1209
1219 1210 PREFIX = 'mphase'
1220 1211
1221 1212 def __init__(self):
1222 1213
1223 1214 self.timerange = 24*60*60
1224 1215 self.isConfig = False
1225 1216 self.__nsubplots = 1
1226 1217 self.counter_imagwr = 0
1227 1218 self.WIDTH = 600
1228 1219 self.HEIGHT = 300
1229 1220 self.WIDTHPROF = 120
1230 1221 self.HEIGHTPROF = 0
1231 1222 self.xdata = None
1232 1223 self.ydata = None
1233 1224
1234 1225 self.PLOT_CODE = MPHASE_CODE
1235 1226
1236 1227 self.FTP_WEI = None
1237 1228 self.EXP_CODE = None
1238 1229 self.SUB_EXP_CODE = None
1239 1230 self.PLOT_POS = None
1240 1231
1241 1232
1242 1233 self.filename_phase = None
1243 1234
1244 1235 self.figfile = None
1245 1236
1246 1237 def getSubplots(self):
1247 1238
1248 1239 ncol = 1
1249 1240 nrow = 1
1250 1241
1251 1242 return nrow, ncol
1252 1243
1253 1244 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1254 1245
1255 1246 self.__showprofile = showprofile
1256 1247 self.nplots = nplots
1257 1248
1258 1249 ncolspan = 7
1259 1250 colspan = 6
1260 1251 self.__nsubplots = 2
1261 1252
1262 1253 self.createFigure(id = id,
1263 1254 wintitle = wintitle,
1264 1255 widthplot = self.WIDTH+self.WIDTHPROF,
1265 1256 heightplot = self.HEIGHT+self.HEIGHTPROF,
1266 1257 show=show)
1267 1258
1268 1259 nrow, ncol = self.getSubplots()
1269 1260
1270 1261 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1271 1262
1272 1263
1273 1264 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1274 1265 xmin=None, xmax=None, ymin=None, ymax=None,
1275 1266 timerange=None,
1276 1267 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1277 1268 server=None, folder=None, username=None, password=None,
1278 1269 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1279 1270
1280 1271
1281 1272 tmin = None
1282 1273 tmax = None
1283 x = dataOut.getTimeRange1()
1274 x = dataOut.getTimeRange1(dataOut.outputInterval)
1284 1275 y = dataOut.getHeiRange()
1285 1276
1286 1277
1287 1278 #thisDatetime = dataOut.datatime
1288 1279 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1289 1280 title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1290 1281 xlabel = "Local Time"
1291 1282 ylabel = "Phase"
1292 1283
1293 1284
1294 1285 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1295 1286 phase_beacon = dataOut.data_output
1296 1287 update_figfile = False
1297 1288
1298 1289 if not self.isConfig:
1299 1290
1300 1291 self.nplots = phase_beacon.size
1301 1292
1302 1293 self.setup(id=id,
1303 1294 nplots=self.nplots,
1304 1295 wintitle=wintitle,
1305 1296 showprofile=showprofile,
1306 1297 show=show)
1307 1298
1308 1299 if timerange is not None:
1309 1300 self.timerange = timerange
1310 1301
1311 1302 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1312 1303
1313 1304 if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0
1314 1305 if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0
1315 1306
1316 1307 self.FTP_WEI = ftp_wei
1317 1308 self.EXP_CODE = exp_code
1318 1309 self.SUB_EXP_CODE = sub_exp_code
1319 1310 self.PLOT_POS = plot_pos
1320 1311
1321 1312 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1322 1313 self.isConfig = True
1323 1314 self.figfile = figfile
1324 1315 self.xdata = numpy.array([])
1325 1316 self.ydata = numpy.array([])
1326 1317
1327 1318 #open file beacon phase
1328 1319 path = '%s%03d' %(self.PREFIX, self.id)
1329 1320 beacon_file = os.path.join(path,'%s.txt'%self.name)
1330 1321 self.filename_phase = os.path.join(figpath,beacon_file)
1331 1322 update_figfile = True
1332 1323
1333 1324
1334 1325 #store data beacon phase
1335 1326 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1336 1327
1337 1328 self.setWinTitle(title)
1338 1329
1339 1330
1340 1331 title = "Phase Offset %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1341 1332
1342 1333 legendlabels = ["phase %d"%(chan) for chan in numpy.arange(self.nplots)]
1343 1334
1344 1335 axes = self.axesList[0]
1345 1336
1346 1337 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1347 1338
1348 1339 if len(self.ydata)==0:
1349 1340 self.ydata = phase_beacon.reshape(-1,1)
1350 1341 else:
1351 1342 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1352 1343
1353 1344
1354 1345 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1355 1346 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1356 1347 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1357 1348 XAxisAsTime=True, grid='both'
1358 1349 )
1359 1350
1360 1351 self.draw()
1361 1352
1362 1353 if dataOut.ltctime >= self.xmax:
1363 1354 self.counter_imagwr = wr_period
1364 1355 self.isConfig = False
1365 1356 update_figfile = True
1366 1357
1367 1358 self.save(figpath=figpath,
1368 1359 figfile=figfile,
1369 1360 save=save,
1370 1361 ftp=ftp,
1371 1362 wr_period=wr_period,
1372 1363 thisDatetime=thisDatetime,
1373 1364 update_figfile=update_figfile)
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (728 lines changed) Show them Hide them
@@ -1,1026 +1,1054
1 1 import numpy
2 2 import time
3 3 import os
4 4 import h5py
5 5 import re
6 import datetime
6 7
7 8 from schainpy.model.data.jrodata import *
8 9 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
10 # from jroIO_base import *
9 11 from schainpy.model.io.jroIO_base import *
12 import schainpy
10 13
11 14
12 15 class HDF5Reader(ProcessingUnit):
16 '''
17 Reads HDF5 format files
18
19 path
20
21 startDate
22
23 endDate
24
25 startTime
26
27 endTime
28 '''
13 29
14 30 ext = ".hdf5"
15 31
16 32 optchar = "D"
17 33
18 34 timezone = None
19 35
20 secStart = None
36 startTime = None
21 37
22 secEnd = None
38 endTime = None
23 39
24 40 fileIndex = None
25 41
26 blockIndex = None
42 utcList = None #To select data in the utctime list
27 43
28 blocksPerFile = None
44 blockList = None #List to blocks to be read from the file
45
46 blocksPerFile = None #Number of blocks to be read
47
48 blockIndex = None
29 49
30 50 path = None
31 51
32 52 #List of Files
33 53
34 54 filenameList = None
35 55
36 56 datetimeList = None
37 57
38 58 #Hdf5 File
39 59
40 fpMetadata = None
41
42 pathMeta = None
43
44 60 listMetaname = None
45 61
46 62 listMeta = None
47 63
48 64 listDataname = None
49 65
50 66 listData = None
51 67
52 68 listShapes = None
53 69
54 70 fp = None
55 71
56 72 #dataOut reconstruction
57 73
58 74 dataOut = None
59 75
60 nRecords = None
61
62 76
63 77 def __init__(self):
64 self.dataOut = self.__createObjByDefault()
78 self.dataOut = Parameters()
65 79 return
66
67 def __createObjByDefault(self):
68
69 dataObj = Parameters()
70
71 return dataObj
72
73 def setup(self,path=None,
74 startDate=None,
75 endDate=None,
76 startTime=datetime.time(0,0,0),
77 endTime=datetime.time(23,59,59),
78 walk=True,
79 timezone='ut',
80 all=0,
81 online=False,
82 ext=None):
83
84 if ext==None:
85 ext = self.ext
86 self.timezone = timezone
87 # self.all = all
88 # self.online = online
89 self.path = path
90
91 startDateTime = datetime.datetime.combine(startDate,startTime)
92 endDateTime = datetime.datetime.combine(endDate,endTime)
93 secStart = (startDateTime-datetime.datetime(1970,1,1)).total_seconds()
94 secEnd = (endDateTime-datetime.datetime(1970,1,1)).total_seconds()
95
96 self.secStart = secStart
97 self.secEnd = secEnd
98
99 if not(online):
100 #Busqueda de archivos offline
101 self.__searchFilesOffline(path, startDate, endDate, ext, startTime, endTime, secStart, secEnd, walk)
102 else:
103 self.__searchFilesOnline(path, walk)
104 80
105 if not(self.filenameList):
81 def setup(self, **kwargs):
82
83 path = kwargs['path']
84 startDate = kwargs['startDate']
85 endDate = kwargs['endDate']
86 startTime = kwargs['startTime']
87 endTime = kwargs['endTime']
88 walk = kwargs['walk']
89 if kwargs.has_key('ext'):
90 ext = kwargs['ext']
91 else:
92 ext = '.hdf5'
93
94 print "[Reading] Searching files in offline mode ..."
95 pathList, filenameList = self.__searchFilesOffLine(path, startDate=startDate, endDate=endDate,
96 startTime=startTime, endTime=endTime,
97 ext=ext, walk=walk)
98
99 if not(filenameList):
106 100 print "There is no files into the folder: %s"%(path)
107 101 sys.exit(-1)
108 102
109 # self.__getExpParameters()
110
111 103 self.fileIndex = -1
112
113 self.__setNextFileOffline()
104 self.startTime = startTime
105 self.endTime = endTime
114 106
115 107 self.__readMetadata()
116 108
117 self.blockIndex = 0
109 self.__setNextFileOffline()
118 110
119 111 return
120
121 def __searchFilesOffline(self,
112
113 def __searchFilesOffLine(self,
122 114 path,
123 startDate,
124 endDate,
125 ext,
115 startDate=None,
116 endDate=None,
126 117 startTime=datetime.time(0,0,0),
127 118 endTime=datetime.time(23,59,59),
128 secStart = 0,
129 secEnd = numpy.inf,
119 ext='.hdf5',
130 120 walk=True):
131 121
132 # self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
133 #
134 # self.__checkPath()
135 #
136 # self.__findDataForDates()
137 #
138 # self.__selectDataForTimes()
139 #
140 # for i in range(len(self.filenameList)):
141 # print "%s" %(self.filenameList[i])
122 expLabel = ''
123 self.filenameList = []
124 self.datetimeList = []
142 125
143 126 pathList = []
144 127
145 if not walk:
146 #pathList.append(path)
147 multi_path = path.split(',')
148 for single_path in multi_path:
149 pathList.append(single_path)
150
151 else:
152 #dirList = []
153 multi_path = path.split(',')
154 for single_path in multi_path:
155 dirList = []
156 for thisPath in os.listdir(single_path):
157 if not os.path.isdir(os.path.join(single_path,thisPath)):
158 continue
159 if not isDoyFolder(thisPath):
160 continue
128 JRODataObj = JRODataReader()
129 dateList, pathList = JRODataObj.findDatafiles(path, startDate, endDate, expLabel, ext, walk, include_path=True)
130
131 if dateList == []:
132 print "[Reading] No *%s files in %s from %s to %s)"%(ext, path,
133 datetime.datetime.combine(startDate,startTime).ctime(),
134 datetime.datetime.combine(endDate,endTime).ctime())
161 135
162 dirList.append(thisPath)
163
164 if not(dirList):
165 return None, None
166
167 thisDate = startDate
168
169 while(thisDate <= endDate):
170 year = thisDate.timetuple().tm_year
171 doy = thisDate.timetuple().tm_yday
172
173 matchlist = fnmatch.filter(dirList, '?' + '%4.4d%3.3d' % (year,doy) + '*')
174 if len(matchlist) == 0:
175 thisDate += datetime.timedelta(1)
176 continue
177 for match in matchlist:
178 pathList.append(os.path.join(single_path,match))
179
180 thisDate += datetime.timedelta(1)
181
182 if pathList == []:
183 print "Any folder was found for the date range: %s-%s" %(startDate, endDate)
184 136 return None, None
185 137
186 print "%d folder(s) was(were) found for the date range: %s - %s" %(len(pathList), startDate, endDate)
187
138 if len(dateList) > 1:
139 print "[Reading] %d days were found in date range: %s - %s" %(len(dateList), startDate, endDate)
140 else:
141 print "[Reading] data was found for the date %s" %(dateList[0])
142
188 143 filenameList = []
189 144 datetimeList = []
190 pathDict = {}
191 filenameList_to_sort = []
192
193 for i in range(len(pathList)):
194
195 thisPath = pathList[i]
196
197 fileList = glob.glob1(thisPath, "*%s" %ext)
198 fileList.sort()
199 pathDict.setdefault(fileList[0])
200 pathDict[fileList[0]] = i
201 filenameList_to_sort.append(fileList[0])
202 145
203 filenameList_to_sort.sort()
146 #----------------------------------------------------------------------------------
204 147
205 for file in filenameList_to_sort:
206 thisPath = pathList[pathDict[file]]
148 for thisPath in pathList:
149 # thisPath = pathList[pathDict[file]]
207 150
208 151 fileList = glob.glob1(thisPath, "*%s" %ext)
209 152 fileList.sort()
210 153
211 154 for file in fileList:
212 155
213 156 filename = os.path.join(thisPath,file)
214 thisDatetime = self.__isFileinThisTime(filename, secStart, secEnd)
157
158 if not isFileInDateRange(filename, startDate, endDate):
159 continue
160
161 thisDatetime = self.__isFileInTimeRange(filename, startDate, endDate, startTime, endTime)
215 162
216 163 if not(thisDatetime):
217 164 continue
218 165
219 166 filenameList.append(filename)
220 167 datetimeList.append(thisDatetime)
221 168
222 169 if not(filenameList):
223 print "Any file was found for the time range %s - %s" %(startTime, endTime)
170 print "[Reading] Any file was found int time range %s - %s" %(datetime.datetime.combine(startDate,startTime).ctime(), datetime.datetime.combine(endDate,endTime).ctime())
224 171 return None, None
225 172
226 print "%d file(s) was(were) found for the time range: %s - %s" %(len(filenameList), startTime, endTime)
173 print "[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime)
227 174 print
228 175
229 176 for i in range(len(filenameList)):
230 print "%s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
177 print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
231 178
232 179 self.filenameList = filenameList
233 180 self.datetimeList = datetimeList
234
181
235 182 return pathList, filenameList
236
237 def __isFileinThisTime(self, filename, startSeconds, endSeconds):
183
184 def __isFileInTimeRange(self,filename, startDate, endDate, startTime, endTime):
185
238 186 """
239 187 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
240 188
241 189 Inputs:
242 190 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
243 191
192 startDate : fecha inicial del rango seleccionado en formato datetime.date
193
194 endDate : fecha final del rango seleccionado en formato datetime.date
195
244 196 startTime : tiempo inicial del rango seleccionado en formato datetime.time
245 197
246 198 endTime : tiempo final del rango seleccionado en formato datetime.time
247 199
248 200 Return:
249 201 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
250 202 fecha especificado, de lo contrario retorna False.
251 203
252 204 Excepciones:
253 205 Si el archivo no existe o no puede ser abierto
254 206 Si la cabecera no puede ser leida.
255 207
256 208 """
257
209
258 210 try:
259 fp = fp = h5py.File(filename,'r')
211 fp = h5py.File(filename,'r')
212 grp1 = fp['Data']
213
260 214 except IOError:
261 215 traceback.print_exc()
262 216 raise IOError, "The file %s can't be opened" %(filename)
217 #chino rata
218 #In case has utctime attribute
219 grp2 = grp1['utctime']
220 # thisUtcTime = grp2.value[0] - 5*3600 #To convert to local time
221 thisUtcTime = grp2.value[0]
222
223 fp.close()
263 224
264 grp = fp['Data']
265 timeAux = grp['time']
266 time0 = timeAux[:][0].astype(numpy.float) #Time Vector
225 thisDatetime = datetime.datetime.fromtimestamp(thisUtcTime[0])
226 thisDate = thisDatetime.date()
227 thisTime = thisDatetime.time()
267 228
268 fp.close()
229 startUtcTime = (datetime.datetime.combine(thisDate,startTime)- datetime.datetime(1970, 1, 1)).total_seconds()
230 endUtcTime = (datetime.datetime.combine(thisDate,endTime)- datetime.datetime(1970, 1, 1)).total_seconds()
269 231
270 if self.timezone == 'lt':
271 time0 -= 5*3600
272
273 boolTimer = numpy.logical_and(time0 >= startSeconds,time0 < endSeconds)
274
275 if not (numpy.any(boolTimer)):
232 #General case
233 # o>>>>>>>>>>>>>><<<<<<<<<<<<<<o
234 #-----------o----------------------------o-----------
235 # startTime endTime
236
237 if endTime >= startTime:
238 thisUtcLog = numpy.logical_and(thisUtcTime > startUtcTime, thisUtcTime < endUtcTime)
239 if numpy.any(thisUtcLog): #If there is one block between the hours mentioned
240 return thisDatetime
241 return None
242
243 #If endTime < startTime then endTime belongs to the next day
244 #<<<<<<<<<<<o o>>>>>>>>>>>
245 #-----------o----------------------------o-----------
246 # endTime startTime
247
248 if (thisDate == startDate) and numpy.all(thisUtcTime < startUtcTime):
249 return None
250
251 if (thisDate == endDate) and numpy.all(thisUtcTime > endUtcTime):
252 return None
253
254 if numpy.all(thisUtcTime < startUtcTime) and numpy.all(thisUtcTime > endUtcTime):
276 255 return None
277 256
278 thisDatetime = datetime.datetime.utcfromtimestamp(time0[0])
279 257 return thisDatetime
280
281 def __checkPath(self):
282 if os.path.exists(self.path):
283 self.status = 1
284 else:
285 self.status = 0
286 print 'Path:%s does not exists'%self.path
287
288 return
289
258
290 259 def __setNextFileOffline(self):
291 idFile = self.fileIndex
292 idFile += 1
293 260
261 self.fileIndex += 1
262 idFile = self.fileIndex
263
294 264 if not(idFile < len(self.filenameList)):
295 265 print "No more Files"
296 266 return 0
297 267
298 268 filename = self.filenameList[idFile]
299 269
300 270 filePointer = h5py.File(filename,'r')
301
302 self.flagIsNewFile = 1
303 self.fileIndex = idFile
271
304 272 self.filename = filename
305 273
306 274 self.fp = filePointer
307 275
308 276 print "Setting the file: %s"%self.filename
309 277
310 self.__readMetadata()
278 # self.__readMetadata()
311 279 self.__setBlockList()
280 self.__readData()
312 281 # self.nRecords = self.fp['Data'].attrs['blocksPerFile']
313 self.nRecords = self.fp['Data'].attrs['nRecords']
282 # self.nRecords = self.fp['Data'].attrs['nRecords']
314 283 self.blockIndex = 0
315 284 return 1
316 285
317 286 def __setBlockList(self):
318 287 '''
288 Selects the data within the times defined
289
319 290 self.fp
320 self.startDateTime
321 self.endDateTime
291 self.startTime
292 self.endTime
322 293
323 294 self.blockList
324 295 self.blocksPerFile
325 296
326 297 '''
327 filePointer = self.fp
328 secStart = self.secStart
329 secEnd = self.secEnd
298 fp = self.fp
299 startTime = self.startTime
300 endTime = self.endTime
330 301
331 grp = filePointer['Data']
332 timeVector = grp['time'].value.astype(numpy.float)[0]
302 grp = fp['Data']
303 thisUtcTime = grp['utctime'].value.astype(numpy.float)[0]
333 304
334 305 if self.timezone == 'lt':
335 timeVector -= 5*3600
306 thisUtcTime -= 5*3600
307
308 thisDatetime = datetime.datetime.fromtimestamp(thisUtcTime[0])
309 thisDate = thisDatetime.date()
310 thisTime = thisDatetime.time()
336 311
337 ind = numpy.where(numpy.logical_and(timeVector >= secStart , timeVector < secEnd))[0]
312 startUtcTime = (datetime.datetime.combine(thisDate,startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
313 endUtcTime = (datetime.datetime.combine(thisDate,endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
314
315 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
338 316
339 317 self.blockList = ind
340 318 self.blocksPerFile = len(ind)
341 319
342 320 return
343 321
344 322 def __readMetadata(self):
345 323 '''
324 Reads Metadata
325
346 326 self.pathMeta
347 327
348 328 self.listShapes
349 329 self.listMetaname
350 330 self.listMeta
351 331
352 332 '''
353 333
354 grp = self.fp['Data']
355 pathMeta = os.path.join(self.path, grp.attrs['metadata'])
356
357 if pathMeta == self.pathMeta:
358 return
359 else:
360 self.pathMeta = pathMeta
361
362 filePointer = h5py.File(self.pathMeta,'r')
363 groupPointer = filePointer['Metadata']
334 # grp = self.fp['Data']
335 # pathMeta = os.path.join(self.path, grp.attrs['metadata'])
336 #
337 # if pathMeta == self.pathMeta:
338 # return
339 # else:
340 # self.pathMeta = pathMeta
341 #
342 # filePointer = h5py.File(self.pathMeta,'r')
343 # groupPointer = filePointer['Metadata']
344
345 filename = self.filenameList[0]
346
347 fp = h5py.File(filename,'r')
348
349 gp = fp['Metadata']
364 350
365 351 listMetaname = []
366 352 listMetadata = []
367 for item in groupPointer.items():
353 for item in gp.items():
368 354 name = item[0]
369 355
370 356 if name=='array dimensions':
371 table = groupPointer[name][:]
357 table = gp[name][:]
372 358 listShapes = {}
373 359 for shapes in table:
374 listShapes[shapes[0]] = numpy.array([shapes[1],shapes[2],shapes[3],shapes[4]])
360 listShapes[shapes[0]] = numpy.array([shapes[1],shapes[2],shapes[3],shapes[4],shapes[5]])
375 361 else:
376 data = groupPointer[name].value
362 data = gp[name].value
377 363 listMetaname.append(name)
378 364 listMetadata.append(data)
379 365
380 if name=='type':
381 self.__initDataOut(data)
382
383 filePointer.close()
384
366 # if name=='type':
367 # self.__initDataOut(data)
368
385 369 self.listShapes = listShapes
386 370 self.listMetaname = listMetaname
387 371 self.listMeta = listMetadata
388 372
373 fp.close()
389 374 return
390 375
391 376 def __readData(self):
392 377 grp = self.fp['Data']
393 378 listdataname = []
394 379 listdata = []
395 380
396 381 for item in grp.items():
397 382 name = item[0]
398
399 if name == 'time':
400 listdataname.append('utctime')
401 timeAux = grp[name].value.astype(numpy.float)[0]
402 listdata.append(timeAux)
403 continue
404
405 383 listdataname.append(name)
406 array = self.__setDataArray(self.nRecords, grp[name],self.listShapes[name])
384
385 array = self.__setDataArray(grp[name],self.listShapes[name])
407 386 listdata.append(array)
408 387
409 388 self.listDataname = listdataname
410 389 self.listData = listdata
411 390 return
412 391
413 def __setDataArray(self, nRecords, dataset, shapes):
392 def __setDataArray(self, dataset, shapes):
393
394 nDims = shapes[0]
414 395
415 nChannels = shapes[0] #Dimension 0
396 nDim2 = shapes[1] #Dimension 0
416 397
417 nPoints = shapes[1] #Dimension 1, number of Points or Parameters
398 nDim1 = shapes[2] #Dimension 1, number of Points or Parameters
418 399
419 nSamples = shapes[2] #Dimension 2, number of samples or ranges
400 nDim0 = shapes[3] #Dimension 2, number of samples or ranges
420 401
421 mode = shapes[3]
402 mode = shapes[4] #Mode of storing
422 403
423 # if nPoints>1:
424 # arrayData = numpy.zeros((nRecords,nChannels,nPoints,nSamples))
425 # else:
426 # arrayData = numpy.zeros((nRecords,nChannels,nSamples))
427 #
428 # chn = 'channel'
429 #
430 # for i in range(nChannels):
431 #
432 # data = dataset[chn + str(i)].value
433 #
434 # if nPoints>1:
435 # data = numpy.rollaxis(data,2)
436 #
437 # arrayData[:,i,:] = data
438
439 arrayData = numpy.zeros((nRecords,nChannels,nPoints,nSamples))
440 doSqueeze = False
441 if mode == 0:
442 strds = 'channel'
443 nDatas = nChannels
444 newShapes = (nRecords,nPoints,nSamples)
445 if nPoints == 1:
446 doSqueeze = True
447 axisSqueeze = 2
448 else:
404 blockList = self.blockList
405
406 blocksPerFile = self.blocksPerFile
407
408 #Depending on what mode the data was stored
409 # if mode == 0: #Divided in channels
410 # strds = 'channel'
411 # nDatas = nDim2
412 # newShapes = (blocksPerFile,nDim1,nDim0)
413 if mode == 1: #Divided in parameter
449 414 strds = 'param'
450 nDatas = nPoints
451 newShapes = (nRecords,nChannels,nSamples)
452 if nChannels == 1:
453 doSqueeze = True
454 axisSqueeze = 1
455
456 for i in range(nDatas):
415 nDatas = nDim1
416 newShapes = (blocksPerFile,nDim2,nDim0)
417 elif mode==2: #Concatenated in a table
418 strds = 'table0'
419 arrayData = dataset[strds].value
420 #Selecting part of the dataset
421 utctime = arrayData[:,0]
422 u, indices = numpy.unique(utctime, return_index=True)
457 423
458 data = dataset[strds + str(i)].value
459 data = data.reshape(newShapes)
460
461 if mode == 0:
462 arrayData[:,i,:,:] = data
463 else:
424 if blockList.size != indices.size:
425 indMin = indices[blockList[0]]
426 indMax = indices[blockList[-1] + 1]
427 arrayData = arrayData[indMin:indMax,:]
428 return arrayData
429
430 #------- One dimension ---------------
431 if nDims == 1:
432 arrayData = dataset.value.astype(numpy.float)[0][blockList]
433
434 #------- Two dimensions -----------
435 elif nDims == 2:
436 arrayData = numpy.zeros((blocksPerFile,nDim1,nDim0))
437 newShapes = (blocksPerFile,nDim0)
438 nDatas = nDim1
439
440 for i in range(nDatas):
441 data = dataset[strds + str(i)].value
442 arrayData[:,i,:] = data[blockList,:]
443
444 #------- Three dimensions ---------
445 else:
446 arrayData = numpy.zeros((blocksPerFile,nDim2,nDim1,nDim0))
447 for i in range(nDatas):
448
449 data = dataset[strds + str(i)].value
450 data = data[blockList,:,:]
451 data = data.reshape(newShapes)
452 # if mode == 0:
453 # arrayData[:,i,:,:] = data
454 # else:
464 455 arrayData[:,:,i,:] = data
465
466 if doSqueeze:
467 arrayData = numpy.squeeze(arrayData, axis=axisSqueeze)
468
456
469 457 return arrayData
470
471 def __initDataOut(self, type):
472 458
473 # if type =='Parameters':
474 # self.dataOut = Parameters()
475 # elif type =='Spectra':
476 # self.dataOut = Spectra()
477 # elif type =='Voltage':
478 # self.dataOut = Voltage()
479 # elif type =='Correlation':
480 # self.dataOut = Correlation()
481
482 return
483
484 459 def __setDataOut(self):
485 460 listMeta = self.listMeta
486 461 listMetaname = self.listMetaname
487 462 listDataname = self.listDataname
488 463 listData = self.listData
464 listShapes = self.listShapes
489 465
490 466 blockIndex = self.blockIndex
491 blockList = self.blockList
467 # blockList = self.blockList
492 468
493 469 for i in range(len(listMeta)):
494 470 setattr(self.dataOut,listMetaname[i],listMeta[i])
495 471
496 472 for j in range(len(listData)):
497 if listDataname[j]=='utctime':
498 # setattr(self.dataOut,listDataname[j],listData[j][blockList[blockIndex]])
499 setattr(self.dataOut,'utctimeInit',listData[j][blockList[blockIndex]])
500 continue
501
502 setattr(self.dataOut,listDataname[j],listData[j][blockList[blockIndex],:])
473 nShapes = listShapes[listDataname[j]][0]
474 mode = listShapes[listDataname[j]][4]
475 if nShapes == 1:
476 setattr(self.dataOut,listDataname[j],listData[j][blockIndex])
477 elif nShapes > 1:
478 setattr(self.dataOut,listDataname[j],listData[j][blockIndex,:])
479 #Mode Meteors
480 elif mode ==2:
481 selectedData = self.__selectDataMode2(listData[j], blockIndex)
482 setattr(self.dataOut, listDataname[j], selectedData)
483 return
484
485 def __selectDataMode2(self, data, blockIndex):
486 utctime = data[:,0]
487 aux, indices = numpy.unique(utctime, return_inverse=True)
488 selInd = numpy.where(indices == blockIndex)[0]
489 selData = data[selInd,:]
503 490
504 return self.dataOut.data_param
491 return selData
505 492
506 493 def getData(self):
507 494
508 495 # if self.flagNoMoreFiles:
509 496 # self.dataOut.flagNoData = True
510 497 # print 'Process finished'
511 498 # return 0
512 499 #
513 500 if self.blockIndex==self.blocksPerFile:
514 501 if not( self.__setNextFileOffline() ):
515 502 self.dataOut.flagNoData = True
516 503 return 0
517
518 #
504
519 505 # if self.datablock == None: # setear esta condicion cuando no hayan datos por leers
520 506 # self.dataOut.flagNoData = True
521 507 # return 0
522
523 self.__readData()
508 # self.__readData()
524 509 self.__setDataOut()
525 510 self.dataOut.flagNoData = False
526 511
527 512 self.blockIndex += 1
528 513
529 514 return
530 515
531 516 def run(self, **kwargs):
532 517
533 518 if not(self.isConfig):
534 519 self.setup(**kwargs)
535 520 # self.setObjProperties()
536 521 self.isConfig = True
537 522
538 523 self.getData()
539 524
540 525 return
541 526
542 527 class HDF5Writer(Operation):
528 '''
529 HDF5 Writer, stores parameters data in HDF5 format files
530
531 path: path where the files will be stored
532
533 blocksPerFile: number of blocks that will be saved in per HDF5 format file
534
535 mode: selects the data stacking mode: '0' channels, '1' parameters, '3' table (for meteors)
536
537 metadataList: list of attributes that will be stored as metadata
538
539 dataList: list of attributes that will be stores as data
540
541 '''
542
543 543
544 544 ext = ".hdf5"
545 545
546 546 optchar = "D"
547 547
548 548 metaoptchar = "M"
549 549
550 550 metaFile = None
551 551
552 552 filename = None
553 553
554 554 path = None
555 555
556 556 setFile = None
557 557
558 558 fp = None
559 559
560 560 grp = None
561 561
562 562 ds = None
563 563
564 564 firsttime = True
565 565
566 566 #Configurations
567 567
568 568 blocksPerFile = None
569 569
570 570 blockIndex = None
571 571
572 572 dataOut = None
573 573
574 574 #Data Arrays
575 575
576 576 dataList = None
577 577
578 578 metadataList = None
579 579
580 580 arrayDim = None
581 581
582 582 tableDim = None
583 583
584 584 # dtype = [('arrayName', 'S20'),('nChannels', 'i'), ('nPoints', 'i'), ('nSamples', 'i'),('mode', 'b')]
585 585
586 586 dtype = [('arrayName', 'S20'),('nDimensions', 'i'), ('dim2', 'i'), ('dim1', 'i'),('dim0', 'i'),('mode', 'b')]
587 587
588 588 mode = None
589 589
590 590 nDatas = None #Number of datasets to be stored per array
591 591
592 592 nDims = None #Number Dimensions in each dataset
593 593
594 594 nDimsForDs = None
595 595
596 596 currentDay = None
597 597
598 598 def __init__(self):
599 599
600 600 Operation.__init__(self)
601 601 self.isConfig = False
602 602 return
603 603
604 604 def setup(self, dataOut, **kwargs):
605 605
606 606 self.path = kwargs['path']
607 607
608 if kwargs.has_key('ext'):
609 self.ext = kwargs['ext']
610
611 608 if kwargs.has_key('blocksPerFile'):
612 609 self.blocksPerFile = kwargs['blocksPerFile']
613 610 else:
614 611 self.blocksPerFile = 10
615 612
616 613 self.metadataList = kwargs['metadataList']
617 614
618 615 self.dataList = kwargs['dataList']
619 616
620 617 self.dataOut = dataOut
621 618
622 619 if kwargs.has_key('mode'):
623 620 mode = kwargs['mode']
624 621
625 622 if type(mode) == int:
626 623 mode = numpy.zeros(len(self.dataList)) + mode
627 624 else:
628 mode = numpy.zeros(len(self.dataList))
625 mode = numpy.ones(len(self.dataList))
629 626
630 627 self.mode = mode
631 628
632 629 arrayDim = numpy.zeros((len(self.dataList),5))
633 630
634 631 #Table dimensions
635 632
636 633 dtype0 = self.dtype
637 634
638 635 tableList = []
639 636
640 637 for i in range(len(self.dataList)):
641 638
642 639 dataAux = getattr(self.dataOut, self.dataList[i])
643 640
644 if type(dataAux)==float or type(dataAux)==int:
641 #--------------------- Conditionals ------------------------
642 #There is no data
643 if dataAux == None:
644 return 0
645
646 #Not array, just a number
647 if type(dataAux)==float or type(dataAux)==int:
645 648 arrayDim[i,0] = 1
649 mode[i] = 0
650
651 #Mode meteors
652 elif mode[i] == 2:
653 arrayDim[i,3] = dataAux.shape[-1]
654 arrayDim[i,4] = mode[i] #Mode the data was stored
655
656 #All the rest
646 657 else:
658 arrayDim0 = dataAux.shape #Data dimensions
659 arrayDim[i,0] = len(arrayDim0) #Number of array dimensions
660 arrayDim[i,4] = mode[i] #Mode the data was stored
647 661
648 if dataAux == None:
649 return 0
650
651 arrayDim0 = dataAux.shape
652 arrayDim[i,0] = len(arrayDim0)
653 arrayDim[i,4] = mode[i]
654
662 # Three-dimension arrays
655 663 if len(arrayDim0) == 3:
656 664 arrayDim[i,1:-1] = numpy.array(arrayDim0)
665
666 # Two-dimension arrays
657 667 elif len(arrayDim0) == 2:
658 arrayDim[i,2:-1] = numpy.array(arrayDim0) #nHeights
668 arrayDim[i,2:-1] = numpy.array(arrayDim0)
669
670 # One-dimension arrays
659 671 elif len(arrayDim0) == 1:
660 672 arrayDim[i,3] = arrayDim0
673
674 # No array, just a number
661 675 elif len(arrayDim0) == 0:
662 676 arrayDim[i,0] = 1
663 677 arrayDim[i,3] = 1
664 678
665 679 table = numpy.array((self.dataList[i],) + tuple(arrayDim[i,:]),dtype = dtype0)
666 680 tableList.append(table)
667 681
668 682 self.arrayDim = arrayDim
669 683 self.tableDim = numpy.array(tableList, dtype = dtype0)
670 684 self.blockIndex = 0
671 685
672 686 timeTuple = time.localtime(dataOut.utctime)
673 687 self.currentDay = timeTuple.tm_yday
674 688 return 1
675 689
676 690 def putMetadata(self):
677 691
678 692 fp = self.createMetadataFile()
679 693 self.writeMetadata(fp)
680 694 fp.close()
681 695 return
682 696
683 697 def createMetadataFile(self):
684 698 ext = self.ext
685 699 path = self.path
686 700 setFile = self.setFile
687 701
688 702 timeTuple = time.localtime(self.dataOut.utctime)
689 703
690 704 subfolder = ''
691 705 fullpath = os.path.join( path, subfolder )
692 706
693 707 if not( os.path.exists(fullpath) ):
694 708 os.mkdir(fullpath)
695 709 setFile = -1 #inicializo mi contador de seteo
696 710
697 711 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
698 712 fullpath = os.path.join( path, subfolder )
699 713
700 714 if not( os.path.exists(fullpath) ):
701 715 os.mkdir(fullpath)
702 716 setFile = -1 #inicializo mi contador de seteo
703 717
704 718 else:
705 719 filesList = os.listdir( fullpath )
706 720 filesList = sorted( filesList, key=str.lower )
707 721 if len( filesList ) > 0:
708 722 filesList = [k for k in filesList if 'M' in k]
709 723 filen = filesList[-1]
710 724 # el filename debera tener el siguiente formato
711 725 # 0 1234 567 89A BCDE (hex)
712 726 # x YYYY DDD SSS .ext
713 727 if isNumber( filen[8:11] ):
714 728 setFile = int( filen[8:11] ) #inicializo mi contador de seteo al seteo del ultimo file
715 729 else:
716 730 setFile = -1
717 731 else:
718 732 setFile = -1 #inicializo mi contador de seteo
719 733
720 734 setFile += 1
721 735
722 736 file = '%s%4.4d%3.3d%3.3d%s' % (self.metaoptchar,
723 737 timeTuple.tm_year,
724 738 timeTuple.tm_yday,
725 739 setFile,
726 740 ext )
727 741
728 742 filename = os.path.join( path, subfolder, file )
729 743 self.metaFile = file
730 744 #Setting HDF5 File
731 745 fp = h5py.File(filename,'w')
732 746
733 747 return fp
734 748
735 749 def writeMetadata(self, fp):
736 750
737 751 grp = fp.create_group("Metadata")
738 752 grp.create_dataset('array dimensions', data = self.tableDim, dtype = self.dtype)
739 753
740 754 for i in range(len(self.metadataList)):
741 755 grp.create_dataset(self.metadataList[i], data=getattr(self.dataOut, self.metadataList[i]))
742 756 return
743 757
744 758 def dateFlag(self):
745 759
746 760 timeTuple = time.localtime(self.dataOut.utctime)
747 761 dataDay = timeTuple.tm_yday
748 762
749 763 if dataDay == self.currentDay:
750 764 return False
751 765
752 766 self.currentDay = dataDay
753 767 return True
754 768
755 769 def setNextFile(self):
756 770
757 771 ext = self.ext
758 772 path = self.path
759 773 setFile = self.setFile
760 774 mode = self.mode
761 775
762 776 timeTuple = time.localtime(self.dataOut.utctime)
763 777 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
764 778
765 779 fullpath = os.path.join( path, subfolder )
766 780
767 781 if os.path.exists(fullpath):
768 782 filesList = os.listdir( fullpath )
769 783 filesList = [k for k in filesList if 'D' in k]
770 784 if len( filesList ) > 0:
771 785 filesList = sorted( filesList, key=str.lower )
772 786 filen = filesList[-1]
773 787 # el filename debera tener el siguiente formato
774 788 # 0 1234 567 89A BCDE (hex)
775 789 # x YYYY DDD SSS .ext
776 790 if isNumber( filen[8:11] ):
777 791 setFile = int( filen[8:11] ) #inicializo mi contador de seteo al seteo del ultimo file
778 792 else:
779 793 setFile = -1
780 794 else:
781 795 setFile = -1 #inicializo mi contador de seteo
782 796 else:
783 797 os.mkdir(fullpath)
784 798 setFile = -1 #inicializo mi contador de seteo
785 799
786 800 setFile += 1
787 801
788 802 file = '%s%4.4d%3.3d%3.3d%s' % (self.optchar,
789 803 timeTuple.tm_year,
790 804 timeTuple.tm_yday,
791 805 setFile,
792 806 ext )
793 807
794 808 filename = os.path.join( path, subfolder, file )
795 809
796 810 #Setting HDF5 File
797 811 fp = h5py.File(filename,'w')
798 812
799 813 #writemetadata
800 814 self.writeMetadata(fp)
801 815
802 816 grp = fp.create_group("Data")
803 817 # grp.attrs['metadata'] = self.metaFile
804 818
805 819 # grp.attrs['blocksPerFile'] = 0
806 820
807 821 ds = []
808 822 data = []
809 823 nDimsForDs = []
810 824
811 825 nDatas = numpy.zeros(len(self.dataList))
812 826 nDims = self.arrayDim[:,0]
813 827
814 828 nDim1 = self.arrayDim[:,2]
815 829 nDim0 = self.arrayDim[:,3]
816 830
817 831 for i in range(len(self.dataList)):
818 832
833 #One-dimension data
819 834 if nDims[i]==1:
820 835 # ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype='S20')
821 836 ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype=numpy.float64)
822 837 ds.append(ds0)
823 838 data.append([])
824 839 nDimsForDs.append(nDims[i])
825 840 else:
826 841
827 if mode[i]==0:
828 strMode = "channel"
829 nDatas[i] = self.arrayDim[i,1]
830 else:
842 #Channel mode
843 # if mode[i] == 0:
844 # strMode = "channel"
845 #
846 # #nDatas is the number of arrays per variable
847 # if nDims[i] == 1:
848 # nDatas[i] = self.arrayDim[i,1]
849 # elif nDims[i] == 2:
850 # nDatas[i] = self.arrayDim[i,2]
851
852 #Parameters mode
853 if mode[i] == 1:
831 854 strMode = "param"
832 855 nDatas[i] = self.arrayDim[i,2]
833
834 if nDims[i]==2:
835 nDatas[i] = self.arrayDim[i,2]
836
856
857 #Meteors mode
858 elif mode[i] == 2:
859 strMode = "table"
860 nDatas[i] = 1
861
837 862 grp0 = grp.create_group(self.dataList[i])
838 863
839 864 for j in range(int(nDatas[i])):
840 865 tableName = strMode + str(j)
841 866
842 867 if nDims[i] == 3:
843 868 ds0 = grp0.create_dataset(tableName, (nDim1[i],nDim0[i],1) , data = numpy.zeros((nDim1[i],nDim0[i],1)) ,maxshape=(None,nDim0[i],None), chunks=True)
869
844 870 else:
845 871 ds0 = grp0.create_dataset(tableName, (1,nDim0[i]), data = numpy.zeros((1,nDim0[i])) , maxshape=(None,nDim0[i]), chunks=True)
846 872
847 873 ds.append(ds0)
848 874 data.append([])
849 875 nDimsForDs.append(nDims[i])
876
877 fp.flush()
878 fp.close()
879
850 880 self.nDatas = nDatas
851 881 self.nDims = nDims
852 882 self.nDimsForDs = nDimsForDs
853 883 #Saving variables
854 884 print 'Writing the file: %s'%filename
855 885 self.filename = filename
856 self.fp = fp
857 self.grp = grp
858 self.grp.attrs.modify('nRecords', 1)
886 # self.fp = fp
887 # self.grp = grp
888 # self.grp.attrs.modify('nRecords', 1)
859 889 self.ds = ds
860 890 self.data = data
861
862 self.setFile = setFile
891 #
892 # self.setFile = setFile
863 893 self.firsttime = True
864 894 self.blockIndex = 0
865 895 return
866 896
867 897 def putData(self):
868
869 if not self.firsttime:
870 self.readBlock()
871 898
872 899 if self.blockIndex == self.blocksPerFile or self.dateFlag():
873
874 900 self.setNextFile()
875 901
876 self.setBlock()
877 self.writeBlock()
878
879 self.fp.flush()
880 self.fp.close()
902 # if not self.firsttime:
903 self.readBlock()
904 self.setBlock() #Prepare data to be written
905 self.writeBlock() #Write data
881 906
882 907 return
883 908
884 909 def readBlock(self):
885 910
886 911 '''
887 912 data Array configured
888 913
889 914
890 915 self.data
891 916 '''
892 917 ds = self.ds
893 918 #Setting HDF5 File
894 919 fp = h5py.File(self.filename,'r+')
895 920 grp = fp["Data"]
896 921 ind = 0
897 922
898 923 # grp.attrs['blocksPerFile'] = 0
899 924 for i in range(len(self.dataList)):
900 925
901 926 if self.nDims[i]==1:
902 927 ds0 = grp[self.dataList[i]]
903 928 ds[ind] = ds0
904 929 ind += 1
905 930 else:
906 if self.mode[i]==0:
907 strMode = "channel"
908 else:
931 # if self.mode[i] == 0:
932 # strMode = "channel"
933 if self.mode[i] == 1:
909 934 strMode = "param"
910
935 elif self.mode[i] == 2:
936 strMode = "table"
937
911 938 grp0 = grp[self.dataList[i]]
912 939
913 940 for j in range(int(self.nDatas[i])):
914 941 tableName = strMode + str(j)
915 942 ds0 = grp0[tableName]
916 943 ds[ind] = ds0
917 944 ind += 1
918
919
945
920 946 self.fp = fp
921 947 self.grp = grp
922 948 self.ds = ds
923 949
924 950 return
925 951
926 952 def setBlock(self):
927 953 '''
928 954 data Array configured
929 955
930 956
931 957 self.data
932 958 '''
933 959 #Creating Arrays
934 960 data = self.data
935 961 nDatas = self.nDatas
936 962 nDims = self.nDims
937 963 mode = self.mode
938 964 ind = 0
939 965
940 966 for i in range(len(self.dataList)):
941 967 dataAux = getattr(self.dataOut,self.dataList[i])
942 968
943 if nDims[i] == 1:
944 # data[ind] = numpy.array([str(dataAux)]).reshape((1,1))
969 if nDims[i] == 1 or mode[i] == 2:
945 970 data[ind] = dataAux
946 # if not self.firsttime:
947 # data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind]))
948 ind += 1
949 else:
971 ind += 1
972
973 elif nDims[i] == 2:
950 974 for j in range(int(nDatas[i])):
951 if (mode[i] == 0) or (nDims[i] == 2): #In case division per channel or Dimensions is only 1
952 data[ind] = dataAux[j,:]
953 else:
954 data[ind] = dataAux[:,j,:]
975 data[ind] = dataAux[j,:]
976 ind += 1
955 977
956 # if nDims[i] == 3:
957 # data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1))
958
959 # if not self.firsttime:
960 # data[ind] = numpy.dstack((self.ds[ind][:], data[ind]))
961
978 elif nDims[i] == 3:
979 for j in range(int(nDatas[i])):
980 # Extinct mode 0
981 # if (mode[i] == 0):
982 # data[ind] = dataAux[j,:,:]
962 983 # else:
963 # data[ind] = data[ind].reshape((1,data[ind].shape[0]))
964
965 # if not self.firsttime:
966 # data[ind] = numpy.vstack((self.ds[ind][:], data[ind]))
984 data[ind] = dataAux[:,j,:]
967 985 ind += 1
968
986
969 987 self.data = data
970 988 return
971 989
972 990 def writeBlock(self):
973 991 '''
974 992 Saves the block in the HDF5 file
975 993 '''
976 994 for i in range(len(self.ds)):
995
996 # First time
977 997 if self.firsttime:
978 998 # self.ds[i].resize(self.data[i].shape)
979 999 # self.ds[i][self.blockIndex,:] = self.data[i]
980 1000 if type(self.data[i]) == numpy.ndarray:
981 1001 nDims1 = len(self.ds[i].shape)
982 1002
983 1003 if nDims1 == 3:
984 1004 self.data[i] = self.data[i].reshape((self.data[i].shape[0],self.data[i].shape[1],1))
985
986 self.ds[i].resize(self.data[i].shape)
987 self.ds[i][:] = self.data[i]
1005
1006 self.ds[i].resize(self.data[i].shape)
1007
1008 self.ds[i][:] = self.data[i]
988 1009 else:
989 if self.nDimsForDs[i] == 1:
1010
1011 # From second time
1012 # Meteors!
1013 if self.mode[i] == 2:
1014 dataShape = self.data[i].shape
1015 dsShape = self.ds[i].shape
1016 self.ds[i].resize((self.ds[i].shape[0] + dataShape[0],self.ds[i].shape[1]))
1017 self.ds[i][dsShape[0]:,:] = self.data[i]
1018 # One dimension
1019 elif self.nDimsForDs[i] == 1:
990 1020 self.ds[i].resize((self.ds[i].shape[0], self.ds[i].shape[1] + 1))
991 1021 self.ds[i][0,-1] = self.data[i]
1022 # Two dimension
992 1023 elif self.nDimsForDs[i] == 2:
993 1024 self.ds[i].resize((self.ds[i].shape[0] + 1,self.ds[i].shape[1]))
994 1025 self.ds[i][self.blockIndex,:] = self.data[i]
1026 # Three dimensions
995 1027 elif self.nDimsForDs[i] == 3:
996
997 dataShape = self.data[i].shape
998 dsShape = self.ds[i].shape
999
1000 if dataShape[0]==dsShape[0]:
1001 self.ds[i].resize((self.ds[i].shape[0],self.ds[i].shape[1],self.ds[i].shape[2]+1))
1002 self.ds[i][:,:,-1] = self.data[i]
1003 else:
1004 self.ds[i].resize((self.ds[i].shape[0] + dataShape[0],self.ds[i].shape[1],self.ds[i].shape[2]))
1005 self.ds[i][dsShape[0]:,:,0] = self.data[i]
1028 self.ds[i].resize((self.ds[i].shape[0],self.ds[i].shape[1],self.ds[i].shape[2]+1))
1029 self.ds[i][:,:,-1] = self.data[i]
1006 1030
1031 self.firsttime = False
1007 1032 self.blockIndex += 1
1008 self.firsttime = False
1033
1034 #Close to save changes
1035 self.fp.flush()
1036 self.fp.close()
1009 1037 return
1010 1038
1011 1039 def run(self, dataOut, **kwargs):
1012 1040
1013 1041 if not(self.isConfig):
1014 1042 flagdata = self.setup(dataOut, **kwargs)
1015 1043
1016 1044 if not(flagdata):
1017 1045 return
1018 1046
1019 1047 self.isConfig = True
1020 1048 # self.putMetadata()
1021 1049 self.setNextFile()
1022 1050
1023 1051 self.putData()
1024 1052 return
1025 1053
1026 1054
Show file before | Show file after | Hide comments
@@ -1,2154 +1,2169
1 1 import numpy
2 2 import math
3 3 from scipy import optimize
4 4 from scipy import interpolate
5 5 from scipy import signal
6 6 from scipy import stats
7 7 import re
8 8 import datetime
9 9 import copy
10 10 import sys
11 11 import importlib
12 12 import itertools
13 13
14 14 from jroproc_base import ProcessingUnit, Operation
15 15 from schainpy.model.data.jrodata import Parameters
16 16
17 17
18 18 class ParametersProc(ProcessingUnit):
19 19
20 20 nSeconds = None
21 21
22 22 def __init__(self):
23 23 ProcessingUnit.__init__(self)
24 24
25 25 # self.objectDict = {}
26 26 self.buffer = None
27 27 self.firstdatatime = None
28 28 self.profIndex = 0
29 29 self.dataOut = Parameters()
30 30
31 31 def __updateObjFromInput(self):
32 32
33 33 self.dataOut.inputUnit = self.dataIn.type
34 34
35 35 self.dataOut.timeZone = self.dataIn.timeZone
36 36 self.dataOut.dstFlag = self.dataIn.dstFlag
37 37 self.dataOut.errorCount = self.dataIn.errorCount
38 38 self.dataOut.useLocalTime = self.dataIn.useLocalTime
39 39
40 40 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
41 41 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
42 42 self.dataOut.channelList = self.dataIn.channelList
43 43 self.dataOut.heightList = self.dataIn.heightList
44 44 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
45 45 # self.dataOut.nHeights = self.dataIn.nHeights
46 46 # self.dataOut.nChannels = self.dataIn.nChannels
47 47 self.dataOut.nBaud = self.dataIn.nBaud
48 48 self.dataOut.nCode = self.dataIn.nCode
49 49 self.dataOut.code = self.dataIn.code
50 50 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
51 51 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
52 52 self.dataOut.utctime = self.firstdatatime
53 53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
54 54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
55 55 # self.dataOut.nCohInt = self.dataIn.nCohInt
56 56 # self.dataOut.nIncohInt = 1
57 57 self.dataOut.ippSeconds = self.dataIn.ippSeconds
58 58 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59 self.dataOut.timeInterval = self.dataIn.timeInterval
60 60 self.dataOut.heightList = self.dataIn.getHeiRange()
61 61 self.dataOut.frequency = self.dataIn.frequency
62 62
63 63 def run(self, nSeconds = 100, nProfiles = None):
64 64
65 65
66 66
67 67 if self.firstdatatime == None:
68 68 self.firstdatatime = self.dataIn.utctime
69 69
70 70 #---------------------- Voltage Data ---------------------------
71 71
72 72 if self.dataIn.type == "Voltage":
73 73 self.dataOut.flagNoData = True
74 74
75 75
76 76 if self.buffer == None:
77 77 self.nSeconds = nSeconds
78 78 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
79 79
80 80 self.buffer = numpy.zeros((self.dataIn.nChannels,
81 81 self.nProfiles,
82 82 self.dataIn.nHeights),
83 83 dtype='complex')
84 84
85 85 if self.profIndex == 7990:
86 86 a = 1
87 87
88 88 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
89 89 self.profIndex += 1
90 90
91 91 if self.profIndex == self.nProfiles:
92 92
93 93 self.__updateObjFromInput()
94 94 self.dataOut.data_pre = self.buffer.copy()
95 95 self.dataOut.paramInterval = nSeconds
96 96 self.dataOut.flagNoData = False
97 97
98 98 self.buffer = None
99 99 self.firstdatatime = None
100 100 self.profIndex = 0
101 101 return
102 102
103 103 #---------------------- Spectra Data ---------------------------
104 104
105 105 if self.dataIn.type == "Spectra":
106 106 self.dataOut.data_pre = self.dataIn.data_spc.copy()
107 107 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
108 108 self.dataOut.noise = self.dataIn.getNoise()
109 109 self.dataOut.normFactor = self.dataIn.normFactor
110 110 self.dataOut.groupList = self.dataIn.pairsList
111 111 self.dataOut.flagNoData = False
112 112
113 113 #---------------------- Correlation Data ---------------------------
114 114
115 115 if self.dataIn.type == "Correlation":
116 116 lagRRange = self.dataIn.lagR
117 117 indR = numpy.where(lagRRange == 0)[0][0]
118 118
119 119 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
120 120 self.dataOut.abscissaList = self.dataIn.getLagTRange(1)
121 121 self.dataOut.noise = self.dataIn.noise
122 122 self.dataOut.normFactor = self.dataIn.normFactor
123 123 self.dataOut.data_SNR = self.dataIn.SNR
124 124 self.dataOut.groupList = self.dataIn.pairsList
125 125 self.dataOut.flagNoData = False
126 126
127 127 #---------------------- Correlation Data ---------------------------
128 128
129 129 if self.dataIn.type == "Parameters":
130 130 self.dataOut.copy(self.dataIn)
131 131 self.dataOut.flagNoData = False
132 132
133 133 return True
134 134
135 135 self.__updateObjFromInput()
136 136 self.firstdatatime = None
137 137 self.dataOut.utctimeInit = self.dataIn.utctime
138 138 self.dataOut.paramInterval = self.dataIn.timeInterval
139 139
140 140 #------------------- Get Moments ----------------------------------
141 141 def GetMoments(self, channelList = None):
142 142 '''
143 143 Function GetMoments()
144 144
145 145 Input:
146 146 channelList : simple channel list to select e.g. [2,3,7]
147 147 self.dataOut.data_pre
148 148 self.dataOut.abscissaList
149 149 self.dataOut.noise
150 150
151 151 Affected:
152 152 self.dataOut.data_param
153 153 self.dataOut.data_SNR
154 154
155 155 '''
156 156 data = self.dataOut.data_pre
157 157 absc = self.dataOut.abscissaList[:-1]
158 158 noise = self.dataOut.noise
159 159
160 160 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
161 161
162 162 if channelList== None:
163 163 channelList = self.dataIn.channelList
164 164 self.dataOut.channelList = channelList
165 165
166 166 for ind in channelList:
167 167 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
168 168
169 169 self.dataOut.data_param = data_param[:,1:,:]
170 170 self.dataOut.data_SNR = data_param[:,0]
171 171 return
172 172
173 173 def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
174 174
175 175 if (nicoh == None): nicoh = 1
176 176 if (graph == None): graph = 0
177 177 if (smooth == None): smooth = 0
178 178 elif (self.smooth < 3): smooth = 0
179 179
180 180 if (type1 == None): type1 = 0
181 181 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
182 182 if (snrth == None): snrth = -3
183 183 if (dc == None): dc = 0
184 184 if (aliasing == None): aliasing = 0
185 185 if (oldfd == None): oldfd = 0
186 186 if (wwauto == None): wwauto = 0
187 187
188 188 if (n0 < 1.e-20): n0 = 1.e-20
189 189
190 190 freq = oldfreq
191 191 vec_power = numpy.zeros(oldspec.shape[1])
192 192 vec_fd = numpy.zeros(oldspec.shape[1])
193 193 vec_w = numpy.zeros(oldspec.shape[1])
194 194 vec_snr = numpy.zeros(oldspec.shape[1])
195 195
196 196 for ind in range(oldspec.shape[1]):
197 197
198 198 spec = oldspec[:,ind]
199 199 aux = spec*fwindow
200 200 max_spec = aux.max()
201 201 m = list(aux).index(max_spec)
202 202
203 203 #Smooth
204 204 if (smooth == 0): spec2 = spec
205 205 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
206 206
207 207 # Calculo de Momentos
208 208 bb = spec2[range(m,spec2.size)]
209 209 bb = (bb<n0).nonzero()
210 210 bb = bb[0]
211 211
212 212 ss = spec2[range(0,m + 1)]
213 213 ss = (ss<n0).nonzero()
214 214 ss = ss[0]
215 215
216 216 if (bb.size == 0):
217 217 bb0 = spec.size - 1 - m
218 218 else:
219 219 bb0 = bb[0] - 1
220 220 if (bb0 < 0):
221 221 bb0 = 0
222 222
223 223 if (ss.size == 0): ss1 = 1
224 224 else: ss1 = max(ss) + 1
225 225
226 226 if (ss1 > m): ss1 = m
227 227
228 228 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
229 229 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
230 230 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
231 231 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
232 232 snr = (spec2.mean()-n0)/n0
233 233
234 234 if (snr < 1.e-20) :
235 235 snr = 1.e-20
236 236
237 237 vec_power[ind] = power
238 238 vec_fd[ind] = fd
239 239 vec_w[ind] = w
240 240 vec_snr[ind] = snr
241 241
242 242 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
243 243 return moments
244 244
245 245 #------------------ Get SA Parameters --------------------------
246 246
247 247 def GetSAParameters(self):
248 248 pairslist = self.dataOut.groupList
249 249 num_pairs = len(pairslist)
250 250
251 251 vel = self.dataOut.abscissaList
252 252 spectra = self.dataOut.data_pre
253 253 cspectra = self.dataIn.data_cspc
254 254 delta_v = vel[1] - vel[0]
255 255
256 256 #Calculating the power spectrum
257 257 spc_pow = numpy.sum(spectra, 3)*delta_v
258 258 #Normalizing Spectra
259 259 norm_spectra = spectra/spc_pow
260 260 #Calculating the norm_spectra at peak
261 261 max_spectra = numpy.max(norm_spectra, 3)
262 262
263 263 #Normalizing Cross Spectra
264 264 norm_cspectra = numpy.zeros(cspectra.shape)
265 265
266 266 for i in range(num_chan):
267 267 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
268 268
269 269 max_cspectra = numpy.max(norm_cspectra,2)
270 270 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
271 271
272 272 for i in range(num_pairs):
273 273 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
274 274 #------------------- Get Lags ----------------------------------
275 275
276 276 def GetLags(self):
277 277 '''
278 278 Function GetMoments()
279 279
280 280 Input:
281 281 self.dataOut.data_pre
282 282 self.dataOut.abscissaList
283 283 self.dataOut.noise
284 284 self.dataOut.normFactor
285 285 self.dataOut.data_SNR
286 286 self.dataOut.groupList
287 287 self.dataOut.nChannels
288 288
289 289 Affected:
290 290 self.dataOut.data_param
291 291
292 292 '''
293 293
294 294 data = self.dataOut.data_pre
295 295 normFactor = self.dataOut.normFactor
296 296 nHeights = self.dataOut.nHeights
297 297 absc = self.dataOut.abscissaList[:-1]
298 298 noise = self.dataOut.noise
299 299 SNR = self.dataOut.data_SNR
300 300 pairsList = self.dataOut.groupList
301 301 nChannels = self.dataOut.nChannels
302 302 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
303 303 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
304 304
305 305 dataNorm = numpy.abs(data)
306 306 for l in range(len(pairsList)):
307 307 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
308 308
309 309 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
310 310 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
311 311 return
312 312
313 313 def __getPairsAutoCorr(self, pairsList, nChannels):
314 314
315 315 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
316 316
317 317 for l in range(len(pairsList)):
318 318 firstChannel = pairsList[l][0]
319 319 secondChannel = pairsList[l][1]
320 320
321 321 #Obteniendo pares de Autocorrelacion
322 322 if firstChannel == secondChannel:
323 323 pairsAutoCorr[firstChannel] = int(l)
324 324
325 325 pairsAutoCorr = pairsAutoCorr.astype(int)
326 326
327 327 pairsCrossCorr = range(len(pairsList))
328 328 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
329 329
330 330 return pairsAutoCorr, pairsCrossCorr
331 331
332 332 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
333 333
334 334 Pt0 = data.shape[1]/2
335 335 #Funcion de Autocorrelacion
336 336 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
337 337
338 338 #Obtencion Indice de TauCross
339 339 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
340 340 #Obtencion Indice de TauAuto
341 341 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
342 342 CCValue = data[pairsCrossCorr,Pt0,:]
343 343 for i in range(pairsCrossCorr.size):
344 344 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
345 345
346 346 #Obtencion de TauCross y TauAuto
347 347 tauCross = lagTRange[indCross]
348 348 tauAuto = lagTRange[indAuto]
349 349
350 350 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
351 351
352 352 tauCross[Nan1,Nan2] = numpy.nan
353 353 tauAuto[Nan1,Nan2] = numpy.nan
354 354 tau = numpy.vstack((tauCross,tauAuto))
355 355
356 356 return tau
357 357
358 358 def __calculateLag1Phase(self, data, pairs, lagTRange):
359 359 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
360 360 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
361 361
362 362 phase = numpy.angle(data1[lag1,:])
363 363
364 364 return phase
365 365 #------------------- Detect Meteors ------------------------------
366 366
367 367 def MeteorDetection(self, hei_ref = None, tauindex = 0,
368 predefinedPhaseShifts = None, centerReceiverIndex = 2, saveAll = False,
368 predefinedPhaseShifts = None, centerReceiverIndex = 2,
369 369 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
370 370 noise_timeStep = 4, noise_multiple = 4,
371 371 multDet_timeLimit = 1, multDet_rangeLimit = 3,
372 372 phaseThresh = 20, SNRThresh = 8,
373 373 hmin = 70, hmax=110, azimuth = 0) :
374 374
375 375 '''
376 376 Function DetectMeteors()
377 377 Project developed with paper:
378 378 HOLDSWORTH ET AL. 2004
379 379
380 380 Input:
381 381 self.dataOut.data_pre
382 382
383 383 centerReceiverIndex: From the channels, which is the center receiver
384 384
385 385 hei_ref: Height reference for the Beacon signal extraction
386 386 tauindex:
387 387 predefinedPhaseShifts: Predefined phase offset for the voltge signals
388 388
389 389 cohDetection: Whether to user Coherent detection or not
390 390 cohDet_timeStep: Coherent Detection calculation time step
391 391 cohDet_thresh: Coherent Detection phase threshold to correct phases
392 392
393 393 noise_timeStep: Noise calculation time step
394 394 noise_multiple: Noise multiple to define signal threshold
395 395
396 396 multDet_timeLimit: Multiple Detection Removal time limit in seconds
397 397 multDet_rangeLimit: Multiple Detection Removal range limit in km
398 398
399 399 phaseThresh: Maximum phase difference between receiver to be consider a meteor
400 400 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
401 401
402 402 hmin: Minimum Height of the meteor to use it in the further wind estimations
403 403 hmax: Maximum Height of the meteor to use it in the further wind estimations
404 404 azimuth: Azimuth angle correction
405 405
406 406 Affected:
407 407 self.dataOut.data_param
408 408
409 409 Rejection Criteria (Errors):
410 410 0: No error; analysis OK
411 411 1: SNR < SNR threshold
412 412 2: angle of arrival (AOA) ambiguously determined
413 413 3: AOA estimate not feasible
414 414 4: Large difference in AOAs obtained from different antenna baselines
415 415 5: echo at start or end of time series
416 416 6: echo less than 5 examples long; too short for analysis
417 417 7: echo rise exceeds 0.3s
418 418 8: echo decay time less than twice rise time
419 419 9: large power level before echo
420 420 10: large power level after echo
421 421 11: poor fit to amplitude for estimation of decay time
422 422 12: poor fit to CCF phase variation for estimation of radial drift velocity
423 423 13: height unresolvable echo: not valid height within 70 to 110 km
424 424 14: height ambiguous echo: more then one possible height within 70 to 110 km
425 425 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
426 426 16: oscilatory echo, indicating event most likely not an underdense echo
427 427
428 428 17: phase difference in meteor Reestimation
429 429
430 430 Data Storage:
431 431 Meteors for Wind Estimation (8):
432 432 Day Hour | Range Height
433 433 Azimuth Zenith errorCosDir
434 434 VelRad errorVelRad
435 435 TypeError
436 436
437 437 '''
438 438 #Get Beacon signal
439 439 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
440 440
441 441 if hei_ref != None:
442 442 newheis = numpy.where(self.dataOut.heightList>hei_ref)
443 443
444 444 heiRang = self.dataOut.getHeiRange()
445 445 #Pairs List
446 446 pairslist = []
447 447 nChannel = self.dataOut.nChannels
448 448 for i in range(nChannel):
449 449 if i != centerReceiverIndex:
450 450 pairslist.append((centerReceiverIndex,i))
451 451
452 452 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
453 453 # see if the user put in pre defined phase shifts
454 454 voltsPShift = self.dataOut.data_pre.copy()
455 455
456 456 if predefinedPhaseShifts != None:
457 457 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
458 458
459 elif beaconPhaseShifts:
460 #get hardware phase shifts using beacon signal
461 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
462 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
459 # elif beaconPhaseShifts:
460 # #get hardware phase shifts using beacon signal
461 # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
462 # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
463 463
464 464 else:
465 465 hardwarePhaseShifts = numpy.zeros(5)
466 466
467 467
468 468 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
469 469 for i in range(self.dataOut.data_pre.shape[0]):
470 470 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
471 471
472 472
473 473 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
474 474
475 475 #Remove DC
476 476 voltsDC = numpy.mean(voltsPShift,1)
477 477 voltsDC = numpy.mean(voltsDC,1)
478 478 for i in range(voltsDC.shape[0]):
479 479 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
480 480
481 481 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
482 482 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
483 483
484 484 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
485 485 #Coherent Detection
486 486 if cohDetection:
487 487 #use coherent detection to get the net power
488 488 cohDet_thresh = cohDet_thresh*numpy.pi/180
489 489 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
490 490
491 491 #Non-coherent detection!
492 492 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
493 493 #********** END OF COH/NON-COH POWER CALCULATION**********************
494 494
495 495 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
496 496 #Get noise
497 497 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
498 498 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
499 499 #Get signal threshold
500 500 signalThresh = noise_multiple*noise
501 501 #Meteor echoes detection
502 502 listMeteors = self.__findMeteors(powerNet, signalThresh)
503 503 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
504 504
505 505 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
506 506 #Parameters
507 507 heiRange = self.dataOut.getHeiRange()
508 508 rangeInterval = heiRange[1] - heiRange[0]
509 509 rangeLimit = multDet_rangeLimit/rangeInterval
510 510 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
511 511 #Multiple detection removals
512 512 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
513 513 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
514 514
515 515 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
516 516 #Parameters
517 517 phaseThresh = phaseThresh*numpy.pi/180
518 518 thresh = [phaseThresh, noise_multiple, SNRThresh]
519 519 #Meteor reestimation (Errors N 1, 6, 12, 17)
520 520 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
521 521 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
522 522 #Estimation of decay times (Errors N 7, 8, 11)
523 523 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
524 524 #******************* END OF METEOR REESTIMATION *******************
525 525
526 526 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
527 527 #Calculating Radial Velocity (Error N 15)
528 528 radialStdThresh = 10
529 529 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
530 530
531 531 if len(listMeteors4) > 0:
532 #Setting New Array
533 date = self.dataOut.utctime
534 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
532 535
533 536 pairsList = []
534 537 pairx = (0,3)
535 538 pairy = (1,2)
536 539 pairsList.append(pairx)
537 540 pairsList.append(pairy)
538 541
539 #Setting New Array
540 date = repr(self.dataOut.datatime)
541 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
542
543 542 meteorOps = MeteorOperations()
544 543 jph = numpy.array([0,0,0,0])
545 544 h = (hmin,hmax)
546 545 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph)
547 546
548 547 # #Calculate AOA (Error N 3, 4)
549 548 # #JONES ET AL. 1998
550 549 # error = arrayParameters[:,-1]
551 550 # AOAthresh = numpy.pi/8
552 551 # phases = -arrayParameters[:,9:13]
553 552 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
554 553 #
555 554 # #Calculate Heights (Error N 13 and 14)
556 555 # error = arrayParameters[:,-1]
557 556 # Ranges = arrayParameters[:,2]
558 557 # zenith = arrayParameters[:,5]
559 558 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
560 559 # error = arrayParameters[:,-1]
561 560 #********************* END OF PARAMETERS CALCULATION **************************
562 561
563 562 #***************************+ PASS DATA TO NEXT STEP **********************
564 arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
565 self.dataOut.data_param = arrayFinal
563 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
564 self.dataOut.data_param = arrayParameters
566 565
567 if arrayFinal == None:
566 if arrayParameters == None:
568 567 self.dataOut.flagNoData = True
569 568
570 569 return
571 570
571 def correctMeteorPhases(self):
572
573 arrayParameters = self.dataOut.data_param
574 pairsList = []
575 pairx = (0,3)
576 pairy = (1,2)
577 pairsList.append(pairx)
578 pairsList.append(pairy)
579
580 meteorOps = MeteorOperations()
581 jph = numpy.array([0,0,0,0])
582 h = (hmin,hmax)
583 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph)
584 self.dataOut.data_param = arrayParameters
585 return
586
572 587 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
573 588
574 589 minIndex = min(newheis[0])
575 590 maxIndex = max(newheis[0])
576 591
577 592 voltage = voltage0[:,:,minIndex:maxIndex+1]
578 593 nLength = voltage.shape[1]/n
579 594 nMin = 0
580 595 nMax = 0
581 596 phaseOffset = numpy.zeros((len(pairslist),n))
582 597
583 598 for i in range(n):
584 599 nMax += nLength
585 600 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
586 601 phaseCCF = numpy.mean(phaseCCF, axis = 2)
587 602 phaseOffset[:,i] = phaseCCF.transpose()
588 603 nMin = nMax
589 604 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
590 605
591 606 #Remove Outliers
592 607 factor = 2
593 608 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
594 609 dw = numpy.std(wt,axis = 1)
595 610 dw = dw.reshape((dw.size,1))
596 611 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
597 612 phaseOffset[ind] = numpy.nan
598 613 phaseOffset = stats.nanmean(phaseOffset, axis=1)
599 614
600 615 return phaseOffset
601 616
602 617 def __shiftPhase(self, data, phaseShift):
603 618 #this will shift the phase of a complex number
604 619 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
605 620 return dataShifted
606 621
607 622 def __estimatePhaseDifference(self, array, pairslist):
608 623 nChannel = array.shape[0]
609 624 nHeights = array.shape[2]
610 625 numPairs = len(pairslist)
611 626 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
612 627 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
613 628
614 629 #Correct phases
615 630 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
616 631 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
617 632
618 633 if indDer[0].shape[0] > 0:
619 634 for i in range(indDer[0].shape[0]):
620 635 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
621 636 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
622 637
623 638 # for j in range(numSides):
624 639 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
625 640 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
626 641 #
627 642 #Linear
628 643 phaseInt = numpy.zeros((numPairs,1))
629 644 angAllCCF = phaseCCF[:,[0,1,3,4],0]
630 645 for j in range(numPairs):
631 646 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
632 647 phaseInt[j] = fit[1]
633 648 #Phase Differences
634 649 phaseDiff = phaseInt - phaseCCF[:,2,:]
635 650 phaseArrival = phaseInt.reshape(phaseInt.size)
636 651
637 652 #Dealias
638 653 indAlias = numpy.where(phaseArrival > numpy.pi)
639 654 phaseArrival[indAlias] -= 2*numpy.pi
640 655 indAlias = numpy.where(phaseArrival < -numpy.pi)
641 656 phaseArrival[indAlias] += 2*numpy.pi
642 657
643 658 return phaseDiff, phaseArrival
644 659
645 660 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
646 661 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
647 662 #find the phase shifts of each channel over 1 second intervals
648 663 #only look at ranges below the beacon signal
649 664 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
650 665 numBlocks = int(volts.shape[1]/numProfPerBlock)
651 666 numHeights = volts.shape[2]
652 667 nChannel = volts.shape[0]
653 668 voltsCohDet = volts.copy()
654 669
655 670 pairsarray = numpy.array(pairslist)
656 671 indSides = pairsarray[:,1]
657 672 # indSides = numpy.array(range(nChannel))
658 673 # indSides = numpy.delete(indSides, indCenter)
659 674 #
660 675 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
661 676 listBlocks = numpy.array_split(volts, numBlocks, 1)
662 677
663 678 startInd = 0
664 679 endInd = 0
665 680
666 681 for i in range(numBlocks):
667 682 startInd = endInd
668 683 endInd = endInd + listBlocks[i].shape[1]
669 684
670 685 arrayBlock = listBlocks[i]
671 686 # arrayBlockCenter = listCenter[i]
672 687
673 688 #Estimate the Phase Difference
674 689 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
675 690 #Phase Difference RMS
676 691 arrayPhaseRMS = numpy.abs(phaseDiff)
677 692 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
678 693 indPhase = numpy.where(phaseRMSaux==4)
679 694 #Shifting
680 695 if indPhase[0].shape[0] > 0:
681 696 for j in range(indSides.size):
682 697 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
683 698 voltsCohDet[:,startInd:endInd,:] = arrayBlock
684 699
685 700 return voltsCohDet
686 701
687 702 def __calculateCCF(self, volts, pairslist ,laglist):
688 703
689 704 nHeights = volts.shape[2]
690 705 nPoints = volts.shape[1]
691 706 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
692 707
693 708 for i in range(len(pairslist)):
694 709 volts1 = volts[pairslist[i][0]]
695 710 volts2 = volts[pairslist[i][1]]
696 711
697 712 for t in range(len(laglist)):
698 713 idxT = laglist[t]
699 714 if idxT >= 0:
700 715 vStacked = numpy.vstack((volts2[idxT:,:],
701 716 numpy.zeros((idxT, nHeights),dtype='complex')))
702 717 else:
703 718 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
704 719 volts2[:(nPoints + idxT),:]))
705 720 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
706 721
707 722 vStacked = None
708 723 return voltsCCF
709 724
710 725 def __getNoise(self, power, timeSegment, timeInterval):
711 726 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
712 727 numBlocks = int(power.shape[0]/numProfPerBlock)
713 728 numHeights = power.shape[1]
714 729
715 730 listPower = numpy.array_split(power, numBlocks, 0)
716 731 noise = numpy.zeros((power.shape[0], power.shape[1]))
717 732 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
718 733
719 734 startInd = 0
720 735 endInd = 0
721 736
722 737 for i in range(numBlocks): #split por canal
723 738 startInd = endInd
724 739 endInd = endInd + listPower[i].shape[0]
725 740
726 741 arrayBlock = listPower[i]
727 742 noiseAux = numpy.mean(arrayBlock, 0)
728 743 # noiseAux = numpy.median(noiseAux)
729 744 # noiseAux = numpy.mean(arrayBlock)
730 745 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
731 746
732 747 noiseAux1 = numpy.mean(arrayBlock)
733 748 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
734 749
735 750 return noise, noise1
736 751
737 752 def __findMeteors(self, power, thresh):
738 753 nProf = power.shape[0]
739 754 nHeights = power.shape[1]
740 755 listMeteors = []
741 756
742 757 for i in range(nHeights):
743 758 powerAux = power[:,i]
744 759 threshAux = thresh[:,i]
745 760
746 761 indUPthresh = numpy.where(powerAux > threshAux)[0]
747 762 indDNthresh = numpy.where(powerAux <= threshAux)[0]
748 763
749 764 j = 0
750 765
751 766 while (j < indUPthresh.size - 2):
752 767 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
753 768 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
754 769 indDNthresh = indDNthresh[indDNAux]
755 770
756 771 if (indDNthresh.size > 0):
757 772 indEnd = indDNthresh[0] - 1
758 773 indInit = indUPthresh[j]
759 774
760 775 meteor = powerAux[indInit:indEnd + 1]
761 776 indPeak = meteor.argmax() + indInit
762 777 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
763 778
764 779 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
765 780 j = numpy.where(indUPthresh == indEnd)[0] + 1
766 781 else: j+=1
767 782 else: j+=1
768 783
769 784 return listMeteors
770 785
771 786 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
772 787
773 788 arrayMeteors = numpy.asarray(listMeteors)
774 789 listMeteors1 = []
775 790
776 791 while arrayMeteors.shape[0] > 0:
777 792 FLAs = arrayMeteors[:,4]
778 793 maxFLA = FLAs.argmax()
779 794 listMeteors1.append(arrayMeteors[maxFLA,:])
780 795
781 796 MeteorInitTime = arrayMeteors[maxFLA,1]
782 797 MeteorEndTime = arrayMeteors[maxFLA,3]
783 798 MeteorHeight = arrayMeteors[maxFLA,0]
784 799
785 800 #Check neighborhood
786 801 maxHeightIndex = MeteorHeight + rangeLimit
787 802 minHeightIndex = MeteorHeight - rangeLimit
788 803 minTimeIndex = MeteorInitTime - timeLimit
789 804 maxTimeIndex = MeteorEndTime + timeLimit
790 805
791 806 #Check Heights
792 807 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
793 808 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
794 809 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
795 810
796 811 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
797 812
798 813 return listMeteors1
799 814
800 815 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
801 816 numHeights = volts.shape[2]
802 817 nChannel = volts.shape[0]
803 818
804 819 thresholdPhase = thresh[0]
805 820 thresholdNoise = thresh[1]
806 821 thresholdDB = float(thresh[2])
807 822
808 823 thresholdDB1 = 10**(thresholdDB/10)
809 824 pairsarray = numpy.array(pairslist)
810 825 indSides = pairsarray[:,1]
811 826
812 827 pairslist1 = list(pairslist)
813 828 pairslist1.append((0,1))
814 829 pairslist1.append((3,4))
815 830
816 831 listMeteors1 = []
817 832 listPowerSeries = []
818 833 listVoltageSeries = []
819 834 #volts has the war data
820 835
821 836 if frequency == 30e6:
822 837 timeLag = 45*10**-3
823 838 else:
824 839 timeLag = 15*10**-3
825 840 lag = numpy.ceil(timeLag/timeInterval)
826 841
827 842 for i in range(len(listMeteors)):
828 843
829 844 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
830 845 meteorAux = numpy.zeros(16)
831 846
832 847 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
833 848 mHeight = listMeteors[i][0]
834 849 mStart = listMeteors[i][1]
835 850 mPeak = listMeteors[i][2]
836 851 mEnd = listMeteors[i][3]
837 852
838 853 #get the volt data between the start and end times of the meteor
839 854 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
840 855 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
841 856
842 857 #3.6. Phase Difference estimation
843 858 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
844 859
845 860 #3.7. Phase difference removal & meteor start, peak and end times reestimated
846 861 #meteorVolts0.- all Channels, all Profiles
847 862 meteorVolts0 = volts[:,:,mHeight]
848 863 meteorThresh = noise[:,mHeight]*thresholdNoise
849 864 meteorNoise = noise[:,mHeight]
850 865 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
851 866 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
852 867
853 868 #Times reestimation
854 869 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
855 870 if mStart1.size > 0:
856 871 mStart1 = mStart1[-1] + 1
857 872
858 873 else:
859 874 mStart1 = mPeak
860 875
861 876 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
862 877 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
863 878 if mEndDecayTime1.size == 0:
864 879 mEndDecayTime1 = powerNet0.size
865 880 else:
866 881 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
867 882 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
868 883
869 884 #meteorVolts1.- all Channels, from start to end
870 885 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
871 886 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
872 887 if meteorVolts2.shape[1] == 0:
873 888 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
874 889 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
875 890 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
876 891 ##################### END PARAMETERS REESTIMATION #########################
877 892
878 893 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
879 894 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
880 895 if meteorVolts2.shape[1] > 0:
881 896 #Phase Difference re-estimation
882 897 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
883 898 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
884 899 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
885 900 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
886 901 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
887 902
888 903 #Phase Difference RMS
889 904 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
890 905 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
891 906 #Data from Meteor
892 907 mPeak1 = powerNet1.argmax() + mStart1
893 908 mPeakPower1 = powerNet1.max()
894 909 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
895 910 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
896 911 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
897 912 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
898 913 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
899 914 #Vectorize
900 915 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
901 916 meteorAux[7:11] = phaseDiffint[0:4]
902 917
903 918 #Rejection Criterions
904 919 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
905 920 meteorAux[-1] = 17
906 921 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
907 922 meteorAux[-1] = 1
908 923
909 924
910 925 else:
911 926 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
912 927 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
913 928 PowerSeries = 0
914 929
915 930 listMeteors1.append(meteorAux)
916 931 listPowerSeries.append(PowerSeries)
917 932 listVoltageSeries.append(meteorVolts1)
918 933
919 934 return listMeteors1, listPowerSeries, listVoltageSeries
920 935
921 936 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
922 937
923 938 threshError = 10
924 939 #Depending if it is 30 or 50 MHz
925 940 if frequency == 30e6:
926 941 timeLag = 45*10**-3
927 942 else:
928 943 timeLag = 15*10**-3
929 944 lag = numpy.ceil(timeLag/timeInterval)
930 945
931 946 listMeteors1 = []
932 947
933 948 for i in range(len(listMeteors)):
934 949 meteorPower = listPower[i]
935 950 meteorAux = listMeteors[i]
936 951
937 952 if meteorAux[-1] == 0:
938 953
939 954 try:
940 955 indmax = meteorPower.argmax()
941 956 indlag = indmax + lag
942 957
943 958 y = meteorPower[indlag:]
944 959 x = numpy.arange(0, y.size)*timeLag
945 960
946 961 #first guess
947 962 a = y[0]
948 963 tau = timeLag
949 964 #exponential fit
950 965 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
951 966 y1 = self.__exponential_function(x, *popt)
952 967 #error estimation
953 968 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
954 969
955 970 decayTime = popt[1]
956 971 riseTime = indmax*timeInterval
957 972 meteorAux[11:13] = [decayTime, error]
958 973
959 974 #Table items 7, 8 and 11
960 975 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
961 976 meteorAux[-1] = 7
962 977 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
963 978 meteorAux[-1] = 8
964 979 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
965 980 meteorAux[-1] = 11
966 981
967 982
968 983 except:
969 984 meteorAux[-1] = 11
970 985
971 986
972 987 listMeteors1.append(meteorAux)
973 988
974 989 return listMeteors1
975 990
976 991 #Exponential Function
977 992
978 993 def __exponential_function(self, x, a, tau):
979 994 y = a*numpy.exp(-x/tau)
980 995 return y
981 996
982 997 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
983 998
984 999 pairslist1 = list(pairslist)
985 1000 pairslist1.append((0,1))
986 1001 pairslist1.append((3,4))
987 1002 numPairs = len(pairslist1)
988 1003 #Time Lag
989 1004 timeLag = 45*10**-3
990 1005 c = 3e8
991 1006 lag = numpy.ceil(timeLag/timeInterval)
992 1007 freq = 30e6
993 1008
994 1009 listMeteors1 = []
995 1010
996 1011 for i in range(len(listMeteors)):
997 1012 meteorAux = listMeteors[i]
998 1013 if meteorAux[-1] == 0:
999 1014 mStart = listMeteors[i][1]
1000 1015 mPeak = listMeteors[i][2]
1001 1016 mLag = mPeak - mStart + lag
1002 1017
1003 1018 #get the volt data between the start and end times of the meteor
1004 1019 meteorVolts = listVolts[i]
1005 1020 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
1006 1021
1007 1022 #Get CCF
1008 1023 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
1009 1024
1010 1025 #Method 2
1011 1026 slopes = numpy.zeros(numPairs)
1012 1027 time = numpy.array([-2,-1,1,2])*timeInterval
1013 1028 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
1014 1029
1015 1030 #Correct phases
1016 1031 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
1017 1032 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
1018 1033
1019 1034 if indDer[0].shape[0] > 0:
1020 1035 for i in range(indDer[0].shape[0]):
1021 1036 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
1022 1037 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
1023 1038
1024 1039 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
1025 1040 for j in range(numPairs):
1026 1041 fit = stats.linregress(time, angAllCCF[j,:])
1027 1042 slopes[j] = fit[0]
1028 1043
1029 1044 #Remove Outlier
1030 1045 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
1031 1046 # slopes = numpy.delete(slopes,indOut)
1032 1047 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
1033 1048 # slopes = numpy.delete(slopes,indOut)
1034 1049
1035 1050 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
1036 1051 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
1037 1052 meteorAux[-2] = radialError
1038 1053 meteorAux[-3] = radialVelocity
1039 1054
1040 1055 #Setting Error
1041 1056 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
1042 1057 if numpy.abs(radialVelocity) > 200:
1043 1058 meteorAux[-1] = 15
1044 1059 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
1045 1060 elif radialError > radialStdThresh:
1046 1061 meteorAux[-1] = 12
1047 1062
1048 1063 listMeteors1.append(meteorAux)
1049 1064 return listMeteors1
1050 1065
1051 1066 def __setNewArrays(self, listMeteors, date, heiRang):
1052 1067
1053 1068 #New arrays
1054 1069 arrayMeteors = numpy.array(listMeteors)
1055 1070 arrayParameters = numpy.zeros((len(listMeteors), 14))
1056 1071
1057 1072 #Date inclusion
1058 date = re.findall(r'\((.*?)\)', date)
1059 date = date[0].split(',')
1060 date = map(int, date)
1061
1062 if len(date)<6:
1063 date.append(0)
1064
1065 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
1066 arrayDate = numpy.tile(date, (len(listMeteors), 1))
1073 # date = re.findall(r'\((.*?)\)', date)
1074 # date = date[0].split(',')
1075 # date = map(int, date)
1076 #
1077 # if len(date)<6:
1078 # date.append(0)
1079 #
1080 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
1081 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
1082 arrayDate = numpy.tile(date, (len(listMeteors)))
1067 1083
1068 1084 #Meteor array
1069 1085 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1070 1086 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1071 1087
1072 1088 #Parameters Array
1073 arrayParameters[:,:2] = arrayDate #Date
1074 arrayParameters[:,2] = heiRang[arrayMeteors[:,0].astype(int)] #Range
1075 arrayParameters[:,7:9] = arrayMeteors[:,-3:-1] #Radial velocity and its error
1076 arrayParameters[:,9:13] = arrayMeteors[:,7:11] #Phases
1089 arrayParameters[:,0] = arrayDate #Date
1090 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
1091 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
1092 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
1077 1093 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
1078 1094
1079 1095
1080 1096 return arrayParameters
1081 1097
1082 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1083
1084 arrayAOA = numpy.zeros((phases.shape[0],3))
1085 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1086
1087 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1088 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1089 arrayAOA[:,2] = cosDirError
1090
1091 azimuthAngle = arrayAOA[:,0]
1092 zenithAngle = arrayAOA[:,1]
1093
1094 #Setting Error
1095 #Number 3: AOA not fesible
1096 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1097 error[indInvalid] = 3
1098 #Number 4: Large difference in AOAs obtained from different antenna baselines
1099 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1100 error[indInvalid] = 4
1101 return arrayAOA, error
1102
1103 def __getDirectionCosines(self, arrayPhase, pairsList):
1104
1105 #Initializing some variables
1106 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1107 ang_aux = ang_aux.reshape(1,ang_aux.size)
1108
1109 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1110 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1111
1112
1113 for i in range(2):
1114 #First Estimation
1115 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1116 #Dealias
1117 indcsi = numpy.where(phi0_aux > numpy.pi)
1118 phi0_aux[indcsi] -= 2*numpy.pi
1119 indcsi = numpy.where(phi0_aux < -numpy.pi)
1120 phi0_aux[indcsi] += 2*numpy.pi
1121 #Direction Cosine 0
1122 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1123
1124 #Most-Accurate Second Estimation
1125 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1126 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1127 #Direction Cosine 1
1128 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1129
1130 #Searching the correct Direction Cosine
1131 cosdir0_aux = cosdir0[:,i]
1132 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1133 #Minimum Distance
1134 cosDiff = (cosdir1 - cosdir0_aux)**2
1135 indcos = cosDiff.argmin(axis = 1)
1136 #Saving Value obtained
1137 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1138
1139 return cosdir0, cosdir
1140
1141 def __calculateAOA(self, cosdir, azimuth):
1142 cosdirX = cosdir[:,0]
1143 cosdirY = cosdir[:,1]
1144
1145 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1146 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1147 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1148
1149 return angles
1150
1151 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1152
1153 Ramb = 375 #Ramb = c/(2*PRF)
1154 Re = 6371 #Earth Radius
1155 heights = numpy.zeros(Ranges.shape)
1156
1157 R_aux = numpy.array([0,1,2])*Ramb
1158 R_aux = R_aux.reshape(1,R_aux.size)
1159
1160 Ranges = Ranges.reshape(Ranges.size,1)
1161
1162 Ri = Ranges + R_aux
1163 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1164
1165 #Check if there is a height between 70 and 110 km
1166 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1167 ind_h = numpy.where(h_bool == 1)[0]
1168
1169 hCorr = hi[ind_h, :]
1170 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1171
1172 hCorr = hi[ind_hCorr]
1173 heights[ind_h] = hCorr
1174
1175 #Setting Error
1176 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1177 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1178
1179 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1180 error[indInvalid2] = 14
1181 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1182 error[indInvalid1] = 13
1183
1184 return heights, error
1098 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1099 #
1100 # arrayAOA = numpy.zeros((phases.shape[0],3))
1101 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1102 #
1103 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1104 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1105 # arrayAOA[:,2] = cosDirError
1106 #
1107 # azimuthAngle = arrayAOA[:,0]
1108 # zenithAngle = arrayAOA[:,1]
1109 #
1110 # #Setting Error
1111 # #Number 3: AOA not fesible
1112 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1113 # error[indInvalid] = 3
1114 # #Number 4: Large difference in AOAs obtained from different antenna baselines
1115 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1116 # error[indInvalid] = 4
1117 # return arrayAOA, error
1118 #
1119 # def __getDirectionCosines(self, arrayPhase, pairsList):
1120 #
1121 # #Initializing some variables
1122 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1123 # ang_aux = ang_aux.reshape(1,ang_aux.size)
1124 #
1125 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
1126 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1127 #
1128 #
1129 # for i in range(2):
1130 # #First Estimation
1131 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1132 # #Dealias
1133 # indcsi = numpy.where(phi0_aux > numpy.pi)
1134 # phi0_aux[indcsi] -= 2*numpy.pi
1135 # indcsi = numpy.where(phi0_aux < -numpy.pi)
1136 # phi0_aux[indcsi] += 2*numpy.pi
1137 # #Direction Cosine 0
1138 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1139 #
1140 # #Most-Accurate Second Estimation
1141 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1142 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1143 # #Direction Cosine 1
1144 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1145 #
1146 # #Searching the correct Direction Cosine
1147 # cosdir0_aux = cosdir0[:,i]
1148 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1149 # #Minimum Distance
1150 # cosDiff = (cosdir1 - cosdir0_aux)**2
1151 # indcos = cosDiff.argmin(axis = 1)
1152 # #Saving Value obtained
1153 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1154 #
1155 # return cosdir0, cosdir
1156 #
1157 # def __calculateAOA(self, cosdir, azimuth):
1158 # cosdirX = cosdir[:,0]
1159 # cosdirY = cosdir[:,1]
1160 #
1161 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1162 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1163 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1164 #
1165 # return angles
1166 #
1167 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1168 #
1169 # Ramb = 375 #Ramb = c/(2*PRF)
1170 # Re = 6371 #Earth Radius
1171 # heights = numpy.zeros(Ranges.shape)
1172 #
1173 # R_aux = numpy.array([0,1,2])*Ramb
1174 # R_aux = R_aux.reshape(1,R_aux.size)
1175 #
1176 # Ranges = Ranges.reshape(Ranges.size,1)
1177 #
1178 # Ri = Ranges + R_aux
1179 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1180 #
1181 # #Check if there is a height between 70 and 110 km
1182 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1183 # ind_h = numpy.where(h_bool == 1)[0]
1184 #
1185 # hCorr = hi[ind_h, :]
1186 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1187 #
1188 # hCorr = hi[ind_hCorr]
1189 # heights[ind_h] = hCorr
1190 #
1191 # #Setting Error
1192 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1193 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1194 #
1195 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1196 # error[indInvalid2] = 14
1197 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1198 # error[indInvalid1] = 13
1199 #
1200 # return heights, error
1185 1201
1186 1202 def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
1187 1203
1188 1204 '''
1189 1205 Function GetMoments()
1190 1206
1191 1207 Input:
1192 1208 Output:
1193 1209 Variables modified:
1194 1210 '''
1195 1211 if path != None:
1196 1212 sys.path.append(path)
1197 1213 self.dataOut.library = importlib.import_module(file)
1198 1214
1199 1215 #To be inserted as a parameter
1200 1216 groupArray = numpy.array(groupList)
1201 1217 # groupArray = numpy.array([[0,1],[2,3]])
1202 1218 self.dataOut.groupList = groupArray
1203 1219
1204 1220 nGroups = groupArray.shape[0]
1205 1221 nChannels = self.dataIn.nChannels
1206 1222 nHeights=self.dataIn.heightList.size
1207 1223
1208 1224 #Parameters Array
1209 1225 self.dataOut.data_param = None
1210 1226
1211 1227 #Set constants
1212 1228 constants = self.dataOut.library.setConstants(self.dataIn)
1213 1229 self.dataOut.constants = constants
1214 1230 M = self.dataIn.normFactor
1215 1231 N = self.dataIn.nFFTPoints
1216 1232 ippSeconds = self.dataIn.ippSeconds
1217 1233 K = self.dataIn.nIncohInt
1218 1234 pairsArray = numpy.array(self.dataIn.pairsList)
1219 1235
1220 1236 #List of possible combinations
1221 1237 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1222 1238 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1223 1239
1224 1240 if getSNR:
1225 1241 listChannels = groupArray.reshape((groupArray.size))
1226 1242 listChannels.sort()
1227 1243 noise = self.dataIn.getNoise()
1228 1244 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1229 1245
1230 1246 for i in range(nGroups):
1231 1247 coord = groupArray[i,:]
1232 1248
1233 1249 #Input data array
1234 1250 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1235 1251 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1236 1252
1237 1253 #Cross Spectra data array for Covariance Matrixes
1238 1254 ind = 0
1239 1255 for pairs in listComb:
1240 1256 pairsSel = numpy.array([coord[x],coord[y]])
1241 1257 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1242 1258 ind += 1
1243 1259 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1244 1260 dataCross = dataCross**2/K
1245 1261
1246 1262 for h in range(nHeights):
1247 1263 # print self.dataOut.heightList[h]
1248 1264
1249 1265 #Input
1250 1266 d = data[:,h]
1251 1267
1252 1268 #Covariance Matrix
1253 1269 D = numpy.diag(d**2/K)
1254 1270 ind = 0
1255 1271 for pairs in listComb:
1256 1272 #Coordinates in Covariance Matrix
1257 1273 x = pairs[0]
1258 1274 y = pairs[1]
1259 1275 #Channel Index
1260 1276 S12 = dataCross[ind,:,h]
1261 1277 D12 = numpy.diag(S12)
1262 1278 #Completing Covariance Matrix with Cross Spectras
1263 1279 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1264 1280 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1265 1281 ind += 1
1266 1282 Dinv=numpy.linalg.inv(D)
1267 1283 L=numpy.linalg.cholesky(Dinv)
1268 1284 LT=L.T
1269 1285
1270 1286 dp = numpy.dot(LT,d)
1271 1287
1272 1288 #Initial values
1273 1289 data_spc = self.dataIn.data_spc[coord,:,h]
1274 1290
1275 1291 if (h>0)and(error1[3]<5):
1276 1292 p0 = self.dataOut.data_param[i,:,h-1]
1277 1293 else:
1278 1294 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1279 1295
1280 1296 try:
1281 1297 #Least Squares
1282 1298 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1283 1299 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1284 1300 #Chi square error
1285 1301 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1286 1302 #Error with Jacobian
1287 1303 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1288 1304 except:
1289 1305 minp = p0*numpy.nan
1290 1306 error0 = numpy.nan
1291 1307 error1 = p0*numpy.nan
1292 1308
1293 1309 #Save
1294 1310 if self.dataOut.data_param == None:
1295 1311 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1296 1312 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1297 1313
1298 1314 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1299 1315 self.dataOut.data_param[i,:,h] = minp
1300 1316 return
1301 1317
1302 1318 def __residFunction(self, p, dp, LT, constants):
1303 1319
1304 1320 fm = self.dataOut.library.modelFunction(p, constants)
1305 1321 fmp=numpy.dot(LT,fm)
1306 1322
1307 1323 return dp-fmp
1308 1324
1309 1325 def __getSNR(self, z, noise):
1310 1326
1311 1327 avg = numpy.average(z, axis=1)
1312 1328 SNR = (avg.T-noise)/noise
1313 1329 SNR = SNR.T
1314 1330 return SNR
1315 1331
1316 1332 def __chisq(p,chindex,hindex):
1317 1333 #similar to Resid but calculates CHI**2
1318 1334 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1319 1335 dp=numpy.dot(LT,d)
1320 1336 fmp=numpy.dot(LT,fm)
1321 1337 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1322 1338 return chisq
1323 1339
1324
1325 1340
1326 1341 class WindProfiler(Operation):
1327 1342
1328 1343 __isConfig = False
1329 1344
1330 1345 __initime = None
1331 1346 __lastdatatime = None
1332 1347 __integrationtime = None
1333 1348
1334 1349 __buffer = None
1335 1350
1336 1351 __dataReady = False
1337 1352
1338 1353 __firstdata = None
1339 1354
1340 1355 n = None
1341 1356
1342 1357 def __init__(self):
1343 1358 Operation.__init__(self)
1344 1359
1345 1360 def __calculateCosDir(self, elev, azim):
1346 1361 zen = (90 - elev)*numpy.pi/180
1347 1362 azim = azim*numpy.pi/180
1348 1363 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1349 1364 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1350 1365
1351 1366 signX = numpy.sign(numpy.cos(azim))
1352 1367 signY = numpy.sign(numpy.sin(azim))
1353 1368
1354 1369 cosDirX = numpy.copysign(cosDirX, signX)
1355 1370 cosDirY = numpy.copysign(cosDirY, signY)
1356 1371 return cosDirX, cosDirY
1357 1372
1358 1373 def __calculateAngles(self, theta_x, theta_y, azimuth):
1359 1374
1360 1375 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1361 1376 zenith_arr = numpy.arccos(dir_cosw)
1362 1377 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1363 1378
1364 1379 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1365 1380 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1366 1381
1367 1382 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1368 1383
1369 1384 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1370 1385
1371 1386 #
1372 1387 if horOnly:
1373 1388 A = numpy.c_[dir_cosu,dir_cosv]
1374 1389 else:
1375 1390 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1376 1391 A = numpy.asmatrix(A)
1377 1392 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1378 1393
1379 1394 return A1
1380 1395
1381 1396 def __correctValues(self, heiRang, phi, velRadial, SNR):
1382 1397 listPhi = phi.tolist()
1383 1398 maxid = listPhi.index(max(listPhi))
1384 1399 minid = listPhi.index(min(listPhi))
1385 1400
1386 1401 rango = range(len(phi))
1387 1402 # rango = numpy.delete(rango,maxid)
1388 1403
1389 1404 heiRang1 = heiRang*math.cos(phi[maxid])
1390 1405 heiRangAux = heiRang*math.cos(phi[minid])
1391 1406 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1392 1407 heiRang1 = numpy.delete(heiRang1,indOut)
1393 1408
1394 1409 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1395 1410 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1396 1411
1397 1412 for i in rango:
1398 1413 x = heiRang*math.cos(phi[i])
1399 1414 y1 = velRadial[i,:]
1400 1415 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1401 1416
1402 1417 x1 = heiRang1
1403 1418 y11 = f1(x1)
1404 1419
1405 1420 y2 = SNR[i,:]
1406 1421 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1407 1422 y21 = f2(x1)
1408 1423
1409 1424 velRadial1[i,:] = y11
1410 1425 SNR1[i,:] = y21
1411 1426
1412 1427 return heiRang1, velRadial1, SNR1
1413 1428
1414 1429 def __calculateVelUVW(self, A, velRadial):
1415 1430
1416 1431 #Operacion Matricial
1417 1432 # velUVW = numpy.zeros((velRadial.shape[1],3))
1418 1433 # for ind in range(velRadial.shape[1]):
1419 1434 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1420 1435 # velUVW = velUVW.transpose()
1421 1436 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1422 1437 velUVW[:,:] = numpy.dot(A,velRadial)
1423 1438
1424 1439
1425 1440 return velUVW
1426 1441
1427 1442 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1428 1443 """
1429 1444 Function that implements Doppler Beam Swinging (DBS) technique.
1430 1445
1431 1446 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1432 1447 Direction correction (if necessary), Ranges and SNR
1433 1448
1434 1449 Output: Winds estimation (Zonal, Meridional and Vertical)
1435 1450
1436 1451 Parameters affected: Winds, height range, SNR
1437 1452 """
1438 1453 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1439 1454 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1440 1455 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1441 1456
1442 1457 #Calculo de Componentes de la velocidad con DBS
1443 1458 winds = self.__calculateVelUVW(A,velRadial1)
1444 1459
1445 1460 return winds, heiRang1, SNR1
1446 1461
1447 1462 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1448 1463
1449 1464 posx = numpy.asarray(posx)
1450 1465 posy = numpy.asarray(posy)
1451 1466
1452 1467 #Rotacion Inversa para alinear con el azimuth
1453 1468 if azimuth!= None:
1454 1469 azimuth = azimuth*math.pi/180
1455 1470 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1456 1471 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1457 1472 else:
1458 1473 posx1 = posx
1459 1474 posy1 = posy
1460 1475
1461 1476 #Calculo de Distancias
1462 1477 distx = numpy.zeros(pairsCrossCorr.size)
1463 1478 disty = numpy.zeros(pairsCrossCorr.size)
1464 1479 dist = numpy.zeros(pairsCrossCorr.size)
1465 1480 ang = numpy.zeros(pairsCrossCorr.size)
1466 1481
1467 1482 for i in range(pairsCrossCorr.size):
1468 1483 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1469 1484 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1470 1485 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1471 1486 ang[i] = numpy.arctan2(disty[i],distx[i])
1472 1487 #Calculo de Matrices
1473 1488 nPairs = len(pairs)
1474 1489 ang1 = numpy.zeros((nPairs, 2, 1))
1475 1490 dist1 = numpy.zeros((nPairs, 2, 1))
1476 1491
1477 1492 for j in range(nPairs):
1478 1493 dist1[j,0,0] = dist[pairs[j][0]]
1479 1494 dist1[j,1,0] = dist[pairs[j][1]]
1480 1495 ang1[j,0,0] = ang[pairs[j][0]]
1481 1496 ang1[j,1,0] = ang[pairs[j][1]]
1482 1497
1483 1498 return distx,disty, dist1,ang1
1484 1499
1485 1500 def __calculateVelVer(self, phase, lagTRange, _lambda):
1486 1501
1487 1502 Ts = lagTRange[1] - lagTRange[0]
1488 1503 velW = -_lambda*phase/(4*math.pi*Ts)
1489 1504
1490 1505 return velW
1491 1506
1492 1507 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1493 1508 nPairs = tau1.shape[0]
1494 1509 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1495 1510
1496 1511 angCos = numpy.cos(ang)
1497 1512 angSin = numpy.sin(ang)
1498 1513
1499 1514 vel0 = dist*tau1/(2*tau2**2)
1500 1515 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1501 1516 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1502 1517
1503 1518 ind = numpy.where(numpy.isinf(vel))
1504 1519 vel[ind] = numpy.nan
1505 1520
1506 1521 return vel
1507 1522
1508 1523 def __getPairsAutoCorr(self, pairsList, nChannels):
1509 1524
1510 1525 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1511 1526
1512 1527 for l in range(len(pairsList)):
1513 1528 firstChannel = pairsList[l][0]
1514 1529 secondChannel = pairsList[l][1]
1515 1530
1516 1531 #Obteniendo pares de Autocorrelacion
1517 1532 if firstChannel == secondChannel:
1518 1533 pairsAutoCorr[firstChannel] = int(l)
1519 1534
1520 1535 pairsAutoCorr = pairsAutoCorr.astype(int)
1521 1536
1522 1537 pairsCrossCorr = range(len(pairsList))
1523 1538 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1524 1539
1525 1540 return pairsAutoCorr, pairsCrossCorr
1526 1541
1527 1542 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1528 1543 """
1529 1544 Function that implements Spaced Antenna (SA) technique.
1530 1545
1531 1546 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1532 1547 Direction correction (if necessary), Ranges and SNR
1533 1548
1534 1549 Output: Winds estimation (Zonal, Meridional and Vertical)
1535 1550
1536 1551 Parameters affected: Winds
1537 1552 """
1538 1553 #Cross Correlation pairs obtained
1539 1554 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1540 1555 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1541 1556 pairsSelArray = numpy.array(pairsSelected)
1542 1557 pairs = []
1543 1558
1544 1559 #Wind estimation pairs obtained
1545 1560 for i in range(pairsSelArray.shape[0]/2):
1546 1561 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1547 1562 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1548 1563 pairs.append((ind1,ind2))
1549 1564
1550 1565 indtau = tau.shape[0]/2
1551 1566 tau1 = tau[:indtau,:]
1552 1567 tau2 = tau[indtau:-1,:]
1553 1568 tau1 = tau1[pairs,:]
1554 1569 tau2 = tau2[pairs,:]
1555 1570 phase1 = tau[-1,:]
1556 1571
1557 1572 #---------------------------------------------------------------------
1558 1573 #Metodo Directo
1559 1574 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1560 1575 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1561 1576 winds = stats.nanmean(winds, axis=0)
1562 1577 #---------------------------------------------------------------------
1563 1578 #Metodo General
1564 1579 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1565 1580 # #Calculo Coeficientes de Funcion de Correlacion
1566 1581 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1567 1582 # #Calculo de Velocidades
1568 1583 # winds = self.calculateVelUV(F,G,A,B,H)
1569 1584
1570 1585 #---------------------------------------------------------------------
1571 1586 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1572 1587 winds = correctFactor*winds
1573 1588 return winds
1574 1589
1575 1590 def __checkTime(self, currentTime, paramInterval, outputInterval):
1576 1591
1577 1592 dataTime = currentTime + paramInterval
1578 1593 deltaTime = dataTime - self.__initime
1579 1594
1580 1595 if deltaTime >= outputInterval or deltaTime < 0:
1581 1596 self.__dataReady = True
1582 1597 return
1583 1598
1584 1599 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1585 1600 '''
1586 1601 Function that implements winds estimation technique with detected meteors.
1587 1602
1588 1603 Input: Detected meteors, Minimum meteor quantity to wind estimation
1589 1604
1590 1605 Output: Winds estimation (Zonal and Meridional)
1591 1606
1592 1607 Parameters affected: Winds
1593 1608 '''
1594 1609 # print arrayMeteor.shape
1595 1610 #Settings
1596 1611 nInt = (heightMax - heightMin)/2
1597 1612 # print nInt
1598 1613 nInt = int(nInt)
1599 1614 # print nInt
1600 1615 winds = numpy.zeros((2,nInt))*numpy.nan
1601 1616
1602 1617 #Filter errors
1603 error = numpy.where(arrayMeteor[0,:,-1] == 0)[0]
1604 finalMeteor = arrayMeteor[0,error,:]
1618 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1619 finalMeteor = arrayMeteor[error,:]
1605 1620
1606 1621 #Meteor Histogram
1607 finalHeights = finalMeteor[:,3]
1622 finalHeights = finalMeteor[:,2]
1608 1623 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1609 1624 nMeteorsPerI = hist[0]
1610 1625 heightPerI = hist[1]
1611 1626
1612 1627 #Sort of meteors
1613 1628 indSort = finalHeights.argsort()
1614 1629 finalMeteor2 = finalMeteor[indSort,:]
1615 1630
1616 1631 # Calculating winds
1617 1632 ind1 = 0
1618 1633 ind2 = 0
1619 1634
1620 1635 for i in range(nInt):
1621 1636 nMet = nMeteorsPerI[i]
1622 1637 ind1 = ind2
1623 1638 ind2 = ind1 + nMet
1624 1639
1625 1640 meteorAux = finalMeteor2[ind1:ind2,:]
1626 1641
1627 1642 if meteorAux.shape[0] >= meteorThresh:
1628 vel = meteorAux[:, 7]
1629 zen = meteorAux[:, 5]*numpy.pi/180
1630 azim = meteorAux[:, 4]*numpy.pi/180
1643 vel = meteorAux[:, 6]
1644 zen = meteorAux[:, 4]*numpy.pi/180
1645 azim = meteorAux[:, 3]*numpy.pi/180
1631 1646
1632 1647 n = numpy.cos(zen)
1633 1648 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1634 1649 # l = m*numpy.tan(azim)
1635 1650 l = numpy.sin(zen)*numpy.sin(azim)
1636 1651 m = numpy.sin(zen)*numpy.cos(azim)
1637 1652
1638 1653 A = numpy.vstack((l, m)).transpose()
1639 1654 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1640 1655 windsAux = numpy.dot(A1, vel)
1641 1656
1642 1657 winds[0,i] = windsAux[0]
1643 1658 winds[1,i] = windsAux[1]
1644 1659
1645 1660 return winds, heightPerI[:-1]
1646 1661
1647 1662 def run(self, dataOut, technique, **kwargs):
1648 1663
1649 1664 param = dataOut.data_param
1650 1665 if dataOut.abscissaList != None:
1651 1666 absc = dataOut.abscissaList[:-1]
1652 1667 noise = dataOut.noise
1653 1668 heightList = dataOut.heightList
1654 1669 SNR = dataOut.data_SNR
1655 1670
1656 1671 if technique == 'DBS':
1657 1672
1658 1673 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1659 1674 theta_x = numpy.array(kwargs['dirCosx'])
1660 1675 theta_y = numpy.array(kwargs['dirCosy'])
1661 1676 else:
1662 1677 elev = numpy.array(kwargs['elevation'])
1663 1678 azim = numpy.array(kwargs['azimuth'])
1664 1679 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1665 1680 azimuth = kwargs['correctAzimuth']
1666 1681 if kwargs.has_key('horizontalOnly'):
1667 1682 horizontalOnly = kwargs['horizontalOnly']
1668 1683 else: horizontalOnly = False
1669 1684 if kwargs.has_key('correctFactor'):
1670 1685 correctFactor = kwargs['correctFactor']
1671 1686 else: correctFactor = 1
1672 1687 if kwargs.has_key('channelList'):
1673 1688 channelList = kwargs['channelList']
1674 1689 if len(channelList) == 2:
1675 1690 horizontalOnly = True
1676 1691 arrayChannel = numpy.array(channelList)
1677 1692 param = param[arrayChannel,:,:]
1678 1693 theta_x = theta_x[arrayChannel]
1679 1694 theta_y = theta_y[arrayChannel]
1680 1695
1681 1696 velRadial0 = param[:,1,:] #Radial velocity
1682 1697 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function
1683 1698 dataOut.utctimeInit = dataOut.utctime
1684 1699 dataOut.outputInterval = dataOut.timeInterval
1685 1700
1686 1701 elif technique == 'SA':
1687 1702
1688 1703 #Parameters
1689 1704 position_x = kwargs['positionX']
1690 1705 position_y = kwargs['positionY']
1691 1706 azimuth = kwargs['azimuth']
1692 1707
1693 1708 if kwargs.has_key('crosspairsList'):
1694 1709 pairs = kwargs['crosspairsList']
1695 1710 else:
1696 1711 pairs = None
1697 1712
1698 1713 if kwargs.has_key('correctFactor'):
1699 1714 correctFactor = kwargs['correctFactor']
1700 1715 else:
1701 1716 correctFactor = 1
1702 1717
1703 1718 tau = dataOut.data_param
1704 1719 _lambda = dataOut.C/dataOut.frequency
1705 1720 pairsList = dataOut.groupList
1706 1721 nChannels = dataOut.nChannels
1707 1722
1708 1723 dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1709 1724 dataOut.utctimeInit = dataOut.utctime
1710 1725 dataOut.outputInterval = dataOut.timeInterval
1711 1726
1712 1727 elif technique == 'Meteors':
1713 1728 dataOut.flagNoData = True
1714 1729 self.__dataReady = False
1715 1730
1716 1731 if kwargs.has_key('nHours'):
1717 1732 nHours = kwargs['nHours']
1718 1733 else:
1719 1734 nHours = 1
1720 1735
1721 1736 if kwargs.has_key('meteorsPerBin'):
1722 1737 meteorThresh = kwargs['meteorsPerBin']
1723 1738 else:
1724 1739 meteorThresh = 6
1725 1740
1726 1741 if kwargs.has_key('hmin'):
1727 1742 hmin = kwargs['hmin']
1728 1743 else: hmin = 70
1729 1744 if kwargs.has_key('hmax'):
1730 1745 hmax = kwargs['hmax']
1731 1746 else: hmax = 110
1732 1747
1733 1748 dataOut.outputInterval = nHours*3600
1734 1749
1735 1750 if self.__isConfig == False:
1736 1751 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1737 1752 #Get Initial LTC time
1738 1753 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1739 1754 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1740 1755
1741 1756 self.__isConfig = True
1742 1757
1743 1758 if self.__buffer == None:
1744 1759 self.__buffer = dataOut.data_param
1745 1760 self.__firstdata = copy.copy(dataOut)
1746 1761
1747 1762 else:
1748 self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
1763 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1749 1764
1750 1765 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1751 1766
1752 1767 if self.__dataReady:
1753 1768 dataOut.utctimeInit = self.__initime
1754 1769
1755 1770 self.__initime += dataOut.outputInterval #to erase time offset
1756 1771
1757 1772 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1758 1773 dataOut.flagNoData = False
1759 1774 self.__buffer = None
1760 1775
1761 1776 return
1762 1777
1763 1778 class EWDriftsEstimation(Operation):
1764 1779
1765 1780
1766 1781 def __init__(self):
1767 1782 Operation.__init__(self)
1768 1783
1769 1784 def __correctValues(self, heiRang, phi, velRadial, SNR):
1770 1785 listPhi = phi.tolist()
1771 1786 maxid = listPhi.index(max(listPhi))
1772 1787 minid = listPhi.index(min(listPhi))
1773 1788
1774 1789 rango = range(len(phi))
1775 1790 # rango = numpy.delete(rango,maxid)
1776 1791
1777 1792 heiRang1 = heiRang*math.cos(phi[maxid])
1778 1793 heiRangAux = heiRang*math.cos(phi[minid])
1779 1794 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1780 1795 heiRang1 = numpy.delete(heiRang1,indOut)
1781 1796
1782 1797 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1783 1798 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1784 1799
1785 1800 for i in rango:
1786 1801 x = heiRang*math.cos(phi[i])
1787 1802 y1 = velRadial[i,:]
1788 1803 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1789 1804
1790 1805 x1 = heiRang1
1791 1806 y11 = f1(x1)
1792 1807
1793 1808 y2 = SNR[i,:]
1794 1809 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1795 1810 y21 = f2(x1)
1796 1811
1797 1812 velRadial1[i,:] = y11
1798 1813 SNR1[i,:] = y21
1799 1814
1800 1815 return heiRang1, velRadial1, SNR1
1801 1816
1802 1817 def run(self, dataOut, zenith, zenithCorrection):
1803 1818 heiRang = dataOut.heightList
1804 1819 velRadial = dataOut.data_param[:,3,:]
1805 1820 SNR = dataOut.data_SNR
1806 1821
1807 1822 zenith = numpy.array(zenith)
1808 1823 zenith -= zenithCorrection
1809 1824 zenith *= numpy.pi/180
1810 1825
1811 1826 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1812 1827
1813 1828 alp = zenith[0]
1814 1829 bet = zenith[1]
1815 1830
1816 1831 w_w = velRadial1[0,:]
1817 1832 w_e = velRadial1[1,:]
1818 1833
1819 1834 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1820 1835 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1821 1836
1822 1837 winds = numpy.vstack((u,w))
1823 1838
1824 1839 dataOut.heightList = heiRang1
1825 1840 dataOut.data_output = winds
1826 1841 dataOut.data_SNR = SNR1
1827 1842
1828 1843 dataOut.utctimeInit = dataOut.utctime
1829 1844 dataOut.outputInterval = dataOut.timeInterval
1830 1845 return
1831 1846
1832 1847 class PhaseCalibration(Operation):
1833 1848
1834 1849 __buffer = None
1835 1850
1836 1851 __initime = None
1837 1852
1838 1853 __dataReady = False
1839 1854
1840 1855 __isConfig = False
1841 1856
1842 1857 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
1843 1858
1844 1859 dataTime = currentTime + paramInterval
1845 1860 deltaTime = dataTime - initTime
1846 1861
1847 1862 if deltaTime >= outputInterval or deltaTime < 0:
1848 1863 return True
1849 1864
1850 1865 return False
1851 1866
1852 1867 def __getGammas(self, pairs, k, d, phases):
1853 1868 gammas = numpy.zeros(2)
1854 1869
1855 1870 for i in range(len(pairs)):
1856 1871
1857 1872 pairi = pairs[i]
1858 1873
1859 1874 #Calculating gamma
1860 1875 jdcos = phases[:,pairi[1]]/(k*d[pairi[1]])
1861 1876 jgamma = numpy.angle(numpy.exp(1j*(k*d[pairi[0]]*jdcos - phases[:,pairi[0]])))
1862 1877
1863 1878 #Revised distribution
1864 1879 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
1865 1880
1866 1881 #Histogram
1867 1882 nBins = 64.0
1868 1883 rmin = -0.5*numpy.pi
1869 1884 rmax = 0.5*numpy.pi
1870 1885 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
1871 1886
1872 1887 meteorsY = phaseHisto[0]
1873 1888 phasesX = phaseHisto[1][:-1]
1874 1889 width = phasesX[1] - phasesX[0]
1875 1890 phasesX += width/2
1876 1891
1877 1892 #Gaussian aproximation
1878 1893 bpeak = meteorsY.argmax()
1879 1894 peak = meteorsY.max()
1880 1895 jmin = bpeak - 5
1881 1896 jmax = bpeak + 5 + 1
1882 1897
1883 1898 if jmin<0:
1884 1899 jmin = 0
1885 1900 jmax = 6
1886 1901 elif jmax > meteorsY.size:
1887 1902 jmin = meteorsY.size - 6
1888 1903 jmax = meteorsY.size
1889 1904
1890 1905 x0 = numpy.array([peak,bpeak,50])
1891 1906 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
1892 1907
1893 1908 #Gammas
1894 1909 gammas[i] = coeff[0][1]
1895 1910
1896 1911 return gammas
1897 1912
1898 1913 def __residualFunction(self, coeffs, y, t):
1899 1914
1900 1915 return y - self.__gauss_function(t, coeffs)
1901 1916
1902 1917 def __gauss_function(self, t, coeffs):
1903 1918
1904 1919 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
1905 1920
1906 1921 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
1907 1922 meteorOps = MeteorOperations()
1908 1923 nchan = 4
1909 1924 pairx = pairsList[0]
1910 1925 pairy = pairsList[1]
1911 1926 center_xangle = 0
1912 1927 center_yangle = 0
1913 1928 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
1914 1929 ntimes = len(range_angle)
1915 1930
1916 1931 nstepsx = 20.0
1917 1932 nstepsy = 20.0
1918 1933
1919 1934 for iz in range(ntimes):
1920 1935 min_xangle = -range_angle[iz]/2 + center_xangle
1921 1936 max_xangle = range_angle[iz]/2 + center_xangle
1922 1937 min_yangle = -range_angle[iz]/2 + center_yangle
1923 1938 max_yangle = range_angle[iz]/2 + center_yangle
1924 1939
1925 1940 inc_x = (max_xangle-min_xangle)/nstepsx
1926 1941 inc_y = (max_yangle-min_yangle)/nstepsy
1927 1942
1928 1943 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
1929 1944 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
1930 1945 penalty = numpy.zeros((nstepsx,nstepsy))
1931 1946 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
1932 1947 jph = numpy.zeros(nchan)
1933 1948
1934 1949 # Iterations looking for the offset
1935 1950 for iy in range(int(nstepsy)):
1936 1951 for ix in range(int(nstepsx)):
1937 1952 jph[pairy[1]] = alpha_y[iy]
1938 1953 jph[pairy[0]] = -gammas[1] + alpha_y[iy]*d[pairy[0]]/d[pairy[1]]
1939 1954
1940 1955 jph[pairx[1]] = alpha_x[ix]
1941 1956 jph[pairx[0]] = -gammas[0] + alpha_x[ix]*d[pairx[0]]/d[pairx[1]]
1942 1957
1943 1958 jph_array[:,ix,iy] = jph
1944 1959
1945 1960 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, jph)
1946 1961 error = meteorsArray1[:,-1]
1947 1962 ind1 = numpy.where(error==0)[0]
1948 1963 penalty[ix,iy] = ind1.size
1949 1964
1950 1965 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
1951 1966 phOffset = jph_array[:,i,j]
1952 1967
1953 1968 center_xangle = phOffset[pairx[1]]
1954 1969 center_yangle = phOffset[pairy[1]]
1955 1970
1956 1971 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
1957 1972 phOffset = phOffset*180/numpy.pi
1958 1973 return phOffset
1959 1974
1960 1975
1961 1976 def run(self, dataOut, pairs, distances, hmin, hmax, nHours = None):
1962 1977
1963 1978 dataOut.flagNoData = True
1964 1979 self.__dataReady = False
1965 1980
1966 1981 if nHours == None:
1967 1982 nHours = 1
1968 1983
1969 1984 dataOut.outputInterval = nHours*3600
1970 1985
1971 1986 if self.__isConfig == False:
1972 1987 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1973 1988 #Get Initial LTC time
1974 1989 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1975 1990 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1976 1991
1977 1992 self.__isConfig = True
1978 1993
1979 1994 if self.__buffer == None:
1980 1995 self.__buffer = dataOut.data_param.copy()
1981 1996
1982 1997 else:
1983 1998 self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
1984 1999
1985 2000 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1986 2001
1987 2002 if self.__dataReady:
1988 2003 dataOut.utctimeInit = self.__initime
1989 2004 self.__initime += dataOut.outputInterval #to erase time offset
1990 2005
1991 2006 freq = dataOut.frequency
1992 2007 c = dataOut.C #m/s
1993 2008 lamb = c/freq
1994 2009 k = 2*numpy.pi/lamb
1995 2010 azimuth = 0
1996 2011 h = (hmin, hmax)
1997 2012 pairsList = ((0,3),(1,2))
1998 2013
1999 meteorsArray = self.__buffer[0,:,:]
2014 meteorsArray = self.__buffer
2000 2015 error = meteorsArray[:,-1]
2001 2016 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
2002 2017 ind1 = numpy.where(boolError)[0]
2003 2018 meteorsArray = meteorsArray[ind1,:]
2004 2019 meteorsArray[:,-1] = 0
2005 phases = meteorsArray[:,9:13]
2020 phases = meteorsArray[:,8:12]
2006 2021
2007 2022 #Calculate Gammas
2008 2023 gammas = self.__getGammas(pairs, k, distances, phases)
2009 2024 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2010 2025 #Calculate Phases
2011 2026 phasesOff = self.__getPhases(azimuth, h, pairsList, distances, gammas, meteorsArray)
2012 2027 phasesOff = phasesOff.reshape((1,phasesOff.size))
2013 2028 dataOut.data_output = -phasesOff
2014 2029 dataOut.flagNoData = False
2015 2030 self.__buffer = None
2016 2031
2017 2032
2018 2033 return
2019 2034
2020 2035 class MeteorOperations():
2021 2036
2022 2037 def __init__(self):
2023 2038
2024 2039 return
2025 2040
2026 2041 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, jph):
2027 2042
2028 2043 arrayParameters = arrayParameters0.copy()
2029 2044 hmin = h[0]
2030 2045 hmax = h[1]
2031 2046
2032 2047 #Calculate AOA (Error N 3, 4)
2033 2048 #JONES ET AL. 1998
2034 2049 AOAthresh = numpy.pi/8
2035 2050 error = arrayParameters[:,-1]
2036 phases = -arrayParameters[:,9:13] + jph
2037 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
2051 phases = -arrayParameters[:,8:12] + jph
2052 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
2038 2053
2039 2054 #Calculate Heights (Error N 13 and 14)
2040 2055 error = arrayParameters[:,-1]
2041 Ranges = arrayParameters[:,2]
2042 zenith = arrayParameters[:,5]
2043 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2056 Ranges = arrayParameters[:,1]
2057 zenith = arrayParameters[:,4]
2058 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2044 2059
2045 2060 #----------------------- Get Final data ------------------------------------
2046 2061 # error = arrayParameters[:,-1]
2047 2062 # ind1 = numpy.where(error==0)[0]
2048 2063 # arrayParameters = arrayParameters[ind1,:]
2049 2064
2050 2065 return arrayParameters
2051 2066
2052 2067 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2053 2068
2054 2069 arrayAOA = numpy.zeros((phases.shape[0],3))
2055 2070 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2056 2071
2057 2072 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2058 2073 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2059 2074 arrayAOA[:,2] = cosDirError
2060 2075
2061 2076 azimuthAngle = arrayAOA[:,0]
2062 2077 zenithAngle = arrayAOA[:,1]
2063 2078
2064 2079 #Setting Error
2065 2080 #Number 3: AOA not fesible
2066 2081 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2067 2082 error[indInvalid] = 3
2068 2083 #Number 4: Large difference in AOAs obtained from different antenna baselines
2069 2084 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2070 2085 error[indInvalid] = 4
2071 2086 return arrayAOA, error
2072 2087
2073 2088 def __getDirectionCosines(self, arrayPhase, pairsList):
2074 2089
2075 2090 #Initializing some variables
2076 2091 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2077 2092 ang_aux = ang_aux.reshape(1,ang_aux.size)
2078 2093
2079 2094 cosdir = numpy.zeros((arrayPhase.shape[0],2))
2080 2095 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2081 2096
2082 2097
2083 2098 for i in range(2):
2084 2099 #First Estimation
2085 2100 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2086 2101 #Dealias
2087 2102 indcsi = numpy.where(phi0_aux > numpy.pi)
2088 2103 phi0_aux[indcsi] -= 2*numpy.pi
2089 2104 indcsi = numpy.where(phi0_aux < -numpy.pi)
2090 2105 phi0_aux[indcsi] += 2*numpy.pi
2091 2106 #Direction Cosine 0
2092 2107 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2093 2108
2094 2109 #Most-Accurate Second Estimation
2095 2110 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2096 2111 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2097 2112 #Direction Cosine 1
2098 2113 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2099 2114
2100 2115 #Searching the correct Direction Cosine
2101 2116 cosdir0_aux = cosdir0[:,i]
2102 2117 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2103 2118 #Minimum Distance
2104 2119 cosDiff = (cosdir1 - cosdir0_aux)**2
2105 2120 indcos = cosDiff.argmin(axis = 1)
2106 2121 #Saving Value obtained
2107 2122 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2108 2123
2109 2124 return cosdir0, cosdir
2110 2125
2111 2126 def __calculateAOA(self, cosdir, azimuth):
2112 2127 cosdirX = cosdir[:,0]
2113 2128 cosdirY = cosdir[:,1]
2114 2129
2115 2130 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2116 2131 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2117 2132 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2118 2133
2119 2134 return angles
2120 2135
2121 2136 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2122 2137
2123 2138 Ramb = 375 #Ramb = c/(2*PRF)
2124 2139 Re = 6371 #Earth Radius
2125 2140 heights = numpy.zeros(Ranges.shape)
2126 2141
2127 2142 R_aux = numpy.array([0,1,2])*Ramb
2128 2143 R_aux = R_aux.reshape(1,R_aux.size)
2129 2144
2130 2145 Ranges = Ranges.reshape(Ranges.size,1)
2131 2146
2132 2147 Ri = Ranges + R_aux
2133 2148 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2134 2149
2135 2150 #Check if there is a height between 70 and 110 km
2136 2151 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2137 2152 ind_h = numpy.where(h_bool == 1)[0]
2138 2153
2139 2154 hCorr = hi[ind_h, :]
2140 2155 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2141 2156
2142 2157 hCorr = hi[ind_hCorr]
2143 2158 heights[ind_h] = hCorr
2144 2159
2145 2160 #Setting Error
2146 2161 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2147 2162 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2148 2163
2149 2164 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2150 2165 error[indInvalid2] = 14
2151 2166 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2152 2167 error[indInvalid1] = 13
2153 2168
2154 2169 return heights, error No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now