##// END OF EJS Templates
schain
RSS
cambios para amisr ISR getNoise
joabAM -
r1468:1fb12846b034
parent child
Show More
Side by Side Unified
Context file:
r1468:1fb12846b034 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1076 +1,1076
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Definition of diferent Data objects for different types of data
6 6
7 7 Here you will find the diferent data objects for the different types
8 8 of data, this data objects must be used as dataIn or dataOut objects in
9 9 processing units and operations. Currently the supported data objects are:
10 10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 11 """
12 12
13 13 import copy
14 14 import numpy
15 15 import datetime
16 16 import json
17 17
18 18 import schainpy.admin
19 19 from schainpy.utils import log
20 20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 21 from schainpy.model.data import _noise
22 22
23 23
24 24 def getNumpyDtype(dataTypeCode):
25 25
26 26 if dataTypeCode == 0:
27 27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 28 elif dataTypeCode == 1:
29 29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 30 elif dataTypeCode == 2:
31 31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 32 elif dataTypeCode == 3:
33 33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 34 elif dataTypeCode == 4:
35 35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 36 elif dataTypeCode == 5:
37 37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 38 else:
39 39 raise ValueError('dataTypeCode was not defined')
40 40
41 41 return numpyDtype
42 42
43 43
44 44 def getDataTypeCode(numpyDtype):
45 45
46 46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 47 datatype = 0
48 48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 49 datatype = 1
50 50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 51 datatype = 2
52 52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 53 datatype = 3
54 54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 55 datatype = 4
56 56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 57 datatype = 5
58 58 else:
59 59 datatype = None
60 60
61 61 return datatype
62 62
63 63
64 64 def hildebrand_sekhon(data, navg):
65 65 """
66 66 This method is for the objective determination of the noise level in Doppler spectra. This
67 67 implementation technique is based on the fact that the standard deviation of the spectral
68 68 densities is equal to the mean spectral density for white Gaussian noise
69 69
70 70 Inputs:
71 71 Data : heights
72 72 navg : numbers of averages
73 73
74 74 Return:
75 75 mean : noise's level
76 76 """
77 77
78 78 sortdata = numpy.sort(data, axis=None)
79 79 '''
80 80 lenOfData = len(sortdata)
81 81 nums_min = lenOfData*0.2
82 82
83 83 if nums_min <= 5:
84 84
85 85 nums_min = 5
86 86
87 87 sump = 0.
88 88 sumq = 0.
89 89
90 90 j = 0
91 91 cont = 1
92 92
93 93 while((cont == 1)and(j < lenOfData)):
94 94
95 95 sump += sortdata[j]
96 96 sumq += sortdata[j]**2
97 97
98 98 if j > nums_min:
99 99 rtest = float(j)/(j-1) + 1.0/navg
100 100 if ((sumq*j) > (rtest*sump**2)):
101 101 j = j - 1
102 102 sump = sump - sortdata[j]
103 103 sumq = sumq - sortdata[j]**2
104 104 cont = 0
105 105
106 106 j += 1
107 107
108 108 lnoise = sump / j
109 109 '''
110 110 return _noise.hildebrand_sekhon(sortdata, navg)
111 111
112 112
113 113 class Beam:
114 114
115 115 def __init__(self):
116 116 self.codeList = []
117 117 self.azimuthList = []
118 118 self.zenithList = []
119 119
120 120
121 121
122 122 class GenericData(object):
123 123
124 124 flagNoData = True
125 125
126 126 def copy(self, inputObj=None):
127 127
128 128 if inputObj == None:
129 129 return copy.deepcopy(self)
130 130
131 131 for key in list(inputObj.__dict__.keys()):
132 132
133 133 attribute = inputObj.__dict__[key]
134 134
135 135 # If this attribute is a tuple or list
136 136 if type(inputObj.__dict__[key]) in (tuple, list):
137 137 self.__dict__[key] = attribute[:]
138 138 continue
139 139
140 140 # If this attribute is another object or instance
141 141 if hasattr(attribute, '__dict__'):
142 142 self.__dict__[key] = attribute.copy()
143 143 continue
144 144
145 145 self.__dict__[key] = inputObj.__dict__[key]
146 146
147 147 def deepcopy(self):
148 148
149 149 return copy.deepcopy(self)
150 150
151 151 def isEmpty(self):
152 152
153 153 return self.flagNoData
154 154
155 155 def isReady(self):
156 156
157 157 return not self.flagNoData
158 158
159 159
160 160 class JROData(GenericData):
161 161
162 162 systemHeaderObj = SystemHeader()
163 163 radarControllerHeaderObj = RadarControllerHeader()
164 164 type = None
165 165 datatype = None # dtype but in string
166 166 nProfiles = None
167 167 heightList = None
168 168 channelList = None
169 169 flagDiscontinuousBlock = False
170 170 useLocalTime = False
171 171 utctime = None
172 172 timeZone = None
173 173 dstFlag = None
174 174 errorCount = None
175 175 blocksize = None
176 176 flagDecodeData = False # asumo q la data no esta decodificada
177 177 flagDeflipData = False # asumo q la data no esta sin flip
178 178 flagShiftFFT = False
179 179 nCohInt = None
180 180 windowOfFilter = 1
181 181 C = 3e8
182 182 frequency = 49.92e6
183 183 realtime = False
184 184 beacon_heiIndexList = None
185 185 last_block = None
186 186 blocknow = None
187 187 azimuth = None
188 188 zenith = None
189 189 beam = Beam()
190 190 profileIndex = None
191 191 error = None
192 192 data = None
193 193 nmodes = None
194 194 metadata_list = ['heightList', 'timeZone', 'type']
195 195 codeList = None
196 196 azimuthList = None
197 197 elevationList = None
198 198
199 199 def __str__(self):
200 200
201 201 return '{} - {}'.format(self.type, self.datatime())
202 202
203 203 def getNoise(self):
204 204
205 205 raise NotImplementedError
206 206
207 207 @property
208 208 def nChannels(self):
209 209
210 210 return len(self.channelList)
211 211
212 212 @property
213 213 def channelIndexList(self):
214 214
215 215 return list(range(self.nChannels))
216 216
217 217 @property
218 218 def nHeights(self):
219 219
220 220 return len(self.heightList)
221 221
222 222 def getDeltaH(self):
223 223
224 224 return self.heightList[1] - self.heightList[0]
225 225
226 226 @property
227 227 def ltctime(self):
228 228
229 229 if self.useLocalTime:
230 230 return self.utctime - self.timeZone * 60
231 231
232 232 return self.utctime
233 233
234 234 @property
235 235 def datatime(self):
236 236
237 237 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
238 238 return datatimeValue
239 239
240 240 def getTimeRange(self):
241 241
242 242 datatime = []
243 243
244 244 datatime.append(self.ltctime)
245 245 datatime.append(self.ltctime + self.timeInterval + 1)
246 246
247 247 datatime = numpy.array(datatime)
248 248
249 249 return datatime
250 250
251 251 def getFmaxTimeResponse(self):
252 252
253 253 period = (10**-6) * self.getDeltaH() / (0.15)
254 254
255 255 PRF = 1. / (period * self.nCohInt)
256 256
257 257 fmax = PRF
258 258
259 259 return fmax
260 260
261 261 def getFmax(self):
262 262 PRF = 1. / (self.ippSeconds * self.nCohInt)
263 263
264 264 fmax = PRF
265 265 return fmax
266 266
267 267 def getVmax(self):
268 268
269 269 _lambda = self.C / self.frequency
270 270
271 271 vmax = self.getFmax() * _lambda / 2
272 272
273 273 return vmax
274 274
275 275 @property
276 276 def ippSeconds(self):
277 277 '''
278 278 '''
279 279 return self.radarControllerHeaderObj.ippSeconds
280 280
281 281 @ippSeconds.setter
282 282 def ippSeconds(self, ippSeconds):
283 283 '''
284 284 '''
285 285 self.radarControllerHeaderObj.ippSeconds = ippSeconds
286 286
287 287 @property
288 288 def code(self):
289 289 '''
290 290 '''
291 291 return self.radarControllerHeaderObj.code
292 292
293 293 @code.setter
294 294 def code(self, code):
295 295 '''
296 296 '''
297 297 self.radarControllerHeaderObj.code = code
298 298
299 299 @property
300 300 def nCode(self):
301 301 '''
302 302 '''
303 303 return self.radarControllerHeaderObj.nCode
304 304
305 305 @nCode.setter
306 306 def nCode(self, ncode):
307 307 '''
308 308 '''
309 309 self.radarControllerHeaderObj.nCode = ncode
310 310
311 311 @property
312 312 def nBaud(self):
313 313 '''
314 314 '''
315 315 return self.radarControllerHeaderObj.nBaud
316 316
317 317 @nBaud.setter
318 318 def nBaud(self, nbaud):
319 319 '''
320 320 '''
321 321 self.radarControllerHeaderObj.nBaud = nbaud
322 322
323 323 @property
324 324 def ipp(self):
325 325 '''
326 326 '''
327 327 return self.radarControllerHeaderObj.ipp
328 328
329 329 @ipp.setter
330 330 def ipp(self, ipp):
331 331 '''
332 332 '''
333 333 self.radarControllerHeaderObj.ipp = ipp
334 334
335 335 @property
336 336 def metadata(self):
337 337 '''
338 338 '''
339 339
340 340 return {attr: getattr(self, attr) for attr in self.metadata_list}
341 341
342 342
343 343 class Voltage(JROData):
344 344
345 345 dataPP_POW = None
346 346 dataPP_DOP = None
347 347 dataPP_WIDTH = None
348 348 dataPP_SNR = None
349 349
350 350 def __init__(self):
351 351 '''
352 352 Constructor
353 353 '''
354 354
355 355 self.useLocalTime = True
356 356 self.radarControllerHeaderObj = RadarControllerHeader()
357 357 self.systemHeaderObj = SystemHeader()
358 358 self.type = "Voltage"
359 359 self.data = None
360 360 self.nProfiles = None
361 361 self.heightList = None
362 362 self.channelList = None
363 363 self.flagNoData = True
364 364 self.flagDiscontinuousBlock = False
365 365 self.utctime = None
366 366 self.timeZone = 0
367 367 self.dstFlag = None
368 368 self.errorCount = None
369 369 self.nCohInt = None
370 370 self.blocksize = None
371 371 self.flagCohInt = False
372 372 self.flagDecodeData = False # asumo q la data no esta decodificada
373 373 self.flagDeflipData = False # asumo q la data no esta sin flip
374 374 self.flagShiftFFT = False
375 375 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
376 376 self.profileIndex = 0
377 377 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
378 378 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
379 379
380 380 def getNoisebyHildebrand(self, channel=None):
381 381 """
382 382 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
383 383
384 384 Return:
385 385 noiselevel
386 386 """
387 387
388 388 if channel != None:
389 389 data = self.data[channel]
390 390 nChannels = 1
391 391 else:
392 392 data = self.data
393 393 nChannels = self.nChannels
394 394
395 395 noise = numpy.zeros(nChannels)
396 396 power = data * numpy.conjugate(data)
397 397
398 398 for thisChannel in range(nChannels):
399 399 if nChannels == 1:
400 400 daux = power[:].real
401 401 else:
402 402 daux = power[thisChannel, :].real
403 403 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
404 404
405 405 return noise
406 406
407 407 def getNoise(self, type=1, channel=None):
408 408
409 409 if type == 1:
410 410 noise = self.getNoisebyHildebrand(channel)
411 411
412 412 return noise
413 413
414 414 def getPower(self, channel=None):
415 415
416 416 if channel != None:
417 417 data = self.data[channel]
418 418 else:
419 419 data = self.data
420 420
421 421 power = data * numpy.conjugate(data)
422 422 powerdB = 10 * numpy.log10(power.real)
423 423 powerdB = numpy.squeeze(powerdB)
424 424
425 425 return powerdB
426 426
427 427 @property
428 428 def timeInterval(self):
429 429
430 430 return self.ippSeconds * self.nCohInt
431 431
432 432 noise = property(getNoise, "I'm the 'nHeights' property.")
433 433
434 434
435 435 class Spectra(JROData):
436 436
437 437 def __init__(self):
438 438 '''
439 439 Constructor
440 440 '''
441 441
442 442 self.data_dc = None
443 443 self.data_spc = None
444 444 self.data_cspc = None
445 445 self.useLocalTime = True
446 446 self.radarControllerHeaderObj = RadarControllerHeader()
447 447 self.systemHeaderObj = SystemHeader()
448 448 self.type = "Spectra"
449 449 self.timeZone = 0
450 450 self.nProfiles = None
451 451 self.heightList = None
452 452 self.channelList = None
453 453 self.pairsList = None
454 454 self.flagNoData = True
455 455 self.flagDiscontinuousBlock = False
456 456 self.utctime = None
457 457 self.nCohInt = None
458 458 self.nIncohInt = None
459 459 self.blocksize = None
460 460 self.nFFTPoints = None
461 461 self.wavelength = None
462 462 self.flagDecodeData = False # asumo q la data no esta decodificada
463 463 self.flagDeflipData = False # asumo q la data no esta sin flip
464 464 self.flagShiftFFT = False
465 465 self.ippFactor = 1
466 466 self.beacon_heiIndexList = []
467 467 self.noise_estimation = None
468 468 self.codeList = []
469 469 self.azimuthList = []
470 470 self.elevationList = []
471 471 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
472 472 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
473 473
474 474
475 475 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
476 476 """
477 477 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
478 478
479 479 Return:
480 480 noiselevel
481 481 """
482 482
483 483 noise = numpy.zeros(self.nChannels)
484 484 for channel in range(self.nChannels):
485 485 daux = self.data_spc[channel,
486 486 xmin_index:xmax_index, ymin_index:ymax_index]
487 487 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
488 488
489 489 return noise
490 490
491 491 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
492
492
493 493 if self.noise_estimation is not None:
494 494 # this was estimated by getNoise Operation defined in jroproc_spectra.py
495 495 return self.noise_estimation
496 496 else:
497 497 noise = self.getNoisebyHildebrand(
498 498 xmin_index, xmax_index, ymin_index, ymax_index)
499 499 return noise
500 500
501 501 def getFreqRangeTimeResponse(self, extrapoints=0):
502 502
503 503 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
504 504 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
505 505
506 506 return freqrange
507 507
508 508 def getAcfRange(self, extrapoints=0):
509 509
510 510 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
511 511 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
512 512
513 513 return freqrange
514 514
515 515 def getFreqRange(self, extrapoints=0):
516 516
517 517 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
518 518 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
519 519
520 520 return freqrange
521 521
522 522 def getVelRange(self, extrapoints=0):
523 523
524 524 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
525 525 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
526 526
527 527 if self.nmodes:
528 528 return velrange/self.nmodes
529 529 else:
530 530 return velrange
531 531
532 532 @property
533 533 def nPairs(self):
534 534
535 535 return len(self.pairsList)
536 536
537 537 @property
538 538 def pairsIndexList(self):
539 539
540 540 return list(range(self.nPairs))
541 541
542 542 @property
543 543 def normFactor(self):
544 544
545 545 pwcode = 1
546 546
547 547 if self.flagDecodeData:
548 548 pwcode = numpy.sum(self.code[0]**2)
549 549 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
550 550 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
551 551
552 552 return normFactor
553 553
554 554 @property
555 555 def flag_cspc(self):
556 556
557 557 if self.data_cspc is None:
558 558 return True
559 559
560 560 return False
561 561
562 562 @property
563 563 def flag_dc(self):
564 564
565 565 if self.data_dc is None:
566 566 return True
567 567
568 568 return False
569 569
570 570 @property
571 571 def timeInterval(self):
572 572
573 573 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
574 574 if self.nmodes:
575 575 return self.nmodes*timeInterval
576 576 else:
577 577 return timeInterval
578 578
579 579 def getPower(self):
580 580
581 581 factor = self.normFactor
582 582 z = self.data_spc / factor
583 583 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
584 584 avg = numpy.average(z, axis=1)
585 585
586 586 return 10 * numpy.log10(avg)
587 587
588 588 def getCoherence(self, pairsList=None, phase=False):
589 589
590 590 z = []
591 591 if pairsList is None:
592 592 pairsIndexList = self.pairsIndexList
593 593 else:
594 594 pairsIndexList = []
595 595 for pair in pairsList:
596 596 if pair not in self.pairsList:
597 597 raise ValueError("Pair %s is not in dataOut.pairsList" % (
598 598 pair))
599 599 pairsIndexList.append(self.pairsList.index(pair))
600 600 for i in range(len(pairsIndexList)):
601 601 pair = self.pairsList[pairsIndexList[i]]
602 602 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
603 603 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
604 604 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
605 605 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
606 606 if phase:
607 607 data = numpy.arctan2(avgcoherenceComplex.imag,
608 608 avgcoherenceComplex.real) * 180 / numpy.pi
609 609 else:
610 610 data = numpy.abs(avgcoherenceComplex)
611 611
612 612 z.append(data)
613 613
614 614 return numpy.array(z)
615 615
616 616 def setValue(self, value):
617 617
618 618 print("This property should not be initialized")
619 619
620 620 return
621 621
622 622 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
623 623
624 624
625 625 class SpectraHeis(Spectra):
626 626
627 627 def __init__(self):
628 628
629 629 self.radarControllerHeaderObj = RadarControllerHeader()
630 630 self.systemHeaderObj = SystemHeader()
631 631 self.type = "SpectraHeis"
632 632 self.nProfiles = None
633 633 self.heightList = None
634 634 self.channelList = None
635 635 self.flagNoData = True
636 636 self.flagDiscontinuousBlock = False
637 637 self.utctime = None
638 638 self.blocksize = None
639 639 self.profileIndex = 0
640 640 self.nCohInt = 1
641 641 self.nIncohInt = 1
642 642
643 643 @property
644 644 def normFactor(self):
645 645 pwcode = 1
646 646 if self.flagDecodeData:
647 647 pwcode = numpy.sum(self.code[0]**2)
648 648
649 649 normFactor = self.nIncohInt * self.nCohInt * pwcode
650 650
651 651 return normFactor
652 652
653 653 @property
654 654 def timeInterval(self):
655 655
656 656 return self.ippSeconds * self.nCohInt * self.nIncohInt
657 657
658 658
659 659 class Fits(JROData):
660 660
661 661 def __init__(self):
662 662
663 663 self.type = "Fits"
664 664 self.nProfiles = None
665 665 self.heightList = None
666 666 self.channelList = None
667 667 self.flagNoData = True
668 668 self.utctime = None
669 669 self.nCohInt = 1
670 670 self.nIncohInt = 1
671 671 self.useLocalTime = True
672 672 self.profileIndex = 0
673 673 self.timeZone = 0
674 674
675 675 def getTimeRange(self):
676 676
677 677 datatime = []
678 678
679 679 datatime.append(self.ltctime)
680 680 datatime.append(self.ltctime + self.timeInterval)
681 681
682 682 datatime = numpy.array(datatime)
683 683
684 684 return datatime
685 685
686 686 def getChannelIndexList(self):
687 687
688 688 return list(range(self.nChannels))
689 689
690 690 def getNoise(self, type=1):
691 691
692 692
693 693 if type == 1:
694 694 noise = self.getNoisebyHildebrand()
695 695
696 696 if type == 2:
697 697 noise = self.getNoisebySort()
698 698
699 699 if type == 3:
700 700 noise = self.getNoisebyWindow()
701 701
702 702 return noise
703 703
704 704 @property
705 705 def timeInterval(self):
706 706
707 707 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
708 708
709 709 return timeInterval
710 710
711 711 @property
712 712 def ippSeconds(self):
713 713 '''
714 714 '''
715 715 return self.ipp_sec
716 716
717 717 noise = property(getNoise, "I'm the 'nHeights' property.")
718 718
719 719
720 720 class Correlation(JROData):
721 721
722 722 def __init__(self):
723 723 '''
724 724 Constructor
725 725 '''
726 726 self.radarControllerHeaderObj = RadarControllerHeader()
727 727 self.systemHeaderObj = SystemHeader()
728 728 self.type = "Correlation"
729 729 self.data = None
730 730 self.dtype = None
731 731 self.nProfiles = None
732 732 self.heightList = None
733 733 self.channelList = None
734 734 self.flagNoData = True
735 735 self.flagDiscontinuousBlock = False
736 736 self.utctime = None
737 737 self.timeZone = 0
738 738 self.dstFlag = None
739 739 self.errorCount = None
740 740 self.blocksize = None
741 741 self.flagDecodeData = False # asumo q la data no esta decodificada
742 742 self.flagDeflipData = False # asumo q la data no esta sin flip
743 743 self.pairsList = None
744 744 self.nPoints = None
745 745
746 746 def getPairsList(self):
747 747
748 748 return self.pairsList
749 749
750 750 def getNoise(self, mode=2):
751 751
752 752 indR = numpy.where(self.lagR == 0)[0][0]
753 753 indT = numpy.where(self.lagT == 0)[0][0]
754 754
755 755 jspectra0 = self.data_corr[:, :, indR, :]
756 756 jspectra = copy.copy(jspectra0)
757 757
758 758 num_chan = jspectra.shape[0]
759 759 num_hei = jspectra.shape[2]
760 760
761 761 freq_dc = jspectra.shape[1] / 2
762 762 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
763 763
764 764 if ind_vel[0] < 0:
765 765 ind_vel[list(range(0, 1))] = ind_vel[list(
766 766 range(0, 1))] + self.num_prof
767 767
768 768 if mode == 1:
769 769 jspectra[:, freq_dc, :] = (
770 770 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
771 771
772 772 if mode == 2:
773 773
774 774 vel = numpy.array([-2, -1, 1, 2])
775 775 xx = numpy.zeros([4, 4])
776 776
777 777 for fil in range(4):
778 778 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
779 779
780 780 xx_inv = numpy.linalg.inv(xx)
781 781 xx_aux = xx_inv[0, :]
782 782
783 783 for ich in range(num_chan):
784 784 yy = jspectra[ich, ind_vel, :]
785 785 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
786 786
787 787 junkid = jspectra[ich, freq_dc, :] <= 0
788 788 cjunkid = sum(junkid)
789 789
790 790 if cjunkid.any():
791 791 jspectra[ich, freq_dc, junkid.nonzero()] = (
792 792 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
793 793
794 794 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
795 795
796 796 return noise
797 797
798 798 @property
799 799 def timeInterval(self):
800 800
801 801 return self.ippSeconds * self.nCohInt * self.nProfiles
802 802
803 803 def splitFunctions(self):
804 804
805 805 pairsList = self.pairsList
806 806 ccf_pairs = []
807 807 acf_pairs = []
808 808 ccf_ind = []
809 809 acf_ind = []
810 810 for l in range(len(pairsList)):
811 811 chan0 = pairsList[l][0]
812 812 chan1 = pairsList[l][1]
813 813
814 814 # Obteniendo pares de Autocorrelacion
815 815 if chan0 == chan1:
816 816 acf_pairs.append(chan0)
817 817 acf_ind.append(l)
818 818 else:
819 819 ccf_pairs.append(pairsList[l])
820 820 ccf_ind.append(l)
821 821
822 822 data_acf = self.data_cf[acf_ind]
823 823 data_ccf = self.data_cf[ccf_ind]
824 824
825 825 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
826 826
827 827 @property
828 828 def normFactor(self):
829 829 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
830 830 acf_pairs = numpy.array(acf_pairs)
831 831 normFactor = numpy.zeros((self.nPairs, self.nHeights))
832 832
833 833 for p in range(self.nPairs):
834 834 pair = self.pairsList[p]
835 835
836 836 ch0 = pair[0]
837 837 ch1 = pair[1]
838 838
839 839 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
840 840 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
841 841 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
842 842
843 843 return normFactor
844 844
845 845
846 846 class Parameters(Spectra):
847 847
848 848 groupList = None # List of Pairs, Groups, etc
849 849 data_param = None # Parameters obtained
850 850 data_pre = None # Data Pre Parametrization
851 851 data_SNR = None # Signal to Noise Ratio
852 852 abscissaList = None # Abscissa, can be velocities, lags or time
853 853 utctimeInit = None # Initial UTC time
854 854 paramInterval = None # Time interval to calculate Parameters in seconds
855 855 useLocalTime = True
856 856 # Fitting
857 857 data_error = None # Error of the estimation
858 858 constants = None
859 859 library = None
860 860 # Output signal
861 861 outputInterval = None # Time interval to calculate output signal in seconds
862 862 data_output = None # Out signal
863 863 nAvg = None
864 864 noise_estimation = None
865 865 GauSPC = None # Fit gaussian SPC
866 866
867 867 def __init__(self):
868 868 '''
869 869 Constructor
870 870 '''
871 871 self.radarControllerHeaderObj = RadarControllerHeader()
872 872 self.systemHeaderObj = SystemHeader()
873 873 self.type = "Parameters"
874 874 self.timeZone = 0
875 875
876 876 def getTimeRange1(self, interval):
877 877
878 878 datatime = []
879 879
880 880 if self.useLocalTime:
881 881 time1 = self.utctimeInit - self.timeZone * 60
882 882 else:
883 883 time1 = self.utctimeInit
884 884
885 885 datatime.append(time1)
886 886 datatime.append(time1 + interval)
887 887 datatime = numpy.array(datatime)
888 888
889 889 return datatime
890 890
891 891 @property
892 892 def timeInterval(self):
893 893
894 894 if hasattr(self, 'timeInterval1'):
895 895 return self.timeInterval1
896 896 else:
897 897 return self.paramInterval
898 898
899 899 def setValue(self, value):
900 900
901 901 print("This property should not be initialized")
902 902
903 903 return
904 904
905 905 def getNoise(self):
906 906
907 907 return self.spc_noise
908 908
909 909 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
910 910
911 911
912 912 class PlotterData(object):
913 913 '''
914 914 Object to hold data to be plotted
915 915 '''
916 916
917 917 MAXNUMX = 200
918 918 MAXNUMY = 200
919 919
920 920 def __init__(self, code, exp_code, localtime=True):
921 921
922 922 self.key = code
923 923 self.exp_code = exp_code
924 924 self.ready = False
925 925 self.flagNoData = False
926 926 self.localtime = localtime
927 927 self.data = {}
928 928 self.meta = {}
929 929 self.__heights = []
930 930
931 931 def __str__(self):
932 932 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
933 933 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
934 934
935 935 def __len__(self):
936 936 return len(self.data)
937 937
938 938 def __getitem__(self, key):
939 939 if isinstance(key, int):
940 940 return self.data[self.times[key]]
941 941 elif isinstance(key, str):
942 942 ret = numpy.array([self.data[x][key] for x in self.times])
943 943 if ret.ndim > 1:
944 944 ret = numpy.swapaxes(ret, 0, 1)
945 945 return ret
946 946
947 947 def __contains__(self, key):
948 948 return key in self.data[self.min_time]
949 949
950 950 def setup(self):
951 951 '''
952 952 Configure object
953 953 '''
954 954 self.type = ''
955 955 self.ready = False
956 956 del self.data
957 957 self.data = {}
958 958 self.__heights = []
959 959 self.__all_heights = set()
960 960
961 961 def shape(self, key):
962 962 '''
963 963 Get the shape of the one-element data for the given key
964 964 '''
965 965
966 966 if len(self.data[self.min_time][key]):
967 967 return self.data[self.min_time][key].shape
968 968 return (0,)
969 969
970 970 def update(self, data, tm, meta={}):
971 971 '''
972 972 Update data object with new dataOut
973 973 '''
974 974
975 975 self.data[tm] = data
976 976
977 977 for key, value in meta.items():
978 978 setattr(self, key, value)
979 979
980 980 def normalize_heights(self):
981 981 '''
982 982 Ensure same-dimension of the data for different heighList
983 983 '''
984 984
985 985 H = numpy.array(list(self.__all_heights))
986 986 H.sort()
987 987 for key in self.data:
988 988 shape = self.shape(key)[:-1] + H.shape
989 989 for tm, obj in list(self.data[key].items()):
990 990 h = self.__heights[self.times.tolist().index(tm)]
991 991 if H.size == h.size:
992 992 continue
993 993 index = numpy.where(numpy.in1d(H, h))[0]
994 994 dummy = numpy.zeros(shape) + numpy.nan
995 995 if len(shape) == 2:
996 996 dummy[:, index] = obj
997 997 else:
998 998 dummy[index] = obj
999 999 self.data[key][tm] = dummy
1000 1000
1001 1001 self.__heights = [H for tm in self.times]
1002 1002
1003 1003 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1004 1004 '''
1005 1005 Convert data to json
1006 1006 '''
1007 1007
1008 1008 meta = {}
1009 1009 meta['xrange'] = []
1010 1010 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1011 1011 tmp = self.data[tm][self.key]
1012 1012 shape = tmp.shape
1013 1013 if len(shape) == 2:
1014 1014 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1015 1015 elif len(shape) == 3:
1016 1016 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1017 1017 data = self.roundFloats(
1018 1018 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1019 1019 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1020 1020 else:
1021 1021 data = self.roundFloats(self.data[tm][self.key].tolist())
1022 1022
1023 1023 ret = {
1024 1024 'plot': plot_name,
1025 1025 'code': self.exp_code,
1026 1026 'time': float(tm),
1027 1027 'data': data,
1028 1028 }
1029 1029 meta['type'] = plot_type
1030 1030 meta['interval'] = float(self.interval)
1031 1031 meta['localtime'] = self.localtime
1032 1032 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1033 1033 meta.update(self.meta)
1034 1034 ret['metadata'] = meta
1035 1035 return json.dumps(ret)
1036 1036
1037 1037 @property
1038 1038 def times(self):
1039 1039 '''
1040 1040 Return the list of times of the current data
1041 1041 '''
1042 1042
1043 1043 ret = [t for t in self.data]
1044 1044 ret.sort()
1045 1045 return numpy.array(ret)
1046 1046
1047 1047 @property
1048 1048 def min_time(self):
1049 1049 '''
1050 1050 Return the minimun time value
1051 1051 '''
1052 1052
1053 1053 return self.times[0]
1054 1054
1055 1055 @property
1056 1056 def max_time(self):
1057 1057 '''
1058 1058 Return the maximun time value
1059 1059 '''
1060 1060
1061 1061 return self.times[-1]
1062 1062
1063 1063 # @property
1064 1064 # def heights(self):
1065 1065 # '''
1066 1066 # Return the list of heights of the current data
1067 1067 # '''
1068 1068
1069 1069 # return numpy.array(self.__heights[-1])
1070 1070
1071 1071 @staticmethod
1072 1072 def roundFloats(obj):
1073 1073 if isinstance(obj, list):
1074 1074 return list(map(PlotterData.roundFloats, obj))
1075 1075 elif isinstance(obj, float):
1076 1076 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,961 +1,962
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13 from itertools import combinations
14 14
15 15
16 16 class SpectraPlot(Plot):
17 17 '''
18 18 Plot for Spectra data
19 19 '''
20 20
21 21 CODE = 'spc'
22 22 colormap = 'jet'
23 23 plot_type = 'pcolor'
24 24 buffering = False
25 25 channelList = []
26 26
27 27 def setup(self):
28 28
29 29 self.nplots = len(self.data.channels)
30 30 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
31 31 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
32 32 self.height = 2.6 * self.nrows
33 33
34 34 self.cb_label = 'dB'
35 35 if self.showprofile:
36 36 self.width = 4 * self.ncols
37 37 else:
38 38 self.width = 3.5 * self.ncols
39 39 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
40 40 self.ylabel = 'Range [km]'
41 41
42 42
43 43 def update_list(self,dataOut):
44 44 if len(self.channelList) == 0:
45 45 self.channelList = dataOut.channelList
46 46
47 47 def update(self, dataOut):
48 48
49 49 self.update_list(dataOut)
50 50 data = {}
51 51 meta = {}
52 52 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
53 53 data['spc'] = spc
54 54 data['rti'] = dataOut.getPower()
55 55 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
56 56 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
57 57 if self.CODE == 'spc_moments':
58 58 data['moments'] = dataOut.moments
59 59
60 60 return data, meta
61 61
62 62 def plot(self):
63 63 if self.xaxis == "frequency":
64 64 x = self.data.xrange[0]
65 65 self.xlabel = "Frequency (kHz)"
66 66 elif self.xaxis == "time":
67 67 x = self.data.xrange[1]
68 68 self.xlabel = "Time (ms)"
69 69 else:
70 70 x = self.data.xrange[2]
71 71 self.xlabel = "Velocity (m/s)"
72 72
73 73 if self.CODE == 'spc_moments':
74 74 x = self.data.xrange[2]
75 75 self.xlabel = "Velocity (m/s)"
76 76
77 77 self.titles = []
78 78 y = self.data.yrange
79 79 self.y = y
80 80
81 81 data = self.data[-1]
82 82 z = data['spc']
83 83
84 84 for n, ax in enumerate(self.axes):
85 85 noise = data['noise'][n]
86 86 if self.CODE == 'spc_moments':
87 87 mean = data['moments'][n, 1]
88 88 if ax.firsttime:
89 89 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
90 90 self.xmin = self.xmin if self.xmin else -self.xmax
91 91 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
92 92 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
93 93 ax.plt = ax.pcolormesh(x, y, z[n].T,
94 94 vmin=self.zmin,
95 95 vmax=self.zmax,
96 96 cmap=plt.get_cmap(self.colormap)
97 97 )
98 98
99 99 if self.showprofile:
100 100 ax.plt_profile = self.pf_axes[n].plot(
101 101 data['rti'][n], y)[0]
102 102 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
103 103 color="k", linestyle="dashed", lw=1)[0]
104 104 if self.CODE == 'spc_moments':
105 105 ax.plt_mean = ax.plot(mean, y, color='k')[0]
106 106 else:
107 107 ax.plt.set_array(z[n].T.ravel())
108 108 if self.showprofile:
109 109 ax.plt_profile.set_data(data['rti'][n], y)
110 110 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
111 111 if self.CODE == 'spc_moments':
112 112 ax.plt_mean.set_data(mean, y)
113 113 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
114 114
115 115
116 116 class CrossSpectraPlot(Plot):
117 117
118 118 CODE = 'cspc'
119 119 colormap = 'jet'
120 120 plot_type = 'pcolor'
121 121 zmin_coh = None
122 122 zmax_coh = None
123 123 zmin_phase = None
124 124 zmax_phase = None
125 125 realChannels = None
126 126 crossPairs = None
127 127
128 128 def setup(self):
129 129
130 130 self.ncols = 4
131 131 self.nplots = len(self.data.pairs) * 2
132 132 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
133 133 self.width = 3.1 * self.ncols
134 134 self.height = 2.6 * self.nrows
135 135 self.ylabel = 'Range [km]'
136 136 self.showprofile = False
137 137 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
138 138
139 139 def update(self, dataOut):
140 140
141 141 data = {}
142 142 meta = {}
143 143
144 144 spc = dataOut.data_spc
145 145 cspc = dataOut.data_cspc
146 146 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
147 147 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
148 148 meta['pairs'] = rawPairs
149 149
150 150 if self.crossPairs == None:
151 151 self.crossPairs = dataOut.pairsList
152 152
153 153 tmp = []
154 154
155 155 for n, pair in enumerate(meta['pairs']):
156 156
157 157 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
158 158 coh = numpy.abs(out)
159 159 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
160 160 tmp.append(coh)
161 161 tmp.append(phase)
162 162
163 163 data['cspc'] = numpy.array(tmp)
164 164
165 165 return data, meta
166 166
167 167 def plot(self):
168 168
169 169 if self.xaxis == "frequency":
170 170 x = self.data.xrange[0]
171 171 self.xlabel = "Frequency (kHz)"
172 172 elif self.xaxis == "time":
173 173 x = self.data.xrange[1]
174 174 self.xlabel = "Time (ms)"
175 175 else:
176 176 x = self.data.xrange[2]
177 177 self.xlabel = "Velocity (m/s)"
178 178
179 179 self.titles = []
180 180
181 181 y = self.data.yrange
182 182 self.y = y
183 183
184 184 data = self.data[-1]
185 185 cspc = data['cspc']
186 186
187 187 for n in range(len(self.data.pairs)):
188 188
189 189 pair = self.crossPairs[n]
190 190
191 191 coh = cspc[n*2]
192 192 phase = cspc[n*2+1]
193 193 ax = self.axes[2 * n]
194 194
195 195 if ax.firsttime:
196 196 ax.plt = ax.pcolormesh(x, y, coh.T,
197 197 vmin=self.zmin_coh,
198 198 vmax=self.zmax_coh,
199 199 cmap=plt.get_cmap(self.colormap_coh)
200 200 )
201 201 else:
202 202 ax.plt.set_array(coh.T.ravel())
203 203 self.titles.append(
204 204 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
205 205
206 206 ax = self.axes[2 * n + 1]
207 207 if ax.firsttime:
208 208 ax.plt = ax.pcolormesh(x, y, phase.T,
209 209 vmin=-180,
210 210 vmax=180,
211 211 cmap=plt.get_cmap(self.colormap_phase)
212 212 )
213 213 else:
214 214 ax.plt.set_array(phase.T.ravel())
215 215
216 216 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
217 217
218 218
219 219 class RTIPlot(Plot):
220 220 '''
221 221 Plot for RTI data
222 222 '''
223 223
224 224 CODE = 'rti'
225 225 colormap = 'jet'
226 226 plot_type = 'pcolorbuffer'
227 227 titles = None
228 228 channelList = []
229 229
230 230 def setup(self):
231 231 self.xaxis = 'time'
232 232 self.ncols = 1
233 233 #print("dataChannels ",self.data.channels)
234 234 self.nrows = len(self.data.channels)
235 235 self.nplots = len(self.data.channels)
236 236 self.ylabel = 'Range [km]'
237 237 self.xlabel = 'Time'
238 238 self.cb_label = 'dB'
239 239 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
240 240 self.titles = ['{} Channel {}'.format(
241 241 self.CODE.upper(), x) for x in range(self.nplots)]
242 242
243 243 def update_list(self,dataOut):
244 244
245 245 self.channelList = dataOut.channelList
246 246
247 247
248 248 def update(self, dataOut):
249 249 if len(self.channelList) == 0:
250 250 self.update_list(dataOut)
251 251 data = {}
252 252 meta = {}
253 253 data['rti'] = dataOut.getPower()
254 254 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
255 255 return data, meta
256 256
257 257 def plot(self):
258 258
259 259 self.x = self.data.times
260 260 self.y = self.data.yrange
261 261 self.z = self.data[self.CODE]
262 262 self.z = numpy.array(self.z, dtype=float)
263 263 self.z = numpy.ma.masked_invalid(self.z)
264 264
265 265 try:
266 266 if self.channelList != None:
267 267 self.titles = ['{} Channel {}'.format(
268 268 self.CODE.upper(), x) for x in self.channelList]
269 269 except:
270 270 if self.channelList.any() != None:
271 271 self.titles = ['{} Channel {}'.format(
272 272 self.CODE.upper(), x) for x in self.channelList]
273 273 if self.decimation is None:
274 274 x, y, z = self.fill_gaps(self.x, self.y, self.z)
275 275 else:
276 276 x, y, z = self.fill_gaps(*self.decimate())
277 277 dummy_var = self.axes #ExtraΓ±amente esto actualiza el valor axes
278 278 for n, ax in enumerate(self.axes):
279 279 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
280 280 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
281 281 data = self.data[-1]
282 282 if ax.firsttime:
283 283 ax.plt = ax.pcolormesh(x, y, z[n].T,
284 284 vmin=self.zmin,
285 285 vmax=self.zmax,
286 286 cmap=plt.get_cmap(self.colormap)
287 287 )
288 288 if self.showprofile:
289 289 ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0]
290 290
291 291 if "noise" in self.data:
292 292 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
293 293 color="k", linestyle="dashed", lw=1)[0]
294 294 else:
295 295 ax.collections.remove(ax.collections[0])
296 296 ax.plt = ax.pcolormesh(x, y, z[n].T,
297 297 vmin=self.zmin,
298 298 vmax=self.zmax,
299 299 cmap=plt.get_cmap(self.colormap)
300 300 )
301 301 if self.showprofile:
302 302 ax.plot_profile.set_data(data[self.CODE][n], self.y)
303 303 if "noise" in self.data:
304 304 ax.plot_noise.set_data(numpy.repeat(
305 305 data['noise'][n], len(self.y)), self.y)
306 306
307 307
308 308 class CoherencePlot(RTIPlot):
309 309 '''
310 310 Plot for Coherence data
311 311 '''
312 312
313 313 CODE = 'coh'
314 314
315 315 def setup(self):
316 316 self.xaxis = 'time'
317 317 self.ncols = 1
318 318 self.nrows = len(self.data.pairs)
319 319 self.nplots = len(self.data.pairs)
320 320 self.ylabel = 'Range [km]'
321 321 self.xlabel = 'Time'
322 322 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
323 323 if self.CODE == 'coh':
324 324 self.cb_label = ''
325 325 self.titles = [
326 326 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
327 327 else:
328 328 self.cb_label = 'Degrees'
329 329 self.titles = [
330 330 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
331 331
332 332 def update(self, dataOut):
333 333 self.update_list(dataOut)
334 334 data = {}
335 335 meta = {}
336 336 data['coh'] = dataOut.getCoherence()
337 337 meta['pairs'] = dataOut.pairsList
338 338
339 339
340 340 return data, meta
341 341
342 342 class PhasePlot(CoherencePlot):
343 343 '''
344 344 Plot for Phase map data
345 345 '''
346 346
347 347 CODE = 'phase'
348 348 colormap = 'seismic'
349 349
350 350 def update(self, dataOut):
351 351
352 352 data = {}
353 353 meta = {}
354 354 data['phase'] = dataOut.getCoherence(phase=True)
355 355 meta['pairs'] = dataOut.pairsList
356 356
357 357 return data, meta
358 358
359 359 class NoisePlot(Plot):
360 360 '''
361 361 Plot for noise
362 362 '''
363 363
364 364 CODE = 'noise'
365 365 plot_type = 'scatterbuffer'
366 366
367 367 def setup(self):
368 368 self.xaxis = 'time'
369 369 self.ncols = 1
370 370 self.nrows = 1
371 371 self.nplots = 1
372 372 self.ylabel = 'Intensity [dB]'
373 373 self.xlabel = 'Time'
374 374 self.titles = ['Noise']
375 375 self.colorbar = False
376 376 self.plots_adjust.update({'right': 0.85 })
377 377
378 378 def update(self, dataOut):
379 379
380 380 data = {}
381 381 meta = {}
382 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
382 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
383 data['noise'] = noise
383 384 meta['yrange'] = numpy.array([])
384 385
385 386 return data, meta
386 387
387 388 def plot(self):
388 389
389 390 x = self.data.times
390 391 xmin = self.data.min_time
391 392 xmax = xmin + self.xrange * 60 * 60
392 393 Y = self.data['noise']
393 394
394 395 if self.axes[0].firsttime:
395 self.ymin = numpy.nanmin(Y) - 5
396 self.ymax = numpy.nanmax(Y) + 5
396 if self.ymin is None: self.ymin = numpy.nanmin(Y) - 5
397 if self.ymax is None: self.ymax = numpy.nanmax(Y) + 5
397 398 for ch in self.data.channels:
398 399 y = Y[ch]
399 400 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
400 401 plt.legend(bbox_to_anchor=(1.18, 1.0))
401 402 else:
402 403 for ch in self.data.channels:
403 404 y = Y[ch]
404 405 self.axes[0].lines[ch].set_data(x, y)
405 406
406 407
407 408 class PowerProfilePlot(Plot):
408 409
409 410 CODE = 'pow_profile'
410 411 plot_type = 'scatter'
411 412
412 413 def setup(self):
413 414
414 415 self.ncols = 1
415 416 self.nrows = 1
416 417 self.nplots = 1
417 418 self.height = 4
418 419 self.width = 3
419 420 self.ylabel = 'Range [km]'
420 421 self.xlabel = 'Intensity [dB]'
421 422 self.titles = ['Power Profile']
422 423 self.colorbar = False
423 424
424 425 def update(self, dataOut):
425 426
426 427 data = {}
427 428 meta = {}
428 429 data[self.CODE] = dataOut.getPower()
429 430
430 431 return data, meta
431 432
432 433 def plot(self):
433 434
434 435 y = self.data.yrange
435 436 self.y = y
436 437
437 438 x = self.data[-1][self.CODE]
438 439
439 440 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
440 441 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
441 442
442 443 if self.axes[0].firsttime:
443 444 for ch in self.data.channels:
444 445 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
445 446 plt.legend()
446 447 else:
447 448 for ch in self.data.channels:
448 449 self.axes[0].lines[ch].set_data(x[ch], y)
449 450
450 451
451 452 class SpectraCutPlot(Plot):
452 453
453 454 CODE = 'spc_cut'
454 455 plot_type = 'scatter'
455 456 buffering = False
456 457 heights = []
457 458 channelList = []
458 459 maintitle = "Spectra Cuts"
459 460
460 461 def setup(self):
461 462
462 463 self.nplots = len(self.data.channels)
463 464 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
464 465 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
465 466 self.width = 3.6 * self.ncols + 1.5
466 467 self.height = 3 * self.nrows
467 468 self.ylabel = 'Power [dB]'
468 469 self.colorbar = False
469 470 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
470 471 if self.selectedHeight:
471 472 self.maintitle = "Spectra Cut for %d km " %(int(self.selectedHeight))
472 473
473 474 def update(self, dataOut):
474 475 if len(self.channelList) == 0:
475 476 self.channelList = dataOut.channelList
476 477
477 478 self.heights = dataOut.heightList
478 479 if self.selectedHeight:
479 480 index_list = numpy.where(self.heights >= self.selectedHeight)
480 481 self.height_index = index_list[0]
481 482 self.height_index = self.height_index[0]
482 483 #print(self.height_index)
483 484 data = {}
484 485 meta = {}
485 486 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
486 487 if self.selectedHeight:
487 488 data['spc'] = spc[:,:,self.height_index]
488 489 else:
489 490 data['spc'] = spc
490 491 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
491 492
492 493 return data, meta
493 494
494 495 def plot(self):
495 496 if self.xaxis == "frequency":
496 497 x = self.data.xrange[0][1:]
497 498 self.xlabel = "Frequency (kHz)"
498 499 elif self.xaxis == "time":
499 500 x = self.data.xrange[1]
500 501 self.xlabel = "Time (ms)"
501 502 else:
502 503 x = self.data.xrange[2]
503 504 self.xlabel = "Velocity (m/s)"
504 505
505 506 self.titles = []
506 507
507 508 y = self.data.yrange
508 509 z = self.data[-1]['spc']
509 510 #print(z.shape)
510 511 if self.height_index:
511 512 index = numpy.array(self.height_index)
512 513 else:
513 514 index = numpy.arange(0, len(y), int((len(y))/9))
514 515
515 516 for n, ax in enumerate(self.axes):
516 517 if ax.firsttime:
517 518 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
518 519 self.xmin = self.xmin if self.xmin else -self.xmax
519 520 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
520 521 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
521 522 if self.selectedHeight:
522 523 ax.plt = ax.plot(x, z[n,:])
523 524
524 525 else:
525 526 ax.plt = ax.plot(x, z[n, :, index].T)
526 527 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
527 528 self.figures[0].legend(ax.plt, labels, loc='center right')
528 529 else:
529 530 for i, line in enumerate(ax.plt):
530 531 if self.selectedHeight:
531 532 line.set_data(x, z[n, :])
532 533 else:
533 534 line.set_data(x, z[n, :, index[i]])
534 535 self.titles.append('CH {}'.format(self.channelList[n]))
535 536 plt.suptitle(self.maintitle)
536 537
537 538 class BeaconPhase(Plot):
538 539
539 540 __isConfig = None
540 541 __nsubplots = None
541 542
542 543 PREFIX = 'beacon_phase'
543 544
544 545 def __init__(self):
545 546 Plot.__init__(self)
546 547 self.timerange = 24*60*60
547 548 self.isConfig = False
548 549 self.__nsubplots = 1
549 550 self.counter_imagwr = 0
550 551 self.WIDTH = 800
551 552 self.HEIGHT = 400
552 553 self.WIDTHPROF = 120
553 554 self.HEIGHTPROF = 0
554 555 self.xdata = None
555 556 self.ydata = None
556 557
557 558 self.PLOT_CODE = BEACON_CODE
558 559
559 560 self.FTP_WEI = None
560 561 self.EXP_CODE = None
561 562 self.SUB_EXP_CODE = None
562 563 self.PLOT_POS = None
563 564
564 565 self.filename_phase = None
565 566
566 567 self.figfile = None
567 568
568 569 self.xmin = None
569 570 self.xmax = None
570 571
571 572 def getSubplots(self):
572 573
573 574 ncol = 1
574 575 nrow = 1
575 576
576 577 return nrow, ncol
577 578
578 579 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
579 580
580 581 self.__showprofile = showprofile
581 582 self.nplots = nplots
582 583
583 584 ncolspan = 7
584 585 colspan = 6
585 586 self.__nsubplots = 2
586 587
587 588 self.createFigure(id = id,
588 589 wintitle = wintitle,
589 590 widthplot = self.WIDTH+self.WIDTHPROF,
590 591 heightplot = self.HEIGHT+self.HEIGHTPROF,
591 592 show=show)
592 593
593 594 nrow, ncol = self.getSubplots()
594 595
595 596 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
596 597
597 598 def save_phase(self, filename_phase):
598 599 f = open(filename_phase,'w+')
599 600 f.write('\n\n')
600 601 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
601 602 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
602 603 f.close()
603 604
604 605 def save_data(self, filename_phase, data, data_datetime):
605 606 f=open(filename_phase,'a')
606 607 timetuple_data = data_datetime.timetuple()
607 608 day = str(timetuple_data.tm_mday)
608 609 month = str(timetuple_data.tm_mon)
609 610 year = str(timetuple_data.tm_year)
610 611 hour = str(timetuple_data.tm_hour)
611 612 minute = str(timetuple_data.tm_min)
612 613 second = str(timetuple_data.tm_sec)
613 614 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
614 615 f.close()
615 616
616 617 def plot(self):
617 618 log.warning('TODO: Not yet implemented...')
618 619
619 620 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
620 621 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
621 622 timerange=None,
622 623 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
623 624 server=None, folder=None, username=None, password=None,
624 625 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
625 626
626 627 if dataOut.flagNoData:
627 628 return dataOut
628 629
629 630 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
630 631 return
631 632
632 633 if pairsList == None:
633 634 pairsIndexList = dataOut.pairsIndexList[:10]
634 635 else:
635 636 pairsIndexList = []
636 637 for pair in pairsList:
637 638 if pair not in dataOut.pairsList:
638 639 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
639 640 pairsIndexList.append(dataOut.pairsList.index(pair))
640 641
641 642 if pairsIndexList == []:
642 643 return
643 644
644 645 # if len(pairsIndexList) > 4:
645 646 # pairsIndexList = pairsIndexList[0:4]
646 647
647 648 hmin_index = None
648 649 hmax_index = None
649 650
650 651 if hmin != None and hmax != None:
651 652 indexes = numpy.arange(dataOut.nHeights)
652 653 hmin_list = indexes[dataOut.heightList >= hmin]
653 654 hmax_list = indexes[dataOut.heightList <= hmax]
654 655
655 656 if hmin_list.any():
656 657 hmin_index = hmin_list[0]
657 658
658 659 if hmax_list.any():
659 660 hmax_index = hmax_list[-1]+1
660 661
661 662 x = dataOut.getTimeRange()
662 663
663 664 thisDatetime = dataOut.datatime
664 665
665 666 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
666 667 xlabel = "Local Time"
667 668 ylabel = "Phase (degrees)"
668 669
669 670 update_figfile = False
670 671
671 672 nplots = len(pairsIndexList)
672 673 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
673 674 phase_beacon = numpy.zeros(len(pairsIndexList))
674 675 for i in range(nplots):
675 676 pair = dataOut.pairsList[pairsIndexList[i]]
676 677 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
677 678 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
678 679 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
679 680 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
680 681 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
681 682
682 683 if dataOut.beacon_heiIndexList:
683 684 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
684 685 else:
685 686 phase_beacon[i] = numpy.average(phase)
686 687
687 688 if not self.isConfig:
688 689
689 690 nplots = len(pairsIndexList)
690 691
691 692 self.setup(id=id,
692 693 nplots=nplots,
693 694 wintitle=wintitle,
694 695 showprofile=showprofile,
695 696 show=show)
696 697
697 698 if timerange != None:
698 699 self.timerange = timerange
699 700
700 701 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
701 702
702 703 if ymin == None: ymin = 0
703 704 if ymax == None: ymax = 360
704 705
705 706 self.FTP_WEI = ftp_wei
706 707 self.EXP_CODE = exp_code
707 708 self.SUB_EXP_CODE = sub_exp_code
708 709 self.PLOT_POS = plot_pos
709 710
710 711 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
711 712 self.isConfig = True
712 713 self.figfile = figfile
713 714 self.xdata = numpy.array([])
714 715 self.ydata = numpy.array([])
715 716
716 717 update_figfile = True
717 718
718 719 #open file beacon phase
719 720 path = '%s%03d' %(self.PREFIX, self.id)
720 721 beacon_file = os.path.join(path,'%s.txt'%self.name)
721 722 self.filename_phase = os.path.join(figpath,beacon_file)
722 723 #self.save_phase(self.filename_phase)
723 724
724 725
725 726 #store data beacon phase
726 727 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
727 728
728 729 self.setWinTitle(title)
729 730
730 731
731 732 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
732 733
733 734 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
734 735
735 736 axes = self.axesList[0]
736 737
737 738 self.xdata = numpy.hstack((self.xdata, x[0:1]))
738 739
739 740 if len(self.ydata)==0:
740 741 self.ydata = phase_beacon.reshape(-1,1)
741 742 else:
742 743 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
743 744
744 745
745 746 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
746 747 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
747 748 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
748 749 XAxisAsTime=True, grid='both'
749 750 )
750 751
751 752 self.draw()
752 753
753 754 if dataOut.ltctime >= self.xmax:
754 755 self.counter_imagwr = wr_period
755 756 self.isConfig = False
756 757 update_figfile = True
757 758
758 759 self.save(figpath=figpath,
759 760 figfile=figfile,
760 761 save=save,
761 762 ftp=ftp,
762 763 wr_period=wr_period,
763 764 thisDatetime=thisDatetime,
764 765 update_figfile=update_figfile)
765 766
766 767 return dataOut
767 768
768 769 class NoiselessSpectraPlot(Plot):
769 770 '''
770 771 Plot for Spectra data, subtracting
771 772 the noise in all channels, using for
772 773 amisr-14 data
773 774 '''
774 775
775 776 CODE = 'noiseless_spc'
776 777 colormap = 'nipy_spectral'
777 778 plot_type = 'pcolor'
778 779 buffering = False
779 780 channelList = []
780 781
781 782 def setup(self):
782 783
783 784 self.nplots = len(self.data.channels)
784 785 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
785 786 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
786 787 self.height = 2.6 * self.nrows
787 788
788 789 self.cb_label = 'dB'
789 790 if self.showprofile:
790 791 self.width = 4 * self.ncols
791 792 else:
792 793 self.width = 3.5 * self.ncols
793 794 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
794 795 self.ylabel = 'Range [km]'
795 796
796 797
797 798 def update_list(self,dataOut):
798 799 if len(self.channelList) == 0:
799 800 self.channelList = dataOut.channelList
800 801
801 802 def update(self, dataOut):
802 803
803 804 self.update_list(dataOut)
804 805 data = {}
805 806 meta = {}
806 807 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
807 808 (nch, nff, nh) = dataOut.data_spc.shape
808 809 n1 = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
809 810 noise = numpy.repeat(n1,nff, axis=1).reshape((nch,nff,nh))
810 811 #print(noise.shape, "noise", noise)
811 812
812 813 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - noise
813 814
814 815 data['spc'] = spc
815 816 data['rti'] = dataOut.getPower() - n1
816 817
817 818 data['noise'] = n0
818 819 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
819 820
820 821 return data, meta
821 822
822 823 def plot(self):
823 824 if self.xaxis == "frequency":
824 825 x = self.data.xrange[0]
825 826 self.xlabel = "Frequency (kHz)"
826 827 elif self.xaxis == "time":
827 828 x = self.data.xrange[1]
828 829 self.xlabel = "Time (ms)"
829 830 else:
830 831 x = self.data.xrange[2]
831 832 self.xlabel = "Velocity (m/s)"
832 833
833 834 self.titles = []
834 835 y = self.data.yrange
835 836 self.y = y
836 837
837 838 data = self.data[-1]
838 839 z = data['spc']
839 840
840 841 for n, ax in enumerate(self.axes):
841 842 noise = data['noise'][n]
842 843
843 844 if ax.firsttime:
844 845 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
845 846 self.xmin = self.xmin if self.xmin else -self.xmax
846 847 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
847 848 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
848 849 ax.plt = ax.pcolormesh(x, y, z[n].T,
849 850 vmin=self.zmin,
850 851 vmax=self.zmax,
851 852 cmap=plt.get_cmap(self.colormap)
852 853 )
853 854
854 855 if self.showprofile:
855 856 ax.plt_profile = self.pf_axes[n].plot(
856 857 data['rti'][n], y)[0]
857 858 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
858 859 color="k", linestyle="dashed", lw=1)[0]
859 860
860 861 else:
861 862 ax.plt.set_array(z[n].T.ravel())
862 863 if self.showprofile:
863 864 ax.plt_profile.set_data(data['rti'][n], y)
864 865 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
865 866
866 867 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
867 868
868 869
869 870 class NoiselessRTIPlot(Plot):
870 871 '''
871 872 Plot for RTI data
872 873 '''
873 874
874 875 CODE = 'noiseless_rti'
875 876 colormap = 'jet'
876 877 plot_type = 'pcolorbuffer'
877 878 titles = None
878 879 channelList = []
879 880
880 881 def setup(self):
881 882 self.xaxis = 'time'
882 883 self.ncols = 1
883 884 #print("dataChannels ",self.data.channels)
884 885 self.nrows = len(self.data.channels)
885 886 self.nplots = len(self.data.channels)
886 887 self.ylabel = 'Range [km]'
887 888 self.xlabel = 'Time'
888 889 self.cb_label = 'dB'
889 890 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
890 891 self.titles = ['{} Channel {}'.format(
891 892 self.CODE.upper(), x) for x in range(self.nplots)]
892 893
893 894 def update_list(self,dataOut):
894 895
895 896 self.channelList = dataOut.channelList
896 897
897 898
898 899 def update(self, dataOut):
899 900 if len(self.channelList) == 0:
900 901 self.update_list(dataOut)
901 902 data = {}
902 903 meta = {}
903 904
904 905 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
905 906 (nch, nff, nh) = dataOut.data_spc.shape
906 907 noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
907 908
908 909
909 910 data['noiseless_rti'] = dataOut.getPower() - noise
910 911 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
911 912 return data, meta
912 913
913 914 def plot(self):
914 915
915 916 self.x = self.data.times
916 917 self.y = self.data.yrange
917 918 self.z = self.data['noiseless_rti']
918 919 self.z = numpy.array(self.z, dtype=float)
919 920 self.z = numpy.ma.masked_invalid(self.z)
920 921
921 922 try:
922 923 if self.channelList != None:
923 924 self.titles = ['{} Channel {}'.format(
924 925 self.CODE.upper(), x) for x in self.channelList]
925 926 except:
926 927 if self.channelList.any() != None:
927 928 self.titles = ['{} Channel {}'.format(
928 929 self.CODE.upper(), x) for x in self.channelList]
929 930 if self.decimation is None:
930 931 x, y, z = self.fill_gaps(self.x, self.y, self.z)
931 932 else:
932 933 x, y, z = self.fill_gaps(*self.decimate())
933 934 dummy_var = self.axes #ExtraΓ±amente esto actualiza el valor axes
934 935 for n, ax in enumerate(self.axes):
935 936 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
936 937 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
937 938 data = self.data[-1]
938 939 if ax.firsttime:
939 940 ax.plt = ax.pcolormesh(x, y, z[n].T,
940 941 vmin=self.zmin,
941 942 vmax=self.zmax,
942 943 cmap=plt.get_cmap(self.colormap)
943 944 )
944 945 if self.showprofile:
945 946 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
946 947
947 948 if "noise" in self.data:
948 949 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
949 950 color="k", linestyle="dashed", lw=1)[0]
950 951 else:
951 952 ax.collections.remove(ax.collections[0])
952 953 ax.plt = ax.pcolormesh(x, y, z[n].T,
953 954 vmin=self.zmin,
954 955 vmax=self.zmax,
955 956 cmap=plt.get_cmap(self.colormap)
956 957 )
957 958 if self.showprofile:
958 959 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
959 960 if "noise" in self.data:
960 961 ax.plot_noise.set_data(numpy.repeat(
961 962 data['noise'][n], len(self.y)), self.y)
Show file before | Show file after | Hide comments
@@ -1,1688 +1,1814
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15 import math
16 16
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 18 from schainpy.model.data.jrodata import Spectra
19 19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 20 from schainpy.utils import log
21 21
22 22 from scipy.optimize import curve_fit
23 23
24 24 class SpectraProc(ProcessingUnit):
25 25
26 26 def __init__(self):
27 27
28 28 ProcessingUnit.__init__(self)
29 29
30 30 self.buffer = None
31 31 self.firstdatatime = None
32 32 self.profIndex = 0
33 33 self.dataOut = Spectra()
34 34 self.id_min = None
35 35 self.id_max = None
36 36 self.setupReq = False #Agregar a todas las unidades de proc
37 37
38 38 def __updateSpecFromVoltage(self):
39 39
40 40 self.dataOut.timeZone = self.dataIn.timeZone
41 41 self.dataOut.dstFlag = self.dataIn.dstFlag
42 42 self.dataOut.errorCount = self.dataIn.errorCount
43 43 self.dataOut.useLocalTime = self.dataIn.useLocalTime
44 44 try:
45 45 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
46 46 except:
47 47 pass
48 48 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
49 49 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
50 50 self.dataOut.channelList = self.dataIn.channelList
51 51 self.dataOut.heightList = self.dataIn.heightList
52 52 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
53 53 self.dataOut.nProfiles = self.dataOut.nFFTPoints
54 54 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
55 55 self.dataOut.utctime = self.firstdatatime
56 56 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
57 57 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
58 58 self.dataOut.flagShiftFFT = False
59 59 self.dataOut.nCohInt = self.dataIn.nCohInt
60 60 self.dataOut.nIncohInt = 1
61 61 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
62 62 self.dataOut.frequency = self.dataIn.frequency
63 63 self.dataOut.realtime = self.dataIn.realtime
64 64 self.dataOut.azimuth = self.dataIn.azimuth
65 65 self.dataOut.zenith = self.dataIn.zenith
66 66 self.dataOut.codeList = self.dataIn.codeList
67 67 self.dataOut.azimuthList = self.dataIn.azimuthList
68 68 self.dataOut.elevationList = self.dataIn.elevationList
69 69
70 70
71 71 def __getFft(self):
72 72 """
73 73 Convierte valores de Voltaje a Spectra
74 74
75 75 Affected:
76 76 self.dataOut.data_spc
77 77 self.dataOut.data_cspc
78 78 self.dataOut.data_dc
79 79 self.dataOut.heightList
80 80 self.profIndex
81 81 self.buffer
82 82 self.dataOut.flagNoData
83 83 """
84 84 fft_volt = numpy.fft.fft(
85 85 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
86 86 fft_volt = fft_volt.astype(numpy.dtype('complex'))
87 87 dc = fft_volt[:, 0, :]
88 88
89 89 # calculo de self-spectra
90 90 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
91 91 spc = fft_volt * numpy.conjugate(fft_volt)
92 92 spc = spc.real
93 93
94 94 blocksize = 0
95 95 blocksize += dc.size
96 96 blocksize += spc.size
97 97
98 98 cspc = None
99 99 pairIndex = 0
100 100 if self.dataOut.pairsList != None:
101 101 # calculo de cross-spectra
102 102 cspc = numpy.zeros(
103 103 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
104 104 for pair in self.dataOut.pairsList:
105 105 if pair[0] not in self.dataOut.channelList:
106 106 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
107 107 str(pair), str(self.dataOut.channelList)))
108 108 if pair[1] not in self.dataOut.channelList:
109 109 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
110 110 str(pair), str(self.dataOut.channelList)))
111 111
112 112 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
113 113 numpy.conjugate(fft_volt[pair[1], :, :])
114 114 pairIndex += 1
115 115 blocksize += cspc.size
116 116
117 117 self.dataOut.data_spc = spc
118 118 self.dataOut.data_cspc = cspc
119 119 self.dataOut.data_dc = dc
120 120 self.dataOut.blockSize = blocksize
121 121 self.dataOut.flagShiftFFT = False
122 122
123 123 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
124 124
125 125 if self.dataIn.type == "Spectra":
126 126
127 127 try:
128 128 self.dataOut.copy(self.dataIn)
129 129
130 130 except Exception as e:
131 131 print(e)
132 132
133 133 if shift_fft:
134 134 #desplaza a la derecha en el eje 2 determinadas posiciones
135 135 shift = int(self.dataOut.nFFTPoints/2)
136 136 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
137 137
138 138 if self.dataOut.data_cspc is not None:
139 139 #desplaza a la derecha en el eje 2 determinadas posiciones
140 140 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
141 141 if pairsList:
142 142 self.__selectPairs(pairsList)
143 143
144 144
145 145 elif self.dataIn.type == "Voltage":
146 146
147 147 self.dataOut.flagNoData = True
148 148
149 149 if nFFTPoints == None:
150 150 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
151 151
152 152 if nProfiles == None:
153 153 nProfiles = nFFTPoints
154 154
155 155 if ippFactor == None:
156 156 self.dataOut.ippFactor = 1
157 157
158 158 self.dataOut.nFFTPoints = nFFTPoints
159 159
160 160 if self.buffer is None:
161 161 self.buffer = numpy.zeros((self.dataIn.nChannels,
162 162 nProfiles,
163 163 self.dataIn.nHeights),
164 164 dtype='complex')
165 165
166 166 if self.dataIn.flagDataAsBlock:
167 167 nVoltProfiles = self.dataIn.data.shape[1]
168 168
169 169 if nVoltProfiles == nProfiles:
170 170 self.buffer = self.dataIn.data.copy()
171 171 self.profIndex = nVoltProfiles
172 172
173 173 elif nVoltProfiles < nProfiles:
174 174
175 175 if self.profIndex == 0:
176 176 self.id_min = 0
177 177 self.id_max = nVoltProfiles
178 178
179 179 self.buffer[:, self.id_min:self.id_max,
180 180 :] = self.dataIn.data
181 181 self.profIndex += nVoltProfiles
182 182 self.id_min += nVoltProfiles
183 183 self.id_max += nVoltProfiles
184 184 else:
185 185 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
186 186 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
187 187 self.dataOut.flagNoData = True
188 188 else:
189 189 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
190 190 self.profIndex += 1
191 191
192 192 if self.firstdatatime == None:
193 193 self.firstdatatime = self.dataIn.utctime
194 194
195 195 if self.profIndex == nProfiles:
196 196 self.__updateSpecFromVoltage()
197 197 if pairsList == None:
198 198 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
199 199 else:
200 200 self.dataOut.pairsList = pairsList
201 201 self.__getFft()
202 202 self.dataOut.flagNoData = False
203 203 self.firstdatatime = None
204 204 self.profIndex = 0
205 self.dataOut.noise_estimation = None
205 206 else:
206 207 raise ValueError("The type of input object '%s' is not valid".format(
207 208 self.dataIn.type))
208 209
209 210 def __selectPairs(self, pairsList):
210 211
211 212 if not pairsList:
212 213 return
213 214
214 215 pairs = []
215 216 pairsIndex = []
216 217
217 218 for pair in pairsList:
218 219 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
219 220 continue
220 221 pairs.append(pair)
221 222 pairsIndex.append(pairs.index(pair))
222 223
223 224 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
224 225 self.dataOut.pairsList = pairs
225 226
226 227 return
227 228
228 229 def selectFFTs(self, minFFT, maxFFT ):
229 230 """
230 231 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
231 232 minFFT<= FFT <= maxFFT
232 233 """
233 234
234 235 if (minFFT > maxFFT):
235 236 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
236 237
237 238 if (minFFT < self.dataOut.getFreqRange()[0]):
238 239 minFFT = self.dataOut.getFreqRange()[0]
239 240
240 241 if (maxFFT > self.dataOut.getFreqRange()[-1]):
241 242 maxFFT = self.dataOut.getFreqRange()[-1]
242 243
243 244 minIndex = 0
244 245 maxIndex = 0
245 246 FFTs = self.dataOut.getFreqRange()
246 247
247 248 inda = numpy.where(FFTs >= minFFT)
248 249 indb = numpy.where(FFTs <= maxFFT)
249 250
250 251 try:
251 252 minIndex = inda[0][0]
252 253 except:
253 254 minIndex = 0
254 255
255 256 try:
256 257 maxIndex = indb[0][-1]
257 258 except:
258 259 maxIndex = len(FFTs)
259 260
260 261 self.selectFFTsByIndex(minIndex, maxIndex)
261 262
262 263 return 1
263 264
264 265 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
265 266 newheis = numpy.where(
266 267 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
267 268
268 269 if hei_ref != None:
269 270 newheis = numpy.where(self.dataOut.heightList > hei_ref)
270 271
271 272 minIndex = min(newheis[0])
272 273 maxIndex = max(newheis[0])
273 274 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
274 275 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
275 276
276 277 # determina indices
277 278 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
278 279 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
279 280 avg_dB = 10 * \
280 281 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
281 282 beacon_dB = numpy.sort(avg_dB)[-nheis:]
282 283 beacon_heiIndexList = []
283 284 for val in avg_dB.tolist():
284 285 if val >= beacon_dB[0]:
285 286 beacon_heiIndexList.append(avg_dB.tolist().index(val))
286 287
287 288 #data_spc = data_spc[:,:,beacon_heiIndexList]
288 289 data_cspc = None
289 290 if self.dataOut.data_cspc is not None:
290 291 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
291 292 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
292 293
293 294 data_dc = None
294 295 if self.dataOut.data_dc is not None:
295 296 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
296 297 #data_dc = data_dc[:,beacon_heiIndexList]
297 298
298 299 self.dataOut.data_spc = data_spc
299 300 self.dataOut.data_cspc = data_cspc
300 301 self.dataOut.data_dc = data_dc
301 302 self.dataOut.heightList = heightList
302 303 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
303 304
304 305 return 1
305 306
306 307 def selectFFTsByIndex(self, minIndex, maxIndex):
307 308 """
308 309
309 310 """
310 311
311 312 if (minIndex < 0) or (minIndex > maxIndex):
312 313 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
313 314
314 315 if (maxIndex >= self.dataOut.nProfiles):
315 316 maxIndex = self.dataOut.nProfiles-1
316 317
317 318 #Spectra
318 319 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
319 320
320 321 data_cspc = None
321 322 if self.dataOut.data_cspc is not None:
322 323 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
323 324
324 325 data_dc = None
325 326 if self.dataOut.data_dc is not None:
326 327 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
327 328
328 329 self.dataOut.data_spc = data_spc
329 330 self.dataOut.data_cspc = data_cspc
330 331 self.dataOut.data_dc = data_dc
331 332
332 333 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
333 334 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
334 335 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
335 336
336 337 return 1
337 338
338 339 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
339 340 # validacion de rango
340 341 if minHei == None:
341 342 minHei = self.dataOut.heightList[0]
342 343
343 344 if maxHei == None:
344 345 maxHei = self.dataOut.heightList[-1]
345 346
346 347 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
347 348 print('minHei: %.2f is out of the heights range' % (minHei))
348 349 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
349 350 minHei = self.dataOut.heightList[0]
350 351
351 352 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
352 353 print('maxHei: %.2f is out of the heights range' % (maxHei))
353 354 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
354 355 maxHei = self.dataOut.heightList[-1]
355 356
356 357 # validacion de velocidades
357 358 velrange = self.dataOut.getVelRange(1)
358 359
359 360 if minVel == None:
360 361 minVel = velrange[0]
361 362
362 363 if maxVel == None:
363 364 maxVel = velrange[-1]
364 365
365 366 if (minVel < velrange[0]) or (minVel > maxVel):
366 367 print('minVel: %.2f is out of the velocity range' % (minVel))
367 368 print('minVel is setting to %.2f' % (velrange[0]))
368 369 minVel = velrange[0]
369 370
370 371 if (maxVel > velrange[-1]) or (maxVel < minVel):
371 372 print('maxVel: %.2f is out of the velocity range' % (maxVel))
372 373 print('maxVel is setting to %.2f' % (velrange[-1]))
373 374 maxVel = velrange[-1]
374 375
375 376 # seleccion de indices para rango
376 377 minIndex = 0
377 378 maxIndex = 0
378 379 heights = self.dataOut.heightList
379 380
380 381 inda = numpy.where(heights >= minHei)
381 382 indb = numpy.where(heights <= maxHei)
382 383
383 384 try:
384 385 minIndex = inda[0][0]
385 386 except:
386 387 minIndex = 0
387 388
388 389 try:
389 390 maxIndex = indb[0][-1]
390 391 except:
391 392 maxIndex = len(heights)
392 393
393 394 if (minIndex < 0) or (minIndex > maxIndex):
394 395 raise ValueError("some value in (%d,%d) is not valid" % (
395 396 minIndex, maxIndex))
396 397
397 398 if (maxIndex >= self.dataOut.nHeights):
398 399 maxIndex = self.dataOut.nHeights - 1
399 400
400 401 # seleccion de indices para velocidades
401 402 indminvel = numpy.where(velrange >= minVel)
402 403 indmaxvel = numpy.where(velrange <= maxVel)
403 404 try:
404 405 minIndexVel = indminvel[0][0]
405 406 except:
406 407 minIndexVel = 0
407 408
408 409 try:
409 410 maxIndexVel = indmaxvel[0][-1]
410 411 except:
411 412 maxIndexVel = len(velrange)
412 413
413 414 # seleccion del espectro
414 415 data_spc = self.dataOut.data_spc[:,
415 416 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
416 417 # estimacion de ruido
417 418 noise = numpy.zeros(self.dataOut.nChannels)
418 419
419 420 for channel in range(self.dataOut.nChannels):
420 421 daux = data_spc[channel, :, :]
421 422 sortdata = numpy.sort(daux, axis=None)
422 423 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
423 424
424 425 self.dataOut.noise_estimation = noise.copy()
425 426
426 427 return 1
427 428
428 429 class removeDC(Operation):
429 430
430 431 def run(self, dataOut, mode=2):
431 432 self.dataOut = dataOut
432 433 jspectra = self.dataOut.data_spc
433 434 jcspectra = self.dataOut.data_cspc
434 435
435 436 num_chan = jspectra.shape[0]
436 437 num_hei = jspectra.shape[2]
437 438
438 439 if jcspectra is not None:
439 440 jcspectraExist = True
440 441 num_pairs = jcspectra.shape[0]
441 442 else:
442 443 jcspectraExist = False
443 444
444 445 freq_dc = int(jspectra.shape[1] / 2)
445 446 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
446 447 ind_vel = ind_vel.astype(int)
447 448
448 449 if ind_vel[0] < 0:
449 450 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
450 451
451 452 if mode == 1:
452 453 jspectra[:, freq_dc, :] = (
453 454 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
454 455
455 456 if jcspectraExist:
456 457 jcspectra[:, freq_dc, :] = (
457 458 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
458 459
459 460 if mode == 2:
460 461
461 462 vel = numpy.array([-2, -1, 1, 2])
462 463 xx = numpy.zeros([4, 4])
463 464
464 465 for fil in range(4):
465 466 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
466 467
467 468 xx_inv = numpy.linalg.inv(xx)
468 469 xx_aux = xx_inv[0, :]
469 470
470 471 for ich in range(num_chan):
471 472 yy = jspectra[ich, ind_vel, :]
472 473 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
473 474
474 475 junkid = jspectra[ich, freq_dc, :] <= 0
475 476 cjunkid = sum(junkid)
476 477
477 478 if cjunkid.any():
478 479 jspectra[ich, freq_dc, junkid.nonzero()] = (
479 480 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
480 481
481 482 if jcspectraExist:
482 483 for ip in range(num_pairs):
483 484 yy = jcspectra[ip, ind_vel, :]
484 485 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
485 486
486 487 self.dataOut.data_spc = jspectra
487 488 self.dataOut.data_cspc = jcspectra
488 489
489 490 return self.dataOut
490 491
492 class getNoise(Operation):
493 def __init__(self):
494
495 Operation.__init__(self)
496
497 def run(self, dataOut, minHei=None, maxHei=None, minVel=None, maxVel=None, minFreq= None, maxFreq=None,):
498 self.dataOut = dataOut.copy()
499 print("1: ",dataOut.noise_estimation, dataOut.normFactor)
500
501 if minHei == None:
502 minHei = self.dataOut.heightList[0]
503
504 if maxHei == None:
505 maxHei = self.dataOut.heightList[-1]
506
507 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
508 print('minHei: %.2f is out of the heights range' % (minHei))
509 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
510 minHei = self.dataOut.heightList[0]
511
512 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
513 print('maxHei: %.2f is out of the heights range' % (maxHei))
514 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
515 maxHei = self.dataOut.heightList[-1]
516
517
518 #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
519 minIndexFFT = 0
520 maxIndexFFT = 0
521 # validacion de velocidades
522 indminPoint = None
523 indmaxPoint = None
524
525 if minVel == None and maxVel == None:
526
527 freqrange = self.dataOut.getFreqRange(1)
528
529 if minFreq == None:
530 minFreq = freqrange[0]
531
532 if maxFreq == None:
533 maxFreq = freqrange[-1]
534
535 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
536 print('minFreq: %.2f is out of the frequency range' % (minFreq))
537 print('minFreq is setting to %.2f' % (freqrange[0]))
538 minFreq = freqrange[0]
539
540 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
541 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
542 print('maxFreq is setting to %.2f' % (freqrange[-1]))
543 maxFreq = freqrange[-1]
544
545 indminPoint = numpy.where(freqrange >= minFreq)
546 indmaxPoint = numpy.where(freqrange <= maxFreq)
547
548 else:
549 velrange = self.dataOut.getVelRange(1)
550
551 if minVel == None:
552 minVel = velrange[0]
553
554 if maxVel == None:
555 maxVel = velrange[-1]
556
557 if (minVel < velrange[0]) or (minVel > maxVel):
558 print('minVel: %.2f is out of the velocity range' % (minVel))
559 print('minVel is setting to %.2f' % (velrange[0]))
560 minVel = velrange[0]
561
562 if (maxVel > velrange[-1]) or (maxVel < minVel):
563 print('maxVel: %.2f is out of the velocity range' % (maxVel))
564 print('maxVel is setting to %.2f' % (velrange[-1]))
565 maxVel = velrange[-1]
566
567 indminPoint = numpy.where(velrange >= minVel)
568 indmaxPoint = numpy.where(velrange <= maxVel)
569
570
571 # seleccion de indices para rango
572 minIndex = 0
573 maxIndex = 0
574 heights = self.dataOut.heightList
575
576 inda = numpy.where(heights >= minHei)
577 indb = numpy.where(heights <= maxHei)
578
579 try:
580 minIndex = inda[0][0]
581 except:
582 minIndex = 0
583
584 try:
585 maxIndex = indb[0][-1]
586 except:
587 maxIndex = len(heights)
588
589 if (minIndex < 0) or (minIndex > maxIndex):
590 raise ValueError("some value in (%d,%d) is not valid" % (
591 minIndex, maxIndex))
592
593 if (maxIndex >= self.dataOut.nHeights):
594 maxIndex = self.dataOut.nHeights - 1
595 #############################################################3
596 # seleccion de indices para velocidades
597
598 try:
599 minIndexFFT = indminPoint[0][0]
600 except:
601 minIndexFFT = 0
602
603 try:
604 maxIndexFFT = indmaxPoint[0][-1]
605 except:
606 maxIndexFFT = len( self.dataOut.getFreqRange(1))
607
608 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
609 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
610
611 self.dataOut.noise_estimation = noise.copy()
612 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
613 return self.dataOut
614
615
616
491 617 # import matplotlib.pyplot as plt
492 618
493 619 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
494 620 z = (x - a1) / a2
495 621 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
496 622 return y
497 623
498 624
499 625 class CleanRayleigh(Operation):
500 626
501 627 def __init__(self):
502 628
503 629 Operation.__init__(self)
504 630 self.i=0
505 631 self.isConfig = False
506 632 self.__dataReady = False
507 633 self.__profIndex = 0
508 634 self.byTime = False
509 635 self.byProfiles = False
510 636
511 637 self.bloques = None
512 638 self.bloque0 = None
513 639
514 640 self.index = 0
515 641
516 642 self.buffer = 0
517 643 self.buffer2 = 0
518 644 self.buffer3 = 0
519 645
520 646
521 647 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
522 648
523 649 self.nChannels = dataOut.nChannels
524 650 self.nProf = dataOut.nProfiles
525 651 self.nPairs = dataOut.data_cspc.shape[0]
526 652 self.pairsArray = numpy.array(dataOut.pairsList)
527 653 self.spectra = dataOut.data_spc
528 654 self.cspectra = dataOut.data_cspc
529 655 self.heights = dataOut.heightList #alturas totales
530 656 self.nHeights = len(self.heights)
531 657 self.min_hei = min_hei
532 658 self.max_hei = max_hei
533 659 if (self.min_hei == None):
534 660 self.min_hei = 0
535 661 if (self.max_hei == None):
536 662 self.max_hei = dataOut.heightList[-1]
537 663 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
538 664 self.heightsClean = self.heights[self.hval] #alturas filtradas
539 665 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
540 666 self.nHeightsClean = len(self.heightsClean)
541 667 self.channels = dataOut.channelList
542 668 self.nChan = len(self.channels)
543 669 self.nIncohInt = dataOut.nIncohInt
544 670 self.__initime = dataOut.utctime
545 671 self.maxAltInd = self.hval[-1]+1
546 672 self.minAltInd = self.hval[0]
547 673
548 674 self.crosspairs = dataOut.pairsList
549 675 self.nPairs = len(self.crosspairs)
550 676 self.normFactor = dataOut.normFactor
551 677 self.nFFTPoints = dataOut.nFFTPoints
552 678 self.ippSeconds = dataOut.ippSeconds
553 679 self.currentTime = self.__initime
554 680 self.pairsArray = numpy.array(dataOut.pairsList)
555 681 self.factor_stdv = factor_stdv
556 682
557 683 if n != None :
558 684 self.byProfiles = True
559 685 self.nIntProfiles = n
560 686 else:
561 687 self.__integrationtime = timeInterval
562 688
563 689 self.__dataReady = False
564 690 self.isConfig = True
565 691
566 692
567 693
568 694 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
569 695
570 696 if not self.isConfig :
571 697
572 698 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
573 699
574 700 tini=dataOut.utctime
575 701
576 702 if self.byProfiles:
577 703 if self.__profIndex == self.nIntProfiles:
578 704 self.__dataReady = True
579 705 else:
580 706 if (tini - self.__initime) >= self.__integrationtime:
581 707
582 708 self.__dataReady = True
583 709 self.__initime = tini
584 710
585 711 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
586 712
587 713 if self.__dataReady:
588 714
589 715 self.__profIndex = 0
590 716 jspc = self.buffer
591 717 jcspc = self.buffer2
592 718 #jnoise = self.buffer3
593 719 self.buffer = dataOut.data_spc
594 720 self.buffer2 = dataOut.data_cspc
595 721 #self.buffer3 = dataOut.noise
596 722 self.currentTime = dataOut.utctime
597 723 if numpy.any(jspc) :
598 724 #print( jspc.shape, jcspc.shape)
599 725 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
600 726 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
601 727 self.__dataReady = False
602 728 #print( jspc.shape, jcspc.shape)
603 729 dataOut.flagNoData = False
604 730 else:
605 731 dataOut.flagNoData = True
606 732 self.__dataReady = False
607 733 return dataOut
608 734 else:
609 735 #print( len(self.buffer))
610 736 if numpy.any(self.buffer):
611 737 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
612 738 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
613 739 self.buffer3 += dataOut.data_dc
614 740 else:
615 741 self.buffer = dataOut.data_spc
616 742 self.buffer2 = dataOut.data_cspc
617 743 self.buffer3 = dataOut.data_dc
618 744 #print self.index, self.fint
619 745 #print self.buffer2.shape
620 746 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
621 747 self.__profIndex += 1
622 748 return dataOut ## NOTE: REV
623 749
624 750
625 751 #index = tini.tm_hour*12+tini.tm_min/5
626 752 '''REVISAR'''
627 753 # jspc = jspc/self.nFFTPoints/self.normFactor
628 754 # jcspc = jcspc/self.nFFTPoints/self.normFactor
629 755
630 756
631 757
632 758 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
633 759 dataOut.data_spc = tmp_spectra
634 760 dataOut.data_cspc = tmp_cspectra
635 761
636 762 #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
637 763
638 764 dataOut.data_dc = self.buffer3
639 765 dataOut.nIncohInt *= self.nIntProfiles
640 766 dataOut.utctime = self.currentTime #tiempo promediado
641 767 #print("Time: ",time.localtime(dataOut.utctime))
642 768 # dataOut.data_spc = sat_spectra
643 769 # dataOut.data_cspc = sat_cspectra
644 770 self.buffer = 0
645 771 self.buffer2 = 0
646 772 self.buffer3 = 0
647 773
648 774 return dataOut
649 775
650 776 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
651 777 #print("OP cleanRayleigh")
652 778 #import matplotlib.pyplot as plt
653 779 #for k in range(149):
654 780 #channelsProcssd = []
655 781 #channelA_ok = False
656 782 #rfunc = cspectra.copy() #self.bloques
657 783 rfunc = spectra.copy()
658 784 #rfunc = cspectra
659 785 #val_spc = spectra*0.0 #self.bloque0*0.0
660 786 #val_cspc = cspectra*0.0 #self.bloques*0.0
661 787 #in_sat_spectra = spectra.copy() #self.bloque0
662 788 #in_sat_cspectra = cspectra.copy() #self.bloques
663 789
664 790
665 791 ###ONLY FOR TEST:
666 792 raxs = math.ceil(math.sqrt(self.nPairs))
667 793 caxs = math.ceil(self.nPairs/raxs)
668 794 if self.nPairs <4:
669 795 raxs = 2
670 796 caxs = 2
671 797 #print(raxs, caxs)
672 798 fft_rev = 14 #nFFT to plot
673 799 hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
674 800 hei_rev = hei_rev[0]
675 801 #print(hei_rev)
676 802
677 803 #print numpy.absolute(rfunc[:,0,0,14])
678 804
679 805 gauss_fit, covariance = None, None
680 806 for ih in range(self.minAltInd,self.maxAltInd):
681 807 for ifreq in range(self.nFFTPoints):
682 808 '''
683 809 ###ONLY FOR TEST:
684 810 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
685 811 fig, axs = plt.subplots(raxs, caxs)
686 812 fig2, axs2 = plt.subplots(raxs, caxs)
687 813 col_ax = 0
688 814 row_ax = 0
689 815 '''
690 816 #print(self.nPairs)
691 817 for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
692 818 # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
693 819 # continue
694 820 # if not self.crosspairs[ii][0] in channelsProcssd:
695 821 # channelA_ok = True
696 822 #print("pair: ",self.crosspairs[ii])
697 823 '''
698 824 ###ONLY FOR TEST:
699 825 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
700 826 col_ax = 0
701 827 row_ax += 1
702 828 '''
703 829 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
704 830 #print(func2clean.shape)
705 831 val = (numpy.isfinite(func2clean)==True).nonzero()
706 832
707 833 if len(val)>0: #limitador
708 834 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
709 835 if min_val <= -40 :
710 836 min_val = -40
711 837 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
712 838 if max_val >= 200 :
713 839 max_val = 200
714 840 #print min_val, max_val
715 841 step = 1
716 842 #print("Getting bins and the histogram")
717 843 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
718 844 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
719 845 #print(len(y_dist),len(binstep[:-1]))
720 846 #print(row_ax,col_ax, " ..")
721 847 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
722 848 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
723 849 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
724 850 parg = [numpy.amax(y_dist),mean,sigma]
725 851
726 852 newY = None
727 853
728 854 try :
729 855 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
730 856 mode = gauss_fit[1]
731 857 stdv = gauss_fit[2]
732 858 #print(" FIT OK",gauss_fit)
733 859 '''
734 860 ###ONLY FOR TEST:
735 861 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
736 862 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
737 863 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
738 864 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
739 865 axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
740 866 '''
741 867 except:
742 868 mode = mean
743 869 stdv = sigma
744 870 #print("FIT FAIL")
745 871 #continue
746 872
747 873
748 874 #print(mode,stdv)
749 875 #Removing echoes greater than mode + std_factor*stdv
750 876 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
751 877 #noval tiene los indices que se van a remover
752 878 #print("Chan ",ii," novals: ",len(noval[0]))
753 879 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
754 880 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
755 881 #print(novall)
756 882 #print(" ",self.pairsArray[ii])
757 883 #cross_pairs = self.pairsArray[ii]
758 884 #Getting coherent echoes which are removed.
759 885 # if len(novall[0]) > 0:
760 886 #
761 887 # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
762 888 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
763 889 # val_cspc[novall[0],ii,ifreq,ih] = 1
764 890 #print("OUT NOVALL 1")
765 891 try:
766 892 pair = (self.channels[ii],self.channels[ii + 1])
767 893 except:
768 894 pair = (99,99)
769 895 #print("par ", pair)
770 896 if ( pair in self.crosspairs):
771 897 q = self.crosspairs.index(pair)
772 898 #print("estΓ‘ aqui: ", q, (ii,ii + 1))
773 899 new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
774 900 cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
775 901
776 902 #if channelA_ok:
777 903 #chA = self.channels.index(cross_pairs[0])
778 904 new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
779 905 spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
780 906 #channelA_ok = False
781 907
782 908 # chB = self.channels.index(cross_pairs[1])
783 909 # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
784 910 # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
785 911 #
786 912 # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
787 913 # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
788 914 '''
789 915 ###ONLY FOR TEST:
790 916 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
791 917 func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
792 918 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
793 919 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
794 920 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
795 921 axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
796 922 '''
797 923 '''
798 924 ###ONLY FOR TEST:
799 925 col_ax += 1 #contador de ploteo columnas
800 926 ##print(col_ax)
801 927 ###ONLY FOR TEST:
802 928 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
803 929 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
804 930 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
805 931 fig.suptitle(title)
806 932 fig2.suptitle(title2)
807 933 plt.show()
808 934 '''
809 935 ##################################################################################################
810 936
811 937 #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
812 938 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
813 939 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
814 940 for ih in range(self.nHeights):
815 941 for ifreq in range(self.nFFTPoints):
816 942 for ich in range(self.nChan):
817 943 tmp = spectra[:,ich,ifreq,ih]
818 944 valid = (numpy.isfinite(tmp[:])==True).nonzero()
819 945
820 946 if len(valid[0]) >0 :
821 947 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
822 948
823 949 for icr in range(self.nPairs):
824 950 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
825 951 valid = (numpy.isfinite(tmp)==True).nonzero()
826 952 if len(valid[0]) > 0:
827 953 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
828 954
829 955 return out_spectra, out_cspectra
830 956
831 957 def REM_ISOLATED_POINTS(self,array,rth):
832 958 # import matplotlib.pyplot as plt
833 959 if rth == None :
834 960 rth = 4
835 961 #print("REM ISO")
836 962 num_prof = len(array[0,:,0])
837 963 num_hei = len(array[0,0,:])
838 964 n2d = len(array[:,0,0])
839 965
840 966 for ii in range(n2d) :
841 967 #print ii,n2d
842 968 tmp = array[ii,:,:]
843 969 #print tmp.shape, array[ii,101,:],array[ii,102,:]
844 970
845 971 # fig = plt.figure(figsize=(6,5))
846 972 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
847 973 # ax = fig.add_axes([left, bottom, width, height])
848 974 # x = range(num_prof)
849 975 # y = range(num_hei)
850 976 # cp = ax.contour(y,x,tmp)
851 977 # ax.clabel(cp, inline=True,fontsize=10)
852 978 # plt.show()
853 979
854 980 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
855 981 tmp = numpy.reshape(tmp,num_prof*num_hei)
856 982 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
857 983 indxs2 = (tmp > 0).nonzero()
858 984
859 985 indxs1 = (indxs1[0])
860 986 indxs2 = indxs2[0]
861 987 #indxs1 = numpy.array(indxs1[0])
862 988 #indxs2 = numpy.array(indxs2[0])
863 989 indxs = None
864 990 #print indxs1 , indxs2
865 991 for iv in range(len(indxs2)):
866 992 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
867 993 #print len(indxs2), indv
868 994 if len(indv[0]) > 0 :
869 995 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
870 996 # print indxs
871 997 indxs = indxs[1:]
872 998 #print(indxs, len(indxs))
873 999 if len(indxs) < 4 :
874 1000 array[ii,:,:] = 0.
875 1001 return
876 1002
877 1003 xpos = numpy.mod(indxs ,num_hei)
878 1004 ypos = (indxs / num_hei)
879 1005 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
880 1006 #print sx
881 1007 xpos = xpos[sx]
882 1008 ypos = ypos[sx]
883 1009
884 1010 # *********************************** Cleaning isolated points **********************************
885 1011 ic = 0
886 1012 while True :
887 1013 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
888 1014 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
889 1015 #plt.plot(r)
890 1016 #plt.show()
891 1017 no_coh1 = (numpy.isfinite(r)==True).nonzero()
892 1018 no_coh2 = (r <= rth).nonzero()
893 1019 #print r, no_coh1, no_coh2
894 1020 no_coh1 = numpy.array(no_coh1[0])
895 1021 no_coh2 = numpy.array(no_coh2[0])
896 1022 no_coh = None
897 1023 #print valid1 , valid2
898 1024 for iv in range(len(no_coh2)):
899 1025 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
900 1026 if len(indv[0]) > 0 :
901 1027 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
902 1028 no_coh = no_coh[1:]
903 1029 #print len(no_coh), no_coh
904 1030 if len(no_coh) < 4 :
905 1031 #print xpos[ic], ypos[ic], ic
906 1032 # plt.plot(r)
907 1033 # plt.show()
908 1034 xpos[ic] = numpy.nan
909 1035 ypos[ic] = numpy.nan
910 1036
911 1037 ic = ic + 1
912 1038 if (ic == len(indxs)) :
913 1039 break
914 1040 #print( xpos, ypos)
915 1041
916 1042 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
917 1043 #print indxs[0]
918 1044 if len(indxs[0]) < 4 :
919 1045 array[ii,:,:] = 0.
920 1046 return
921 1047
922 1048 xpos = xpos[indxs[0]]
923 1049 ypos = ypos[indxs[0]]
924 1050 for i in range(0,len(ypos)):
925 1051 ypos[i]=int(ypos[i])
926 1052 junk = tmp
927 1053 tmp = junk*0.0
928 1054
929 1055 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
930 1056 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
931 1057
932 1058 #print array.shape
933 1059 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
934 1060 #print tmp.shape
935 1061
936 1062 # fig = plt.figure(figsize=(6,5))
937 1063 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
938 1064 # ax = fig.add_axes([left, bottom, width, height])
939 1065 # x = range(num_prof)
940 1066 # y = range(num_hei)
941 1067 # cp = ax.contour(y,x,array[ii,:,:])
942 1068 # ax.clabel(cp, inline=True,fontsize=10)
943 1069 # plt.show()
944 1070 return array
945 1071
946 1072
947 1073 class IntegrationFaradaySpectra(Operation):
948 1074
949 1075 __profIndex = 0
950 1076 __withOverapping = False
951 1077
952 1078 __byTime = False
953 1079 __initime = None
954 1080 __lastdatatime = None
955 1081 __integrationtime = None
956 1082
957 1083 __buffer_spc = None
958 1084 __buffer_cspc = None
959 1085 __buffer_dc = None
960 1086
961 1087 __dataReady = False
962 1088
963 1089 __timeInterval = None
964 1090
965 1091 n = None
966 1092
967 1093 def __init__(self):
968 1094
969 1095 Operation.__init__(self)
970 1096
971 1097 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None):
972 1098 """
973 1099 Set the parameters of the integration class.
974 1100
975 1101 Inputs:
976 1102
977 1103 n : Number of coherent integrations
978 1104 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
979 1105 overlapping :
980 1106
981 1107 """
982 1108
983 1109 self.__initime = None
984 1110 self.__lastdatatime = 0
985 1111
986 1112 self.__buffer_spc = []
987 1113 self.__buffer_cspc = []
988 1114 self.__buffer_dc = 0
989 1115
990 1116 self.__profIndex = 0
991 1117 self.__dataReady = False
992 1118 self.__byTime = False
993 1119
994 1120 #self.ByLags = dataOut.ByLags ###REDEFINIR
995 1121 self.ByLags = False
996 1122
997 1123 if DPL != None:
998 1124 self.DPL=DPL
999 1125 else:
1000 1126 #self.DPL=dataOut.DPL ###REDEFINIR
1001 1127 self.DPL=0
1002 1128
1003 1129 if n is None and timeInterval is None:
1004 1130 raise ValueError("n or timeInterval should be specified ...")
1005 1131
1006 1132 if n is not None:
1007 1133 self.n = int(n)
1008 1134 else:
1009 1135
1010 1136 self.__integrationtime = int(timeInterval)
1011 1137 self.n = None
1012 1138 self.__byTime = True
1013 1139
1014 1140 def putData(self, data_spc, data_cspc, data_dc):
1015 1141 """
1016 1142 Add a profile to the __buffer_spc and increase in one the __profileIndex
1017 1143
1018 1144 """
1019 1145
1020 1146 self.__buffer_spc.append(data_spc)
1021 1147
1022 1148 if data_cspc is None:
1023 1149 self.__buffer_cspc = None
1024 1150 else:
1025 1151 self.__buffer_cspc.append(data_cspc)
1026 1152
1027 1153 if data_dc is None:
1028 1154 self.__buffer_dc = None
1029 1155 else:
1030 1156 self.__buffer_dc += data_dc
1031 1157
1032 1158 self.__profIndex += 1
1033 1159
1034 1160 return
1035 1161
1036 1162 def hildebrand_sekhon_Integration(self,data,navg):
1037 1163
1038 1164 sortdata = numpy.sort(data, axis=None)
1039 1165 sortID=data.argsort()
1040 1166 lenOfData = len(sortdata)
1041 1167 nums_min = lenOfData*0.75
1042 1168 if nums_min <= 5:
1043 1169 nums_min = 5
1044 1170 sump = 0.
1045 1171 sumq = 0.
1046 1172 j = 0
1047 1173 cont = 1
1048 1174 while((cont == 1)and(j < lenOfData)):
1049 1175 sump += sortdata[j]
1050 1176 sumq += sortdata[j]**2
1051 1177 if j > nums_min:
1052 1178 rtest = float(j)/(j-1) + 1.0/navg
1053 1179 if ((sumq*j) > (rtest*sump**2)):
1054 1180 j = j - 1
1055 1181 sump = sump - sortdata[j]
1056 1182 sumq = sumq - sortdata[j]**2
1057 1183 cont = 0
1058 1184 j += 1
1059 1185 #lnoise = sump / j
1060 1186
1061 1187 return j,sortID
1062 1188
1063 1189 def pushData(self):
1064 1190 """
1065 1191 Return the sum of the last profiles and the profiles used in the sum.
1066 1192
1067 1193 Affected:
1068 1194
1069 1195 self.__profileIndex
1070 1196
1071 1197 """
1072 1198 bufferH=None
1073 1199 buffer=None
1074 1200 buffer1=None
1075 1201 buffer_cspc=None
1076 1202 self.__buffer_spc=numpy.array(self.__buffer_spc)
1077 1203 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1078 1204 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1079 1205 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1080 1206 for k in range(7,self.nHeights):
1081 1207 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1082 1208 outliers_IDs_cspc=[]
1083 1209 cspc_outliers_exist=False
1084 1210 for i in range(self.nChannels):#dataOut.nChannels):
1085 1211
1086 1212 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1087 1213 indexes=[]
1088 1214 #sortIDs=[]
1089 1215 outliers_IDs=[]
1090 1216
1091 1217 for j in range(self.nProfiles):
1092 1218 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1093 1219 # continue
1094 1220 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1095 1221 # continue
1096 1222 buffer=buffer1[:,j]
1097 1223 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1098 1224
1099 1225 indexes.append(index)
1100 1226 #sortIDs.append(sortID)
1101 1227 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1102 1228
1103 1229 outliers_IDs=numpy.array(outliers_IDs)
1104 1230 outliers_IDs=outliers_IDs.ravel()
1105 1231 outliers_IDs=numpy.unique(outliers_IDs)
1106 1232 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1107 1233 indexes=numpy.array(indexes)
1108 1234 indexmin=numpy.min(indexes)
1109 1235
1110 1236 if indexmin != buffer1.shape[0]:
1111 1237 cspc_outliers_exist=True
1112 1238 ###sortdata=numpy.sort(buffer1,axis=0)
1113 1239 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
1114 1240 lt=outliers_IDs
1115 1241 avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1116 1242
1117 1243 for p in list(outliers_IDs):
1118 1244 buffer1[p,:]=avg
1119 1245
1120 1246 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1121 1247 ###cspc IDs
1122 1248 #indexmin_cspc+=indexmin_cspc
1123 1249 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1124 1250
1125 1251 #if not breakFlag:
1126 1252 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1127 1253 if cspc_outliers_exist:
1128 1254 #sortdata=numpy.sort(buffer_cspc,axis=0)
1129 1255 #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0)
1130 1256 lt=outliers_IDs_cspc
1131 1257
1132 1258 avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1133 1259 for p in list(outliers_IDs_cspc):
1134 1260 buffer_cspc[p,:]=avg
1135 1261
1136 1262 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1137 1263 #else:
1138 1264 #break
1139 1265
1140 1266
1141 1267
1142 1268
1143 1269 buffer=None
1144 1270 bufferH=None
1145 1271 buffer1=None
1146 1272 buffer_cspc=None
1147 1273
1148 1274 #print("cpsc",self.__buffer_cspc[:,0,0,0,0])
1149 1275 #print(self.__profIndex)
1150 1276 #exit()
1151 1277
1152 1278 buffer=None
1153 1279 #print(self.__buffer_spc[:,1,3,20,0])
1154 1280 #print(self.__buffer_spc[:,1,5,37,0])
1155 1281 data_spc = numpy.sum(self.__buffer_spc,axis=0)
1156 1282 data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1157 1283
1158 1284 #print(numpy.shape(data_spc))
1159 1285 #data_spc[1,4,20,0]=numpy.nan
1160 1286
1161 1287 #data_cspc = self.__buffer_cspc
1162 1288 data_dc = self.__buffer_dc
1163 1289 n = self.__profIndex
1164 1290
1165 1291 self.__buffer_spc = []
1166 1292 self.__buffer_cspc = []
1167 1293 self.__buffer_dc = 0
1168 1294 self.__profIndex = 0
1169 1295
1170 1296 return data_spc, data_cspc, data_dc, n
1171 1297
1172 1298 def byProfiles(self, *args):
1173 1299
1174 1300 self.__dataReady = False
1175 1301 avgdata_spc = None
1176 1302 avgdata_cspc = None
1177 1303 avgdata_dc = None
1178 1304
1179 1305 self.putData(*args)
1180 1306
1181 1307 if self.__profIndex == self.n:
1182 1308
1183 1309 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1184 1310 self.n = n
1185 1311 self.__dataReady = True
1186 1312
1187 1313 return avgdata_spc, avgdata_cspc, avgdata_dc
1188 1314
1189 1315 def byTime(self, datatime, *args):
1190 1316
1191 1317 self.__dataReady = False
1192 1318 avgdata_spc = None
1193 1319 avgdata_cspc = None
1194 1320 avgdata_dc = None
1195 1321
1196 1322 self.putData(*args)
1197 1323
1198 1324 if (datatime - self.__initime) >= self.__integrationtime:
1199 1325 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1200 1326 self.n = n
1201 1327 self.__dataReady = True
1202 1328
1203 1329 return avgdata_spc, avgdata_cspc, avgdata_dc
1204 1330
1205 1331 def integrate(self, datatime, *args):
1206 1332
1207 1333 if self.__profIndex == 0:
1208 1334 self.__initime = datatime
1209 1335
1210 1336 if self.__byTime:
1211 1337 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1212 1338 datatime, *args)
1213 1339 else:
1214 1340 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1215 1341
1216 1342 if not self.__dataReady:
1217 1343 return None, None, None, None
1218 1344
1219 1345 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1220 1346
1221 1347 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False):
1222 1348 if n == 1:
1223 1349 return dataOut
1224 1350
1225 1351 dataOut.flagNoData = True
1226 1352
1227 1353 if not self.isConfig:
1228 1354 self.setup(dataOut, n, timeInterval, overlapping,DPL )
1229 1355 self.isConfig = True
1230 1356
1231 1357 if not self.ByLags:
1232 1358 self.nProfiles=dataOut.nProfiles
1233 1359 self.nChannels=dataOut.nChannels
1234 1360 self.nHeights=dataOut.nHeights
1235 1361 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1236 1362 dataOut.data_spc,
1237 1363 dataOut.data_cspc,
1238 1364 dataOut.data_dc)
1239 1365 else:
1240 1366 self.nProfiles=dataOut.nProfiles
1241 1367 self.nChannels=dataOut.nChannels
1242 1368 self.nHeights=dataOut.nHeights
1243 1369 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1244 1370 dataOut.dataLag_spc,
1245 1371 dataOut.dataLag_cspc,
1246 1372 dataOut.dataLag_dc)
1247 1373
1248 1374 if self.__dataReady:
1249 1375
1250 1376 if not self.ByLags:
1251 1377
1252 1378 dataOut.data_spc = numpy.squeeze(avgdata_spc)
1253 1379 dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1254 1380 dataOut.data_dc = avgdata_dc
1255 1381 else:
1256 1382 dataOut.dataLag_spc = avgdata_spc
1257 1383 dataOut.dataLag_cspc = avgdata_cspc
1258 1384 dataOut.dataLag_dc = avgdata_dc
1259 1385
1260 1386 dataOut.data_spc=dataOut.dataLag_spc[:,:,:,dataOut.LagPlot]
1261 1387 dataOut.data_cspc=dataOut.dataLag_cspc[:,:,:,dataOut.LagPlot]
1262 1388 dataOut.data_dc=dataOut.dataLag_dc[:,:,dataOut.LagPlot]
1263 1389
1264 1390
1265 1391 dataOut.nIncohInt *= self.n
1266 1392 dataOut.utctime = avgdatatime
1267 1393 dataOut.flagNoData = False
1268 1394
1269 1395 return dataOut
1270 1396
1271 1397 class removeInterference(Operation):
1272 1398
1273 1399 def removeInterference2(self):
1274 1400
1275 1401 cspc = self.dataOut.data_cspc
1276 1402 spc = self.dataOut.data_spc
1277 1403 Heights = numpy.arange(cspc.shape[2])
1278 1404 realCspc = numpy.abs(cspc)
1279 1405
1280 1406 for i in range(cspc.shape[0]):
1281 1407 LinePower= numpy.sum(realCspc[i], axis=0)
1282 1408 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1283 1409 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1284 1410 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1285 1411 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1286 1412 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1287 1413
1288 1414
1289 1415 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1290 1416 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1291 1417 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1292 1418 cspc[i,InterferenceRange,:] = numpy.NaN
1293 1419
1294 1420 self.dataOut.data_cspc = cspc
1295 1421
1296 1422 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1297 1423
1298 1424 jspectra = self.dataOut.data_spc
1299 1425 jcspectra = self.dataOut.data_cspc
1300 1426 jnoise = self.dataOut.getNoise()
1301 1427 num_incoh = self.dataOut.nIncohInt
1302 1428
1303 1429 num_channel = jspectra.shape[0]
1304 1430 num_prof = jspectra.shape[1]
1305 1431 num_hei = jspectra.shape[2]
1306 1432
1307 1433 # hei_interf
1308 1434 if hei_interf is None:
1309 1435 count_hei = int(num_hei / 2)
1310 1436 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1311 1437 hei_interf = numpy.asarray(hei_interf)[0]
1312 1438 # nhei_interf
1313 1439 if (nhei_interf == None):
1314 1440 nhei_interf = 5
1315 1441 if (nhei_interf < 1):
1316 1442 nhei_interf = 1
1317 1443 if (nhei_interf > count_hei):
1318 1444 nhei_interf = count_hei
1319 1445 if (offhei_interf == None):
1320 1446 offhei_interf = 0
1321 1447
1322 1448 ind_hei = list(range(num_hei))
1323 1449 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1324 1450 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1325 1451 mask_prof = numpy.asarray(list(range(num_prof)))
1326 1452 num_mask_prof = mask_prof.size
1327 1453 comp_mask_prof = [0, num_prof / 2]
1328 1454
1329 1455 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1330 1456 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1331 1457 jnoise = numpy.nan
1332 1458 noise_exist = jnoise[0] < numpy.Inf
1333 1459
1334 1460 # Subrutina de Remocion de la Interferencia
1335 1461 for ich in range(num_channel):
1336 1462 # Se ordena los espectros segun su potencia (menor a mayor)
1337 1463 power = jspectra[ich, mask_prof, :]
1338 1464 power = power[:, hei_interf]
1339 1465 power = power.sum(axis=0)
1340 1466 psort = power.ravel().argsort()
1341 1467
1342 1468 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1343 1469 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1344 1470 offhei_interf, nhei_interf + offhei_interf))]]]
1345 1471
1346 1472 if noise_exist:
1347 1473 # tmp_noise = jnoise[ich] / num_prof
1348 1474 tmp_noise = jnoise[ich]
1349 1475 junkspc_interf = junkspc_interf - tmp_noise
1350 1476 #junkspc_interf[:,comp_mask_prof] = 0
1351 1477
1352 1478 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1353 1479 jspc_interf = jspc_interf.transpose()
1354 1480 # Calculando el espectro de interferencia promedio
1355 1481 noiseid = numpy.where(
1356 1482 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1357 1483 noiseid = noiseid[0]
1358 1484 cnoiseid = noiseid.size
1359 1485 interfid = numpy.where(
1360 1486 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1361 1487 interfid = interfid[0]
1362 1488 cinterfid = interfid.size
1363 1489
1364 1490 if (cnoiseid > 0):
1365 1491 jspc_interf[noiseid] = 0
1366 1492
1367 1493 # Expandiendo los perfiles a limpiar
1368 1494 if (cinterfid > 0):
1369 1495 new_interfid = (
1370 1496 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1371 1497 new_interfid = numpy.asarray(new_interfid)
1372 1498 new_interfid = {x for x in new_interfid}
1373 1499 new_interfid = numpy.array(list(new_interfid))
1374 1500 new_cinterfid = new_interfid.size
1375 1501 else:
1376 1502 new_cinterfid = 0
1377 1503
1378 1504 for ip in range(new_cinterfid):
1379 1505 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1380 1506 jspc_interf[new_interfid[ip]
1381 1507 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1382 1508
1383 1509 jspectra[ich, :, ind_hei] = jspectra[ich, :,
1384 1510 ind_hei] - jspc_interf # Corregir indices
1385 1511
1386 1512 # Removiendo la interferencia del punto de mayor interferencia
1387 1513 ListAux = jspc_interf[mask_prof].tolist()
1388 1514 maxid = ListAux.index(max(ListAux))
1389 1515
1390 1516 if cinterfid > 0:
1391 1517 for ip in range(cinterfid * (interf == 2) - 1):
1392 1518 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1393 1519 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1394 1520 cind = len(ind)
1395 1521
1396 1522 if (cind > 0):
1397 1523 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1398 1524 (1 + (numpy.random.uniform(cind) - 0.5) /
1399 1525 numpy.sqrt(num_incoh))
1400 1526
1401 1527 ind = numpy.array([-2, -1, 1, 2])
1402 1528 xx = numpy.zeros([4, 4])
1403 1529
1404 1530 for id1 in range(4):
1405 1531 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1406 1532
1407 1533 xx_inv = numpy.linalg.inv(xx)
1408 1534 xx = xx_inv[:, 0]
1409 1535 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1410 1536 yy = jspectra[ich, mask_prof[ind], :]
1411 1537 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
1412 1538 yy.transpose(), xx)
1413 1539
1414 1540 indAux = (jspectra[ich, :, :] < tmp_noise *
1415 1541 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1416 1542 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1417 1543 (1 - 1 / numpy.sqrt(num_incoh))
1418 1544
1419 1545 # Remocion de Interferencia en el Cross Spectra
1420 1546 if jcspectra is None:
1421 1547 return jspectra, jcspectra
1422 1548 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1423 1549 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1424 1550
1425 1551 for ip in range(num_pairs):
1426 1552
1427 1553 #-------------------------------------------
1428 1554
1429 1555 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1430 1556 cspower = cspower[:, hei_interf]
1431 1557 cspower = cspower.sum(axis=0)
1432 1558
1433 1559 cspsort = cspower.ravel().argsort()
1434 1560 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1435 1561 offhei_interf, nhei_interf + offhei_interf))]]]
1436 1562 junkcspc_interf = junkcspc_interf.transpose()
1437 1563 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1438 1564
1439 1565 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1440 1566
1441 1567 median_real = int(numpy.median(numpy.real(
1442 1568 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1443 1569 median_imag = int(numpy.median(numpy.imag(
1444 1570 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1445 1571 comp_mask_prof = [int(e) for e in comp_mask_prof]
1446 1572 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1447 1573 median_real, median_imag)
1448 1574
1449 1575 for iprof in range(num_prof):
1450 1576 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1451 1577 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1452 1578
1453 1579 # Removiendo la Interferencia
1454 1580 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1455 1581 :, ind_hei] - jcspc_interf
1456 1582
1457 1583 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1458 1584 maxid = ListAux.index(max(ListAux))
1459 1585
1460 1586 ind = numpy.array([-2, -1, 1, 2])
1461 1587 xx = numpy.zeros([4, 4])
1462 1588
1463 1589 for id1 in range(4):
1464 1590 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1465 1591
1466 1592 xx_inv = numpy.linalg.inv(xx)
1467 1593 xx = xx_inv[:, 0]
1468 1594
1469 1595 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1470 1596 yy = jcspectra[ip, mask_prof[ind], :]
1471 1597 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1472 1598
1473 1599 # Guardar Resultados
1474 1600 self.dataOut.data_spc = jspectra
1475 1601 self.dataOut.data_cspc = jcspectra
1476 1602
1477 1603 return 1
1478 1604
1479 1605 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
1480 1606
1481 1607 self.dataOut = dataOut
1482 1608
1483 1609 if mode == 1:
1484 1610 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1485 1611 elif mode == 2:
1486 1612 self.removeInterference2()
1487 1613
1488 1614 return self.dataOut
1489 1615
1490 1616
1491 1617 class IncohInt(Operation):
1492 1618
1493 1619 __profIndex = 0
1494 1620 __withOverapping = False
1495 1621
1496 1622 __byTime = False
1497 1623 __initime = None
1498 1624 __lastdatatime = None
1499 1625 __integrationtime = None
1500 1626
1501 1627 __buffer_spc = None
1502 1628 __buffer_cspc = None
1503 1629 __buffer_dc = None
1504 1630
1505 1631 __dataReady = False
1506 1632
1507 1633 __timeInterval = None
1508 1634
1509 1635 n = None
1510 1636
1511 1637 def __init__(self):
1512 1638
1513 1639 Operation.__init__(self)
1514 1640
1515 1641 def setup(self, n=None, timeInterval=None, overlapping=False):
1516 1642 """
1517 1643 Set the parameters of the integration class.
1518 1644
1519 1645 Inputs:
1520 1646
1521 1647 n : Number of coherent integrations
1522 1648 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1523 1649 overlapping :
1524 1650
1525 1651 """
1526 1652
1527 1653 self.__initime = None
1528 1654 self.__lastdatatime = 0
1529 1655
1530 1656 self.__buffer_spc = 0
1531 1657 self.__buffer_cspc = 0
1532 1658 self.__buffer_dc = 0
1533 1659
1534 1660 self.__profIndex = 0
1535 1661 self.__dataReady = False
1536 1662 self.__byTime = False
1537 1663
1538 1664 if n is None and timeInterval is None:
1539 1665 raise ValueError("n or timeInterval should be specified ...")
1540 1666
1541 1667 if n is not None:
1542 1668 self.n = int(n)
1543 1669 else:
1544 1670
1545 1671 self.__integrationtime = int(timeInterval)
1546 1672 self.n = None
1547 1673 self.__byTime = True
1548 1674
1549 1675 def putData(self, data_spc, data_cspc, data_dc):
1550 1676 """
1551 1677 Add a profile to the __buffer_spc and increase in one the __profileIndex
1552 1678
1553 1679 """
1554 1680
1555 1681 self.__buffer_spc += data_spc
1556 1682
1557 1683 if data_cspc is None:
1558 1684 self.__buffer_cspc = None
1559 1685 else:
1560 1686 self.__buffer_cspc += data_cspc
1561 1687
1562 1688 if data_dc is None:
1563 1689 self.__buffer_dc = None
1564 1690 else:
1565 1691 self.__buffer_dc += data_dc
1566 1692
1567 1693 self.__profIndex += 1
1568 1694
1569 1695 return
1570 1696
1571 1697 def pushData(self):
1572 1698 """
1573 1699 Return the sum of the last profiles and the profiles used in the sum.
1574 1700
1575 1701 Affected:
1576 1702
1577 1703 self.__profileIndex
1578 1704
1579 1705 """
1580 1706
1581 1707 data_spc = self.__buffer_spc
1582 1708 data_cspc = self.__buffer_cspc
1583 1709 data_dc = self.__buffer_dc
1584 1710 n = self.__profIndex
1585 1711
1586 1712 self.__buffer_spc = 0
1587 1713 self.__buffer_cspc = 0
1588 1714 self.__buffer_dc = 0
1589 1715 self.__profIndex = 0
1590 1716
1591 1717 return data_spc, data_cspc, data_dc, n
1592 1718
1593 1719 def byProfiles(self, *args):
1594 1720
1595 1721 self.__dataReady = False
1596 1722 avgdata_spc = None
1597 1723 avgdata_cspc = None
1598 1724 avgdata_dc = None
1599 1725
1600 1726 self.putData(*args)
1601 1727
1602 1728 if self.__profIndex == self.n:
1603 1729
1604 1730 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1605 1731 self.n = n
1606 1732 self.__dataReady = True
1607 1733
1608 1734 return avgdata_spc, avgdata_cspc, avgdata_dc
1609 1735
1610 1736 def byTime(self, datatime, *args):
1611 1737
1612 1738 self.__dataReady = False
1613 1739 avgdata_spc = None
1614 1740 avgdata_cspc = None
1615 1741 avgdata_dc = None
1616 1742
1617 1743 self.putData(*args)
1618 1744
1619 1745 if (datatime - self.__initime) >= self.__integrationtime:
1620 1746 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1621 1747 self.n = n
1622 1748 self.__dataReady = True
1623 1749
1624 1750 return avgdata_spc, avgdata_cspc, avgdata_dc
1625 1751
1626 1752 def integrate(self, datatime, *args):
1627 1753
1628 1754 if self.__profIndex == 0:
1629 1755 self.__initime = datatime
1630 1756
1631 1757 if self.__byTime:
1632 1758 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1633 1759 datatime, *args)
1634 1760 else:
1635 1761 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1636 1762
1637 1763 if not self.__dataReady:
1638 1764 return None, None, None, None
1639 1765
1640 1766 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1641 1767
1642 1768 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1643 1769 if n == 1:
1644 1770 return dataOut
1645 1771
1646 1772 dataOut.flagNoData = True
1647 1773
1648 1774 if not self.isConfig:
1649 1775 self.setup(n, timeInterval, overlapping)
1650 1776 self.isConfig = True
1651 1777
1652 1778 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1653 1779 dataOut.data_spc,
1654 1780 dataOut.data_cspc,
1655 1781 dataOut.data_dc)
1656 1782
1657 1783 if self.__dataReady:
1658 1784
1659 1785 dataOut.data_spc = avgdata_spc
1660 1786 dataOut.data_cspc = avgdata_cspc
1661 1787 dataOut.data_dc = avgdata_dc
1662 1788 dataOut.nIncohInt *= self.n
1663 1789 dataOut.utctime = avgdatatime
1664 1790 dataOut.flagNoData = False
1665 1791
1666 1792 return dataOut
1667 1793
1668 1794 class dopplerFlip(Operation):
1669 1795
1670 1796 def run(self, dataOut):
1671 1797 # arreglo 1: (num_chan, num_profiles, num_heights)
1672 1798 self.dataOut = dataOut
1673 1799 # JULIA-oblicua, indice 2
1674 1800 # arreglo 2: (num_profiles, num_heights)
1675 1801 jspectra = self.dataOut.data_spc[2]
1676 1802 jspectra_tmp = numpy.zeros(jspectra.shape)
1677 1803 num_profiles = jspectra.shape[0]
1678 1804 freq_dc = int(num_profiles / 2)
1679 1805 # Flip con for
1680 1806 for j in range(num_profiles):
1681 1807 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1682 1808 # Intercambio perfil de DC con perfil inmediato anterior
1683 1809 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1684 1810 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1685 1811 # canal modificado es re-escrito en el arreglo de canales
1686 1812 self.dataOut.data_spc[2] = jspectra_tmp
1687 1813
1688 1814 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now