##// END OF EJS Templates
schain-jc
RSS
Fix all PlotData, add SpectraMean, CrossSpectra plots, now Parameters extends Spectra fix bugs in ParametersProc
jespinoza -
r922:d680543828ae
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,1216 +1,1218
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12 from schainpy import cSchain
13 13
14 14
15 15 def getNumpyDtype(dataTypeCode):
16 16
17 17 if dataTypeCode == 0:
18 18 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
19 19 elif dataTypeCode == 1:
20 20 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
21 21 elif dataTypeCode == 2:
22 22 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
23 23 elif dataTypeCode == 3:
24 24 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
25 25 elif dataTypeCode == 4:
26 26 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
27 27 elif dataTypeCode == 5:
28 28 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
29 29 else:
30 30 raise ValueError, 'dataTypeCode was not defined'
31 31
32 32 return numpyDtype
33 33
34 34 def getDataTypeCode(numpyDtype):
35 35
36 36 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
37 37 datatype = 0
38 38 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
39 39 datatype = 1
40 40 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
41 41 datatype = 2
42 42 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
43 43 datatype = 3
44 44 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
45 45 datatype = 4
46 46 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
47 47 datatype = 5
48 48 else:
49 49 datatype = None
50 50
51 51 return datatype
52 52
53 53 def hildebrand_sekhon(data, navg):
54 54 """
55 55 This method is for the objective determination of the noise level in Doppler spectra. This
56 56 implementation technique is based on the fact that the standard deviation of the spectral
57 57 densities is equal to the mean spectral density for white Gaussian noise
58 58
59 59 Inputs:
60 60 Data : heights
61 61 navg : numbers of averages
62 62
63 63 Return:
64 64 -1 : any error
65 65 anoise : noise's level
66 66 """
67 67
68 68 sortdata = numpy.sort(data,axis=None)
69 69 # lenOfData = len(sortdata)
70 70 # nums_min = lenOfData*0.2
71 71 #
72 72 # if nums_min <= 5:
73 73 # nums_min = 5
74 74 #
75 75 # sump = 0.
76 76 #
77 77 # sumq = 0.
78 78 #
79 79 # j = 0
80 80 #
81 81 # cont = 1
82 82 #
83 83 # while((cont==1)and(j<lenOfData)):
84 84 #
85 85 # sump += sortdata[j]
86 86 #
87 87 # sumq += sortdata[j]**2
88 88 #
89 89 # if j > nums_min:
90 90 # rtest = float(j)/(j-1) + 1.0/navg
91 91 # if ((sumq*j) > (rtest*sump**2)):
92 92 # j = j - 1
93 93 # sump = sump - sortdata[j]
94 94 # sumq = sumq - sortdata[j]**2
95 95 # cont = 0
96 96 #
97 97 # j += 1
98 98 #
99 99 # lnoise = sump /j
100 100 #
101 101 # return lnoise
102 102
103 103 return cSchain.hildebrand_sekhon(sortdata, navg)
104 104
105 105
106 106 class Beam:
107 107
108 108 def __init__(self):
109 109 self.codeList = []
110 110 self.azimuthList = []
111 111 self.zenithList = []
112 112
113 113 class GenericData(object):
114 114
115 115 flagNoData = True
116 116
117 117 def __init__(self):
118 118
119 119 raise NotImplementedError
120 120
121 121 def copy(self, inputObj=None):
122 122
123 123 if inputObj == None:
124 124 return copy.deepcopy(self)
125 125
126 126 for key in inputObj.__dict__.keys():
127 127
128 128 attribute = inputObj.__dict__[key]
129 129
130 130 #If this attribute is a tuple or list
131 131 if type(inputObj.__dict__[key]) in (tuple, list):
132 132 self.__dict__[key] = attribute[:]
133 133 continue
134 134
135 135 #If this attribute is another object or instance
136 136 if hasattr(attribute, '__dict__'):
137 137 self.__dict__[key] = attribute.copy()
138 138 continue
139 139
140 140 self.__dict__[key] = inputObj.__dict__[key]
141 141
142 142 def deepcopy(self):
143 143
144 144 return copy.deepcopy(self)
145 145
146 146 def isEmpty(self):
147 147
148 148 return self.flagNoData
149 149
150 150 class JROData(GenericData):
151 151
152 152 # m_BasicHeader = BasicHeader()
153 153 # m_ProcessingHeader = ProcessingHeader()
154 154
155 155 systemHeaderObj = SystemHeader()
156 156
157 157 radarControllerHeaderObj = RadarControllerHeader()
158 158
159 159 # data = None
160 160
161 161 type = None
162 162
163 163 datatype = None #dtype but in string
164 164
165 165 # dtype = None
166 166
167 167 # nChannels = None
168 168
169 169 # nHeights = None
170 170
171 171 nProfiles = None
172 172
173 173 heightList = None
174 174
175 175 channelList = None
176 176
177 177 flagDiscontinuousBlock = False
178 178
179 179 useLocalTime = False
180 180
181 181 utctime = None
182 182
183 183 timeZone = None
184 184
185 185 dstFlag = None
186 186
187 187 errorCount = None
188 188
189 189 blocksize = None
190 190
191 191 # nCode = None
192 192 #
193 193 # nBaud = None
194 194 #
195 195 # code = None
196 196
197 197 flagDecodeData = False #asumo q la data no esta decodificada
198 198
199 199 flagDeflipData = False #asumo q la data no esta sin flip
200 200
201 201 flagShiftFFT = False
202 202
203 203 # ippSeconds = None
204 204
205 205 # timeInterval = None
206 206
207 207 nCohInt = None
208 208
209 209 # noise = None
210 210
211 211 windowOfFilter = 1
212 212
213 213 #Speed of ligth
214 214 C = 3e8
215 215
216 216 frequency = 49.92e6
217 217
218 218 realtime = False
219 219
220 220 beacon_heiIndexList = None
221 221
222 222 last_block = None
223 223
224 224 blocknow = None
225 225
226 226 azimuth = None
227 227
228 228 zenith = None
229 229
230 230 beam = Beam()
231 231
232 232 profileIndex = None
233 233
234 234 def __init__(self):
235 235
236 236 raise NotImplementedError
237 237
238 238 def getNoise(self):
239 239
240 240 raise NotImplementedError
241 241
242 242 def getNChannels(self):
243 243
244 244 return len(self.channelList)
245 245
246 246 def getChannelIndexList(self):
247 247
248 248 return range(self.nChannels)
249 249
250 250 def getNHeights(self):
251 251
252 252 return len(self.heightList)
253 253
254 254 def getHeiRange(self, extrapoints=0):
255 255
256 256 heis = self.heightList
257 257 # deltah = self.heightList[1] - self.heightList[0]
258 258 #
259 259 # heis.append(self.heightList[-1])
260 260
261 261 return heis
262 262
263 263 def getDeltaH(self):
264 264
265 265 delta = self.heightList[1] - self.heightList[0]
266 266
267 267 return delta
268 268
269 269 def getltctime(self):
270 270
271 271 if self.useLocalTime:
272 272 return self.utctime - self.timeZone*60
273 273
274 274 return self.utctime
275 275
276 276 def getDatatime(self):
277 277
278 278 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
279 279 return datatimeValue
280 280
281 281 def getTimeRange(self):
282 282
283 283 datatime = []
284 284
285 285 datatime.append(self.ltctime)
286 286 datatime.append(self.ltctime + self.timeInterval+1)
287 287
288 288 datatime = numpy.array(datatime)
289 289
290 290 return datatime
291 291
292 292 def getFmaxTimeResponse(self):
293 293
294 294 period = (10**-6)*self.getDeltaH()/(0.15)
295 295
296 296 PRF = 1./(period * self.nCohInt)
297 297
298 298 fmax = PRF
299 299
300 300 return fmax
301 301
302 302 def getFmax(self):
303 303
304 304 PRF = 1./(self.ippSeconds * self.nCohInt)
305 305
306 306 fmax = PRF
307 307
308 308 return fmax
309 309
310 310 def getVmax(self):
311 311
312 312 _lambda = self.C/self.frequency
313 313
314 314 vmax = self.getFmax() * _lambda/2
315 315
316 316 return vmax
317 317
318 318 def get_ippSeconds(self):
319 319 '''
320 320 '''
321 321 return self.radarControllerHeaderObj.ippSeconds
322 322
323 323 def set_ippSeconds(self, ippSeconds):
324 324 '''
325 325 '''
326 326
327 327 self.radarControllerHeaderObj.ippSeconds = ippSeconds
328 328
329 329 return
330 330
331 331 def get_dtype(self):
332 332 '''
333 333 '''
334 334 return getNumpyDtype(self.datatype)
335 335
336 336 def set_dtype(self, numpyDtype):
337 337 '''
338 338 '''
339 339
340 340 self.datatype = getDataTypeCode(numpyDtype)
341 341
342 342 def get_code(self):
343 343 '''
344 344 '''
345 345 return self.radarControllerHeaderObj.code
346 346
347 347 def set_code(self, code):
348 348 '''
349 349 '''
350 350 self.radarControllerHeaderObj.code = code
351 351
352 352 return
353 353
354 354 def get_ncode(self):
355 355 '''
356 356 '''
357 357 return self.radarControllerHeaderObj.nCode
358 358
359 359 def set_ncode(self, nCode):
360 360 '''
361 361 '''
362 362 self.radarControllerHeaderObj.nCode = nCode
363 363
364 364 return
365 365
366 366 def get_nbaud(self):
367 367 '''
368 368 '''
369 369 return self.radarControllerHeaderObj.nBaud
370 370
371 371 def set_nbaud(self, nBaud):
372 372 '''
373 373 '''
374 374 self.radarControllerHeaderObj.nBaud = nBaud
375 375
376 376 return
377 377
378 378 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
379 379 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
380 380 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
381 381 #noise = property(getNoise, "I'm the 'nHeights' property.")
382 382 datatime = property(getDatatime, "I'm the 'datatime' property")
383 383 ltctime = property(getltctime, "I'm the 'ltctime' property")
384 384 ippSeconds = property(get_ippSeconds, set_ippSeconds)
385 385 dtype = property(get_dtype, set_dtype)
386 386 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
387 387 code = property(get_code, set_code)
388 388 nCode = property(get_ncode, set_ncode)
389 389 nBaud = property(get_nbaud, set_nbaud)
390 390
391 391 class Voltage(JROData):
392 392
393 393 #data es un numpy array de 2 dmensiones (canales, alturas)
394 394 data = None
395 395
396 396 def __init__(self):
397 397 '''
398 398 Constructor
399 399 '''
400 400
401 401 self.useLocalTime = True
402 402
403 403 self.radarControllerHeaderObj = RadarControllerHeader()
404 404
405 405 self.systemHeaderObj = SystemHeader()
406 406
407 407 self.type = "Voltage"
408 408
409 409 self.data = None
410 410
411 411 # self.dtype = None
412 412
413 413 # self.nChannels = 0
414 414
415 415 # self.nHeights = 0
416 416
417 417 self.nProfiles = None
418 418
419 419 self.heightList = None
420 420
421 421 self.channelList = None
422 422
423 423 # self.channelIndexList = None
424 424
425 425 self.flagNoData = True
426 426
427 427 self.flagDiscontinuousBlock = False
428 428
429 429 self.utctime = None
430 430
431 431 self.timeZone = None
432 432
433 433 self.dstFlag = None
434 434
435 435 self.errorCount = None
436 436
437 437 self.nCohInt = None
438 438
439 439 self.blocksize = None
440 440
441 441 self.flagDecodeData = False #asumo q la data no esta decodificada
442 442
443 443 self.flagDeflipData = False #asumo q la data no esta sin flip
444 444
445 445 self.flagShiftFFT = False
446 446
447 447 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
448 448
449 449 self.profileIndex = 0
450 450
451 451 def getNoisebyHildebrand(self, channel = None):
452 452 """
453 453 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
454 454
455 455 Return:
456 456 noiselevel
457 457 """
458 458
459 459 if channel != None:
460 460 data = self.data[channel]
461 461 nChannels = 1
462 462 else:
463 463 data = self.data
464 464 nChannels = self.nChannels
465 465
466 466 noise = numpy.zeros(nChannels)
467 467 power = data * numpy.conjugate(data)
468 468
469 469 for thisChannel in range(nChannels):
470 470 if nChannels == 1:
471 471 daux = power[:].real
472 472 else:
473 473 daux = power[thisChannel,:].real
474 474 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
475 475
476 476 return noise
477 477
478 478 def getNoise(self, type = 1, channel = None):
479 479
480 480 if type == 1:
481 481 noise = self.getNoisebyHildebrand(channel)
482 482
483 483 return noise
484 484
485 485 def getPower(self, channel = None):
486 486
487 487 if channel != None:
488 488 data = self.data[channel]
489 489 else:
490 490 data = self.data
491 491
492 492 power = data * numpy.conjugate(data)
493 493 powerdB = 10*numpy.log10(power.real)
494 494 powerdB = numpy.squeeze(powerdB)
495 495
496 496 return powerdB
497 497
498 498 def getTimeInterval(self):
499 499
500 500 timeInterval = self.ippSeconds * self.nCohInt
501 501
502 502 return timeInterval
503 503
504 504 noise = property(getNoise, "I'm the 'nHeights' property.")
505 505 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
506 506
507 507 class Spectra(JROData):
508 508
509 509 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
510 510 data_spc = None
511 511
512 512 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
513 513 data_cspc = None
514 514
515 515 #data dc es un numpy array de 2 dmensiones (canales, alturas)
516 516 data_dc = None
517 517
518 518 #data power
519 519 data_pwr = None
520 520
521 521 nFFTPoints = None
522 522
523 523 # nPairs = None
524 524
525 525 pairsList = None
526 526
527 527 nIncohInt = None
528 528
529 529 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
530 530
531 531 nCohInt = None #se requiere para determinar el valor de timeInterval
532 532
533 533 ippFactor = None
534 534
535 535 profileIndex = 0
536 536
537 537 plotting = "spectra"
538 538
539 539 def __init__(self):
540 540 '''
541 541 Constructor
542 542 '''
543 543
544 544 self.useLocalTime = True
545 545
546 546 self.radarControllerHeaderObj = RadarControllerHeader()
547 547
548 548 self.systemHeaderObj = SystemHeader()
549 549
550 550 self.type = "Spectra"
551 551
552 552 # self.data = None
553 553
554 554 # self.dtype = None
555 555
556 556 # self.nChannels = 0
557 557
558 558 # self.nHeights = 0
559 559
560 560 self.nProfiles = None
561 561
562 562 self.heightList = None
563 563
564 564 self.channelList = None
565 565
566 566 # self.channelIndexList = None
567 567
568 568 self.pairsList = None
569 569
570 570 self.flagNoData = True
571 571
572 572 self.flagDiscontinuousBlock = False
573 573
574 574 self.utctime = None
575 575
576 576 self.nCohInt = None
577 577
578 578 self.nIncohInt = None
579 579
580 580 self.blocksize = None
581 581
582 582 self.nFFTPoints = None
583 583
584 584 self.wavelength = None
585 585
586 586 self.flagDecodeData = False #asumo q la data no esta decodificada
587 587
588 588 self.flagDeflipData = False #asumo q la data no esta sin flip
589 589
590 590 self.flagShiftFFT = False
591 591
592 592 self.ippFactor = 1
593 593
594 594 #self.noise = None
595 595
596 596 self.beacon_heiIndexList = []
597 597
598 598 self.noise_estimation = None
599 599
600 600
601 601 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
602 602 """
603 603 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
604 604
605 605 Return:
606 606 noiselevel
607 607 """
608 608
609 609 noise = numpy.zeros(self.nChannels)
610 610
611 611 for channel in range(self.nChannels):
612 612 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
613 613 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
614 614
615 615 return noise
616 616
617 617 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
618 618
619 619 if self.noise_estimation is not None:
620 620 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
621 621 else:
622 622 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
623 623 return noise
624 624
625 625 def getFreqRangeTimeResponse(self, extrapoints=0):
626 626
627 627 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
628 628 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
629 629
630 630 return freqrange
631 631
632 632 def getAcfRange(self, extrapoints=0):
633 633
634 634 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
635 635 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
636 636
637 637 return freqrange
638 638
639 639 def getFreqRange(self, extrapoints=0):
640 640
641 641 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
642 642 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
643 643
644 644 return freqrange
645 645
646 646 def getVelRange(self, extrapoints=0):
647 647
648 648 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
649 649 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
650 650
651 651 return velrange
652 652
653 653 def getNPairs(self):
654 654
655 655 return len(self.pairsList)
656 656
657 657 def getPairsIndexList(self):
658 658
659 659 return range(self.nPairs)
660 660
661 661 def getNormFactor(self):
662 662
663 663 pwcode = 1
664 664
665 665 if self.flagDecodeData:
666 666 pwcode = numpy.sum(self.code[0]**2)
667 667 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
668 668 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
669 669
670 670 return normFactor
671 671
672 672 def getFlagCspc(self):
673 673
674 674 if self.data_cspc is None:
675 675 return True
676 676
677 677 return False
678 678
679 679 def getFlagDc(self):
680 680
681 681 if self.data_dc is None:
682 682 return True
683 683
684 684 return False
685 685
686 686 def getTimeInterval(self):
687 687
688 688 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
689 689
690 690 return timeInterval
691 691
692 692 def getPower(self):
693 693
694 694 factor = self.normFactor
695 695 z = self.data_spc/factor
696 696 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
697 697 avg = numpy.average(z, axis=1)
698 698
699 699 return 10*numpy.log10(avg)
700 700
701 701 def getCoherence(self, pairsList=None, phase=False):
702 702
703 703 z = []
704 704 if pairsList is None:
705 705 pairsIndexList = self.pairsIndexList
706 706 else:
707 707 pairsIndexList = []
708 708 for pair in pairsList:
709 709 if pair not in self.pairsList:
710 710 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
711 711 pairsIndexList.append(self.pairsList.index(pair))
712 712 for i in range(len(pairsIndexList)):
713 713 pair = self.pairsList[pairsIndexList[i]]
714 714 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
715 715 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
716 716 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
717 717 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
718 718 if phase:
719 719 data = numpy.arctan2(avgcoherenceComplex.imag,
720 720 avgcoherenceComplex.real)*180/numpy.pi
721 721 else:
722 722 data = numpy.abs(avgcoherenceComplex)
723 723
724 724 z.append(data)
725 725
726 726 return numpy.array(z)
727 727
728 728 def setValue(self, value):
729 729
730 730 print "This property should not be initialized"
731 731
732 732 return
733 733
734 734 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
735 735 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
736 736 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
737 737 flag_cspc = property(getFlagCspc, setValue)
738 738 flag_dc = property(getFlagDc, setValue)
739 739 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
740 740 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
741 741
742 742 class SpectraHeis(Spectra):
743 743
744 744 data_spc = None
745 745
746 746 data_cspc = None
747 747
748 748 data_dc = None
749 749
750 750 nFFTPoints = None
751 751
752 752 # nPairs = None
753 753
754 754 pairsList = None
755 755
756 756 nCohInt = None
757 757
758 758 nIncohInt = None
759 759
760 760 def __init__(self):
761 761
762 762 self.radarControllerHeaderObj = RadarControllerHeader()
763 763
764 764 self.systemHeaderObj = SystemHeader()
765 765
766 766 self.type = "SpectraHeis"
767 767
768 768 # self.dtype = None
769 769
770 770 # self.nChannels = 0
771 771
772 772 # self.nHeights = 0
773 773
774 774 self.nProfiles = None
775 775
776 776 self.heightList = None
777 777
778 778 self.channelList = None
779 779
780 780 # self.channelIndexList = None
781 781
782 782 self.flagNoData = True
783 783
784 784 self.flagDiscontinuousBlock = False
785 785
786 786 # self.nPairs = 0
787 787
788 788 self.utctime = None
789 789
790 790 self.blocksize = None
791 791
792 792 self.profileIndex = 0
793 793
794 794 self.nCohInt = 1
795 795
796 796 self.nIncohInt = 1
797 797
798 798 def getNormFactor(self):
799 799 pwcode = 1
800 800 if self.flagDecodeData:
801 801 pwcode = numpy.sum(self.code[0]**2)
802 802
803 803 normFactor = self.nIncohInt*self.nCohInt*pwcode
804 804
805 805 return normFactor
806 806
807 807 def getTimeInterval(self):
808 808
809 809 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
810 810
811 811 return timeInterval
812 812
813 813 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
814 814 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
815 815
816 816 class Fits(JROData):
817 817
818 818 heightList = None
819 819
820 820 channelList = None
821 821
822 822 flagNoData = True
823 823
824 824 flagDiscontinuousBlock = False
825 825
826 826 useLocalTime = False
827 827
828 828 utctime = None
829 829
830 830 timeZone = None
831 831
832 832 # ippSeconds = None
833 833
834 834 # timeInterval = None
835 835
836 836 nCohInt = None
837 837
838 838 nIncohInt = None
839 839
840 840 noise = None
841 841
842 842 windowOfFilter = 1
843 843
844 844 #Speed of ligth
845 845 C = 3e8
846 846
847 847 frequency = 49.92e6
848 848
849 849 realtime = False
850 850
851 851
852 852 def __init__(self):
853 853
854 854 self.type = "Fits"
855 855
856 856 self.nProfiles = None
857 857
858 858 self.heightList = None
859 859
860 860 self.channelList = None
861 861
862 862 # self.channelIndexList = None
863 863
864 864 self.flagNoData = True
865 865
866 866 self.utctime = None
867 867
868 868 self.nCohInt = 1
869 869
870 870 self.nIncohInt = 1
871 871
872 872 self.useLocalTime = True
873 873
874 874 self.profileIndex = 0
875 875
876 876 # self.utctime = None
877 877 # self.timeZone = None
878 878 # self.ltctime = None
879 879 # self.timeInterval = None
880 880 # self.header = None
881 881 # self.data_header = None
882 882 # self.data = None
883 883 # self.datatime = None
884 884 # self.flagNoData = False
885 885 # self.expName = ''
886 886 # self.nChannels = None
887 887 # self.nSamples = None
888 888 # self.dataBlocksPerFile = None
889 889 # self.comments = ''
890 890 #
891 891
892 892
893 893 def getltctime(self):
894 894
895 895 if self.useLocalTime:
896 896 return self.utctime - self.timeZone*60
897 897
898 898 return self.utctime
899 899
900 900 def getDatatime(self):
901 901
902 902 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
903 903 return datatime
904 904
905 905 def getTimeRange(self):
906 906
907 907 datatime = []
908 908
909 909 datatime.append(self.ltctime)
910 910 datatime.append(self.ltctime + self.timeInterval)
911 911
912 912 datatime = numpy.array(datatime)
913 913
914 914 return datatime
915 915
916 916 def getHeiRange(self):
917 917
918 918 heis = self.heightList
919 919
920 920 return heis
921 921
922 922 def getNHeights(self):
923 923
924 924 return len(self.heightList)
925 925
926 926 def getNChannels(self):
927 927
928 928 return len(self.channelList)
929 929
930 930 def getChannelIndexList(self):
931 931
932 932 return range(self.nChannels)
933 933
934 934 def getNoise(self, type = 1):
935 935
936 936 #noise = numpy.zeros(self.nChannels)
937 937
938 938 if type == 1:
939 939 noise = self.getNoisebyHildebrand()
940 940
941 941 if type == 2:
942 942 noise = self.getNoisebySort()
943 943
944 944 if type == 3:
945 945 noise = self.getNoisebyWindow()
946 946
947 947 return noise
948 948
949 949 def getTimeInterval(self):
950 950
951 951 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
952 952
953 953 return timeInterval
954 954
955 955 datatime = property(getDatatime, "I'm the 'datatime' property")
956 956 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
957 957 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
958 958 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
959 959 noise = property(getNoise, "I'm the 'nHeights' property.")
960 960
961 961 ltctime = property(getltctime, "I'm the 'ltctime' property")
962 962 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
963 963
964 964
965 965 class Correlation(JROData):
966 966
967 967 noise = None
968 968
969 969 SNR = None
970 970
971 971 #--------------------------------------------------
972 972
973 973 mode = None
974 974
975 975 split = False
976 976
977 977 data_cf = None
978 978
979 979 lags = None
980 980
981 981 lagRange = None
982 982
983 983 pairsList = None
984 984
985 985 normFactor = None
986 986
987 987 #--------------------------------------------------
988 988
989 989 # calculateVelocity = None
990 990
991 991 nLags = None
992 992
993 993 nPairs = None
994 994
995 995 nAvg = None
996 996
997 997
998 998 def __init__(self):
999 999 '''
1000 1000 Constructor
1001 1001 '''
1002 1002 self.radarControllerHeaderObj = RadarControllerHeader()
1003 1003
1004 1004 self.systemHeaderObj = SystemHeader()
1005 1005
1006 1006 self.type = "Correlation"
1007 1007
1008 1008 self.data = None
1009 1009
1010 1010 self.dtype = None
1011 1011
1012 1012 self.nProfiles = None
1013 1013
1014 1014 self.heightList = None
1015 1015
1016 1016 self.channelList = None
1017 1017
1018 1018 self.flagNoData = True
1019 1019
1020 1020 self.flagDiscontinuousBlock = False
1021 1021
1022 1022 self.utctime = None
1023 1023
1024 1024 self.timeZone = None
1025 1025
1026 1026 self.dstFlag = None
1027 1027
1028 1028 self.errorCount = None
1029 1029
1030 1030 self.blocksize = None
1031 1031
1032 1032 self.flagDecodeData = False #asumo q la data no esta decodificada
1033 1033
1034 1034 self.flagDeflipData = False #asumo q la data no esta sin flip
1035 1035
1036 1036 self.pairsList = None
1037 1037
1038 1038 self.nPoints = None
1039 1039
1040 1040 def getPairsList(self):
1041 1041
1042 1042 return self.pairsList
1043 1043
1044 1044 def getNoise(self, mode = 2):
1045 1045
1046 1046 indR = numpy.where(self.lagR == 0)[0][0]
1047 1047 indT = numpy.where(self.lagT == 0)[0][0]
1048 1048
1049 1049 jspectra0 = self.data_corr[:,:,indR,:]
1050 1050 jspectra = copy.copy(jspectra0)
1051 1051
1052 1052 num_chan = jspectra.shape[0]
1053 1053 num_hei = jspectra.shape[2]
1054 1054
1055 1055 freq_dc = jspectra.shape[1]/2
1056 1056 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1057 1057
1058 1058 if ind_vel[0]<0:
1059 1059 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1060 1060
1061 1061 if mode == 1:
1062 1062 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1063 1063
1064 1064 if mode == 2:
1065 1065
1066 1066 vel = numpy.array([-2,-1,1,2])
1067 1067 xx = numpy.zeros([4,4])
1068 1068
1069 1069 for fil in range(4):
1070 1070 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1071 1071
1072 1072 xx_inv = numpy.linalg.inv(xx)
1073 1073 xx_aux = xx_inv[0,:]
1074 1074
1075 1075 for ich in range(num_chan):
1076 1076 yy = jspectra[ich,ind_vel,:]
1077 1077 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1078 1078
1079 1079 junkid = jspectra[ich,freq_dc,:]<=0
1080 1080 cjunkid = sum(junkid)
1081 1081
1082 1082 if cjunkid.any():
1083 1083 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1084 1084
1085 1085 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1086 1086
1087 1087 return noise
1088 1088
1089 1089 def getTimeInterval(self):
1090 1090
1091 1091 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1092 1092
1093 1093 return timeInterval
1094 1094
1095 1095 def splitFunctions(self):
1096 1096
1097 1097 pairsList = self.pairsList
1098 1098 ccf_pairs = []
1099 1099 acf_pairs = []
1100 1100 ccf_ind = []
1101 1101 acf_ind = []
1102 1102 for l in range(len(pairsList)):
1103 1103 chan0 = pairsList[l][0]
1104 1104 chan1 = pairsList[l][1]
1105 1105
1106 1106 #Obteniendo pares de Autocorrelacion
1107 1107 if chan0 == chan1:
1108 1108 acf_pairs.append(chan0)
1109 1109 acf_ind.append(l)
1110 1110 else:
1111 1111 ccf_pairs.append(pairsList[l])
1112 1112 ccf_ind.append(l)
1113 1113
1114 1114 data_acf = self.data_cf[acf_ind]
1115 1115 data_ccf = self.data_cf[ccf_ind]
1116 1116
1117 1117 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1118 1118
1119 1119 def getNormFactor(self):
1120 1120 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1121 1121 acf_pairs = numpy.array(acf_pairs)
1122 1122 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1123 1123
1124 1124 for p in range(self.nPairs):
1125 1125 pair = self.pairsList[p]
1126 1126
1127 1127 ch0 = pair[0]
1128 1128 ch1 = pair[1]
1129 1129
1130 1130 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1131 1131 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1132 1132 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1133 1133
1134 1134 return normFactor
1135 1135
1136 1136 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1137 1137 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1138 1138
1139 class Parameters(JROData):
1139 class Parameters(Spectra):
1140 1140
1141 1141 experimentInfo = None #Information about the experiment
1142 1142
1143 1143 #Information from previous data
1144 1144
1145 1145 inputUnit = None #Type of data to be processed
1146 1146
1147 1147 operation = None #Type of operation to parametrize
1148 1148
1149 normFactor = None #Normalization Factor
1149 #normFactor = None #Normalization Factor
1150 1150
1151 1151 groupList = None #List of Pairs, Groups, etc
1152 1152
1153 1153 #Parameters
1154 1154
1155 1155 data_param = None #Parameters obtained
1156 1156
1157 1157 data_pre = None #Data Pre Parametrization
1158 1158
1159 1159 data_SNR = None #Signal to Noise Ratio
1160 1160
1161 1161 # heightRange = None #Heights
1162 1162
1163 1163 abscissaList = None #Abscissa, can be velocities, lags or time
1164 1164
1165 noise = None #Noise Potency
1165 #noise = None #Noise Potency
1166 1166
1167 1167 utctimeInit = None #Initial UTC time
1168 1168
1169 1169 paramInterval = None #Time interval to calculate Parameters in seconds
1170 1170
1171 1171 useLocalTime = True
1172 1172
1173 1173 #Fitting
1174 1174
1175 1175 data_error = None #Error of the estimation
1176 1176
1177 1177 constants = None
1178 1178
1179 1179 library = None
1180 1180
1181 1181 #Output signal
1182 1182
1183 1183 outputInterval = None #Time interval to calculate output signal in seconds
1184 1184
1185 1185 data_output = None #Out signal
1186 1186
1187 1187 nAvg = None
1188 1188
1189 noise_estimation = None
1190
1189 1191
1190 1192 def __init__(self):
1191 1193 '''
1192 1194 Constructor
1193 1195 '''
1194 1196 self.radarControllerHeaderObj = RadarControllerHeader()
1195 1197
1196 1198 self.systemHeaderObj = SystemHeader()
1197 1199
1198 1200 self.type = "Parameters"
1199 1201
1200 1202 def getTimeRange1(self, interval):
1201 1203
1202 1204 datatime = []
1203 1205
1204 1206 if self.useLocalTime:
1205 1207 time1 = self.utctimeInit - self.timeZone*60
1206 1208 else:
1207 1209 time1 = self.utctimeInit
1208 1210
1209 1211 # datatime.append(self.utctimeInit)
1210 1212 # datatime.append(self.utctimeInit + self.outputInterval)
1211 1213 datatime.append(time1)
1212 1214 datatime.append(time1 + interval)
1213 1215
1214 1216 datatime = numpy.array(datatime)
1215 1217
1216 return datatime
1218 return datatime
Show file before | Show file after | Hide comments
@@ -1,426 +1,604
1 1
2 2 import os
3 3 import zmq
4 4 import time
5 5 import numpy
6 6 import datetime
7 7 import numpy as np
8 8 import matplotlib.pyplot as plt
9 9 from mpl_toolkits.axes_grid1 import make_axes_locatable
10 10 from matplotlib.ticker import FuncFormatter, LinearLocator
11 11 from multiprocessing import Process
12 12
13 13 from schainpy.model.proc.jroproc_base import Operation
14 14
15 15 #plt.ion()
16 16
17 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime('%H:%M'))
17 func = lambda x, pos: ('%s') %(datetime.datetime.fromtimestamp(x).strftime('%H:%M'))
18 18
19 19 d1970 = datetime.datetime(1970,1,1)
20 20
21 21 class PlotData(Operation, Process):
22 22
23 23 CODE = 'Figure'
24 colormap = 'jet'
24 colormap = 'jro'
25 25 CONFLATE = True
26 26 __MAXNUMX = 80
27 27 __MAXNUMY = 80
28 28 __missing = 1E30
29 29
30 30 def __init__(self, **kwargs):
31 31
32 32 Operation.__init__(self, plot=True, **kwargs)
33 33 Process.__init__(self)
34 34 self.kwargs['code'] = self.CODE
35 35 self.mp = False
36 36 self.dataOut = None
37 37 self.isConfig = False
38 38 self.figure = None
39 39 self.axes = []
40 40 self.localtime = kwargs.pop('localtime', True)
41 41 self.show = kwargs.get('show', True)
42 42 self.save = kwargs.get('save', False)
43 43 self.colormap = kwargs.get('colormap', self.colormap)
44 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
45 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
44 46 self.showprofile = kwargs.get('showprofile', True)
45 47 self.title = kwargs.get('wintitle', '')
46 self.xaxis = kwargs.get('xaxis', 'time')
48 self.xaxis = kwargs.get('xaxis', 'frequency')
47 49 self.zmin = kwargs.get('zmin', None)
48 50 self.zmax = kwargs.get('zmax', None)
49 51 self.xmin = kwargs.get('xmin', None)
50 52 self.xmax = kwargs.get('xmax', None)
51 53 self.xrange = kwargs.get('xrange', 24)
52 54 self.ymin = kwargs.get('ymin', None)
53 55 self.ymax = kwargs.get('ymax', None)
54 56 self.throttle_value = 5
57
55 58 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
56 59
57 60 if x_buffer.shape[0] < 2:
58 61 return x_buffer, y_buffer, z_buffer
59 62
60 63 deltas = x_buffer[1:] - x_buffer[0:-1]
61 64 x_median = np.median(deltas)
62 65
63 66 index = np.where(deltas > 5*x_median)
64 67
65 68 if len(index[0]) != 0:
66 69 z_buffer[::, index[0], ::] = self.__missing
67 70 z_buffer = np.ma.masked_inside(z_buffer,
68 71 0.99*self.__missing,
69 72 1.01*self.__missing)
70 73
71 74 return x_buffer, y_buffer, z_buffer
72 75
73 76 def decimate(self):
74 77
75 78 # dx = int(len(self.x)/self.__MAXNUMX) + 1
76 79 dy = int(len(self.y)/self.__MAXNUMY) + 1
77 80
78 81 # x = self.x[::dx]
79 82 x = self.x
80 83 y = self.y[::dy]
81 84 z = self.z[::, ::, ::dy]
82 85
83 86 return x, y, z
84 87
85 88 def __plot(self):
86 89
87 90 print 'plotting...{}'.format(self.CODE)
88 91
89 92 if self.show:
90 93 self.figure.show()
91 94
92 95 self.plot()
93 self.figure.suptitle('{} {} - Date:{}'.format(self.title, self.CODE.upper(),
94 datetime.datetime.utcfromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')))
96 plt.tight_layout()
97 self.figure.canvas.manager.set_window_title('{} {} - Date:{}'.format(self.title, self.CODE.upper(),
98 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')))
95 99
96 100 if self.save:
97 101 figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,
98 datetime.datetime.utcfromtimestamp(self.times[0]).strftime('%y%m%d_%H%M%S')))
102 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
99 103 print 'Saving figure: {}'.format(figname)
100 104 self.figure.savefig(figname)
101 105
102 106 self.figure.canvas.draw()
103 107
104 108 def plot(self):
105 109
106 110 print 'plotting...{}'.format(self.CODE.upper())
107 111 return
108 112
109 113 def run(self):
110 114
111 115 print '[Starting] {}'.format(self.name)
112 116 context = zmq.Context()
113 117 receiver = context.socket(zmq.SUB)
114 118 receiver.setsockopt(zmq.SUBSCRIBE, '')
115 119 receiver.setsockopt(zmq.CONFLATE, self.CONFLATE)
116 120 receiver.connect("ipc:///tmp/zmq.plots")
117 121
118 122 while True:
119 123 try:
120 #if True:
121 124 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)
122 125 self.dataOut = self.data['dataOut']
123 126 self.times = self.data['times']
124 127 self.times.sort()
125 128 self.throttle_value = self.data['throttle']
126 129 self.min_time = self.times[0]
127 130 self.max_time = self.times[-1]
128 131
129 132 if self.isConfig is False:
130 133 self.setup()
131 134 self.isConfig = True
132 135 self.__plot()
133 136
134 137 if self.data['ENDED'] is True:
135 # self.__plot()
136 138 self.isConfig = False
137 139
138 140 except zmq.Again as e:
139 141 print 'Waiting for data...'
140 142 plt.pause(self.throttle_value)
141 # time.sleep(3)
142 143
143 144 def close(self):
144 145 if self.dataOut:
145 self._plot()
146 self.__plot()
146 147
147 148
148 149 class PlotSpectraData(PlotData):
149 150
150 151 CODE = 'spc'
151 152 colormap = 'jro'
152 153 CONFLATE = False
154
153 155 def setup(self):
154 156
155 157 ncolspan = 1
156 158 colspan = 1
157 159 self.ncols = int(numpy.sqrt(self.dataOut.nChannels)+0.9)
158 160 self.nrows = int(self.dataOut.nChannels*1./self.ncols + 0.9)
159 161 self.width = 3.6*self.ncols
160 162 self.height = 3.2*self.nrows
161 163 if self.showprofile:
162 164 ncolspan = 3
163 165 colspan = 2
164 166 self.width += 1.2*self.ncols
165 167
166 168 self.ylabel = 'Range [Km]'
167 169 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
168 170
169 171 if self.figure is None:
170 172 self.figure = plt.figure(figsize=(self.width, self.height),
171 173 edgecolor='k',
172 174 facecolor='w')
173 175 else:
174 176 self.figure.clf()
175 177
176 178 n = 0
177 179 for y in range(self.nrows):
178 180 for x in range(self.ncols):
179 181 if n >= self.dataOut.nChannels:
180 182 break
181 183 ax = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan), 1, colspan)
182 184 if self.showprofile:
183 185 ax.ax_profile = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan+colspan), 1, 1)
184 186
185 187 ax.firsttime = True
186 188 self.axes.append(ax)
187 189 n += 1
188 190
189 self.figure.subplots_adjust(left=0.1, right=0.95, bottom=0.15, top=0.85, wspace=0.9, hspace=0.5)
190
191 191 def plot(self):
192 192
193 193 if self.xaxis == "frequency":
194 194 x = self.dataOut.getFreqRange(1)/1000.
195 195 xlabel = "Frequency (kHz)"
196 196 elif self.xaxis == "time":
197 197 x = self.dataOut.getAcfRange(1)
198 198 xlabel = "Time (ms)"
199 199 else:
200 200 x = self.dataOut.getVelRange(1)
201 201 xlabel = "Velocity (m/s)"
202 202
203 203 y = self.dataOut.getHeiRange()
204 204 z = self.data[self.CODE]
205 205
206 206 for n, ax in enumerate(self.axes):
207 207
208 208 if ax.firsttime:
209 209 self.xmax = self.xmax if self.xmax else np.nanmax(x)
210 210 self.xmin = self.xmin if self.xmin else -self.xmax
211 211 self.ymin = self.ymin if self.ymin else np.nanmin(y)
212 212 self.ymax = self.ymax if self.ymax else np.nanmax(y)
213 213 self.zmin = self.zmin if self.zmin else np.nanmin(z)
214 214 self.zmax = self.zmax if self.zmax else np.nanmax(z)
215 215 ax.plot = ax.pcolormesh(x, y, z[n].T,
216 216 vmin=self.zmin,
217 217 vmax=self.zmax,
218 218 cmap=plt.get_cmap(self.colormap)
219 219 )
220 220 divider = make_axes_locatable(ax)
221 221 cax = divider.new_horizontal(size='3%', pad=0.05)
222 222 self.figure.add_axes(cax)
223 223 plt.colorbar(ax.plot, cax)
224 224
225 225 ax.set_xlim(self.xmin, self.xmax)
226 226 ax.set_ylim(self.ymin, self.ymax)
227 227
228 ax.xaxis.set_major_locator(LinearLocator(5))
229 #ax.yaxis.set_major_locator(LinearLocator(4))
230
231 228 ax.set_ylabel(self.ylabel)
232 229 ax.set_xlabel(xlabel)
233 230
234 231 ax.firsttime = False
235 232
236 233 if self.showprofile:
237 234 ax.plot_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
238 235 ax.ax_profile.set_xlim(self.zmin, self.zmax)
239 236 ax.ax_profile.set_ylim(self.ymin, self.ymax)
240 237 ax.ax_profile.set_xlabel('dB')
241 238 ax.ax_profile.grid(b=True, axis='x')
242 239 ax.plot_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
243 240 color="k", linestyle="dashed", lw=2)[0]
244 241 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
245 242 else:
246 243 ax.plot.set_array(z[n].T.ravel())
247 244 if self.showprofile:
248 245 ax.plot_profile.set_data(self.data['rti'][self.max_time][n], y)
249 246 ax.plot_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
250 247
251 248 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
252 249 size=8)
250 self.saveTime = self.max_time
251
252
253 class PlotCrossSpectraData(PlotData):
254
255 CODE = 'cspc'
256 zmin_coh = None
257 zmax_coh = None
258 zmin_phase = None
259 zmax_phase = None
260 CONFLATE = False
261
262 def setup(self):
263
264 ncolspan = 1
265 colspan = 1
266 self.ncols = 2
267 self.nrows = self.dataOut.nPairs
268 self.width = 3.6*self.ncols
269 self.height = 3.2*self.nrows
270
271 self.ylabel = 'Range [Km]'
272 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
273
274 if self.figure is None:
275 self.figure = plt.figure(figsize=(self.width, self.height),
276 edgecolor='k',
277 facecolor='w')
278 else:
279 self.figure.clf()
280
281 for y in range(self.nrows):
282 for x in range(self.ncols):
283 ax = plt.subplot2grid((self.nrows, self.ncols), (y, x), 1, 1)
284 ax.firsttime = True
285 self.axes.append(ax)
286
287 def plot(self):
288
289 if self.xaxis == "frequency":
290 x = self.dataOut.getFreqRange(1)/1000.
291 xlabel = "Frequency (kHz)"
292 elif self.xaxis == "time":
293 x = self.dataOut.getAcfRange(1)
294 xlabel = "Time (ms)"
295 else:
296 x = self.dataOut.getVelRange(1)
297 xlabel = "Velocity (m/s)"
298
299 y = self.dataOut.getHeiRange()
300 z_coh = self.data['cspc_coh']
301 z_phase = self.data['cspc_phase']
302
303 for n in range(self.nrows):
304 ax = self.axes[2*n]
305 ax1 = self.axes[2*n+1]
306 if ax.firsttime:
307 self.xmax = self.xmax if self.xmax else np.nanmax(x)
308 self.xmin = self.xmin if self.xmin else -self.xmax
309 self.ymin = self.ymin if self.ymin else np.nanmin(y)
310 self.ymax = self.ymax if self.ymax else np.nanmax(y)
311 self.zmin_coh = self.zmin_coh if self.zmin_coh else 0.0
312 self.zmax_coh = self.zmax_coh if self.zmax_coh else 1.0
313 self.zmin_phase = self.zmin_phase if self.zmin_phase else -180
314 self.zmax_phase = self.zmax_phase if self.zmax_phase else 180
315
316 ax.plot = ax.pcolormesh(x, y, z_coh[n].T,
317 vmin=self.zmin_coh,
318 vmax=self.zmax_coh,
319 cmap=plt.get_cmap(self.colormap_coh)
320 )
321 divider = make_axes_locatable(ax)
322 cax = divider.new_horizontal(size='3%', pad=0.05)
323 self.figure.add_axes(cax)
324 plt.colorbar(ax.plot, cax)
325
326 ax.set_xlim(self.xmin, self.xmax)
327 ax.set_ylim(self.ymin, self.ymax)
328
329 ax.set_ylabel(self.ylabel)
330 ax.set_xlabel(xlabel)
331 ax.firsttime = False
332
333 ax1.plot = ax1.pcolormesh(x, y, z_phase[n].T,
334 vmin=self.zmin_phase,
335 vmax=self.zmax_phase,
336 cmap=plt.get_cmap(self.colormap_phase)
337 )
338 divider = make_axes_locatable(ax1)
339 cax = divider.new_horizontal(size='3%', pad=0.05)
340 self.figure.add_axes(cax)
341 plt.colorbar(ax1.plot, cax)
342
343 ax1.set_xlim(self.xmin, self.xmax)
344 ax1.set_ylim(self.ymin, self.ymax)
345
346 ax1.set_ylabel(self.ylabel)
347 ax1.set_xlabel(xlabel)
348 ax1.firsttime = False
349 else:
350 ax.plot.set_array(z_coh[n].T.ravel())
351 ax1.plot.set_array(z_phase[n].T.ravel())
352
353 ax.set_title('Coherence Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
354 ax1.set_title('Phase Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
355 self.saveTime = self.max_time
356
357
358 class PlotSpectraMeanData(PlotSpectraData):
359
360 CODE = 'spc_mean'
361 colormap = 'jet'
362
363 def plot(self):
364
365 if self.xaxis == "frequency":
366 x = self.dataOut.getFreqRange(1)/1000.
367 xlabel = "Frequency (kHz)"
368 elif self.xaxis == "time":
369 x = self.dataOut.getAcfRange(1)
370 xlabel = "Time (ms)"
371 else:
372 x = self.dataOut.getVelRange(1)
373 xlabel = "Velocity (m/s)"
374
375 y = self.dataOut.getHeiRange()
376 z = self.data['spc']
377 mean = self.data['mean'][self.max_time]
378
379 for n, ax in enumerate(self.axes):
380
381 if ax.firsttime:
382 self.xmax = self.xmax if self.xmax else np.nanmax(x)
383 self.xmin = self.xmin if self.xmin else -self.xmax
384 self.ymin = self.ymin if self.ymin else np.nanmin(y)
385 self.ymax = self.ymax if self.ymax else np.nanmax(y)
386 self.zmin = self.zmin if self.zmin else np.nanmin(z)
387 self.zmax = self.zmax if self.zmax else np.nanmax(z)
388 ax.plt = ax.pcolormesh(x, y, z[n].T,
389 vmin=self.zmin,
390 vmax=self.zmax,
391 cmap=plt.get_cmap(self.colormap)
392 )
393 ax.plt_dop = ax.plot(mean[n], y,
394 color='k')[0]
395
396 divider = make_axes_locatable(ax)
397 cax = divider.new_horizontal(size='3%', pad=0.05)
398 self.figure.add_axes(cax)
399 plt.colorbar(ax.plt, cax)
400
401 ax.set_xlim(self.xmin, self.xmax)
402 ax.set_ylim(self.ymin, self.ymax)
403
404 ax.set_ylabel(self.ylabel)
405 ax.set_xlabel(xlabel)
406
407 ax.firsttime = False
408
409 if self.showprofile:
410 ax.plt_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
411 ax.ax_profile.set_xlim(self.zmin, self.zmax)
412 ax.ax_profile.set_ylim(self.ymin, self.ymax)
413 ax.ax_profile.set_xlabel('dB')
414 ax.ax_profile.grid(b=True, axis='x')
415 ax.plt_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
416 color="k", linestyle="dashed", lw=2)[0]
417 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
418 else:
419 ax.plt.set_array(z[n].T.ravel())
420 ax.plt_dop.set_data(mean[n], y)
421 if self.showprofile:
422 ax.plt_profile.set_data(self.data['rti'][self.max_time][n], y)
423 ax.plt_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
424
425 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
426 size=8)
427 self.saveTime = self.max_time
428
253 429
254 430 class PlotRTIData(PlotData):
255 431
256 432 CODE = 'rti'
257 433 colormap = 'jro'
258 434
259 435 def setup(self):
260 436 self.ncols = 1
261 437 self.nrows = self.dataOut.nChannels
262 438 self.width = 10
263 self.height = 2.2*self.nrows
439 self.height = 2.2*self.nrows if self.nrows<6 else 12
264 440 if self.nrows==1:
265 441 self.height += 1
266 442 self.ylabel = 'Range [Km]'
267 443 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
268 444
269 445 if self.figure is None:
270 446 self.figure = plt.figure(figsize=(self.width, self.height),
271 447 edgecolor='k',
272 448 facecolor='w')
273 449 else:
274 450 self.figure.clf()
275 451 self.axes = []
276 452
277 453 for n in range(self.nrows):
278 454 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
279 455 ax.firsttime = True
280 456 self.axes.append(ax)
281 self.figure.subplots_adjust(hspace=0.5)
282 457
283 458 def plot(self):
284 459
285 460 self.x = np.array(self.times)
286 461 self.y = self.dataOut.getHeiRange()
287 462 self.z = []
288 463
289 464 for ch in range(self.nrows):
290 465 self.z.append([self.data[self.CODE][t][ch] for t in self.times])
291 466
292 467 self.z = np.array(self.z)
293 468 for n, ax in enumerate(self.axes):
294 469
295 470 x, y, z = self.fill_gaps(*self.decimate())
296 471 xmin = self.min_time
297 472 xmax = xmin+self.xrange*60*60
298 473 if ax.firsttime:
299 474 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
300 475 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
301 476 self.zmin = self.zmin if self.zmin else np.nanmin(self.z)
302 477 self.zmax = self.zmax if self.zmax else np.nanmax(self.z)
303 478 plot = ax.pcolormesh(x, y, z[n].T,
304 479 vmin=self.zmin,
305 480 vmax=self.zmax,
306 481 cmap=plt.get_cmap(self.colormap)
307 482 )
308 483 divider = make_axes_locatable(ax)
309 484 cax = divider.new_horizontal(size='2%', pad=0.05)
310 485 self.figure.add_axes(cax)
311 486 plt.colorbar(plot, cax)
312 487 ax.set_ylim(self.ymin, self.ymax)
313 if self.xaxis == 'time':
314 ax.xaxis.set_major_formatter(FuncFormatter(func))
315 ax.xaxis.set_major_locator(LinearLocator(6))
316 488
317 # ax.yaxis.set_major_locator(LinearLocator(4))
489 ax.xaxis.set_major_formatter(FuncFormatter(func))
490 ax.xaxis.set_major_locator(LinearLocator(6))
318 491
319 492 ax.set_ylabel(self.ylabel)
320 493
321 494 # if self.xmin is None:
322 495 # xmin = self.min_time
323 496 # else:
324 497 # xmin = (datetime.datetime.combine(self.dataOut.datatime.date(),
325 498 # datetime.time(self.xmin, 0, 0))-d1970).total_seconds()
326 499
327 500 ax.set_xlim(xmin, xmax)
328 501 ax.firsttime = False
329 502 else:
330 503 ax.collections.remove(ax.collections[0])
331 504 ax.set_xlim(xmin, xmax)
332 505 plot = ax.pcolormesh(x, y, z[n].T,
333 506 vmin=self.zmin,
334 507 vmax=self.zmax,
335 508 cmap=plt.get_cmap(self.colormap)
336 509 )
337 ax.set_title('{} {}'.format(self.titles[n],
338 datetime.datetime.utcfromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
339 size=8)
510 ax.set_title('{} {}'.format(self.titles[n],
511 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
512 size=8)
513
514 self.saveTime = self.min_time
340 515
341 516
342 517 class PlotCOHData(PlotRTIData):
343 518
344 519 CODE = 'coh'
345 520
346 521 def setup(self):
347 522
348 523 self.ncols = 1
349 524 self.nrows = self.dataOut.nPairs
350 525 self.width = 10
351 self.height = 2.2*self.nrows
526 self.height = 2.2*self.nrows if self.nrows<6 else 12
352 527 if self.nrows==1:
353 528 self.height += 1
354 529 self.ylabel = 'Range [Km]'
355 self.titles = ['Channels {}'.format(x) for x in self.dataOut.pairsList]
530 self.titles = ['{} Ch{} * Ch{}'.format(self.CODE.upper(), x[0], x[1]) for x in self.dataOut.pairsList]
356 531
357 532 if self.figure is None:
358 533 self.figure = plt.figure(figsize=(self.width, self.height),
359 534 edgecolor='k',
360 535 facecolor='w')
361 536 else:
362 537 self.figure.clf()
363 538 self.axes = []
364 539
365 540 for n in range(self.nrows):
366 541 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
367 542 ax.firsttime = True
368 543 self.axes.append(ax)
369 544
370 self.figure.subplots_adjust(hspace=0.5)
371 545
372 546 class PlotNoiseData(PlotData):
373 547 CODE = 'noise'
374 548
375 549 def setup(self):
376 550
377 551 self.ncols = 1
378 552 self.nrows = 1
379 553 self.width = 10
380 554 self.height = 3.2
381 555 self.ylabel = 'Intensity [dB]'
382 556 self.titles = ['Noise']
383 557
384 558 if self.figure is None:
385 559 self.figure = plt.figure(figsize=(self.width, self.height),
386 560 edgecolor='k',
387 561 facecolor='w')
388 562 else:
389 563 self.figure.clf()
390 564 self.axes = []
391 565
392 566 self.ax = self.figure.add_subplot(self.nrows, self.ncols, 1)
393 567 self.ax.firsttime = True
394 568
395 569 def plot(self):
396 570
397 571 x = self.times
398 572 xmin = self.min_time
399 573 xmax = xmin+self.xrange*60*60
400 574 if self.ax.firsttime:
401 575 for ch in self.dataOut.channelList:
402 576 y = [self.data[self.CODE][t][ch] for t in self.times]
403 577 self.ax.plot(x, y, lw=1, label='Ch{}'.format(ch))
404 578 self.ax.firsttime = False
405 579 self.ax.xaxis.set_major_formatter(FuncFormatter(func))
406 580 self.ax.xaxis.set_major_locator(LinearLocator(6))
407 581 self.ax.set_ylabel(self.ylabel)
408 582 plt.legend()
409 583 else:
410 584 for ch in self.dataOut.channelList:
411 585 y = [self.data[self.CODE][t][ch] for t in self.times]
412 586 self.ax.lines[ch].set_data(x, y)
413 587
414 588 self.ax.set_xlim(xmin, xmax)
415 589 self.ax.set_ylim(min(y)-5, max(y)+5)
590 self.saveTime = self.min_time
591
416 592
417 593 class PlotSNRData(PlotRTIData):
418 594 CODE = 'snr'
595 colormap = 'jet'
419 596
420 597 class PlotDOPData(PlotRTIData):
421 598 CODE = 'dop'
422 599 colormap = 'jet'
423 600
601
424 602 class PlotPHASEData(PlotCOHData):
425 603 CODE = 'phase'
426 604 colormap = 'seismic'
Show file before | Show file after | Hide comments
@@ -1,2741 +1,2749
1 1 import numpy
2 2 import math
3 3 from scipy import optimize, interpolate, signal, stats, ndimage
4 4 import re
5 5 import datetime
6 6 import copy
7 7 import sys
8 8 import importlib
9 9 import itertools
10 10
11 11 from jroproc_base import ProcessingUnit, Operation
12 12 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
13 13
14 14
15 15 class ParametersProc(ProcessingUnit):
16 16
17 17 nSeconds = None
18 18
19 19 def __init__(self):
20 20 ProcessingUnit.__init__(self)
21 21
22 22 # self.objectDict = {}
23 23 self.buffer = None
24 24 self.firstdatatime = None
25 25 self.profIndex = 0
26 26 self.dataOut = Parameters()
27 27
28 28 def __updateObjFromInput(self):
29 29
30 30 self.dataOut.inputUnit = self.dataIn.type
31 31
32 32 self.dataOut.timeZone = self.dataIn.timeZone
33 33 self.dataOut.dstFlag = self.dataIn.dstFlag
34 34 self.dataOut.errorCount = self.dataIn.errorCount
35 35 self.dataOut.useLocalTime = self.dataIn.useLocalTime
36 36
37 37 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
38 38 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 self.dataOut.heightList = self.dataIn.heightList
41 41 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
42 42 # self.dataOut.nHeights = self.dataIn.nHeights
43 43 # self.dataOut.nChannels = self.dataIn.nChannels
44 44 self.dataOut.nBaud = self.dataIn.nBaud
45 45 self.dataOut.nCode = self.dataIn.nCode
46 46 self.dataOut.code = self.dataIn.code
47 47 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
48 48 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
49 49 # self.dataOut.utctime = self.firstdatatime
50 50 self.dataOut.utctime = self.dataIn.utctime
51 51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
52 52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
53 53 self.dataOut.nCohInt = self.dataIn.nCohInt
54 54 # self.dataOut.nIncohInt = 1
55 55 self.dataOut.ippSeconds = self.dataIn.ippSeconds
56 56 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
57 self.dataOut.timeInterval = self.dataIn.timeInterval
57 # self.dataOut.timeInterval = self.dataIn.timeInterval
58 58 self.dataOut.heightList = self.dataIn.getHeiRange()
59 59 self.dataOut.frequency = self.dataIn.frequency
60 self.dataOut.noise = self.dataIn.noise
60 #self.dataOut.noise = self.dataIn.noise
61 61
62 62 def run(self):
63 63
64 64 #---------------------- Voltage Data ---------------------------
65 65
66 66 if self.dataIn.type == "Voltage":
67 67
68 68 self.__updateObjFromInput()
69 69 self.dataOut.data_pre = self.dataIn.data.copy()
70 70 self.dataOut.flagNoData = False
71 71 self.dataOut.utctimeInit = self.dataIn.utctime
72 72 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
73 73 return
74 74
75 75 #---------------------- Spectra Data ---------------------------
76 76
77 77 if self.dataIn.type == "Spectra":
78 78
79 self.dataOut.data_pre = (self.dataIn.data_spc,self.dataIn.data_cspc)
80 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
81 # self.dataOut.noise = self.dataIn.getNoise()
82 self.dataOut.normFactor = self.dataIn.normFactor
79 self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)
80 self.dataOut.data_spc = self.dataIn.data_spc
81 self.dataOut.data_cspc = self.dataIn.data_cspc
82 self.dataOut.nProfiles = self.dataIn.nProfiles
83 self.dataOut.nIncohInt = self.dataIn.nIncohInt
84 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
85 self.dataOut.ippFactor = self.dataIn.ippFactor
86 #self.dataOut.normFactor = self.dataIn.getNormFactor()
87 self.dataOut.pairsList = self.dataIn.pairsList
83 88 self.dataOut.groupList = self.dataIn.pairsList
89 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
84 90 self.dataOut.flagNoData = False
85 91
86 92 #---------------------- Correlation Data ---------------------------
87 93
88 94 if self.dataIn.type == "Correlation":
89 95 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
90 96
91 97 self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
92 98 self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
93 99 self.dataOut.groupList = (acf_pairs, ccf_pairs)
94 100
95 101 self.dataOut.abscissaList = self.dataIn.lagRange
96 102 self.dataOut.noise = self.dataIn.noise
97 103 self.dataOut.data_SNR = self.dataIn.SNR
98 104 self.dataOut.flagNoData = False
99 105 self.dataOut.nAvg = self.dataIn.nAvg
100 106
101 107 #---------------------- Parameters Data ---------------------------
102 108
103 109 if self.dataIn.type == "Parameters":
104 110 self.dataOut.copy(self.dataIn)
105 111 self.dataOut.utctimeInit = self.dataIn.utctime
106 112 self.dataOut.flagNoData = False
107 113
108 114 return True
109 115
110 116 self.__updateObjFromInput()
111 117 self.dataOut.utctimeInit = self.dataIn.utctime
112 118 self.dataOut.paramInterval = self.dataIn.timeInterval
113 119
114 120 return
115 121
116 122 class SpectralMoments(Operation):
117 123
118 124 '''
119 125 Function SpectralMoments()
120 126
121 127 Calculates moments (power, mean, standard deviation) and SNR of the signal
122 128
123 129 Type of dataIn: Spectra
124 130
125 131 Configuration Parameters:
126 132
127 133 dirCosx : Cosine director in X axis
128 134 dirCosy : Cosine director in Y axis
129 135
130 136 elevation :
131 137 azimuth :
132 138
133 139 Input:
134 140 channelList : simple channel list to select e.g. [2,3,7]
135 141 self.dataOut.data_pre : Spectral data
136 142 self.dataOut.abscissaList : List of frequencies
137 143 self.dataOut.noise : Noise level per channel
138 144
139 145 Affected:
140 146 self.dataOut.data_param : Parameters per channel
141 147 self.dataOut.data_SNR : SNR per channel
142 148
143 149 '''
144 150
145 151 def run(self, dataOut):
146 152
147 dataOut.data_pre = dataOut.data_pre[0]
148 data = dataOut.data_pre
153 #dataOut.data_pre = dataOut.data_pre[0]
154 data = dataOut.data_pre[0]
149 155 absc = dataOut.abscissaList[:-1]
150 156 noise = dataOut.noise
151 157 nChannel = data.shape[0]
152 158 data_param = numpy.zeros((nChannel, 4, data.shape[2]))
153 159
154 160 for ind in range(nChannel):
155 161 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
156 162
157 163 dataOut.data_param = data_param[:,1:,:]
158 164 dataOut.data_SNR = data_param[:,0]
159 165 dataOut.data_DOP = data_param[:,1]
166 dataOut.data_MEAN = data_param[:,2]
167 dataOut.data_STD = data_param[:,3]
160 168 return
161 169
162 170 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):
163 171
164 172 if (nicoh is None): nicoh = 1
165 173 if (graph is None): graph = 0
166 174 if (smooth is None): smooth = 0
167 175 elif (self.smooth < 3): smooth = 0
168 176
169 177 if (type1 is None): type1 = 0
170 178 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
171 179 if (snrth is None): snrth = -3
172 180 if (dc is None): dc = 0
173 181 if (aliasing is None): aliasing = 0
174 182 if (oldfd is None): oldfd = 0
175 183 if (wwauto is None): wwauto = 0
176 184
177 185 if (n0 < 1.e-20): n0 = 1.e-20
178 186
179 187 freq = oldfreq
180 188 vec_power = numpy.zeros(oldspec.shape[1])
181 189 vec_fd = numpy.zeros(oldspec.shape[1])
182 190 vec_w = numpy.zeros(oldspec.shape[1])
183 191 vec_snr = numpy.zeros(oldspec.shape[1])
184 192
185 193 for ind in range(oldspec.shape[1]):
186 194
187 195 spec = oldspec[:,ind]
188 196 aux = spec*fwindow
189 197 max_spec = aux.max()
190 198 m = list(aux).index(max_spec)
191 199
192 200 #Smooth
193 201 if (smooth == 0): spec2 = spec
194 202 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
195 203
196 204 # Calculo de Momentos
197 205 bb = spec2[range(m,spec2.size)]
198 206 bb = (bb<n0).nonzero()
199 207 bb = bb[0]
200 208
201 209 ss = spec2[range(0,m + 1)]
202 210 ss = (ss<n0).nonzero()
203 211 ss = ss[0]
204 212
205 213 if (bb.size == 0):
206 214 bb0 = spec.size - 1 - m
207 215 else:
208 216 bb0 = bb[0] - 1
209 217 if (bb0 < 0):
210 218 bb0 = 0
211 219
212 220 if (ss.size == 0): ss1 = 1
213 221 else: ss1 = max(ss) + 1
214 222
215 223 if (ss1 > m): ss1 = m
216 224
217 225 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
218 226 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
219 227 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
220 228 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
221 229 snr = (spec2.mean()-n0)/n0
222 230
223 231 if (snr < 1.e-20) :
224 232 snr = 1.e-20
225 233
226 234 vec_power[ind] = power
227 235 vec_fd[ind] = fd
228 236 vec_w[ind] = w
229 237 vec_snr[ind] = snr
230 238
231 239 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
232 240 return moments
233 241
234 242 #------------------ Get SA Parameters --------------------------
235 243
236 244 def GetSAParameters(self):
237 245 #SA en frecuencia
238 246 pairslist = self.dataOut.groupList
239 247 num_pairs = len(pairslist)
240 248
241 249 vel = self.dataOut.abscissaList
242 spectra = self.dataOut.data_pre
243 cspectra = self.dataIn.data_cspc
250 spectra = self.dataOut.data_pre[0]
251 cspectra = self.dataOut.data_pre[1]
244 252 delta_v = vel[1] - vel[0]
245 253
246 254 #Calculating the power spectrum
247 255 spc_pow = numpy.sum(spectra, 3)*delta_v
248 256 #Normalizing Spectra
249 257 norm_spectra = spectra/spc_pow
250 258 #Calculating the norm_spectra at peak
251 259 max_spectra = numpy.max(norm_spectra, 3)
252 260
253 261 #Normalizing Cross Spectra
254 262 norm_cspectra = numpy.zeros(cspectra.shape)
255 263
256 264 for i in range(num_chan):
257 265 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
258 266
259 267 max_cspectra = numpy.max(norm_cspectra,2)
260 268 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
261 269
262 270 for i in range(num_pairs):
263 271 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
264 272 #------------------- Get Lags ----------------------------------
265 273
266 274 class SALags(Operation):
267 275 '''
268 276 Function GetMoments()
269 277
270 278 Input:
271 279 self.dataOut.data_pre
272 280 self.dataOut.abscissaList
273 281 self.dataOut.noise
274 282 self.dataOut.normFactor
275 283 self.dataOut.data_SNR
276 284 self.dataOut.groupList
277 285 self.dataOut.nChannels
278 286
279 287 Affected:
280 288 self.dataOut.data_param
281 289
282 290 '''
283 291 def run(self, dataOut):
284 292 data_acf = dataOut.data_pre[0]
285 293 data_ccf = dataOut.data_pre[1]
286 294 normFactor_acf = dataOut.normFactor[0]
287 295 normFactor_ccf = dataOut.normFactor[1]
288 296 pairs_acf = dataOut.groupList[0]
289 297 pairs_ccf = dataOut.groupList[1]
290 298
291 299 nHeights = dataOut.nHeights
292 300 absc = dataOut.abscissaList
293 301 noise = dataOut.noise
294 302 SNR = dataOut.data_SNR
295 303 nChannels = dataOut.nChannels
296 304 # pairsList = dataOut.groupList
297 305 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
298 306
299 307 for l in range(len(pairs_acf)):
300 308 data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
301 309
302 310 for l in range(len(pairs_ccf)):
303 311 data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
304 312
305 313 dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
306 314 dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
307 315 dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
308 316 return
309 317
310 318 # def __getPairsAutoCorr(self, pairsList, nChannels):
311 319 #
312 320 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
313 321 #
314 322 # for l in range(len(pairsList)):
315 323 # firstChannel = pairsList[l][0]
316 324 # secondChannel = pairsList[l][1]
317 325 #
318 326 # #Obteniendo pares de Autocorrelacion
319 327 # if firstChannel == secondChannel:
320 328 # pairsAutoCorr[firstChannel] = int(l)
321 329 #
322 330 # pairsAutoCorr = pairsAutoCorr.astype(int)
323 331 #
324 332 # pairsCrossCorr = range(len(pairsList))
325 333 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
326 334 #
327 335 # return pairsAutoCorr, pairsCrossCorr
328 336
329 337 def __calculateTaus(self, data_acf, data_ccf, lagRange):
330 338
331 339 lag0 = data_acf.shape[1]/2
332 340 #Funcion de Autocorrelacion
333 341 mean_acf = stats.nanmean(data_acf, axis = 0)
334 342
335 343 #Obtencion Indice de TauCross
336 344 ind_ccf = data_ccf.argmax(axis = 1)
337 345 #Obtencion Indice de TauAuto
338 346 ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
339 347 ccf_lag0 = data_ccf[:,lag0,:]
340 348
341 349 for i in range(ccf_lag0.shape[0]):
342 350 ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
343 351
344 352 #Obtencion de TauCross y TauAuto
345 353 tau_ccf = lagRange[ind_ccf]
346 354 tau_acf = lagRange[ind_acf]
347 355
348 356 Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
349 357
350 358 tau_ccf[Nan1,Nan2] = numpy.nan
351 359 tau_acf[Nan1,Nan2] = numpy.nan
352 360 tau = numpy.vstack((tau_ccf,tau_acf))
353 361
354 362 return tau
355 363
356 364 def __calculateLag1Phase(self, data, lagTRange):
357 365 data1 = stats.nanmean(data, axis = 0)
358 366 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
359 367
360 368 phase = numpy.angle(data1[lag1,:])
361 369
362 370 return phase
363 371
364 372 class SpectralFitting(Operation):
365 373 '''
366 374 Function GetMoments()
367 375
368 376 Input:
369 377 Output:
370 378 Variables modified:
371 379 '''
372 380
373 381 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
374 382
375 383
376 384 if path != None:
377 385 sys.path.append(path)
378 386 self.dataOut.library = importlib.import_module(file)
379 387
380 388 #To be inserted as a parameter
381 389 groupArray = numpy.array(groupList)
382 390 # groupArray = numpy.array([[0,1],[2,3]])
383 391 self.dataOut.groupList = groupArray
384 392
385 393 nGroups = groupArray.shape[0]
386 394 nChannels = self.dataIn.nChannels
387 395 nHeights=self.dataIn.heightList.size
388 396
389 397 #Parameters Array
390 398 self.dataOut.data_param = None
391 399
392 400 #Set constants
393 401 constants = self.dataOut.library.setConstants(self.dataIn)
394 402 self.dataOut.constants = constants
395 403 M = self.dataIn.normFactor
396 404 N = self.dataIn.nFFTPoints
397 405 ippSeconds = self.dataIn.ippSeconds
398 406 K = self.dataIn.nIncohInt
399 407 pairsArray = numpy.array(self.dataIn.pairsList)
400 408
401 409 #List of possible combinations
402 410 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
403 411 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
404 412
405 413 if getSNR:
406 414 listChannels = groupArray.reshape((groupArray.size))
407 415 listChannels.sort()
408 416 noise = self.dataIn.getNoise()
409 417 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
410 418
411 419 for i in range(nGroups):
412 420 coord = groupArray[i,:]
413 421
414 422 #Input data array
415 423 data = self.dataIn.data_spc[coord,:,:]/(M*N)
416 424 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
417 425
418 426 #Cross Spectra data array for Covariance Matrixes
419 427 ind = 0
420 428 for pairs in listComb:
421 429 pairsSel = numpy.array([coord[x],coord[y]])
422 430 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
423 431 ind += 1
424 432 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
425 433 dataCross = dataCross**2/K
426 434
427 435 for h in range(nHeights):
428 436 # print self.dataOut.heightList[h]
429 437
430 438 #Input
431 439 d = data[:,h]
432 440
433 441 #Covariance Matrix
434 442 D = numpy.diag(d**2/K)
435 443 ind = 0
436 444 for pairs in listComb:
437 445 #Coordinates in Covariance Matrix
438 446 x = pairs[0]
439 447 y = pairs[1]
440 448 #Channel Index
441 449 S12 = dataCross[ind,:,h]
442 450 D12 = numpy.diag(S12)
443 451 #Completing Covariance Matrix with Cross Spectras
444 452 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
445 453 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
446 454 ind += 1
447 455 Dinv=numpy.linalg.inv(D)
448 456 L=numpy.linalg.cholesky(Dinv)
449 457 LT=L.T
450 458
451 459 dp = numpy.dot(LT,d)
452 460
453 461 #Initial values
454 462 data_spc = self.dataIn.data_spc[coord,:,h]
455 463
456 464 if (h>0)and(error1[3]<5):
457 465 p0 = self.dataOut.data_param[i,:,h-1]
458 466 else:
459 467 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
460 468
461 469 try:
462 470 #Least Squares
463 471 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
464 472 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
465 473 #Chi square error
466 474 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
467 475 #Error with Jacobian
468 476 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
469 477 except:
470 478 minp = p0*numpy.nan
471 479 error0 = numpy.nan
472 480 error1 = p0*numpy.nan
473 481
474 482 #Save
475 483 if self.dataOut.data_param is None:
476 484 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
477 485 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
478 486
479 487 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
480 488 self.dataOut.data_param[i,:,h] = minp
481 489 return
482 490
483 491 def __residFunction(self, p, dp, LT, constants):
484 492
485 493 fm = self.dataOut.library.modelFunction(p, constants)
486 494 fmp=numpy.dot(LT,fm)
487 495
488 496 return dp-fmp
489 497
490 498 def __getSNR(self, z, noise):
491 499
492 500 avg = numpy.average(z, axis=1)
493 501 SNR = (avg.T-noise)/noise
494 502 SNR = SNR.T
495 503 return SNR
496 504
497 505 def __chisq(p,chindex,hindex):
498 506 #similar to Resid but calculates CHI**2
499 507 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
500 508 dp=numpy.dot(LT,d)
501 509 fmp=numpy.dot(LT,fm)
502 510 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
503 511 return chisq
504 512
505 513 class WindProfiler(Operation):
506 514
507 515 __isConfig = False
508 516
509 517 __initime = None
510 518 __lastdatatime = None
511 519 __integrationtime = None
512 520
513 521 __buffer = None
514 522
515 523 __dataReady = False
516 524
517 525 __firstdata = None
518 526
519 527 n = None
520 528
521 529 def __calculateCosDir(self, elev, azim):
522 530 zen = (90 - elev)*numpy.pi/180
523 531 azim = azim*numpy.pi/180
524 532 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
525 533 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
526 534
527 535 signX = numpy.sign(numpy.cos(azim))
528 536 signY = numpy.sign(numpy.sin(azim))
529 537
530 538 cosDirX = numpy.copysign(cosDirX, signX)
531 539 cosDirY = numpy.copysign(cosDirY, signY)
532 540 return cosDirX, cosDirY
533 541
534 542 def __calculateAngles(self, theta_x, theta_y, azimuth):
535 543
536 544 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
537 545 zenith_arr = numpy.arccos(dir_cosw)
538 546 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
539 547
540 548 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
541 549 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
542 550
543 551 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
544 552
545 553 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
546 554
547 555 #
548 556 if horOnly:
549 557 A = numpy.c_[dir_cosu,dir_cosv]
550 558 else:
551 559 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
552 560 A = numpy.asmatrix(A)
553 561 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
554 562
555 563 return A1
556 564
557 565 def __correctValues(self, heiRang, phi, velRadial, SNR):
558 566 listPhi = phi.tolist()
559 567 maxid = listPhi.index(max(listPhi))
560 568 minid = listPhi.index(min(listPhi))
561 569
562 570 rango = range(len(phi))
563 571 # rango = numpy.delete(rango,maxid)
564 572
565 573 heiRang1 = heiRang*math.cos(phi[maxid])
566 574 heiRangAux = heiRang*math.cos(phi[minid])
567 575 indOut = (heiRang1 < heiRangAux[0]).nonzero()
568 576 heiRang1 = numpy.delete(heiRang1,indOut)
569 577
570 578 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
571 579 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
572 580
573 581 for i in rango:
574 582 x = heiRang*math.cos(phi[i])
575 583 y1 = velRadial[i,:]
576 584 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
577 585
578 586 x1 = heiRang1
579 587 y11 = f1(x1)
580 588
581 589 y2 = SNR[i,:]
582 590 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
583 591 y21 = f2(x1)
584 592
585 593 velRadial1[i,:] = y11
586 594 SNR1[i,:] = y21
587 595
588 596 return heiRang1, velRadial1, SNR1
589 597
590 598 def __calculateVelUVW(self, A, velRadial):
591 599
592 600 #Operacion Matricial
593 601 # velUVW = numpy.zeros((velRadial.shape[1],3))
594 602 # for ind in range(velRadial.shape[1]):
595 603 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
596 604 # velUVW = velUVW.transpose()
597 605 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
598 606 velUVW[:,:] = numpy.dot(A,velRadial)
599 607
600 608
601 609 return velUVW
602 610
603 611 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
604 612
605 613 def techniqueDBS(self, kwargs):
606 614 """
607 615 Function that implements Doppler Beam Swinging (DBS) technique.
608 616
609 617 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
610 618 Direction correction (if necessary), Ranges and SNR
611 619
612 620 Output: Winds estimation (Zonal, Meridional and Vertical)
613 621
614 622 Parameters affected: Winds, height range, SNR
615 623 """
616 624 velRadial0 = kwargs['velRadial']
617 625 heiRang = kwargs['heightList']
618 626 SNR0 = kwargs['SNR']
619 627
620 628 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
621 629 theta_x = numpy.array(kwargs['dirCosx'])
622 630 theta_y = numpy.array(kwargs['dirCosy'])
623 631 else:
624 632 elev = numpy.array(kwargs['elevation'])
625 633 azim = numpy.array(kwargs['azimuth'])
626 634 theta_x, theta_y = self.__calculateCosDir(elev, azim)
627 635 azimuth = kwargs['correctAzimuth']
628 636 if kwargs.has_key('horizontalOnly'):
629 637 horizontalOnly = kwargs['horizontalOnly']
630 638 else: horizontalOnly = False
631 639 if kwargs.has_key('correctFactor'):
632 640 correctFactor = kwargs['correctFactor']
633 641 else: correctFactor = 1
634 642 if kwargs.has_key('channelList'):
635 643 channelList = kwargs['channelList']
636 644 if len(channelList) == 2:
637 645 horizontalOnly = True
638 646 arrayChannel = numpy.array(channelList)
639 647 param = param[arrayChannel,:,:]
640 648 theta_x = theta_x[arrayChannel]
641 649 theta_y = theta_y[arrayChannel]
642 650
643 651 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
644 652 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
645 653 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
646 654
647 655 #Calculo de Componentes de la velocidad con DBS
648 656 winds = self.__calculateVelUVW(A,velRadial1)
649 657
650 658 return winds, heiRang1, SNR1
651 659
652 660 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
653 661
654 662 nPairs = len(pairs_ccf)
655 663 posx = numpy.asarray(posx)
656 664 posy = numpy.asarray(posy)
657 665
658 666 #Rotacion Inversa para alinear con el azimuth
659 667 if azimuth!= None:
660 668 azimuth = azimuth*math.pi/180
661 669 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
662 670 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
663 671 else:
664 672 posx1 = posx
665 673 posy1 = posy
666 674
667 675 #Calculo de Distancias
668 676 distx = numpy.zeros(nPairs)
669 677 disty = numpy.zeros(nPairs)
670 678 dist = numpy.zeros(nPairs)
671 679 ang = numpy.zeros(nPairs)
672 680
673 681 for i in range(nPairs):
674 682 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
675 683 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
676 684 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
677 685 ang[i] = numpy.arctan2(disty[i],distx[i])
678 686
679 687 return distx, disty, dist, ang
680 688 #Calculo de Matrices
681 689 # nPairs = len(pairs)
682 690 # ang1 = numpy.zeros((nPairs, 2, 1))
683 691 # dist1 = numpy.zeros((nPairs, 2, 1))
684 692 #
685 693 # for j in range(nPairs):
686 694 # dist1[j,0,0] = dist[pairs[j][0]]
687 695 # dist1[j,1,0] = dist[pairs[j][1]]
688 696 # ang1[j,0,0] = ang[pairs[j][0]]
689 697 # ang1[j,1,0] = ang[pairs[j][1]]
690 698 #
691 699 # return distx,disty, dist1,ang1
692 700
693 701
694 702 def __calculateVelVer(self, phase, lagTRange, _lambda):
695 703
696 704 Ts = lagTRange[1] - lagTRange[0]
697 705 velW = -_lambda*phase/(4*math.pi*Ts)
698 706
699 707 return velW
700 708
701 709 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
702 710 nPairs = tau1.shape[0]
703 711 nHeights = tau1.shape[1]
704 712 vel = numpy.zeros((nPairs,3,nHeights))
705 713 dist1 = numpy.reshape(dist, (dist.size,1))
706 714
707 715 angCos = numpy.cos(ang)
708 716 angSin = numpy.sin(ang)
709 717
710 718 vel0 = dist1*tau1/(2*tau2**2)
711 719 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
712 720 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
713 721
714 722 ind = numpy.where(numpy.isinf(vel))
715 723 vel[ind] = numpy.nan
716 724
717 725 return vel
718 726
719 727 # def __getPairsAutoCorr(self, pairsList, nChannels):
720 728 #
721 729 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
722 730 #
723 731 # for l in range(len(pairsList)):
724 732 # firstChannel = pairsList[l][0]
725 733 # secondChannel = pairsList[l][1]
726 734 #
727 735 # #Obteniendo pares de Autocorrelacion
728 736 # if firstChannel == secondChannel:
729 737 # pairsAutoCorr[firstChannel] = int(l)
730 738 #
731 739 # pairsAutoCorr = pairsAutoCorr.astype(int)
732 740 #
733 741 # pairsCrossCorr = range(len(pairsList))
734 742 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
735 743 #
736 744 # return pairsAutoCorr, pairsCrossCorr
737 745
738 746 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
739 747 def techniqueSA(self, kwargs):
740 748
741 749 """
742 750 Function that implements Spaced Antenna (SA) technique.
743 751
744 752 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
745 753 Direction correction (if necessary), Ranges and SNR
746 754
747 755 Output: Winds estimation (Zonal, Meridional and Vertical)
748 756
749 757 Parameters affected: Winds
750 758 """
751 759 position_x = kwargs['positionX']
752 760 position_y = kwargs['positionY']
753 761 azimuth = kwargs['azimuth']
754 762
755 763 if kwargs.has_key('correctFactor'):
756 764 correctFactor = kwargs['correctFactor']
757 765 else:
758 766 correctFactor = 1
759 767
760 768 groupList = kwargs['groupList']
761 769 pairs_ccf = groupList[1]
762 770 tau = kwargs['tau']
763 771 _lambda = kwargs['_lambda']
764 772
765 773 #Cross Correlation pairs obtained
766 774 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
767 775 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
768 776 # pairsSelArray = numpy.array(pairsSelected)
769 777 # pairs = []
770 778 #
771 779 # #Wind estimation pairs obtained
772 780 # for i in range(pairsSelArray.shape[0]/2):
773 781 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
774 782 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
775 783 # pairs.append((ind1,ind2))
776 784
777 785 indtau = tau.shape[0]/2
778 786 tau1 = tau[:indtau,:]
779 787 tau2 = tau[indtau:-1,:]
780 788 # tau1 = tau1[pairs,:]
781 789 # tau2 = tau2[pairs,:]
782 790 phase1 = tau[-1,:]
783 791
784 792 #---------------------------------------------------------------------
785 793 #Metodo Directo
786 794 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
787 795 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
788 796 winds = stats.nanmean(winds, axis=0)
789 797 #---------------------------------------------------------------------
790 798 #Metodo General
791 799 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
792 800 # #Calculo Coeficientes de Funcion de Correlacion
793 801 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
794 802 # #Calculo de Velocidades
795 803 # winds = self.calculateVelUV(F,G,A,B,H)
796 804
797 805 #---------------------------------------------------------------------
798 806 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
799 807 winds = correctFactor*winds
800 808 return winds
801 809
802 810 def __checkTime(self, currentTime, paramInterval, outputInterval):
803 811
804 812 dataTime = currentTime + paramInterval
805 813 deltaTime = dataTime - self.__initime
806 814
807 815 if deltaTime >= outputInterval or deltaTime < 0:
808 816 self.__dataReady = True
809 817 return
810 818
811 819 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax, binkm=2):
812 820 '''
813 821 Function that implements winds estimation technique with detected meteors.
814 822
815 823 Input: Detected meteors, Minimum meteor quantity to wind estimation
816 824
817 825 Output: Winds estimation (Zonal and Meridional)
818 826
819 827 Parameters affected: Winds
820 828 '''
821 829 # print arrayMeteor.shape
822 830 #Settings
823 831 nInt = (heightMax - heightMin)/binkm
824 832 # print nInt
825 833 nInt = int(nInt)
826 834 # print nInt
827 835 winds = numpy.zeros((2,nInt))*numpy.nan
828 836
829 837 #Filter errors
830 838 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
831 839 finalMeteor = arrayMeteor[error,:]
832 840
833 841 #Meteor Histogram
834 842 finalHeights = finalMeteor[:,2]
835 843 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
836 844 nMeteorsPerI = hist[0]
837 845 heightPerI = hist[1]
838 846
839 847 #Sort of meteors
840 848 indSort = finalHeights.argsort()
841 849 finalMeteor2 = finalMeteor[indSort,:]
842 850
843 851 # Calculating winds
844 852 ind1 = 0
845 853 ind2 = 0
846 854
847 855 for i in range(nInt):
848 856 nMet = nMeteorsPerI[i]
849 857 ind1 = ind2
850 858 ind2 = ind1 + nMet
851 859
852 860 meteorAux = finalMeteor2[ind1:ind2,:]
853 861
854 862 if meteorAux.shape[0] >= meteorThresh:
855 863 vel = meteorAux[:, 6]
856 864 zen = meteorAux[:, 4]*numpy.pi/180
857 865 azim = meteorAux[:, 3]*numpy.pi/180
858 866
859 867 n = numpy.cos(zen)
860 868 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
861 869 # l = m*numpy.tan(azim)
862 870 l = numpy.sin(zen)*numpy.sin(azim)
863 871 m = numpy.sin(zen)*numpy.cos(azim)
864 872
865 873 A = numpy.vstack((l, m)).transpose()
866 874 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
867 875 windsAux = numpy.dot(A1, vel)
868 876
869 877 winds[0,i] = windsAux[0]
870 878 winds[1,i] = windsAux[1]
871 879
872 880 return winds, heightPerI[:-1]
873 881
874 882 def techniqueNSM_SA(self, **kwargs):
875 883 metArray = kwargs['metArray']
876 884 heightList = kwargs['heightList']
877 885 timeList = kwargs['timeList']
878 886
879 887 rx_location = kwargs['rx_location']
880 888 groupList = kwargs['groupList']
881 889 azimuth = kwargs['azimuth']
882 890 dfactor = kwargs['dfactor']
883 891 k = kwargs['k']
884 892
885 893 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
886 894 d = dist*dfactor
887 895 #Phase calculation
888 896 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
889 897
890 898 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
891 899
892 900 velEst = numpy.zeros((heightList.size,2))*numpy.nan
893 901 azimuth1 = azimuth1*numpy.pi/180
894 902
895 903 for i in range(heightList.size):
896 904 h = heightList[i]
897 905 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
898 906 metHeight = metArray1[indH,:]
899 907 if metHeight.shape[0] >= 2:
900 908 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
901 909 iazim = metHeight[:,1].astype(int)
902 910 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
903 911 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
904 912 A = numpy.asmatrix(A)
905 913 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
906 914 velHor = numpy.dot(A1,velAux)
907 915
908 916 velEst[i,:] = numpy.squeeze(velHor)
909 917 return velEst
910 918
911 919 def __getPhaseSlope(self, metArray, heightList, timeList):
912 920 meteorList = []
913 921 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
914 922 #Putting back together the meteor matrix
915 923 utctime = metArray[:,0]
916 924 uniqueTime = numpy.unique(utctime)
917 925
918 926 phaseDerThresh = 0.5
919 927 ippSeconds = timeList[1] - timeList[0]
920 928 sec = numpy.where(timeList>1)[0][0]
921 929 nPairs = metArray.shape[1] - 6
922 930 nHeights = len(heightList)
923 931
924 932 for t in uniqueTime:
925 933 metArray1 = metArray[utctime==t,:]
926 934 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
927 935 tmet = metArray1[:,1].astype(int)
928 936 hmet = metArray1[:,2].astype(int)
929 937
930 938 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
931 939 metPhase[:,:] = numpy.nan
932 940 metPhase[:,hmet,tmet] = metArray1[:,6:].T
933 941
934 942 #Delete short trails
935 943 metBool = ~numpy.isnan(metPhase[0,:,:])
936 944 heightVect = numpy.sum(metBool, axis = 1)
937 945 metBool[heightVect<sec,:] = False
938 946 metPhase[:,heightVect<sec,:] = numpy.nan
939 947
940 948 #Derivative
941 949 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
942 950 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
943 951 metPhase[phDerAux] = numpy.nan
944 952
945 953 #--------------------------METEOR DETECTION -----------------------------------------
946 954 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
947 955
948 956 for p in numpy.arange(nPairs):
949 957 phase = metPhase[p,:,:]
950 958 phDer = metDer[p,:,:]
951 959
952 960 for h in indMet:
953 961 height = heightList[h]
954 962 phase1 = phase[h,:] #82
955 963 phDer1 = phDer[h,:]
956 964
957 965 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
958 966
959 967 indValid = numpy.where(~numpy.isnan(phase1))[0]
960 968 initMet = indValid[0]
961 969 endMet = 0
962 970
963 971 for i in range(len(indValid)-1):
964 972
965 973 #Time difference
966 974 inow = indValid[i]
967 975 inext = indValid[i+1]
968 976 idiff = inext - inow
969 977 #Phase difference
970 978 phDiff = numpy.abs(phase1[inext] - phase1[inow])
971 979
972 980 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
973 981 sizeTrail = inow - initMet + 1
974 982 if sizeTrail>3*sec: #Too short meteors
975 983 x = numpy.arange(initMet,inow+1)*ippSeconds
976 984 y = phase1[initMet:inow+1]
977 985 ynnan = ~numpy.isnan(y)
978 986 x = x[ynnan]
979 987 y = y[ynnan]
980 988 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
981 989 ylin = x*slope + intercept
982 990 rsq = r_value**2
983 991 if rsq > 0.5:
984 992 vel = slope#*height*1000/(k*d)
985 993 estAux = numpy.array([utctime,p,height, vel, rsq])
986 994 meteorList.append(estAux)
987 995 initMet = inext
988 996 metArray2 = numpy.array(meteorList)
989 997
990 998 return metArray2
991 999
992 1000 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
993 1001
994 1002 azimuth1 = numpy.zeros(len(pairslist))
995 1003 dist = numpy.zeros(len(pairslist))
996 1004
997 1005 for i in range(len(rx_location)):
998 1006 ch0 = pairslist[i][0]
999 1007 ch1 = pairslist[i][1]
1000 1008
1001 1009 diffX = rx_location[ch0][0] - rx_location[ch1][0]
1002 1010 diffY = rx_location[ch0][1] - rx_location[ch1][1]
1003 1011 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
1004 1012 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
1005 1013
1006 1014 azimuth1 -= azimuth0
1007 1015 return azimuth1, dist
1008 1016
1009 1017 def techniqueNSM_DBS(self, **kwargs):
1010 1018 metArray = kwargs['metArray']
1011 1019 heightList = kwargs['heightList']
1012 1020 timeList = kwargs['timeList']
1013 1021 zenithList = kwargs['zenithList']
1014 1022 nChan = numpy.max(cmet) + 1
1015 1023 nHeights = len(heightList)
1016 1024
1017 1025 utctime = metArray[:,0]
1018 1026 cmet = metArray[:,1]
1019 1027 hmet = metArray1[:,3].astype(int)
1020 1028 h1met = heightList[hmet]*zenithList[cmet]
1021 1029 vmet = metArray1[:,5]
1022 1030
1023 1031 for i in range(nHeights - 1):
1024 1032 hmin = heightList[i]
1025 1033 hmax = heightList[i + 1]
1026 1034
1027 1035 vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
1028 1036
1029 1037
1030 1038
1031 1039 return data_output
1032 1040
1033 1041 def run(self, dataOut, technique, **kwargs):
1034 1042
1035 1043 param = dataOut.data_param
1036 1044 if dataOut.abscissaList != None:
1037 1045 absc = dataOut.abscissaList[:-1]
1038 1046 noise = dataOut.noise
1039 1047 heightList = dataOut.heightList
1040 1048 SNR = dataOut.data_SNR
1041 1049
1042 1050 if technique == 'DBS':
1043 1051
1044 1052 kwargs['velRadial'] = param[:,1,:] #Radial velocity
1045 1053 kwargs['heightList'] = heightList
1046 1054 kwargs['SNR'] = SNR
1047 1055
1048 1056 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
1049 1057 dataOut.utctimeInit = dataOut.utctime
1050 1058 dataOut.outputInterval = dataOut.paramInterval
1051 1059
1052 1060 elif technique == 'SA':
1053 1061
1054 1062 #Parameters
1055 1063 # position_x = kwargs['positionX']
1056 1064 # position_y = kwargs['positionY']
1057 1065 # azimuth = kwargs['azimuth']
1058 1066 #
1059 1067 # if kwargs.has_key('crosspairsList'):
1060 1068 # pairs = kwargs['crosspairsList']
1061 1069 # else:
1062 1070 # pairs = None
1063 1071 #
1064 1072 # if kwargs.has_key('correctFactor'):
1065 1073 # correctFactor = kwargs['correctFactor']
1066 1074 # else:
1067 1075 # correctFactor = 1
1068 1076
1069 1077 # tau = dataOut.data_param
1070 1078 # _lambda = dataOut.C/dataOut.frequency
1071 1079 # pairsList = dataOut.groupList
1072 1080 # nChannels = dataOut.nChannels
1073 1081
1074 1082 kwargs['groupList'] = dataOut.groupList
1075 1083 kwargs['tau'] = dataOut.data_param
1076 1084 kwargs['_lambda'] = dataOut.C/dataOut.frequency
1077 1085 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1078 1086 dataOut.data_output = self.techniqueSA(kwargs)
1079 1087 dataOut.utctimeInit = dataOut.utctime
1080 1088 dataOut.outputInterval = dataOut.timeInterval
1081 1089
1082 1090 elif technique == 'Meteors':
1083 1091 dataOut.flagNoData = True
1084 1092 self.__dataReady = False
1085 1093
1086 1094 if kwargs.has_key('nHours'):
1087 1095 nHours = kwargs['nHours']
1088 1096 else:
1089 1097 nHours = 1
1090 1098
1091 1099 if kwargs.has_key('meteorsPerBin'):
1092 1100 meteorThresh = kwargs['meteorsPerBin']
1093 1101 else:
1094 1102 meteorThresh = 6
1095 1103
1096 1104 if kwargs.has_key('hmin'):
1097 1105 hmin = kwargs['hmin']
1098 1106 else: hmin = 70
1099 1107 if kwargs.has_key('hmax'):
1100 1108 hmax = kwargs['hmax']
1101 1109 else: hmax = 110
1102 1110
1103 1111 if kwargs.has_key('BinKm'):
1104 1112 binkm = kwargs['BinKm']
1105 1113 else:
1106 1114 binkm = 2
1107 1115
1108 1116 dataOut.outputInterval = nHours*3600
1109 1117
1110 1118 if self.__isConfig == False:
1111 1119 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1112 1120 #Get Initial LTC time
1113 1121 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1114 1122 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1115 1123
1116 1124 self.__isConfig = True
1117 1125
1118 1126 if self.__buffer is None:
1119 1127 self.__buffer = dataOut.data_param
1120 1128 self.__firstdata = copy.copy(dataOut)
1121 1129
1122 1130 else:
1123 1131 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1124 1132
1125 1133 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1126 1134
1127 1135 if self.__dataReady:
1128 1136 dataOut.utctimeInit = self.__initime
1129 1137
1130 1138 self.__initime += dataOut.outputInterval #to erase time offset
1131 1139
1132 1140 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax, binkm)
1133 1141 dataOut.flagNoData = False
1134 1142 self.__buffer = None
1135 1143
1136 1144 elif technique == 'Meteors1':
1137 1145 dataOut.flagNoData = True
1138 1146 self.__dataReady = False
1139 1147
1140 1148 if kwargs.has_key('nMins'):
1141 1149 nMins = kwargs['nMins']
1142 1150 else: nMins = 20
1143 1151 if kwargs.has_key('rx_location'):
1144 1152 rx_location = kwargs['rx_location']
1145 1153 else: rx_location = [(0,1),(1,1),(1,0)]
1146 1154 if kwargs.has_key('azimuth'):
1147 1155 azimuth = kwargs['azimuth']
1148 1156 else: azimuth = 51
1149 1157 if kwargs.has_key('dfactor'):
1150 1158 dfactor = kwargs['dfactor']
1151 1159 if kwargs.has_key('mode'):
1152 1160 mode = kwargs['mode']
1153 1161 else: mode = 'SA'
1154 1162
1155 1163 #Borrar luego esto
1156 1164 if dataOut.groupList is None:
1157 1165 dataOut.groupList = [(0,1),(0,2),(1,2)]
1158 1166 groupList = dataOut.groupList
1159 1167 C = 3e8
1160 1168 freq = 50e6
1161 1169 lamb = C/freq
1162 1170 k = 2*numpy.pi/lamb
1163 1171
1164 1172 timeList = dataOut.abscissaList
1165 1173 heightList = dataOut.heightList
1166 1174
1167 1175 if self.__isConfig == False:
1168 1176 dataOut.outputInterval = nMins*60
1169 1177 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1170 1178 #Get Initial LTC time
1171 1179 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1172 1180 minuteAux = initime.minute
1173 1181 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
1174 1182 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1175 1183
1176 1184 self.__isConfig = True
1177 1185
1178 1186 if self.__buffer is None:
1179 1187 self.__buffer = dataOut.data_param
1180 1188 self.__firstdata = copy.copy(dataOut)
1181 1189
1182 1190 else:
1183 1191 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1184 1192
1185 1193 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1186 1194
1187 1195 if self.__dataReady:
1188 1196 dataOut.utctimeInit = self.__initime
1189 1197 self.__initime += dataOut.outputInterval #to erase time offset
1190 1198
1191 1199 metArray = self.__buffer
1192 1200 if mode == 'SA':
1193 1201 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
1194 1202 elif mode == 'DBS':
1195 1203 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
1196 1204 dataOut.data_output = dataOut.data_output.T
1197 1205 dataOut.flagNoData = False
1198 1206 self.__buffer = None
1199 1207
1200 1208 return
1201 1209
1202 1210 class EWDriftsEstimation(Operation):
1203 1211
1204 1212
1205 1213 def __correctValues(self, heiRang, phi, velRadial, SNR):
1206 1214 listPhi = phi.tolist()
1207 1215 maxid = listPhi.index(max(listPhi))
1208 1216 minid = listPhi.index(min(listPhi))
1209 1217
1210 1218 rango = range(len(phi))
1211 1219 # rango = numpy.delete(rango,maxid)
1212 1220
1213 1221 heiRang1 = heiRang*math.cos(phi[maxid])
1214 1222 heiRangAux = heiRang*math.cos(phi[minid])
1215 1223 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1216 1224 heiRang1 = numpy.delete(heiRang1,indOut)
1217 1225
1218 1226 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1219 1227 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1220 1228
1221 1229 for i in rango:
1222 1230 x = heiRang*math.cos(phi[i])
1223 1231 y1 = velRadial[i,:]
1224 1232 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1225 1233
1226 1234 x1 = heiRang1
1227 1235 y11 = f1(x1)
1228 1236
1229 1237 y2 = SNR[i,:]
1230 1238 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1231 1239 y21 = f2(x1)
1232 1240
1233 1241 velRadial1[i,:] = y11
1234 1242 SNR1[i,:] = y21
1235 1243
1236 1244 return heiRang1, velRadial1, SNR1
1237 1245
1238 1246 def run(self, dataOut, zenith, zenithCorrection):
1239 1247 heiRang = dataOut.heightList
1240 1248 velRadial = dataOut.data_param[:,3,:]
1241 1249 SNR = dataOut.data_SNR
1242 1250
1243 1251 zenith = numpy.array(zenith)
1244 1252 zenith -= zenithCorrection
1245 1253 zenith *= numpy.pi/180
1246 1254
1247 1255 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1248 1256
1249 1257 alp = zenith[0]
1250 1258 bet = zenith[1]
1251 1259
1252 1260 w_w = velRadial1[0,:]
1253 1261 w_e = velRadial1[1,:]
1254 1262
1255 1263 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1256 1264 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1257 1265
1258 1266 winds = numpy.vstack((u,w))
1259 1267
1260 1268 dataOut.heightList = heiRang1
1261 1269 dataOut.data_output = winds
1262 1270 dataOut.data_SNR = SNR1
1263 1271
1264 1272 dataOut.utctimeInit = dataOut.utctime
1265 1273 dataOut.outputInterval = dataOut.timeInterval
1266 1274 return
1267 1275
1268 1276 #--------------- Non Specular Meteor ----------------
1269 1277
1270 1278 class NonSpecularMeteorDetection(Operation):
1271 1279
1272 1280 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
1273 1281 data_acf = self.dataOut.data_pre[0]
1274 1282 data_ccf = self.dataOut.data_pre[1]
1275 1283
1276 1284 lamb = self.dataOut.C/self.dataOut.frequency
1277 1285 tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
1278 1286 paramInterval = self.dataOut.paramInterval
1279 1287
1280 1288 nChannels = data_acf.shape[0]
1281 1289 nLags = data_acf.shape[1]
1282 1290 nProfiles = data_acf.shape[2]
1283 1291 nHeights = self.dataOut.nHeights
1284 1292 nCohInt = self.dataOut.nCohInt
1285 1293 sec = numpy.round(nProfiles/self.dataOut.paramInterval)
1286 1294 heightList = self.dataOut.heightList
1287 1295 ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
1288 1296 utctime = self.dataOut.utctime
1289 1297
1290 1298 self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
1291 1299
1292 1300 #------------------------ SNR --------------------------------------
1293 1301 power = data_acf[:,0,:,:].real
1294 1302 noise = numpy.zeros(nChannels)
1295 1303 SNR = numpy.zeros(power.shape)
1296 1304 for i in range(nChannels):
1297 1305 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
1298 1306 SNR[i] = (power[i]-noise[i])/noise[i]
1299 1307 SNRm = numpy.nanmean(SNR, axis = 0)
1300 1308 SNRdB = 10*numpy.log10(SNR)
1301 1309
1302 1310 if mode == 'SA':
1303 1311 nPairs = data_ccf.shape[0]
1304 1312 #---------------------- Coherence and Phase --------------------------
1305 1313 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
1306 1314 # phase1 = numpy.copy(phase)
1307 1315 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
1308 1316
1309 1317 for p in range(nPairs):
1310 1318 ch0 = self.dataOut.groupList[p][0]
1311 1319 ch1 = self.dataOut.groupList[p][1]
1312 1320 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
1313 1321 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
1314 1322 # phase1[p,:,:] = numpy.angle(ccf) #median filter
1315 1323 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
1316 1324 # coh1[p,:,:] = numpy.abs(ccf) #median filter
1317 1325 coh = numpy.nanmax(coh1, axis = 0)
1318 1326 # struc = numpy.ones((5,1))
1319 1327 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
1320 1328 #---------------------- Radial Velocity ----------------------------
1321 1329 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
1322 1330 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
1323 1331
1324 1332 if allData:
1325 1333 boolMetFin = ~numpy.isnan(SNRm)
1326 1334 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1327 1335 else:
1328 1336 #------------------------ Meteor mask ---------------------------------
1329 1337 # #SNR mask
1330 1338 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
1331 1339 #
1332 1340 # #Erase small objects
1333 1341 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
1334 1342 #
1335 1343 # auxEEJ = numpy.sum(boolMet1,axis=0)
1336 1344 # indOver = auxEEJ>nProfiles*0.8 #Use this later
1337 1345 # indEEJ = numpy.where(indOver)[0]
1338 1346 # indNEEJ = numpy.where(~indOver)[0]
1339 1347 #
1340 1348 # boolMetFin = boolMet1
1341 1349 #
1342 1350 # if indEEJ.size > 0:
1343 1351 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
1344 1352 #
1345 1353 # boolMet2 = coh > cohThresh
1346 1354 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
1347 1355 #
1348 1356 # #Final Meteor mask
1349 1357 # boolMetFin = boolMet1|boolMet2
1350 1358
1351 1359 #Coherence mask
1352 1360 boolMet1 = coh > 0.75
1353 1361 struc = numpy.ones((30,1))
1354 1362 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
1355 1363
1356 1364 #Derivative mask
1357 1365 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1358 1366 boolMet2 = derPhase < 0.2
1359 1367 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
1360 1368 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
1361 1369 boolMet2 = ndimage.median_filter(boolMet2,size=5)
1362 1370 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
1363 1371 # #Final mask
1364 1372 # boolMetFin = boolMet2
1365 1373 boolMetFin = boolMet1&boolMet2
1366 1374 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
1367 1375 #Creating data_param
1368 1376 coordMet = numpy.where(boolMetFin)
1369 1377
1370 1378 tmet = coordMet[0]
1371 1379 hmet = coordMet[1]
1372 1380
1373 1381 data_param = numpy.zeros((tmet.size, 6 + nPairs))
1374 1382 data_param[:,0] = utctime
1375 1383 data_param[:,1] = tmet
1376 1384 data_param[:,2] = hmet
1377 1385 data_param[:,3] = SNRm[tmet,hmet]
1378 1386 data_param[:,4] = velRad[tmet,hmet]
1379 1387 data_param[:,5] = coh[tmet,hmet]
1380 1388 data_param[:,6:] = phase[:,tmet,hmet].T
1381 1389
1382 1390 elif mode == 'DBS':
1383 1391 self.dataOut.groupList = numpy.arange(nChannels)
1384 1392
1385 1393 #Radial Velocities
1386 1394 # phase = numpy.angle(data_acf[:,1,:,:])
1387 1395 phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
1388 1396 velRad = phase*lamb/(4*numpy.pi*tSamp)
1389 1397
1390 1398 #Spectral width
1391 1399 acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
1392 1400 acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
1393 1401
1394 1402 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
1395 1403 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
1396 1404 if allData:
1397 1405 boolMetFin = ~numpy.isnan(SNRdB)
1398 1406 else:
1399 1407 #SNR
1400 1408 boolMet1 = (SNRdB>SNRthresh) #SNR mask
1401 1409 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
1402 1410
1403 1411 #Radial velocity
1404 1412 boolMet2 = numpy.abs(velRad) < 30
1405 1413 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
1406 1414
1407 1415 #Spectral Width
1408 1416 boolMet3 = spcWidth < 30
1409 1417 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
1410 1418 # boolMetFin = self.__erase_small(boolMet1, 10,5)
1411 1419 boolMetFin = boolMet1&boolMet2&boolMet3
1412 1420
1413 1421 #Creating data_param
1414 1422 coordMet = numpy.where(boolMetFin)
1415 1423
1416 1424 cmet = coordMet[0]
1417 1425 tmet = coordMet[1]
1418 1426 hmet = coordMet[2]
1419 1427
1420 1428 data_param = numpy.zeros((tmet.size, 7))
1421 1429 data_param[:,0] = utctime
1422 1430 data_param[:,1] = cmet
1423 1431 data_param[:,2] = tmet
1424 1432 data_param[:,3] = hmet
1425 1433 data_param[:,4] = SNR[cmet,tmet,hmet].T
1426 1434 data_param[:,5] = velRad[cmet,tmet,hmet].T
1427 1435 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
1428 1436
1429 1437 # self.dataOut.data_param = data_int
1430 1438 if len(data_param) == 0:
1431 1439 self.dataOut.flagNoData = True
1432 1440 else:
1433 1441 self.dataOut.data_param = data_param
1434 1442
1435 1443 def __erase_small(self, binArray, threshX, threshY):
1436 1444 labarray, numfeat = ndimage.measurements.label(binArray)
1437 1445 binArray1 = numpy.copy(binArray)
1438 1446
1439 1447 for i in range(1,numfeat + 1):
1440 1448 auxBin = (labarray==i)
1441 1449 auxSize = auxBin.sum()
1442 1450
1443 1451 x,y = numpy.where(auxBin)
1444 1452 widthX = x.max() - x.min()
1445 1453 widthY = y.max() - y.min()
1446 1454
1447 1455 #width X: 3 seg -> 12.5*3
1448 1456 #width Y:
1449 1457
1450 1458 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
1451 1459 binArray1[auxBin] = False
1452 1460
1453 1461 return binArray1
1454 1462
1455 1463 #--------------- Specular Meteor ----------------
1456 1464
1457 1465 class SMDetection(Operation):
1458 1466 '''
1459 1467 Function DetectMeteors()
1460 1468 Project developed with paper:
1461 1469 HOLDSWORTH ET AL. 2004
1462 1470
1463 1471 Input:
1464 1472 self.dataOut.data_pre
1465 1473
1466 1474 centerReceiverIndex: From the channels, which is the center receiver
1467 1475
1468 1476 hei_ref: Height reference for the Beacon signal extraction
1469 1477 tauindex:
1470 1478 predefinedPhaseShifts: Predefined phase offset for the voltge signals
1471 1479
1472 1480 cohDetection: Whether to user Coherent detection or not
1473 1481 cohDet_timeStep: Coherent Detection calculation time step
1474 1482 cohDet_thresh: Coherent Detection phase threshold to correct phases
1475 1483
1476 1484 noise_timeStep: Noise calculation time step
1477 1485 noise_multiple: Noise multiple to define signal threshold
1478 1486
1479 1487 multDet_timeLimit: Multiple Detection Removal time limit in seconds
1480 1488 multDet_rangeLimit: Multiple Detection Removal range limit in km
1481 1489
1482 1490 phaseThresh: Maximum phase difference between receiver to be consider a meteor
1483 1491 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
1484 1492
1485 1493 hmin: Minimum Height of the meteor to use it in the further wind estimations
1486 1494 hmax: Maximum Height of the meteor to use it in the further wind estimations
1487 1495 azimuth: Azimuth angle correction
1488 1496
1489 1497 Affected:
1490 1498 self.dataOut.data_param
1491 1499
1492 1500 Rejection Criteria (Errors):
1493 1501 0: No error; analysis OK
1494 1502 1: SNR < SNR threshold
1495 1503 2: angle of arrival (AOA) ambiguously determined
1496 1504 3: AOA estimate not feasible
1497 1505 4: Large difference in AOAs obtained from different antenna baselines
1498 1506 5: echo at start or end of time series
1499 1507 6: echo less than 5 examples long; too short for analysis
1500 1508 7: echo rise exceeds 0.3s
1501 1509 8: echo decay time less than twice rise time
1502 1510 9: large power level before echo
1503 1511 10: large power level after echo
1504 1512 11: poor fit to amplitude for estimation of decay time
1505 1513 12: poor fit to CCF phase variation for estimation of radial drift velocity
1506 1514 13: height unresolvable echo: not valid height within 70 to 110 km
1507 1515 14: height ambiguous echo: more then one possible height within 70 to 110 km
1508 1516 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
1509 1517 16: oscilatory echo, indicating event most likely not an underdense echo
1510 1518
1511 1519 17: phase difference in meteor Reestimation
1512 1520
1513 1521 Data Storage:
1514 1522 Meteors for Wind Estimation (8):
1515 1523 Utc Time | Range Height
1516 1524 Azimuth Zenith errorCosDir
1517 1525 VelRad errorVelRad
1518 1526 Phase0 Phase1 Phase2 Phase3
1519 1527 TypeError
1520 1528
1521 1529 '''
1522 1530
1523 1531 def run(self, dataOut, hei_ref = None, tauindex = 0,
1524 1532 phaseOffsets = None,
1525 1533 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
1526 1534 noise_timeStep = 4, noise_multiple = 4,
1527 1535 multDet_timeLimit = 1, multDet_rangeLimit = 3,
1528 1536 phaseThresh = 20, SNRThresh = 5,
1529 1537 hmin = 50, hmax=150, azimuth = 0,
1530 1538 channelPositions = None) :
1531 1539
1532 1540
1533 1541 #Getting Pairslist
1534 1542 if channelPositions is None:
1535 1543 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
1536 1544 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
1537 1545 meteorOps = SMOperations()
1538 1546 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
1539 1547 heiRang = dataOut.getHeiRange()
1540 1548 #Get Beacon signal - No Beacon signal anymore
1541 1549 # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
1542 1550 #
1543 1551 # if hei_ref != None:
1544 1552 # newheis = numpy.where(self.dataOut.heightList>hei_ref)
1545 1553 #
1546 1554
1547 1555
1548 1556 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
1549 1557 # see if the user put in pre defined phase shifts
1550 1558 voltsPShift = dataOut.data_pre.copy()
1551 1559
1552 1560 # if predefinedPhaseShifts != None:
1553 1561 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
1554 1562 #
1555 1563 # # elif beaconPhaseShifts:
1556 1564 # # #get hardware phase shifts using beacon signal
1557 1565 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
1558 1566 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
1559 1567 #
1560 1568 # else:
1561 1569 # hardwarePhaseShifts = numpy.zeros(5)
1562 1570 #
1563 1571 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
1564 1572 # for i in range(self.dataOut.data_pre.shape[0]):
1565 1573 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
1566 1574
1567 1575 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
1568 1576
1569 1577 #Remove DC
1570 1578 voltsDC = numpy.mean(voltsPShift,1)
1571 1579 voltsDC = numpy.mean(voltsDC,1)
1572 1580 for i in range(voltsDC.shape[0]):
1573 1581 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
1574 1582
1575 1583 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
1576 1584 # voltsPShift = voltsPShift[:,:,:newheis[0][0]]
1577 1585
1578 1586 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
1579 1587 #Coherent Detection
1580 1588 if cohDetection:
1581 1589 #use coherent detection to get the net power
1582 1590 cohDet_thresh = cohDet_thresh*numpy.pi/180
1583 1591 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
1584 1592
1585 1593 #Non-coherent detection!
1586 1594 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
1587 1595 #********** END OF COH/NON-COH POWER CALCULATION**********************
1588 1596
1589 1597 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
1590 1598 #Get noise
1591 1599 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
1592 1600 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
1593 1601 #Get signal threshold
1594 1602 signalThresh = noise_multiple*noise
1595 1603 #Meteor echoes detection
1596 1604 listMeteors = self.__findMeteors(powerNet, signalThresh)
1597 1605 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
1598 1606
1599 1607 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
1600 1608 #Parameters
1601 1609 heiRange = dataOut.getHeiRange()
1602 1610 rangeInterval = heiRange[1] - heiRange[0]
1603 1611 rangeLimit = multDet_rangeLimit/rangeInterval
1604 1612 timeLimit = multDet_timeLimit/dataOut.timeInterval
1605 1613 #Multiple detection removals
1606 1614 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
1607 1615 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
1608 1616
1609 1617 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
1610 1618 #Parameters
1611 1619 phaseThresh = phaseThresh*numpy.pi/180
1612 1620 thresh = [phaseThresh, noise_multiple, SNRThresh]
1613 1621 #Meteor reestimation (Errors N 1, 6, 12, 17)
1614 1622 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
1615 1623 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
1616 1624 #Estimation of decay times (Errors N 7, 8, 11)
1617 1625 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
1618 1626 #******************* END OF METEOR REESTIMATION *******************
1619 1627
1620 1628 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
1621 1629 #Calculating Radial Velocity (Error N 15)
1622 1630 radialStdThresh = 10
1623 1631 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
1624 1632
1625 1633 if len(listMeteors4) > 0:
1626 1634 #Setting New Array
1627 1635 date = dataOut.utctime
1628 1636 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
1629 1637
1630 1638 #Correcting phase offset
1631 1639 if phaseOffsets != None:
1632 1640 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
1633 1641 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
1634 1642
1635 1643 #Second Pairslist
1636 1644 pairsList = []
1637 1645 pairx = (0,1)
1638 1646 pairy = (2,3)
1639 1647 pairsList.append(pairx)
1640 1648 pairsList.append(pairy)
1641 1649
1642 1650 jph = numpy.array([0,0,0,0])
1643 1651 h = (hmin,hmax)
1644 1652 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
1645 1653
1646 1654 # #Calculate AOA (Error N 3, 4)
1647 1655 # #JONES ET AL. 1998
1648 1656 # error = arrayParameters[:,-1]
1649 1657 # AOAthresh = numpy.pi/8
1650 1658 # phases = -arrayParameters[:,9:13]
1651 1659 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
1652 1660 #
1653 1661 # #Calculate Heights (Error N 13 and 14)
1654 1662 # error = arrayParameters[:,-1]
1655 1663 # Ranges = arrayParameters[:,2]
1656 1664 # zenith = arrayParameters[:,5]
1657 1665 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
1658 1666 # error = arrayParameters[:,-1]
1659 1667 #********************* END OF PARAMETERS CALCULATION **************************
1660 1668
1661 1669 #***************************+ PASS DATA TO NEXT STEP **********************
1662 1670 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
1663 1671 dataOut.data_param = arrayParameters
1664 1672
1665 1673 if arrayParameters is None:
1666 1674 dataOut.flagNoData = True
1667 1675 else:
1668 1676 dataOut.flagNoData = True
1669 1677
1670 1678 return
1671 1679
1672 1680 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
1673 1681
1674 1682 minIndex = min(newheis[0])
1675 1683 maxIndex = max(newheis[0])
1676 1684
1677 1685 voltage = voltage0[:,:,minIndex:maxIndex+1]
1678 1686 nLength = voltage.shape[1]/n
1679 1687 nMin = 0
1680 1688 nMax = 0
1681 1689 phaseOffset = numpy.zeros((len(pairslist),n))
1682 1690
1683 1691 for i in range(n):
1684 1692 nMax += nLength
1685 1693 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
1686 1694 phaseCCF = numpy.mean(phaseCCF, axis = 2)
1687 1695 phaseOffset[:,i] = phaseCCF.transpose()
1688 1696 nMin = nMax
1689 1697 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
1690 1698
1691 1699 #Remove Outliers
1692 1700 factor = 2
1693 1701 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
1694 1702 dw = numpy.std(wt,axis = 1)
1695 1703 dw = dw.reshape((dw.size,1))
1696 1704 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
1697 1705 phaseOffset[ind] = numpy.nan
1698 1706 phaseOffset = stats.nanmean(phaseOffset, axis=1)
1699 1707
1700 1708 return phaseOffset
1701 1709
1702 1710 def __shiftPhase(self, data, phaseShift):
1703 1711 #this will shift the phase of a complex number
1704 1712 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
1705 1713 return dataShifted
1706 1714
1707 1715 def __estimatePhaseDifference(self, array, pairslist):
1708 1716 nChannel = array.shape[0]
1709 1717 nHeights = array.shape[2]
1710 1718 numPairs = len(pairslist)
1711 1719 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
1712 1720 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
1713 1721
1714 1722 #Correct phases
1715 1723 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
1716 1724 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
1717 1725
1718 1726 if indDer[0].shape[0] > 0:
1719 1727 for i in range(indDer[0].shape[0]):
1720 1728 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
1721 1729 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
1722 1730
1723 1731 # for j in range(numSides):
1724 1732 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
1725 1733 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
1726 1734 #
1727 1735 #Linear
1728 1736 phaseInt = numpy.zeros((numPairs,1))
1729 1737 angAllCCF = phaseCCF[:,[0,1,3,4],0]
1730 1738 for j in range(numPairs):
1731 1739 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
1732 1740 phaseInt[j] = fit[1]
1733 1741 #Phase Differences
1734 1742 phaseDiff = phaseInt - phaseCCF[:,2,:]
1735 1743 phaseArrival = phaseInt.reshape(phaseInt.size)
1736 1744
1737 1745 #Dealias
1738 1746 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
1739 1747 # indAlias = numpy.where(phaseArrival > numpy.pi)
1740 1748 # phaseArrival[indAlias] -= 2*numpy.pi
1741 1749 # indAlias = numpy.where(phaseArrival < -numpy.pi)
1742 1750 # phaseArrival[indAlias] += 2*numpy.pi
1743 1751
1744 1752 return phaseDiff, phaseArrival
1745 1753
1746 1754 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
1747 1755 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
1748 1756 #find the phase shifts of each channel over 1 second intervals
1749 1757 #only look at ranges below the beacon signal
1750 1758 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1751 1759 numBlocks = int(volts.shape[1]/numProfPerBlock)
1752 1760 numHeights = volts.shape[2]
1753 1761 nChannel = volts.shape[0]
1754 1762 voltsCohDet = volts.copy()
1755 1763
1756 1764 pairsarray = numpy.array(pairslist)
1757 1765 indSides = pairsarray[:,1]
1758 1766 # indSides = numpy.array(range(nChannel))
1759 1767 # indSides = numpy.delete(indSides, indCenter)
1760 1768 #
1761 1769 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
1762 1770 listBlocks = numpy.array_split(volts, numBlocks, 1)
1763 1771
1764 1772 startInd = 0
1765 1773 endInd = 0
1766 1774
1767 1775 for i in range(numBlocks):
1768 1776 startInd = endInd
1769 1777 endInd = endInd + listBlocks[i].shape[1]
1770 1778
1771 1779 arrayBlock = listBlocks[i]
1772 1780 # arrayBlockCenter = listCenter[i]
1773 1781
1774 1782 #Estimate the Phase Difference
1775 1783 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
1776 1784 #Phase Difference RMS
1777 1785 arrayPhaseRMS = numpy.abs(phaseDiff)
1778 1786 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
1779 1787 indPhase = numpy.where(phaseRMSaux==4)
1780 1788 #Shifting
1781 1789 if indPhase[0].shape[0] > 0:
1782 1790 for j in range(indSides.size):
1783 1791 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
1784 1792 voltsCohDet[:,startInd:endInd,:] = arrayBlock
1785 1793
1786 1794 return voltsCohDet
1787 1795
1788 1796 def __calculateCCF(self, volts, pairslist ,laglist):
1789 1797
1790 1798 nHeights = volts.shape[2]
1791 1799 nPoints = volts.shape[1]
1792 1800 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
1793 1801
1794 1802 for i in range(len(pairslist)):
1795 1803 volts1 = volts[pairslist[i][0]]
1796 1804 volts2 = volts[pairslist[i][1]]
1797 1805
1798 1806 for t in range(len(laglist)):
1799 1807 idxT = laglist[t]
1800 1808 if idxT >= 0:
1801 1809 vStacked = numpy.vstack((volts2[idxT:,:],
1802 1810 numpy.zeros((idxT, nHeights),dtype='complex')))
1803 1811 else:
1804 1812 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
1805 1813 volts2[:(nPoints + idxT),:]))
1806 1814 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
1807 1815
1808 1816 vStacked = None
1809 1817 return voltsCCF
1810 1818
1811 1819 def __getNoise(self, power, timeSegment, timeInterval):
1812 1820 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1813 1821 numBlocks = int(power.shape[0]/numProfPerBlock)
1814 1822 numHeights = power.shape[1]
1815 1823
1816 1824 listPower = numpy.array_split(power, numBlocks, 0)
1817 1825 noise = numpy.zeros((power.shape[0], power.shape[1]))
1818 1826 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
1819 1827
1820 1828 startInd = 0
1821 1829 endInd = 0
1822 1830
1823 1831 for i in range(numBlocks): #split por canal
1824 1832 startInd = endInd
1825 1833 endInd = endInd + listPower[i].shape[0]
1826 1834
1827 1835 arrayBlock = listPower[i]
1828 1836 noiseAux = numpy.mean(arrayBlock, 0)
1829 1837 # noiseAux = numpy.median(noiseAux)
1830 1838 # noiseAux = numpy.mean(arrayBlock)
1831 1839 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
1832 1840
1833 1841 noiseAux1 = numpy.mean(arrayBlock)
1834 1842 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
1835 1843
1836 1844 return noise, noise1
1837 1845
1838 1846 def __findMeteors(self, power, thresh):
1839 1847 nProf = power.shape[0]
1840 1848 nHeights = power.shape[1]
1841 1849 listMeteors = []
1842 1850
1843 1851 for i in range(nHeights):
1844 1852 powerAux = power[:,i]
1845 1853 threshAux = thresh[:,i]
1846 1854
1847 1855 indUPthresh = numpy.where(powerAux > threshAux)[0]
1848 1856 indDNthresh = numpy.where(powerAux <= threshAux)[0]
1849 1857
1850 1858 j = 0
1851 1859
1852 1860 while (j < indUPthresh.size - 2):
1853 1861 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
1854 1862 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
1855 1863 indDNthresh = indDNthresh[indDNAux]
1856 1864
1857 1865 if (indDNthresh.size > 0):
1858 1866 indEnd = indDNthresh[0] - 1
1859 1867 indInit = indUPthresh[j] if isinstance(indUPthresh[j], (int, float)) else indUPthresh[j][0] ##CHECK!!!!
1860 1868
1861 1869 meteor = powerAux[indInit:indEnd + 1]
1862 1870 indPeak = meteor.argmax() + indInit
1863 1871 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
1864 1872
1865 1873 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
1866 1874 j = numpy.where(indUPthresh == indEnd)[0] + 1
1867 1875 else: j+=1
1868 1876 else: j+=1
1869 1877
1870 1878 return listMeteors
1871 1879
1872 1880 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
1873 1881
1874 1882 arrayMeteors = numpy.asarray(listMeteors)
1875 1883 listMeteors1 = []
1876 1884
1877 1885 while arrayMeteors.shape[0] > 0:
1878 1886 FLAs = arrayMeteors[:,4]
1879 1887 maxFLA = FLAs.argmax()
1880 1888 listMeteors1.append(arrayMeteors[maxFLA,:])
1881 1889
1882 1890 MeteorInitTime = arrayMeteors[maxFLA,1]
1883 1891 MeteorEndTime = arrayMeteors[maxFLA,3]
1884 1892 MeteorHeight = arrayMeteors[maxFLA,0]
1885 1893
1886 1894 #Check neighborhood
1887 1895 maxHeightIndex = MeteorHeight + rangeLimit
1888 1896 minHeightIndex = MeteorHeight - rangeLimit
1889 1897 minTimeIndex = MeteorInitTime - timeLimit
1890 1898 maxTimeIndex = MeteorEndTime + timeLimit
1891 1899
1892 1900 #Check Heights
1893 1901 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
1894 1902 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
1895 1903 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
1896 1904
1897 1905 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
1898 1906
1899 1907 return listMeteors1
1900 1908
1901 1909 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
1902 1910 numHeights = volts.shape[2]
1903 1911 nChannel = volts.shape[0]
1904 1912
1905 1913 thresholdPhase = thresh[0]
1906 1914 thresholdNoise = thresh[1]
1907 1915 thresholdDB = float(thresh[2])
1908 1916
1909 1917 thresholdDB1 = 10**(thresholdDB/10)
1910 1918 pairsarray = numpy.array(pairslist)
1911 1919 indSides = pairsarray[:,1]
1912 1920
1913 1921 pairslist1 = list(pairslist)
1914 1922 pairslist1.append((0,4))
1915 1923 pairslist1.append((1,3))
1916 1924
1917 1925 listMeteors1 = []
1918 1926 listPowerSeries = []
1919 1927 listVoltageSeries = []
1920 1928 #volts has the war data
1921 1929
1922 1930 if frequency == 30.175e6:
1923 1931 timeLag = 45*10**-3
1924 1932 else:
1925 1933 timeLag = 15*10**-3
1926 1934 lag = int(numpy.ceil(timeLag/timeInterval))
1927 1935
1928 1936 for i in range(len(listMeteors)):
1929 1937
1930 1938 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
1931 1939 meteorAux = numpy.zeros(16)
1932 1940
1933 1941 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
1934 1942 mHeight = int(listMeteors[i][0])
1935 1943 mStart = int(listMeteors[i][1])
1936 1944 mPeak = int(listMeteors[i][2])
1937 1945 mEnd = int(listMeteors[i][3])
1938 1946
1939 1947 #get the volt data between the start and end times of the meteor
1940 1948 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
1941 1949 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
1942 1950
1943 1951 #3.6. Phase Difference estimation
1944 1952 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
1945 1953
1946 1954 #3.7. Phase difference removal & meteor start, peak and end times reestimated
1947 1955 #meteorVolts0.- all Channels, all Profiles
1948 1956 meteorVolts0 = volts[:,:,mHeight]
1949 1957 meteorThresh = noise[:,mHeight]*thresholdNoise
1950 1958 meteorNoise = noise[:,mHeight]
1951 1959 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
1952 1960 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
1953 1961
1954 1962 #Times reestimation
1955 1963 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
1956 1964 if mStart1.size > 0:
1957 1965 mStart1 = mStart1[-1] + 1
1958 1966
1959 1967 else:
1960 1968 mStart1 = mPeak
1961 1969
1962 1970 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
1963 1971 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
1964 1972 if mEndDecayTime1.size == 0:
1965 1973 mEndDecayTime1 = powerNet0.size
1966 1974 else:
1967 1975 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
1968 1976 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
1969 1977
1970 1978 #meteorVolts1.- all Channels, from start to end
1971 1979 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
1972 1980 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
1973 1981 if meteorVolts2.shape[1] == 0:
1974 1982 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
1975 1983 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
1976 1984 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
1977 1985 ##################### END PARAMETERS REESTIMATION #########################
1978 1986
1979 1987 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
1980 1988 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
1981 1989 if meteorVolts2.shape[1] > 0:
1982 1990 #Phase Difference re-estimation
1983 1991 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
1984 1992 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
1985 1993 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
1986 1994 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
1987 1995 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
1988 1996
1989 1997 #Phase Difference RMS
1990 1998 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
1991 1999 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
1992 2000 #Data from Meteor
1993 2001 mPeak1 = powerNet1.argmax() + mStart1
1994 2002 mPeakPower1 = powerNet1.max()
1995 2003 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
1996 2004 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
1997 2005 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
1998 2006 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
1999 2007 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
2000 2008 #Vectorize
2001 2009 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
2002 2010 meteorAux[7:11] = phaseDiffint[0:4]
2003 2011
2004 2012 #Rejection Criterions
2005 2013 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
2006 2014 meteorAux[-1] = 17
2007 2015 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
2008 2016 meteorAux[-1] = 1
2009 2017
2010 2018
2011 2019 else:
2012 2020 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
2013 2021 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
2014 2022 PowerSeries = 0
2015 2023
2016 2024 listMeteors1.append(meteorAux)
2017 2025 listPowerSeries.append(PowerSeries)
2018 2026 listVoltageSeries.append(meteorVolts1)
2019 2027
2020 2028 return listMeteors1, listPowerSeries, listVoltageSeries
2021 2029
2022 2030 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
2023 2031
2024 2032 threshError = 10
2025 2033 #Depending if it is 30 or 50 MHz
2026 2034 if frequency == 30.175e6:
2027 2035 timeLag = 45*10**-3
2028 2036 else:
2029 2037 timeLag = 15*10**-3
2030 2038 lag = int(numpy.ceil(timeLag/timeInterval))
2031 2039
2032 2040 listMeteors1 = []
2033 2041
2034 2042 for i in range(len(listMeteors)):
2035 2043 meteorPower = listPower[i]
2036 2044 meteorAux = listMeteors[i]
2037 2045
2038 2046 if meteorAux[-1] == 0:
2039 2047
2040 2048 try:
2041 2049 indmax = meteorPower.argmax()
2042 2050 indlag = indmax + lag
2043 2051
2044 2052 y = meteorPower[indlag:]
2045 2053 x = numpy.arange(0, y.size)*timeLag
2046 2054
2047 2055 #first guess
2048 2056 a = y[0]
2049 2057 tau = timeLag
2050 2058 #exponential fit
2051 2059 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
2052 2060 y1 = self.__exponential_function(x, *popt)
2053 2061 #error estimation
2054 2062 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
2055 2063
2056 2064 decayTime = popt[1]
2057 2065 riseTime = indmax*timeInterval
2058 2066 meteorAux[11:13] = [decayTime, error]
2059 2067
2060 2068 #Table items 7, 8 and 11
2061 2069 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
2062 2070 meteorAux[-1] = 7
2063 2071 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
2064 2072 meteorAux[-1] = 8
2065 2073 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
2066 2074 meteorAux[-1] = 11
2067 2075
2068 2076
2069 2077 except:
2070 2078 meteorAux[-1] = 11
2071 2079
2072 2080
2073 2081 listMeteors1.append(meteorAux)
2074 2082
2075 2083 return listMeteors1
2076 2084
2077 2085 #Exponential Function
2078 2086
2079 2087 def __exponential_function(self, x, a, tau):
2080 2088 y = a*numpy.exp(-x/tau)
2081 2089 return y
2082 2090
2083 2091 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
2084 2092
2085 2093 pairslist1 = list(pairslist)
2086 2094 pairslist1.append((0,4))
2087 2095 pairslist1.append((1,3))
2088 2096 numPairs = len(pairslist1)
2089 2097 #Time Lag
2090 2098 timeLag = 45*10**-3
2091 2099 c = 3e8
2092 2100 lag = numpy.ceil(timeLag/timeInterval)
2093 2101 freq = 30.175e6
2094 2102
2095 2103 listMeteors1 = []
2096 2104
2097 2105 for i in range(len(listMeteors)):
2098 2106 meteorAux = listMeteors[i]
2099 2107 if meteorAux[-1] == 0:
2100 2108 mStart = listMeteors[i][1]
2101 2109 mPeak = listMeteors[i][2]
2102 2110 mLag = mPeak - mStart + lag
2103 2111
2104 2112 #get the volt data between the start and end times of the meteor
2105 2113 meteorVolts = listVolts[i]
2106 2114 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
2107 2115
2108 2116 #Get CCF
2109 2117 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
2110 2118
2111 2119 #Method 2
2112 2120 slopes = numpy.zeros(numPairs)
2113 2121 time = numpy.array([-2,-1,1,2])*timeInterval
2114 2122 angAllCCF = numpy.angle(allCCFs[:,[0,4,2,3],0])
2115 2123
2116 2124 #Correct phases
2117 2125 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
2118 2126 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2119 2127
2120 2128 if indDer[0].shape[0] > 0:
2121 2129 for i in range(indDer[0].shape[0]):
2122 2130 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
2123 2131 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
2124 2132
2125 2133 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
2126 2134 for j in range(numPairs):
2127 2135 fit = stats.linregress(time, angAllCCF[j,:])
2128 2136 slopes[j] = fit[0]
2129 2137
2130 2138 #Remove Outlier
2131 2139 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2132 2140 # slopes = numpy.delete(slopes,indOut)
2133 2141 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2134 2142 # slopes = numpy.delete(slopes,indOut)
2135 2143
2136 2144 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
2137 2145 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
2138 2146 meteorAux[-2] = radialError
2139 2147 meteorAux[-3] = radialVelocity
2140 2148
2141 2149 #Setting Error
2142 2150 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
2143 2151 if numpy.abs(radialVelocity) > 200:
2144 2152 meteorAux[-1] = 15
2145 2153 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
2146 2154 elif radialError > radialStdThresh:
2147 2155 meteorAux[-1] = 12
2148 2156
2149 2157 listMeteors1.append(meteorAux)
2150 2158 return listMeteors1
2151 2159
2152 2160 def __setNewArrays(self, listMeteors, date, heiRang):
2153 2161
2154 2162 #New arrays
2155 2163 arrayMeteors = numpy.array(listMeteors)
2156 2164 arrayParameters = numpy.zeros((len(listMeteors), 13))
2157 2165
2158 2166 #Date inclusion
2159 2167 # date = re.findall(r'\((.*?)\)', date)
2160 2168 # date = date[0].split(',')
2161 2169 # date = map(int, date)
2162 2170 #
2163 2171 # if len(date)<6:
2164 2172 # date.append(0)
2165 2173 #
2166 2174 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
2167 2175 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
2168 2176 arrayDate = numpy.tile(date, (len(listMeteors)))
2169 2177
2170 2178 #Meteor array
2171 2179 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
2172 2180 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
2173 2181
2174 2182 #Parameters Array
2175 2183 arrayParameters[:,0] = arrayDate #Date
2176 2184 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
2177 2185 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
2178 2186 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
2179 2187 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
2180 2188
2181 2189
2182 2190 return arrayParameters
2183 2191
2184 2192 class CorrectSMPhases(Operation):
2185 2193
2186 2194 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
2187 2195
2188 2196 arrayParameters = dataOut.data_param
2189 2197 pairsList = []
2190 2198 pairx = (0,1)
2191 2199 pairy = (2,3)
2192 2200 pairsList.append(pairx)
2193 2201 pairsList.append(pairy)
2194 2202 jph = numpy.zeros(4)
2195 2203
2196 2204 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2197 2205 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2198 2206 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
2199 2207
2200 2208 meteorOps = SMOperations()
2201 2209 if channelPositions is None:
2202 2210 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2203 2211 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2204 2212
2205 2213 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2206 2214 h = (hmin,hmax)
2207 2215
2208 2216 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2209 2217
2210 2218 dataOut.data_param = arrayParameters
2211 2219 return
2212 2220
2213 2221 class SMPhaseCalibration(Operation):
2214 2222
2215 2223 __buffer = None
2216 2224
2217 2225 __initime = None
2218 2226
2219 2227 __dataReady = False
2220 2228
2221 2229 __isConfig = False
2222 2230
2223 2231 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
2224 2232
2225 2233 dataTime = currentTime + paramInterval
2226 2234 deltaTime = dataTime - initTime
2227 2235
2228 2236 if deltaTime >= outputInterval or deltaTime < 0:
2229 2237 return True
2230 2238
2231 2239 return False
2232 2240
2233 2241 def __getGammas(self, pairs, d, phases):
2234 2242 gammas = numpy.zeros(2)
2235 2243
2236 2244 for i in range(len(pairs)):
2237 2245
2238 2246 pairi = pairs[i]
2239 2247
2240 2248 phip3 = phases[:,pairi[1]]
2241 2249 d3 = d[pairi[1]]
2242 2250 phip2 = phases[:,pairi[0]]
2243 2251 d2 = d[pairi[0]]
2244 2252 #Calculating gamma
2245 2253 # jdcos = alp1/(k*d1)
2246 2254 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
2247 2255 jgamma = -phip2*d3/d2 - phip3
2248 2256 jgamma = numpy.angle(numpy.exp(1j*jgamma))
2249 2257 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
2250 2258 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
2251 2259
2252 2260 #Revised distribution
2253 2261 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
2254 2262
2255 2263 #Histogram
2256 2264 nBins = 64.0
2257 2265 rmin = -0.5*numpy.pi
2258 2266 rmax = 0.5*numpy.pi
2259 2267 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
2260 2268
2261 2269 meteorsY = phaseHisto[0]
2262 2270 phasesX = phaseHisto[1][:-1]
2263 2271 width = phasesX[1] - phasesX[0]
2264 2272 phasesX += width/2
2265 2273
2266 2274 #Gaussian aproximation
2267 2275 bpeak = meteorsY.argmax()
2268 2276 peak = meteorsY.max()
2269 2277 jmin = bpeak - 5
2270 2278 jmax = bpeak + 5 + 1
2271 2279
2272 2280 if jmin<0:
2273 2281 jmin = 0
2274 2282 jmax = 6
2275 2283 elif jmax > meteorsY.size:
2276 2284 jmin = meteorsY.size - 6
2277 2285 jmax = meteorsY.size
2278 2286
2279 2287 x0 = numpy.array([peak,bpeak,50])
2280 2288 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
2281 2289
2282 2290 #Gammas
2283 2291 gammas[i] = coeff[0][1]
2284 2292
2285 2293 return gammas
2286 2294
2287 2295 def __residualFunction(self, coeffs, y, t):
2288 2296
2289 2297 return y - self.__gauss_function(t, coeffs)
2290 2298
2291 2299 def __gauss_function(self, t, coeffs):
2292 2300
2293 2301 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
2294 2302
2295 2303 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2296 2304 meteorOps = SMOperations()
2297 2305 nchan = 4
2298 2306 pairx = pairsList[0]
2299 2307 pairy = pairsList[1]
2300 2308 center_xangle = 0
2301 2309 center_yangle = 0
2302 2310 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
2303 2311 ntimes = len(range_angle)
2304 2312
2305 2313 nstepsx = 20.0
2306 2314 nstepsy = 20.0
2307 2315
2308 2316 for iz in range(ntimes):
2309 2317 min_xangle = -range_angle[iz]/2 + center_xangle
2310 2318 max_xangle = range_angle[iz]/2 + center_xangle
2311 2319 min_yangle = -range_angle[iz]/2 + center_yangle
2312 2320 max_yangle = range_angle[iz]/2 + center_yangle
2313 2321
2314 2322 inc_x = (max_xangle-min_xangle)/nstepsx
2315 2323 inc_y = (max_yangle-min_yangle)/nstepsy
2316 2324
2317 2325 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
2318 2326 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
2319 2327 penalty = numpy.zeros((nstepsx,nstepsy))
2320 2328 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
2321 2329 jph = numpy.zeros(nchan)
2322 2330
2323 2331 # Iterations looking for the offset
2324 2332 for iy in range(int(nstepsy)):
2325 2333 for ix in range(int(nstepsx)):
2326 2334 jph[pairy[1]] = alpha_y[iy]
2327 2335 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
2328 2336
2329 2337 jph[pairx[1]] = alpha_x[ix]
2330 2338 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
2331 2339
2332 2340 jph_array[:,ix,iy] = jph
2333 2341
2334 2342 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
2335 2343 error = meteorsArray1[:,-1]
2336 2344 ind1 = numpy.where(error==0)[0]
2337 2345 penalty[ix,iy] = ind1.size
2338 2346
2339 2347 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
2340 2348 phOffset = jph_array[:,i,j]
2341 2349
2342 2350 center_xangle = phOffset[pairx[1]]
2343 2351 center_yangle = phOffset[pairy[1]]
2344 2352
2345 2353 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
2346 2354 phOffset = phOffset*180/numpy.pi
2347 2355 return phOffset
2348 2356
2349 2357
2350 2358 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
2351 2359
2352 2360 dataOut.flagNoData = True
2353 2361 self.__dataReady = False
2354 2362 dataOut.outputInterval = nHours*3600
2355 2363
2356 2364 if self.__isConfig == False:
2357 2365 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2358 2366 #Get Initial LTC time
2359 2367 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2360 2368 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2361 2369
2362 2370 self.__isConfig = True
2363 2371
2364 2372 if self.__buffer is None:
2365 2373 self.__buffer = dataOut.data_param.copy()
2366 2374
2367 2375 else:
2368 2376 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2369 2377
2370 2378 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2371 2379
2372 2380 if self.__dataReady:
2373 2381 dataOut.utctimeInit = self.__initime
2374 2382 self.__initime += dataOut.outputInterval #to erase time offset
2375 2383
2376 2384 freq = dataOut.frequency
2377 2385 c = dataOut.C #m/s
2378 2386 lamb = c/freq
2379 2387 k = 2*numpy.pi/lamb
2380 2388 azimuth = 0
2381 2389 h = (hmin, hmax)
2382 2390 pairs = ((0,1),(2,3))
2383 2391
2384 2392 if channelPositions is None:
2385 2393 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2386 2394 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2387 2395 meteorOps = SMOperations()
2388 2396 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2389 2397
2390 2398 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
2391 2399
2392 2400 meteorsArray = self.__buffer
2393 2401 error = meteorsArray[:,-1]
2394 2402 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
2395 2403 ind1 = numpy.where(boolError)[0]
2396 2404 meteorsArray = meteorsArray[ind1,:]
2397 2405 meteorsArray[:,-1] = 0
2398 2406 phases = meteorsArray[:,8:12]
2399 2407
2400 2408 #Calculate Gammas
2401 2409 gammas = self.__getGammas(pairs, distances, phases)
2402 2410 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2403 2411 #Calculate Phases
2404 2412 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
2405 2413 phasesOff = phasesOff.reshape((1,phasesOff.size))
2406 2414 dataOut.data_output = -phasesOff
2407 2415 dataOut.flagNoData = False
2408 2416 dataOut.channelList = pairslist0
2409 2417 self.__buffer = None
2410 2418
2411 2419
2412 2420 return
2413 2421
2414 2422 class SMOperations():
2415 2423
2416 2424 def __init__(self):
2417 2425
2418 2426 return
2419 2427
2420 2428 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
2421 2429
2422 2430 arrayParameters = arrayParameters0.copy()
2423 2431 hmin = h[0]
2424 2432 hmax = h[1]
2425 2433
2426 2434 #Calculate AOA (Error N 3, 4)
2427 2435 #JONES ET AL. 1998
2428 2436 AOAthresh = numpy.pi/8
2429 2437 error = arrayParameters[:,-1]
2430 2438 phases = -arrayParameters[:,8:12] + jph
2431 2439 # phases = numpy.unwrap(phases)
2432 2440 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
2433 2441
2434 2442 #Calculate Heights (Error N 13 and 14)
2435 2443 error = arrayParameters[:,-1]
2436 2444 Ranges = arrayParameters[:,1]
2437 2445 zenith = arrayParameters[:,4]
2438 2446 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2439 2447
2440 2448 #----------------------- Get Final data ------------------------------------
2441 2449 # error = arrayParameters[:,-1]
2442 2450 # ind1 = numpy.where(error==0)[0]
2443 2451 # arrayParameters = arrayParameters[ind1,:]
2444 2452
2445 2453 return arrayParameters
2446 2454
2447 2455 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
2448 2456
2449 2457 arrayAOA = numpy.zeros((phases.shape[0],3))
2450 2458 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
2451 2459
2452 2460 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2453 2461 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2454 2462 arrayAOA[:,2] = cosDirError
2455 2463
2456 2464 azimuthAngle = arrayAOA[:,0]
2457 2465 zenithAngle = arrayAOA[:,1]
2458 2466
2459 2467 #Setting Error
2460 2468 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
2461 2469 error[indError] = 0
2462 2470 #Number 3: AOA not fesible
2463 2471 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2464 2472 error[indInvalid] = 3
2465 2473 #Number 4: Large difference in AOAs obtained from different antenna baselines
2466 2474 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2467 2475 error[indInvalid] = 4
2468 2476 return arrayAOA, error
2469 2477
2470 2478 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
2471 2479
2472 2480 #Initializing some variables
2473 2481 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2474 2482 ang_aux = ang_aux.reshape(1,ang_aux.size)
2475 2483
2476 2484 cosdir = numpy.zeros((arrayPhase.shape[0],2))
2477 2485 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2478 2486
2479 2487
2480 2488 for i in range(2):
2481 2489 ph0 = arrayPhase[:,pairsList[i][0]]
2482 2490 ph1 = arrayPhase[:,pairsList[i][1]]
2483 2491 d0 = distances[pairsList[i][0]]
2484 2492 d1 = distances[pairsList[i][1]]
2485 2493
2486 2494 ph0_aux = ph0 + ph1
2487 2495 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
2488 2496 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
2489 2497 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
2490 2498 #First Estimation
2491 2499 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
2492 2500
2493 2501 #Most-Accurate Second Estimation
2494 2502 phi1_aux = ph0 - ph1
2495 2503 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2496 2504 #Direction Cosine 1
2497 2505 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
2498 2506
2499 2507 #Searching the correct Direction Cosine
2500 2508 cosdir0_aux = cosdir0[:,i]
2501 2509 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2502 2510 #Minimum Distance
2503 2511 cosDiff = (cosdir1 - cosdir0_aux)**2
2504 2512 indcos = cosDiff.argmin(axis = 1)
2505 2513 #Saving Value obtained
2506 2514 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2507 2515
2508 2516 return cosdir0, cosdir
2509 2517
2510 2518 def __calculateAOA(self, cosdir, azimuth):
2511 2519 cosdirX = cosdir[:,0]
2512 2520 cosdirY = cosdir[:,1]
2513 2521
2514 2522 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2515 2523 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
2516 2524 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2517 2525
2518 2526 return angles
2519 2527
2520 2528 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2521 2529
2522 2530 Ramb = 375 #Ramb = c/(2*PRF)
2523 2531 Re = 6371 #Earth Radius
2524 2532 heights = numpy.zeros(Ranges.shape)
2525 2533
2526 2534 R_aux = numpy.array([0,1,2])*Ramb
2527 2535 R_aux = R_aux.reshape(1,R_aux.size)
2528 2536
2529 2537 Ranges = Ranges.reshape(Ranges.size,1)
2530 2538
2531 2539 Ri = Ranges + R_aux
2532 2540 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2533 2541
2534 2542 #Check if there is a height between 70 and 110 km
2535 2543 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2536 2544 ind_h = numpy.where(h_bool == 1)[0]
2537 2545
2538 2546 hCorr = hi[ind_h, :]
2539 2547 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2540 2548
2541 2549 hCorr = hi[ind_hCorr]
2542 2550 heights[ind_h] = hCorr
2543 2551
2544 2552 #Setting Error
2545 2553 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2546 2554 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2547 2555 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
2548 2556 error[indError] = 0
2549 2557 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2550 2558 error[indInvalid2] = 14
2551 2559 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2552 2560 error[indInvalid1] = 13
2553 2561
2554 2562 return heights, error
2555 2563
2556 2564 def getPhasePairs(self, channelPositions):
2557 2565 chanPos = numpy.array(channelPositions)
2558 2566 listOper = list(itertools.combinations(range(5),2))
2559 2567
2560 2568 distances = numpy.zeros(4)
2561 2569 axisX = []
2562 2570 axisY = []
2563 2571 distX = numpy.zeros(3)
2564 2572 distY = numpy.zeros(3)
2565 2573 ix = 0
2566 2574 iy = 0
2567 2575
2568 2576 pairX = numpy.zeros((2,2))
2569 2577 pairY = numpy.zeros((2,2))
2570 2578
2571 2579 for i in range(len(listOper)):
2572 2580 pairi = listOper[i]
2573 2581
2574 2582 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
2575 2583
2576 2584 if posDif[0] == 0:
2577 2585 axisY.append(pairi)
2578 2586 distY[iy] = posDif[1]
2579 2587 iy += 1
2580 2588 elif posDif[1] == 0:
2581 2589 axisX.append(pairi)
2582 2590 distX[ix] = posDif[0]
2583 2591 ix += 1
2584 2592
2585 2593 for i in range(2):
2586 2594 if i==0:
2587 2595 dist0 = distX
2588 2596 axis0 = axisX
2589 2597 else:
2590 2598 dist0 = distY
2591 2599 axis0 = axisY
2592 2600
2593 2601 side = numpy.argsort(dist0)[:-1]
2594 2602 axis0 = numpy.array(axis0)[side,:]
2595 2603 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
2596 2604 axis1 = numpy.unique(numpy.reshape(axis0,4))
2597 2605 side = axis1[axis1 != chanC]
2598 2606 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
2599 2607 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
2600 2608 if diff1<0:
2601 2609 chan2 = side[0]
2602 2610 d2 = numpy.abs(diff1)
2603 2611 chan1 = side[1]
2604 2612 d1 = numpy.abs(diff2)
2605 2613 else:
2606 2614 chan2 = side[1]
2607 2615 d2 = numpy.abs(diff2)
2608 2616 chan1 = side[0]
2609 2617 d1 = numpy.abs(diff1)
2610 2618
2611 2619 if i==0:
2612 2620 chanCX = chanC
2613 2621 chan1X = chan1
2614 2622 chan2X = chan2
2615 2623 distances[0:2] = numpy.array([d1,d2])
2616 2624 else:
2617 2625 chanCY = chanC
2618 2626 chan1Y = chan1
2619 2627 chan2Y = chan2
2620 2628 distances[2:4] = numpy.array([d1,d2])
2621 2629 # axisXsides = numpy.reshape(axisX[ix,:],4)
2622 2630 #
2623 2631 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
2624 2632 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
2625 2633 #
2626 2634 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
2627 2635 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
2628 2636 # channel25X = int(pairX[0,ind25X])
2629 2637 # channel20X = int(pairX[1,ind20X])
2630 2638 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
2631 2639 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
2632 2640 # channel25Y = int(pairY[0,ind25Y])
2633 2641 # channel20Y = int(pairY[1,ind20Y])
2634 2642
2635 2643 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
2636 2644 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
2637 2645
2638 2646 return pairslist, distances
2639 2647 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2640 2648 #
2641 2649 # arrayAOA = numpy.zeros((phases.shape[0],3))
2642 2650 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2643 2651 #
2644 2652 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2645 2653 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2646 2654 # arrayAOA[:,2] = cosDirError
2647 2655 #
2648 2656 # azimuthAngle = arrayAOA[:,0]
2649 2657 # zenithAngle = arrayAOA[:,1]
2650 2658 #
2651 2659 # #Setting Error
2652 2660 # #Number 3: AOA not fesible
2653 2661 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2654 2662 # error[indInvalid] = 3
2655 2663 # #Number 4: Large difference in AOAs obtained from different antenna baselines
2656 2664 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2657 2665 # error[indInvalid] = 4
2658 2666 # return arrayAOA, error
2659 2667 #
2660 2668 # def __getDirectionCosines(self, arrayPhase, pairsList):
2661 2669 #
2662 2670 # #Initializing some variables
2663 2671 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2664 2672 # ang_aux = ang_aux.reshape(1,ang_aux.size)
2665 2673 #
2666 2674 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
2667 2675 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2668 2676 #
2669 2677 #
2670 2678 # for i in range(2):
2671 2679 # #First Estimation
2672 2680 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2673 2681 # #Dealias
2674 2682 # indcsi = numpy.where(phi0_aux > numpy.pi)
2675 2683 # phi0_aux[indcsi] -= 2*numpy.pi
2676 2684 # indcsi = numpy.where(phi0_aux < -numpy.pi)
2677 2685 # phi0_aux[indcsi] += 2*numpy.pi
2678 2686 # #Direction Cosine 0
2679 2687 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2680 2688 #
2681 2689 # #Most-Accurate Second Estimation
2682 2690 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2683 2691 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2684 2692 # #Direction Cosine 1
2685 2693 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2686 2694 #
2687 2695 # #Searching the correct Direction Cosine
2688 2696 # cosdir0_aux = cosdir0[:,i]
2689 2697 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2690 2698 # #Minimum Distance
2691 2699 # cosDiff = (cosdir1 - cosdir0_aux)**2
2692 2700 # indcos = cosDiff.argmin(axis = 1)
2693 2701 # #Saving Value obtained
2694 2702 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2695 2703 #
2696 2704 # return cosdir0, cosdir
2697 2705 #
2698 2706 # def __calculateAOA(self, cosdir, azimuth):
2699 2707 # cosdirX = cosdir[:,0]
2700 2708 # cosdirY = cosdir[:,1]
2701 2709 #
2702 2710 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2703 2711 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2704 2712 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2705 2713 #
2706 2714 # return angles
2707 2715 #
2708 2716 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2709 2717 #
2710 2718 # Ramb = 375 #Ramb = c/(2*PRF)
2711 2719 # Re = 6371 #Earth Radius
2712 2720 # heights = numpy.zeros(Ranges.shape)
2713 2721 #
2714 2722 # R_aux = numpy.array([0,1,2])*Ramb
2715 2723 # R_aux = R_aux.reshape(1,R_aux.size)
2716 2724 #
2717 2725 # Ranges = Ranges.reshape(Ranges.size,1)
2718 2726 #
2719 2727 # Ri = Ranges + R_aux
2720 2728 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2721 2729 #
2722 2730 # #Check if there is a height between 70 and 110 km
2723 2731 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2724 2732 # ind_h = numpy.where(h_bool == 1)[0]
2725 2733 #
2726 2734 # hCorr = hi[ind_h, :]
2727 2735 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2728 2736 #
2729 2737 # hCorr = hi[ind_hCorr]
2730 2738 # heights[ind_h] = hCorr
2731 2739 #
2732 2740 # #Setting Error
2733 2741 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2734 2742 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2735 2743 #
2736 2744 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2737 2745 # error[indInvalid2] = 14
2738 2746 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2739 2747 # error[indInvalid1] = 13
2740 2748 #
2741 2749 # return heights, error
Show file before | Show file after | Hide comments
@@ -1,425 +1,437
1 1 '''
2 2 @author: Juan C. Espinoza
3 3 '''
4 4
5 5 import time
6 6 import json
7 7 import numpy
8 8 import paho.mqtt.client as mqtt
9 9 import zmq
10 10 import cPickle as pickle
11 11 import datetime
12 12 from zmq.utils.monitor import recv_monitor_message
13 13 from functools import wraps
14 14 from threading import Thread
15 15 from multiprocessing import Process
16 16
17 17 from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit
18 18
19 19 MAXNUMX = 100
20 20 MAXNUMY = 100
21 21
22 22 class PrettyFloat(float):
23 23 def __repr__(self):
24 24 return '%.2f' % self
25 25
26 26 def roundFloats(obj):
27 27 if isinstance(obj, list):
28 28 return map(roundFloats, obj)
29 29 elif isinstance(obj, float):
30 30 return round(obj, 2)
31 31
32 32 def decimate(z):
33 33 # dx = int(len(self.x)/self.__MAXNUMX) + 1
34 34
35 35 dy = int(len(z[0])/MAXNUMY) + 1
36 36
37 37 return z[::, ::dy]
38 38
39 39 class throttle(object):
40 40 """Decorator that prevents a function from being called more than once every
41 41 time period.
42 42 To create a function that cannot be called more than once a minute, but
43 43 will sleep until it can be called:
44 44 @throttle(minutes=1)
45 45 def foo():
46 46 pass
47 47
48 48 for i in range(10):
49 49 foo()
50 50 print "This function has run %s times." % i
51 51 """
52 52
53 53 def __init__(self, seconds=0, minutes=0, hours=0):
54 54 self.throttle_period = datetime.timedelta(
55 55 seconds=seconds, minutes=minutes, hours=hours
56 56 )
57 57
58 58 self.time_of_last_call = datetime.datetime.min
59 59
60 60 def __call__(self, fn):
61 61 @wraps(fn)
62 62 def wrapper(*args, **kwargs):
63 63 now = datetime.datetime.now()
64 64 time_since_last_call = now - self.time_of_last_call
65 65 time_left = self.throttle_period - time_since_last_call
66 66
67 67 if time_left > datetime.timedelta(seconds=0):
68 68 return
69 69
70 70 self.time_of_last_call = datetime.datetime.now()
71 71 return fn(*args, **kwargs)
72 72
73 73 return wrapper
74 74
75 75
76 76 class PublishData(Operation):
77 77 """Clase publish."""
78 78
79 79 def __init__(self, **kwargs):
80 80 """Inicio."""
81 81 Operation.__init__(self, **kwargs)
82 82 self.isConfig = False
83 83 self.client = None
84 84 self.zeromq = None
85 85 self.mqtt = None
86 86
87 87 def on_disconnect(self, client, userdata, rc):
88 88 if rc != 0:
89 89 print("Unexpected disconnection.")
90 90 self.connect()
91 91
92 92 def connect(self):
93 93 print 'trying to connect'
94 94 try:
95 95 self.client.connect(
96 96 host=self.host,
97 97 port=self.port,
98 98 keepalive=60*10,
99 99 bind_address='')
100 100 self.client.loop_start()
101 101 # self.client.publish(
102 102 # self.topic + 'SETUP',
103 103 # json.dumps(setup),
104 104 # retain=True
105 105 # )
106 106 except:
107 107 print "MQTT Conection error."
108 108 self.client = False
109 109
110 110 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, **kwargs):
111 111 self.counter = 0
112 112 self.topic = kwargs.get('topic', 'schain')
113 113 self.delay = kwargs.get('delay', 0)
114 114 self.plottype = kwargs.get('plottype', 'spectra')
115 115 self.host = kwargs.get('host', "10.10.10.82")
116 116 self.port = kwargs.get('port', 3000)
117 117 self.clientId = clientId
118 118 self.cnt = 0
119 119 self.zeromq = zeromq
120 120 self.mqtt = kwargs.get('plottype', 0)
121 121 self.client = None
122 122 setup = []
123 123 if mqtt is 1:
124 124 self.client = mqtt.Client(
125 125 client_id=self.clientId + self.topic + 'SCHAIN',
126 126 clean_session=True)
127 127 self.client.on_disconnect = self.on_disconnect
128 128 self.connect()
129 129 for plot in self.plottype:
130 130 setup.append({
131 131 'plot': plot,
132 132 'topic': self.topic + plot,
133 133 'title': getattr(self, plot + '_' + 'title', False),
134 134 'xlabel': getattr(self, plot + '_' + 'xlabel', False),
135 135 'ylabel': getattr(self, plot + '_' + 'ylabel', False),
136 136 'xrange': getattr(self, plot + '_' + 'xrange', False),
137 137 'yrange': getattr(self, plot + '_' + 'yrange', False),
138 138 'zrange': getattr(self, plot + '_' + 'zrange', False),
139 139 })
140 140 if zeromq is 1:
141 141 context = zmq.Context()
142 142 self.zmq_socket = context.socket(zmq.PUSH)
143 143 server = kwargs.get('server', 'zmq.pipe')
144 144
145 145 if 'tcp://' in server:
146 146 address = server
147 147 else:
148 148 address = 'ipc:///tmp/%s' % server
149 149
150 150 self.zmq_socket.connect(address)
151 151 time.sleep(1)
152 152
153 153 def publish_data(self):
154 154 self.dataOut.finished = False
155 155 if self.mqtt is 1:
156 156 yData = self.dataOut.heightList[:2].tolist()
157 157 if self.plottype == 'spectra':
158 158 data = getattr(self.dataOut, 'data_spc')
159 159 z = data/self.dataOut.normFactor
160 160 zdB = 10*numpy.log10(z)
161 161 xlen, ylen = zdB[0].shape
162 162 dx = int(xlen/MAXNUMX) + 1
163 163 dy = int(ylen/MAXNUMY) + 1
164 164 Z = [0 for i in self.dataOut.channelList]
165 165 for i in self.dataOut.channelList:
166 166 Z[i] = zdB[i][::dx, ::dy].tolist()
167 167 payload = {
168 168 'timestamp': self.dataOut.utctime,
169 169 'data': roundFloats(Z),
170 170 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
171 171 'interval': self.dataOut.getTimeInterval(),
172 172 'type': self.plottype,
173 173 'yData': yData
174 174 }
175 175 # print payload
176 176
177 177 elif self.plottype in ('rti', 'power'):
178 178 data = getattr(self.dataOut, 'data_spc')
179 179 z = data/self.dataOut.normFactor
180 180 avg = numpy.average(z, axis=1)
181 181 avgdB = 10*numpy.log10(avg)
182 182 xlen, ylen = z[0].shape
183 183 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
184 184 AVG = [0 for i in self.dataOut.channelList]
185 185 for i in self.dataOut.channelList:
186 186 AVG[i] = avgdB[i][::dy].tolist()
187 187 payload = {
188 188 'timestamp': self.dataOut.utctime,
189 189 'data': roundFloats(AVG),
190 190 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
191 191 'interval': self.dataOut.getTimeInterval(),
192 192 'type': self.plottype,
193 193 'yData': yData
194 194 }
195 195 elif self.plottype == 'noise':
196 196 noise = self.dataOut.getNoise()/self.dataOut.normFactor
197 197 noisedB = 10*numpy.log10(noise)
198 198 payload = {
199 199 'timestamp': self.dataOut.utctime,
200 200 'data': roundFloats(noisedB.reshape(-1, 1).tolist()),
201 201 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
202 202 'interval': self.dataOut.getTimeInterval(),
203 203 'type': self.plottype,
204 204 'yData': yData
205 205 }
206 206 elif self.plottype == 'snr':
207 207 data = getattr(self.dataOut, 'data_SNR')
208 208 avgdB = 10*numpy.log10(data)
209 209
210 210 ylen = data[0].size
211 211 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
212 212 AVG = [0 for i in self.dataOut.channelList]
213 213 for i in self.dataOut.channelList:
214 214 AVG[i] = avgdB[i][::dy].tolist()
215 215 payload = {
216 216 'timestamp': self.dataOut.utctime,
217 217 'data': roundFloats(AVG),
218 218 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
219 219 'type': self.plottype,
220 220 'yData': yData
221 221 }
222 222 else:
223 223 print "Tipo de grafico invalido"
224 224 payload = {
225 225 'data': 'None',
226 226 'timestamp': 'None',
227 227 'type': None
228 228 }
229 229 # print 'Publishing data to {}'.format(self.host)
230 230 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
231 231
232 232 if self.zeromq is 1:
233 233 print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime)
234 234 self.zmq_socket.send_pyobj(self.dataOut)
235 235
236 236 def run(self, dataOut, **kwargs):
237 237 self.dataOut = dataOut
238 238 if not self.isConfig:
239 239 self.setup(**kwargs)
240 240 self.isConfig = True
241 241
242 242 self.publish_data()
243 243 time.sleep(self.delay)
244 244
245 245 def close(self):
246 246 if self.zeromq is 1:
247 247 self.dataOut.finished = True
248 248 self.zmq_socket.send_pyobj(self.dataOut)
249 249
250 250 if self.client:
251 251 self.client.loop_stop()
252 252 self.client.disconnect()
253 253
254 254
255 255 class ReceiverData(ProcessingUnit, Process):
256 256
257 257 throttle_value = 5
258 258
259 259 def __init__(self, **kwargs):
260 260
261 261 ProcessingUnit.__init__(self, **kwargs)
262 262 Process.__init__(self)
263 263 self.mp = False
264 264 self.isConfig = False
265 265 self.isWebConfig = False
266 266 self.plottypes =[]
267 267 self.connections = 0
268 268 server = kwargs.get('server', 'zmq.pipe')
269 269 plot_server = kwargs.get('plot_server', 'zmq.web')
270 270 if 'tcp://' in server:
271 271 address = server
272 272 else:
273 273 address = 'ipc:///tmp/%s' % server
274 274
275 275 if 'tcp://' in plot_server:
276 276 plot_address = plot_server
277 277 else:
278 278 plot_address = 'ipc:///tmp/%s' % plot_server
279 279
280 280 self.address = address
281 281 self.plot_address = plot_address
282 282 self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')]
283 283 self.realtime = kwargs.get('realtime', False)
284 self.throttle_value = kwargs.get('throttle', 10)
284 self.throttle_value = kwargs.get('throttle', 5)
285 285 self.sendData = self.initThrottle(self.throttle_value)
286 286 self.setup()
287 287
288 288 def setup(self):
289 289
290 290 self.data = {}
291 291 self.data['times'] = []
292 292 for plottype in self.plottypes:
293 293 self.data[plottype] = {}
294 294 self.data['noise'] = {}
295 295 self.data['throttle'] = self.throttle_value
296 296 self.data['ENDED'] = False
297 297 self.isConfig = True
298 298 self.data_web = {}
299 299
300 300 def event_monitor(self, monitor):
301 301
302 302 events = {}
303 303
304 304 for name in dir(zmq):
305 305 if name.startswith('EVENT_'):
306 306 value = getattr(zmq, name)
307 307 events[value] = name
308 308
309 309 while monitor.poll():
310 310 evt = recv_monitor_message(monitor)
311 311 if evt['event'] == 32:
312 312 self.connections += 1
313 313 if evt['event'] == 512:
314 314 pass
315 315 if self.connections == 0 and self.started is True:
316 316 self.ended = True
317 317 # send('ENDED')
318 318 evt.update({'description': events[evt['event']]})
319 319
320 320 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
321 321 break
322 322 monitor.close()
323 323 print("event monitor thread done!")
324 324
325 325 def initThrottle(self, throttle_value):
326 326
327 327 @throttle(seconds=throttle_value)
328 328 def sendDataThrottled(fn_sender, data):
329 329 fn_sender(data)
330 330
331 331 return sendDataThrottled
332 332
333 333 def send(self, data):
334 334 # print '[sending] data=%s size=%s' % (data.keys(), len(data['times']))
335 335 self.sender.send_pyobj(data)
336 336
337 337 def update(self):
338 338 t = self.dataOut.utctime
339 339 self.data['times'].append(t)
340 340 self.data['dataOut'] = self.dataOut
341 341 for plottype in self.plottypes:
342 342 if plottype == 'spc':
343 343 z = self.dataOut.data_spc/self.dataOut.normFactor
344 344 self.data[plottype] = 10*numpy.log10(z)
345 345 self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor)
346 if plottype == 'cspc':
347 jcoherence = self.dataOut.data_cspc/numpy.sqrt(self.dataOut.data_spc*self.dataOut.data_spc)
348 self.data['cspc_coh'] = numpy.abs(jcoherence)
349 self.data['cspc_phase'] = numpy.arctan2(jcoherence.imag, jcoherence.real)*180/numpy.pi
346 350 if plottype == 'rti':
347 351 self.data[plottype][t] = self.dataOut.getPower()
348 352 if plottype == 'snr':
349 353 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR)
350 354 if plottype == 'dop':
351 355 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP)
356 if plottype == 'mean':
357 self.data[plottype][t] = self.dataOut.data_MEAN
358 if plottype == 'std':
359 self.data[plottype][t] = self.dataOut.data_STD
352 360 if plottype == 'coh':
353 361 self.data[plottype][t] = self.dataOut.getCoherence()
354 362 if plottype == 'phase':
355 363 self.data[plottype][t] = self.dataOut.getCoherence(phase=True)
356 364 if self.realtime:
357 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
358 self.data_web['timestamp'] = t
365 self.data_web['timestamp'] = t
359 366 if plottype == 'spc':
360 367 self.data_web[plottype] = roundFloats(decimate(self.data[plottype]).tolist())
368 elif plottype == 'cspc':
369 self.data_web['cspc_coh'] = roundFloats(decimate(self.data['cspc_coh']).tolist())
370 self.data_web['cspc_phase'] = roundFloats(decimate(self.data['cspc_phase']).tolist())
371 elif plottype == 'noise':
372 self.data_web['noise'] = roundFloats(self.data['noise'][t].tolist())
361 373 else:
362 374 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
363 375 self.data_web['interval'] = self.dataOut.getTimeInterval()
364 376 self.data_web['type'] = plottype
365 377
366 378 def run(self):
367 379
368 380 print '[Starting] {} from {}'.format(self.name, self.address)
369 381
370 382 self.context = zmq.Context()
371 383 self.receiver = self.context.socket(zmq.PULL)
372 384 self.receiver.bind(self.address)
373 385 monitor = self.receiver.get_monitor_socket()
374 386 self.sender = self.context.socket(zmq.PUB)
375 387 if self.realtime:
376 388 self.sender_web = self.context.socket(zmq.PUB)
377 389 self.sender_web.connect(self.plot_address)
378 390 time.sleep(1)
379 391 self.sender.bind("ipc:///tmp/zmq.plots")
380 392
381 393 t = Thread(target=self.event_monitor, args=(monitor,))
382 394 t.start()
383 395
384 396 while True:
385 397 self.dataOut = self.receiver.recv_pyobj()
386 398 # print '[Receiving] {} - {}'.format(self.dataOut.type,
387 399 # self.dataOut.datatime.ctime())
388 400
389 401 self.update()
390 402
391 403 if self.dataOut.finished is True:
392 404 self.send(self.data)
393 405 self.connections -= 1
394 406 if self.connections == 0 and self.started:
395 407 self.ended = True
396 408 self.data['ENDED'] = True
397 409 self.send(self.data)
398 410 self.setup()
399 411 else:
400 412 if self.realtime:
401 413 self.send(self.data)
402 414 self.sender_web.send_string(json.dumps(self.data_web))
403 415 else:
404 416 self.sendData(self.send, self.data)
405 417 self.started = True
406 418
407 419 return
408 420
409 421 def sendToWeb(self):
410 422
411 423 if not self.isWebConfig:
412 424 context = zmq.Context()
413 425 sender_web_config = context.socket(zmq.PUB)
414 426 if 'tcp://' in self.plot_address:
415 427 dum, address, port = self.plot_address.split(':')
416 428 conf_address = '{}:{}:{}'.format(dum, address, int(port)+1)
417 429 else:
418 430 conf_address = self.plot_address + '.config'
419 431 sender_web_config.bind(conf_address)
420 432 time.sleep(1)
421 433 for kwargs in self.operationKwargs.values():
422 434 if 'plot' in kwargs:
423 435 print '[Sending] Config data to web for {}'.format(kwargs['code'].upper())
424 436 sender_web_config.send_string(json.dumps(kwargs))
425 437 self.isWebConfig = True
General Comments 0
You need to be logged in to leave comments. Login now