##// END OF EJS Templates
schain
RSS
Interference from YAGI removed, remHeightsIppInterf added
joabAM -
r1579:8a1b75bc0681
parent child
Show More
Side by Side Unified
Context file:
r1579:8a1b75bc0681 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,1154 +1,1159
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,ProcessingHeader
21 21 from schainpy.model.data import _noise
22 22 SPEED_OF_LIGHT = 3e8
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.5
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 return lnoise
110 110 '''
111 111 return _noise.hildebrand_sekhon(sortdata, navg)
112 112
113 113
114 114 class Beam:
115 115
116 116 def __init__(self):
117 117 self.codeList = []
118 118 self.azimuthList = []
119 119 self.zenithList = []
120 120
121 121
122 122
123 123 class GenericData(object):
124 124
125 125 flagNoData = True
126 126
127 127 def copy(self, inputObj=None):
128 128
129 129 if inputObj == None:
130 130 return copy.deepcopy(self)
131 131
132 132 for key in list(inputObj.__dict__.keys()):
133 133
134 134 attribute = inputObj.__dict__[key]
135 135
136 136 # If this attribute is a tuple or list
137 137 if type(inputObj.__dict__[key]) in (tuple, list):
138 138 self.__dict__[key] = attribute[:]
139 139 continue
140 140
141 141 # If this attribute is another object or instance
142 142 if hasattr(attribute, '__dict__'):
143 143 self.__dict__[key] = attribute.copy()
144 144 continue
145 145
146 146 self.__dict__[key] = inputObj.__dict__[key]
147 147
148 148 def deepcopy(self):
149 149
150 150 return copy.deepcopy(self)
151 151
152 152 def isEmpty(self):
153 153
154 154 return self.flagNoData
155 155
156 156 def isReady(self):
157 157
158 158 return not self.flagNoData
159 159
160 160
161 161 class JROData(GenericData):
162 162
163 163 useInputBuffer = False
164 164 buffer_empty = True
165 165
166 166 systemHeaderObj = SystemHeader()
167 167 radarControllerHeaderObj = RadarControllerHeader()
168 168 type = None
169 169 datatype = None # dtype but in string
170 170 nProfiles = None
171 171 heightList = None
172 172 channelList = None
173 173 flagDiscontinuousBlock = False
174 174 useLocalTime = False
175 175 utctime = None
176 176 timeZone = None
177 177 dstFlag = None
178 178 errorCount = None
179 179 blocksize = None
180 180 flagDecodeData = False # asumo q la data no esta decodificada
181 181 flagDeflipData = False # asumo q la data no esta sin flip
182 182 flagShiftFFT = False
183 183 nCohInt = None
184 184 windowOfFilter = 1
185 185 C = 3e8
186 186 frequency = 49.92e6
187 187 realtime = False
188 188 beacon_heiIndexList = None
189 189 last_block = None
190 190 blocknow = None
191 191 azimuth = None
192 192 zenith = None
193 193
194 194 profileIndex = None
195 195 error = None
196 196 data = None
197 197 nmodes = None
198 198 metadata_list = ['heightList', 'timeZone', 'type']
199 199 codeList = []
200 200 azimuthList = []
201 201 elevationList = []
202 202 last_noise = None
203 203 __ipp = None
204 204 __ippSeconds = None
205 205 sampled_heightsFFT = None
206 206 pulseLength_TxA = None
207 207 deltaHeight = None
208 208 __code = None
209 209 __nCode = None
210 210 __nBaud = None
211 211 unitsDescription = "The units of the parameters are according to the International System of units (Seconds, Meter, Hertz, ...), except \
212 212 the parameters related to distances such as heightList, or heightResolution wich are in Km"
213 213
214 214
215 215 def __str__(self):
216 216
217 217 return '{} - {}'.format(self.type, self.datatime())
218 218
219 219 def getNoise(self):
220 220
221 221 raise NotImplementedError
222 222
223 223 @property
224 224 def nChannels(self):
225 225
226 226 return len(self.channelList)
227 227
228 228 @property
229 229 def channelIndexList(self):
230 230
231 231 return list(range(self.nChannels))
232 232
233 233 @property
234 234 def nHeights(self):
235 235
236 236 return len(self.heightList)
237 237
238 238 def getDeltaH(self):
239 239
240 240 return self.heightList[1] - self.heightList[0]
241 241
242 242 @property
243 243 def ltctime(self):
244 244
245 245 if self.useLocalTime:
246 246 return self.utctime - self.timeZone * 60
247 247
248 248 return self.utctime
249 249
250 250 @property
251 251 def datatime(self):
252 252
253 253 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
254 254 return datatimeValue
255 255
256 256 def getTimeRange(self):
257 257
258 258 datatime = []
259 259
260 260 datatime.append(self.ltctime)
261 261 datatime.append(self.ltctime + self.timeInterval + 1)
262 262
263 263 datatime = numpy.array(datatime)
264 264
265 265 return datatime
266 266
267 267 def getFmaxTimeResponse(self):
268 268
269 269 period = (10**-6) * self.getDeltaH() / (0.15)
270 270
271 271 PRF = 1. / (period * self.nCohInt)
272 272
273 273 fmax = PRF
274 274
275 275 return fmax
276 276
277 277 def getFmax(self):
278 278 PRF = 1. / (self.__ippSeconds * self.nCohInt)
279 279
280 280 fmax = PRF
281 281 return fmax
282 282
283 283 def getVmax(self):
284 284
285 285 _lambda = self.C / self.frequency
286 286
287 287 vmax = self.getFmax() * _lambda / 2
288 288
289 289 return vmax
290 290
291 291 @property
292 292 def ippSeconds(self):
293 293 '''
294 294 '''
295 295 #return self.radarControllerHeaderObj.ippSeconds
296 296 return self.__ippSeconds
297 297
298 298 @ippSeconds.setter
299 299 def ippSeconds(self, ippSeconds):
300 300 '''
301 301 '''
302 302 #self.radarControllerHeaderObj.ippSeconds = ippSeconds
303 303 self.__ippSeconds = ippSeconds
304 304 self.__ipp = ippSeconds*SPEED_OF_LIGHT/2000.0
305 305
306 306
307 307 @property
308 308 def code(self):
309 309 '''
310 310 '''
311 311 return self.__code
312 312
313 313 @code.setter
314 314 def code(self, code):
315 315 '''
316 316 '''
317 317 self.__code = code
318 318 #
319 319 @property
320 320 def nCode(self):
321 321 '''
322 322 '''
323 323 return self.__nCode
324 324
325 325 @nCode.setter
326 326 def nCode(self, ncode):
327 327 '''
328 328 '''
329 329 self.__nCode = ncode
330 330
331 331 @property
332 332 def nBaud(self):
333 333 '''
334 334 '''
335 335 return self.__nBaud
336 336
337 337 @nBaud.setter
338 338 def nBaud(self, nbaud):
339 339 '''
340 340 '''
341 341 self.__nBaud = nbaud
342 342
343 343 @property
344 344 def ipp(self):
345 345 '''
346 346 '''
347 347 return self.__ipp
348 348 #return self.radarControllerHeaderObj.ipp
349 349
350 350 @ipp.setter
351 351 def ipp(self, ipp):
352 352 '''
353 353 '''
354 354 self.__ipp = ipp
355 355 #self.radarControllerHeaderObj.ipp = ipp
356 356
357 357 @property
358 358 def metadata(self):
359 359 '''
360 360 '''
361 361
362 362 return {attr: getattr(self, attr) for attr in self.metadata_list}
363 363
364 364
365 365 class Voltage(JROData):
366 366
367 367 dataPP_POW = None
368 368 dataPP_DOP = None
369 369 dataPP_WIDTH = None
370 370 dataPP_SNR = None
371 flagProfilesByRange = False
372 nProfilesByRange = None
371 373
372 374 def __init__(self):
373 375 '''
374 376 Constructor
375 377 '''
376 378
377 379 self.useLocalTime = True
378 380 self.radarControllerHeaderObj = RadarControllerHeader()
379 381 self.systemHeaderObj = SystemHeader()
380 382 self.processingHeaderObj = ProcessingHeader()
381 383 self.type = "Voltage"
382 384 self.data = None
383 385 self.nProfiles = None
384 386 self.heightList = None
385 387 self.channelList = None
386 388 self.flagNoData = True
387 389 self.flagDiscontinuousBlock = False
388 390 self.utctime = None
389 391 self.timeZone = 0
390 392 self.dstFlag = None
391 393 self.errorCount = None
392 394 self.nCohInt = None
393 395 self.blocksize = None
394 396 self.flagCohInt = False
395 397 self.flagDecodeData = False # asumo q la data no esta decodificada
396 398 self.flagDeflipData = False # asumo q la data no esta sin flip
397 399 self.flagShiftFFT = False
398 400 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
399 401 self.profileIndex = 0
400 402 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
401 403 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
402 404
403 405 def getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None):
404 406 """
405 407 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
406 408
407 409 Return:
408 410 noiselevel
409 411 """
410 412
411 413 if channel != None:
412 414 data = self.data[channel,ymin_index:ymax_index]
413 415 nChannels = 1
414 416 else:
415 417 data = self.data[:,ymin_index:ymax_index]
416 418 nChannels = self.nChannels
417 419
418 420 noise = numpy.zeros(nChannels)
419 421 power = data * numpy.conjugate(data)
420 422
421 423 for thisChannel in range(nChannels):
422 424 if nChannels == 1:
423 425 daux = power[:].real
424 426 else:
425 427 daux = power[thisChannel, :].real
426 428 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
427 429
428 430 return noise
429 431
430 432 def getNoise(self, type=1, channel=None,ymin_index=None, ymax_index=None):
431 433
432 434 if type == 1:
433 435 noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index)
434 436
435 437 return noise
436 438
437 439 def getPower(self, channel=None):
438 440
439 441 if channel != None:
440 442 data = self.data[channel]
441 443 else:
442 444 data = self.data
443 445
444 446 power = data * numpy.conjugate(data)
445 447 powerdB = 10 * numpy.log10(power.real)
446 448 powerdB = numpy.squeeze(powerdB)
447 449
448 450 return powerdB
449 451
450 452 @property
451 453 def timeInterval(self):
452 454
453 455 return self.ippSeconds * self.nCohInt
454 456
455 457 noise = property(getNoise, "I'm the 'nHeights' property.")
456 458
457 459
458 460 class Spectra(JROData):
459 461
460 462 data_outlier = None
463 flagProfilesByRange = False
464 nProfilesByRange = None
461 465
462 466 def __init__(self):
463 467 '''
464 468 Constructor
465 469 '''
466 470
467 471 self.data_dc = None
468 472 self.data_spc = None
469 473 self.data_cspc = None
470 474 self.useLocalTime = True
471 475 self.radarControllerHeaderObj = RadarControllerHeader()
472 476 self.systemHeaderObj = SystemHeader()
473 477 self.processingHeaderObj = ProcessingHeader()
474 478 self.type = "Spectra"
475 479 self.timeZone = 0
476 480 self.nProfiles = None
477 481 self.heightList = None
478 482 self.channelList = None
479 483 self.pairsList = None
480 484 self.flagNoData = True
481 485 self.flagDiscontinuousBlock = False
482 486 self.utctime = None
483 487 self.nCohInt = None
484 488 self.nIncohInt = None
485 489 self.blocksize = None
486 490 self.nFFTPoints = None
487 491 self.wavelength = None
488 492 self.flagDecodeData = False # asumo q la data no esta decodificada
489 493 self.flagDeflipData = False # asumo q la data no esta sin flip
490 494 self.flagShiftFFT = False
491 495 self.ippFactor = 1
492 496 self.beacon_heiIndexList = []
493 497 self.noise_estimation = None
494 498 self.codeList = []
495 499 self.azimuthList = []
496 500 self.elevationList = []
497 501 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
498 502 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
499 503
500 504
501 505
502 506
503 507 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
504 508 """
505 509 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
506 510
507 511 Return:
508 512 noiselevel
509 513 """
510 514 # if hasattr(self.nIncohInt, "__len__"): #nIncohInt is a matrix
511 515 #
512 516 # heis = self.data_spc.shape[2]
513 517 #
514 518 # noise = numpy.zeros((self.nChannels, heis))
515 519 # for hei in range(heis):
516 520 # for channel in range(self.nChannels):
517 521 # daux = self.data_spc[channel, xmin_index:xmax_index, hei]
518 522 #
519 523 # noise[channel,hei] = hildebrand_sekhon(daux, self.nIncohInt[channel,hei])
520 524 #
521 525 # else:
522 526 # noise = numpy.zeros(self.nChannels)
523 527 # for channel in range(self.nChannels):
524 528 # daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index]
525 529 #
526 530 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
527 531 noise = numpy.zeros(self.nChannels)
528 532
529 533 for channel in range(self.nChannels):
530 534 daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index]
531 535
532 536 noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt[channel])
533 537
534 538 return noise
535 539
536 540 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
537 541
538 542 if self.noise_estimation is not None:
539 543 # this was estimated by getNoise Operation defined in jroproc_spectra.py
540 544 return self.noise_estimation
541 545 else:
542 546 noise = self.getNoisebyHildebrand(
543 547 xmin_index, xmax_index, ymin_index, ymax_index)
544 548 return noise
545 549
546 550 def getFreqRangeTimeResponse(self, extrapoints=0):
547 551
548 552 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
549 553 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
550 554
551 555 return freqrange
552 556
553 557 def getAcfRange(self, extrapoints=0):
554 558
555 559 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
556 560 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
557 561
558 562 return freqrange
559 563
560 564 def getFreqRange(self, extrapoints=0):
561 565
562 566 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
563 567 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
564 568
565 569 return freqrange
566 570
567 571 def getVelRange(self, extrapoints=0):
568 572
569 573 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
570 574 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
571 575
572 576 if self.nmodes:
573 577 return velrange/self.nmodes
574 578 else:
575 579 return velrange
576 580
577 581 @property
578 582 def nPairs(self):
579 583
580 584 return len(self.pairsList)
581 585
582 586 @property
583 587 def pairsIndexList(self):
584 588
585 589 return list(range(self.nPairs))
586 590
587 591 @property
588 592 def normFactor(self):
589 593
590 594 pwcode = 1
591
592 595 if self.flagDecodeData:
593 596 pwcode = numpy.sum(self.code[0]**2)
594 597 #print(self.flagDecodeData, pwcode)
595 598 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
596 599 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
597
600 if self.flagProfilesByRange:
601 normFactor *= (self.nProfilesByRange/self.nProfilesByRange.max())
598 602
599 603 return normFactor
600 604
601 605 @property
602 606 def flag_cspc(self):
603 607
604 608 if self.data_cspc is None:
605 609 return True
606 610
607 611 return False
608 612
609 613 @property
610 614 def flag_dc(self):
611 615
612 616 if self.data_dc is None:
613 617 return True
614 618
615 619 return False
616 620
617 621 @property
618 622 def timeInterval(self):
619 623
620 624 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
621 625 if self.nmodes:
622 626 return self.nmodes*timeInterval
623 627 else:
624 628 return timeInterval
625 629
626 630 def getPower(self):
627 631
628 632 factor = self.normFactor
629 633 power = numpy.zeros( (self.nChannels,self.nHeights) )
630 634 for ch in range(self.nChannels):
631 635 z = None
632 636 if hasattr(factor,'shape'):
633 637 if factor.ndim > 1:
634 638 z = self.data_spc[ch]/factor[ch]
635 639 else:
636 640 z = self.data_spc[ch]/factor
637 641 else:
638 642 z = self.data_spc[ch]/factor
639 643 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
640 644 avg = numpy.average(z, axis=0)
641 645 power[ch] = 10 * numpy.log10(avg)
642 646 return power
643 647
644 648 @property
645 649 def max_nIncohInt(self):
646 650
647 651 ints = numpy.zeros(self.nChannels)
648 652 for ch in range(self.nChannels):
649 653 if hasattr(self.nIncohInt,'shape'):
650 654 if self.nIncohInt.ndim > 1:
651 655 ints[ch,] = self.nIncohInt[ch].max()
652 656 else:
653 657 ints[ch,] = self.nIncohInt
654 658 self.nIncohInt = int(self.nIncohInt)
655 659 else:
656 660 ints[ch,] = self.nIncohInt
657 661
658 662 return ints
659 663
660 664
661 665 def getCoherence(self, pairsList=None, phase=False):
662 666
663 667 z = []
664 668 if pairsList is None:
665 669 pairsIndexList = self.pairsIndexList
666 670 else:
667 671 pairsIndexList = []
668 672 for pair in pairsList:
669 673 if pair not in self.pairsList:
670 674 raise ValueError("Pair %s is not in dataOut.pairsList" % (
671 675 pair))
672 676 pairsIndexList.append(self.pairsList.index(pair))
673 677 for i in range(len(pairsIndexList)):
674 678 pair = self.pairsList[pairsIndexList[i]]
675 679 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
676 680 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
677 681 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
678 682 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
679 683 if phase:
680 684 data = numpy.arctan2(avgcoherenceComplex.imag,
681 685 avgcoherenceComplex.real) * 180 / numpy.pi
682 686 else:
683 687 data = numpy.abs(avgcoherenceComplex)
684 688
685 689 z.append(data)
686 690
687 691 return numpy.array(z)
688 692
689 693 def setValue(self, value):
690 694
691 695 print("This property should not be initialized", value)
692 696
693 697 return
694 698
695 699 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
696 700
697 701
698 702 class SpectraHeis(Spectra):
699 703
700 704 def __init__(self):
701 705
702 706 self.radarControllerHeaderObj = RadarControllerHeader()
703 707 self.systemHeaderObj = SystemHeader()
704 708 self.type = "SpectraHeis"
705 709 self.nProfiles = None
706 710 self.heightList = None
707 711 self.channelList = None
708 712 self.flagNoData = True
709 713 self.flagDiscontinuousBlock = False
710 714 self.utctime = None
711 715 self.blocksize = None
712 716 self.profileIndex = 0
713 717 self.nCohInt = 1
714 718 self.nIncohInt = 1
715 719
716 720 @property
717 721 def normFactor(self):
718 722 pwcode = 1
719 723 if self.flagDecodeData:
720 724 pwcode = numpy.sum(self.code[0]**2)
721 725
722 726 normFactor = self.nIncohInt * self.nCohInt * pwcode
723 727
724 728 return normFactor
725 729
726 730 @property
727 731 def timeInterval(self):
728 732
729 733 return self.ippSeconds * self.nCohInt * self.nIncohInt
730 734
731 735
732 736 class Fits(JROData):
733 737
734 738 def __init__(self):
735 739
736 740 self.type = "Fits"
737 741 self.nProfiles = None
738 742 self.heightList = None
739 743 self.channelList = None
740 744 self.flagNoData = True
741 745 self.utctime = None
742 746 self.nCohInt = 1
743 747 self.nIncohInt = 1
744 748 self.useLocalTime = True
745 749 self.profileIndex = 0
746 750 self.timeZone = 0
747 751
748 752 def getTimeRange(self):
749 753
750 754 datatime = []
751 755
752 756 datatime.append(self.ltctime)
753 757 datatime.append(self.ltctime + self.timeInterval)
754 758
755 759 datatime = numpy.array(datatime)
756 760
757 761 return datatime
758 762
759 763 def getChannelIndexList(self):
760 764
761 765 return list(range(self.nChannels))
762 766
763 767 def getNoise(self, type=1):
764 768
765 769
766 770 if type == 1:
767 771 noise = self.getNoisebyHildebrand()
768 772
769 773 if type == 2:
770 774 noise = self.getNoisebySort()
771 775
772 776 if type == 3:
773 777 noise = self.getNoisebyWindow()
774 778
775 779 return noise
776 780
777 781 @property
778 782 def timeInterval(self):
779 783
780 784 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
781 785
782 786 return timeInterval
783 787
784 788 @property
785 789 def ippSeconds(self):
786 790 '''
787 791 '''
788 792 return self.ipp_sec
789 793
790 794 noise = property(getNoise, "I'm the 'nHeights' property.")
791 795
792 796
793 797 class Correlation(JROData):
794 798
795 799 def __init__(self):
796 800 '''
797 801 Constructor
798 802 '''
799 803 self.radarControllerHeaderObj = RadarControllerHeader()
800 804 self.systemHeaderObj = SystemHeader()
801 805 self.type = "Correlation"
802 806 self.data = None
803 807 self.dtype = None
804 808 self.nProfiles = None
805 809 self.heightList = None
806 810 self.channelList = None
807 811 self.flagNoData = True
808 812 self.flagDiscontinuousBlock = False
809 813 self.utctime = None
810 814 self.timeZone = 0
811 815 self.dstFlag = None
812 816 self.errorCount = None
813 817 self.blocksize = None
814 818 self.flagDecodeData = False # asumo q la data no esta decodificada
815 819 self.flagDeflipData = False # asumo q la data no esta sin flip
816 820 self.pairsList = None
817 821 self.nPoints = None
818 822
819 823 def getPairsList(self):
820 824
821 825 return self.pairsList
822 826
823 827 def getNoise(self, mode=2):
824 828
825 829 indR = numpy.where(self.lagR == 0)[0][0]
826 830 indT = numpy.where(self.lagT == 0)[0][0]
827 831
828 832 jspectra0 = self.data_corr[:, :, indR, :]
829 833 jspectra = copy.copy(jspectra0)
830 834
831 835 num_chan = jspectra.shape[0]
832 836 num_hei = jspectra.shape[2]
833 837
834 838 freq_dc = jspectra.shape[1] / 2
835 839 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
836 840
837 841 if ind_vel[0] < 0:
838 842 ind_vel[list(range(0, 1))] = ind_vel[list(
839 843 range(0, 1))] + self.num_prof
840 844
841 845 if mode == 1:
842 846 jspectra[:, freq_dc, :] = (
843 847 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
844 848
845 849 if mode == 2:
846 850
847 851 vel = numpy.array([-2, -1, 1, 2])
848 852 xx = numpy.zeros([4, 4])
849 853
850 854 for fil in range(4):
851 855 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
852 856
853 857 xx_inv = numpy.linalg.inv(xx)
854 858 xx_aux = xx_inv[0, :]
855 859
856 860 for ich in range(num_chan):
857 861 yy = jspectra[ich, ind_vel, :]
858 862 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
859 863
860 864 junkid = jspectra[ich, freq_dc, :] <= 0
861 865 cjunkid = sum(junkid)
862 866
863 867 if cjunkid.any():
864 868 jspectra[ich, freq_dc, junkid.nonzero()] = (
865 869 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
866 870
867 871 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
868 872
869 873 return noise
870 874
871 875 @property
872 876 def timeInterval(self):
873 877
874 878 return self.ippSeconds * self.nCohInt * self.nProfiles
875 879
876 880 def splitFunctions(self):
877 881
878 882 pairsList = self.pairsList
879 883 ccf_pairs = []
880 884 acf_pairs = []
881 885 ccf_ind = []
882 886 acf_ind = []
883 887 for l in range(len(pairsList)):
884 888 chan0 = pairsList[l][0]
885 889 chan1 = pairsList[l][1]
886 890
887 891 # Obteniendo pares de Autocorrelacion
888 892 if chan0 == chan1:
889 893 acf_pairs.append(chan0)
890 894 acf_ind.append(l)
891 895 else:
892 896 ccf_pairs.append(pairsList[l])
893 897 ccf_ind.append(l)
894 898
895 899 data_acf = self.data_cf[acf_ind]
896 900 data_ccf = self.data_cf[ccf_ind]
897 901
898 902 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
899 903
900 904 @property
901 905 def normFactor(self):
902 906 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
903 907 acf_pairs = numpy.array(acf_pairs)
904 908 normFactor = numpy.zeros((self.nPairs, self.nHeights))
905 909
906 910 for p in range(self.nPairs):
907 911 pair = self.pairsList[p]
908 912
909 913 ch0 = pair[0]
910 914 ch1 = pair[1]
911 915
912 916 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
913 917 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
914 918 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
915 919
916 920 return normFactor
917 921
918 922
919 923 class Parameters(Spectra):
920 924
921 925 radarControllerHeaderTxt=None #header Controller like text
922 926 groupList = None # List of Pairs, Groups, etc
923 927 data_param = None # Parameters obtained
924 928 data_pre = None # Data Pre Parametrization
925 929 data_SNR = None # Signal to Noise Ratio
926 930 data_outlier = None
927 931 abscissaList = None # Abscissa, can be velocities, lags or time
928 932 utctimeInit = None # Initial UTC time
929 933 paramInterval = None # Time interval to calculate Parameters in seconds
930 934 useLocalTime = True
931 935 # Fitting
932 936 data_error = None # Error of the estimation
933 937 constants = None
934 938 library = None
935 939 # Output signal
936 940 outputInterval = None # Time interval to calculate output signal in seconds
937 941 data_output = None # Out signal
938 942 nAvg = None
939 943 noise_estimation = None
940 944 GauSPC = None # Fit gaussian SPC
941 945 txPower = None
942
946 flagProfilesByRange = False
947 nProfilesByRange = None
943 948
944 949 def __init__(self):
945 950 '''
946 951 Constructor
947 952 '''
948 953 self.radarControllerHeaderObj = RadarControllerHeader()
949 954 self.systemHeaderObj = SystemHeader()
950 955 self.processingHeaderObj = ProcessingHeader()
951 956 self.type = "Parameters"
952 957 self.timeZone = 0
953 958
954 959 def getTimeRange1(self, interval):
955 960
956 961 datatime = []
957 962
958 963 if self.useLocalTime:
959 964 time1 = self.utctimeInit - self.timeZone * 60
960 965 else:
961 966 time1 = self.utctimeInit
962 967
963 968 datatime.append(time1)
964 969 datatime.append(time1 + interval)
965 970 datatime = numpy.array(datatime)
966 971
967 972 return datatime
968 973
969 974 @property
970 975 def timeInterval(self):
971 976
972 977 if hasattr(self, 'timeInterval1'):
973 978 return self.timeInterval1
974 979 else:
975 980 return self.paramInterval
976 981
977 982 def setValue(self, value):
978 983
979 984 print("This property should not be initialized")
980 985
981 986 return
982 987
983 988 def getNoise(self):
984 989
985 990 return self.spc_noise
986 991
987 992 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
988 993
989 994
990 995 class PlotterData(object):
991 996 '''
992 997 Object to hold data to be plotted
993 998 '''
994 999
995 1000 MAXNUMX = 200
996 1001 MAXNUMY = 200
997 1002
998 1003 def __init__(self, code, exp_code, localtime=True):
999 1004
1000 1005 self.key = code
1001 1006 self.exp_code = exp_code
1002 1007 self.ready = False
1003 1008 self.flagNoData = False
1004 1009 self.localtime = localtime
1005 1010 self.data = {}
1006 1011 self.meta = {}
1007 1012 self.__heights = []
1008 1013
1009 1014 def __str__(self):
1010 1015 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1011 1016 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1012 1017
1013 1018 def __len__(self):
1014 1019 return len(self.data)
1015 1020
1016 1021 def __getitem__(self, key):
1017 1022 if isinstance(key, int):
1018 1023 return self.data[self.times[key]]
1019 1024 elif isinstance(key, str):
1020 1025 ret = numpy.array([self.data[x][key] for x in self.times])
1021 1026 if ret.ndim > 1:
1022 1027 ret = numpy.swapaxes(ret, 0, 1)
1023 1028 return ret
1024 1029
1025 1030 def __contains__(self, key):
1026 1031 return key in self.data[self.min_time]
1027 1032
1028 1033 def setup(self):
1029 1034 '''
1030 1035 Configure object
1031 1036 '''
1032 1037 self.type = ''
1033 1038 self.ready = False
1034 1039 del self.data
1035 1040 self.data = {}
1036 1041 self.__heights = []
1037 1042 self.__all_heights = set()
1038 1043
1039 1044 def shape(self, key):
1040 1045 '''
1041 1046 Get the shape of the one-element data for the given key
1042 1047 '''
1043 1048
1044 1049 if len(self.data[self.min_time][key]):
1045 1050 return self.data[self.min_time][key].shape
1046 1051 return (0,)
1047 1052
1048 1053 def update(self, data, tm, meta={}):
1049 1054 '''
1050 1055 Update data object with new dataOut
1051 1056 '''
1052 1057
1053 1058 self.data[tm] = data
1054 1059
1055 1060 for key, value in meta.items():
1056 1061 setattr(self, key, value)
1057 1062
1058 1063 def normalize_heights(self):
1059 1064 '''
1060 1065 Ensure same-dimension of the data for different heighList
1061 1066 '''
1062 1067
1063 1068 H = numpy.array(list(self.__all_heights))
1064 1069 H.sort()
1065 1070 for key in self.data:
1066 1071 shape = self.shape(key)[:-1] + H.shape
1067 1072 for tm, obj in list(self.data[key].items()):
1068 1073 h = self.__heights[self.times.tolist().index(tm)]
1069 1074 if H.size == h.size:
1070 1075 continue
1071 1076 index = numpy.where(numpy.in1d(H, h))[0]
1072 1077 dummy = numpy.zeros(shape) + numpy.nan
1073 1078 if len(shape) == 2:
1074 1079 dummy[:, index] = obj
1075 1080 else:
1076 1081 dummy[index] = obj
1077 1082 self.data[key][tm] = dummy
1078 1083
1079 1084 self.__heights = [H for tm in self.times]
1080 1085
1081 1086 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1082 1087 '''
1083 1088 Convert data to json
1084 1089 '''
1085 1090
1086 1091 meta = {}
1087 1092 meta['xrange'] = []
1088 1093 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1089 1094 tmp = self.data[tm][self.key]
1090 1095 shape = tmp.shape
1091 1096 if len(shape) == 2:
1092 1097 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1093 1098 elif len(shape) == 3:
1094 1099 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1095 1100 data = self.roundFloats(
1096 1101 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1097 1102 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1098 1103 else:
1099 1104 data = self.roundFloats(self.data[tm][self.key].tolist())
1100 1105
1101 1106 ret = {
1102 1107 'plot': plot_name,
1103 1108 'code': self.exp_code,
1104 1109 'time': float(tm),
1105 1110 'data': data,
1106 1111 }
1107 1112 meta['type'] = plot_type
1108 1113 meta['interval'] = float(self.interval)
1109 1114 meta['localtime'] = self.localtime
1110 1115 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1111 1116 meta.update(self.meta)
1112 1117 ret['metadata'] = meta
1113 1118 return json.dumps(ret)
1114 1119
1115 1120 @property
1116 1121 def times(self):
1117 1122 '''
1118 1123 Return the list of times of the current data
1119 1124 '''
1120 1125
1121 1126 ret = [t for t in self.data]
1122 1127 ret.sort()
1123 1128 return numpy.array(ret)
1124 1129
1125 1130 @property
1126 1131 def min_time(self):
1127 1132 '''
1128 1133 Return the minimun time value
1129 1134 '''
1130 1135
1131 1136 return self.times[0]
1132 1137
1133 1138 @property
1134 1139 def max_time(self):
1135 1140 '''
1136 1141 Return the maximun time value
1137 1142 '''
1138 1143
1139 1144 return self.times[-1]
1140 1145
1141 1146 # @property
1142 1147 # def heights(self):
1143 1148 # '''
1144 1149 # Return the list of heights of the current data
1145 1150 # '''
1146 1151
1147 1152 # return numpy.array(self.__heights[-1])
1148 1153
1149 1154 @staticmethod
1150 1155 def roundFloats(obj):
1151 1156 if isinstance(obj, list):
1152 1157 return list(map(PlotterData.roundFloats, obj))
1153 1158 elif isinstance(obj, float):
1154 1159 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,1464 +1,1465
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 from matplotlib.ticker import LinearLocator
15 15
16 16 from matplotlib import __version__ as plt_version
17 17
18 18 if plt_version >='3.3.4':
19 19 EXTRA_POINTS = 0
20 20 else:
21 21 EXTRA_POINTS = 1
22 22
23 23 class SpectraPlot(Plot):
24 24 '''
25 25 Plot for Spectra data
26 26 '''
27 27
28 28 CODE = 'spc'
29 29 colormap = 'jet'
30 30 plot_type = 'pcolor'
31 31 buffering = False
32 32 channelList = []
33 33 elevationList = []
34 34 azimuthList = []
35 35
36 36 def setup(self):
37 37
38 38 self.nplots = len(self.data.channels)
39 39 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
40 40 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
41 41 self.height = 3.4 * self.nrows
42 42
43 43 self.cb_label = 'dB'
44 44 if self.showprofile:
45 45 self.width = 5.2 * self.ncols
46 46 else:
47 47 self.width = 4.2* self.ncols
48 48 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.12})
49 49 self.ylabel = 'Range [km]'
50 50
51 51
52 52 def update_list(self,dataOut):
53 53 if len(self.channelList) == 0:
54 54 self.channelList = dataOut.channelList
55 55 if len(self.elevationList) == 0:
56 56 self.elevationList = dataOut.elevationList
57 57 if len(self.azimuthList) == 0:
58 58 self.azimuthList = dataOut.azimuthList
59 59
60 60 def update(self, dataOut):
61 61
62 62 self.update_list(dataOut)
63 63 data = {}
64 64 meta = {}
65 65
66 66 #data['rti'] = dataOut.getPower()
67 67 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
68 68 noise = 10*numpy.log10(dataOut.getNoise()/norm)
69 69
70 70
71 71 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
72 72 for ch in range(dataOut.nChannels):
73 73 if hasattr(dataOut.normFactor,'ndim'):
74 74 if dataOut.normFactor.ndim > 1:
75 75 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
76 76
77 77 else:
78 78 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
79 79 else:
80 80 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
81 81
82 82
83 83 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
84 84 spc = 10*numpy.log10(z)
85 85
86 86 data['spc'] = spc
87 87 #print(spc[0].shape)
88 88 data['rti'] = spc.mean(axis=1)
89 89 data['noise'] = noise
90 90 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
91 91 if self.CODE == 'spc_moments':
92 92 data['moments'] = dataOut.moments
93 93
94 94 return data, meta
95 95
96 96 def plot(self):
97 97 if self.xaxis == "frequency":
98 98 x = self.data.xrange[0]
99 99 self.xlabel = "Frequency (kHz)"
100 100 elif self.xaxis == "time":
101 101 x = self.data.xrange[1]
102 102 self.xlabel = "Time (ms)"
103 103 else:
104 104 x = self.data.xrange[2]
105 105 self.xlabel = "Velocity (m/s)"
106 106
107 107 if self.CODE == 'spc_moments':
108 108 x = self.data.xrange[2]
109 109 self.xlabel = "Velocity (m/s)"
110 110
111 111 self.titles = []
112 112 y = self.data.yrange
113 113 self.y = y
114 114
115 115 data = self.data[-1]
116 116 z = data['spc']
117 117 #print(z.shape, x.shape, y.shape)
118 118 for n, ax in enumerate(self.axes):
119 119 noise = self.data['noise'][n][0]
120 120 #print(noise)
121 121 if self.CODE == 'spc_moments':
122 122 mean = data['moments'][n, 1]
123 123 if ax.firsttime:
124 124 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
125 125 self.xmin = self.xmin if self.xmin else -self.xmax
126 126 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
127 127 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
128 128 ax.plt = ax.pcolormesh(x, y, z[n].T,
129 129 vmin=self.zmin,
130 130 vmax=self.zmax,
131 131 cmap=plt.get_cmap(self.colormap)
132 132 )
133 133
134 134 if self.showprofile:
135 135 ax.plt_profile = self.pf_axes[n].plot(
136 136 data['rti'][n], y)[0]
137 137 # ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
138 138 # color="k", linestyle="dashed", lw=1)[0]
139 139 if self.CODE == 'spc_moments':
140 140 ax.plt_mean = ax.plot(mean, y, color='k')[0]
141 141 else:
142 142 ax.plt.set_array(z[n].T.ravel())
143 143 if self.showprofile:
144 144 ax.plt_profile.set_data(data['rti'][n], y)
145 145 #ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
146 146 if self.CODE == 'spc_moments':
147 147 ax.plt_mean.set_data(mean, y)
148 148 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
149 149 self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
150 150 else:
151 151 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
152 152
153 153
154 154 class CrossSpectraPlot(Plot):
155 155
156 156 CODE = 'cspc'
157 157 colormap = 'jet'
158 158 plot_type = 'pcolor'
159 159 zmin_coh = None
160 160 zmax_coh = None
161 161 zmin_phase = None
162 162 zmax_phase = None
163 163 realChannels = None
164 164 crossPairs = None
165 165
166 166 def setup(self):
167 167
168 168 self.ncols = 4
169 169 self.nplots = len(self.data.pairs) * 2
170 170 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
171 171 self.width = 3.1 * self.ncols
172 172 self.height = 2.6 * self.nrows
173 173 self.ylabel = 'Range [km]'
174 174 self.showprofile = False
175 175 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
176 176
177 177 def update(self, dataOut):
178 178
179 179 data = {}
180 180 meta = {}
181 181
182 182 spc = dataOut.data_spc
183 183 cspc = dataOut.data_cspc
184 184 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
185 185 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
186 186 meta['pairs'] = rawPairs
187 187
188 188 if self.crossPairs == None:
189 189 self.crossPairs = dataOut.pairsList
190 190
191 191 tmp = []
192 192
193 193 for n, pair in enumerate(meta['pairs']):
194 194
195 195 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
196 196 coh = numpy.abs(out)
197 197 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
198 198 tmp.append(coh)
199 199 tmp.append(phase)
200 200
201 201 data['cspc'] = numpy.array(tmp)
202 202
203 203 return data, meta
204 204
205 205 def plot(self):
206 206
207 207 if self.xaxis == "frequency":
208 208 x = self.data.xrange[0]
209 209 self.xlabel = "Frequency (kHz)"
210 210 elif self.xaxis == "time":
211 211 x = self.data.xrange[1]
212 212 self.xlabel = "Time (ms)"
213 213 else:
214 214 x = self.data.xrange[2]
215 215 self.xlabel = "Velocity (m/s)"
216 216
217 217 self.titles = []
218 218
219 219 y = self.data.yrange
220 220 self.y = y
221 221
222 222 data = self.data[-1]
223 223 cspc = data['cspc']
224 224
225 225 for n in range(len(self.data.pairs)):
226 226
227 227 pair = self.crossPairs[n]
228 228
229 229 coh = cspc[n*2]
230 230 phase = cspc[n*2+1]
231 231 ax = self.axes[2 * n]
232 232
233 233 if ax.firsttime:
234 234 ax.plt = ax.pcolormesh(x, y, coh.T,
235 235 vmin=self.zmin_coh,
236 236 vmax=self.zmax_coh,
237 237 cmap=plt.get_cmap(self.colormap_coh)
238 238 )
239 239 else:
240 240 ax.plt.set_array(coh.T.ravel())
241 241 self.titles.append(
242 242 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
243 243
244 244 ax = self.axes[2 * n + 1]
245 245 if ax.firsttime:
246 246 ax.plt = ax.pcolormesh(x, y, phase.T,
247 247 vmin=-180,
248 248 vmax=180,
249 249 cmap=plt.get_cmap(self.colormap_phase)
250 250 )
251 251 else:
252 252 ax.plt.set_array(phase.T.ravel())
253 253
254 254 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
255 255
256 256
257 257 class RTIPlot(Plot):
258 258 '''
259 259 Plot for RTI data
260 260 '''
261 261
262 262 CODE = 'rti'
263 263 colormap = 'jet'
264 264 plot_type = 'pcolorbuffer'
265 265 titles = None
266 266 channelList = []
267 267 elevationList = []
268 268 azimuthList = []
269 269
270 270 def setup(self):
271 271 self.xaxis = 'time'
272 272 self.ncols = 1
273 273 #print("dataChannels ",self.data.channels)
274 274 self.nrows = len(self.data.channels)
275 275 self.nplots = len(self.data.channels)
276 276 self.ylabel = 'Range [km]'
277 277 self.xlabel = 'Time'
278 278 self.cb_label = 'dB'
279 279 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
280 280 self.titles = ['{} Channel {}'.format(
281 281 self.CODE.upper(), x) for x in range(self.nplots)]
282 282
283 283 def update_list(self,dataOut):
284 284
285 285 if len(self.channelList) == 0:
286 286 self.channelList = dataOut.channelList
287 287 if len(self.elevationList) == 0:
288 288 self.elevationList = dataOut.elevationList
289 289 if len(self.azimuthList) == 0:
290 290 self.azimuthList = dataOut.azimuthList
291 291
292 292
293 293 def update(self, dataOut):
294 294 if len(self.channelList) == 0:
295 295 self.update_list(dataOut)
296 296 data = {}
297 297 meta = {}
298 298
299 299 data['rti'] = dataOut.getPower()
300 300
301 301 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
302 302 noise = 10*numpy.log10(dataOut.getNoise()/norm)
303 303 data['noise'] = noise
304 304
305 305 return data, meta
306 306
307 307 def plot(self):
308 308
309 309 self.x = self.data.times
310 310 self.y = self.data.yrange
311 311 #print(" x, y: ",self.x, self.y)
312 312 self.z = self.data[self.CODE]
313 313 self.z = numpy.array(self.z, dtype=float)
314 314 self.z = numpy.ma.masked_invalid(self.z)
315 315
316 316 try:
317 317 if self.channelList != None:
318 318 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
319 319 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
320 320 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
321 321 else:
322 322 self.titles = ['{} Channel {}'.format(
323 323 self.CODE.upper(), x) for x in self.channelList]
324 324 except:
325 325 if self.channelList.any() != None:
326 326 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
327 327 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
328 328 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
329 329 else:
330 330 self.titles = ['{} Channel {}'.format(
331 331 self.CODE.upper(), x) for x in self.channelList]
332 332
333 333 if self.decimation is None:
334 334 x, y, z = self.fill_gaps(self.x, self.y, self.z)
335 335 else:
336 336 x, y, z = self.fill_gaps(*self.decimate())
337 337
338 338 #dummy_var = self.axes #Extrañamente esto actualiza el valor axes
339 339 for n, ax in enumerate(self.axes):
340 340 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
341 341 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
342 342 data = self.data[-1]
343 343
344 344 if ax.firsttime:
345 345 ax.plt = ax.pcolormesh(x, y, z[n].T,
346 346 vmin=self.zmin,
347 347 vmax=self.zmax,
348 348 cmap=plt.get_cmap(self.colormap)
349 349 )
350 350 if self.showprofile:
351 351 ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0]
352 352 if "noise" in self.data:
353 353
354 354 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
355 355 color="k", linestyle="dashed", lw=1)[0]
356 356 else:
357 357 ax.collections.remove(ax.collections[0])
358 358 ax.plt = ax.pcolormesh(x, y, z[n].T,
359 359 vmin=self.zmin,
360 360 vmax=self.zmax,
361 361 cmap=plt.get_cmap(self.colormap)
362 362 )
363 363 if self.showprofile:
364 364 ax.plot_profile.set_data(data[self.CODE][n], self.y)
365 365 if "noise" in self.data:
366 366 ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
367 367
368 368
369 369 class CoherencePlot(RTIPlot):
370 370 '''
371 371 Plot for Coherence data
372 372 '''
373 373
374 374 CODE = 'coh'
375 375 titles = None
376 376
377 377 def setup(self):
378 378 self.xaxis = 'time'
379 379 self.ncols = 1
380 380 self.nrows = len(self.data.pairs)
381 381 self.nplots = len(self.data.pairs)
382 382 self.ylabel = 'Range [km]'
383 383 self.xlabel = 'Time'
384 384 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
385 385 if self.CODE == 'coh':
386 386 self.cb_label = ''
387 387 self.titles = [
388 388 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
389 389 else:
390 390 self.cb_label = 'Degrees'
391 391 self.titles = [
392 392 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
393 393
394 394
395 395 def update(self, dataOut):
396 396
397 397 data = {}
398 398 meta = {}
399 399 data['coh'] = dataOut.getCoherence()
400 400 meta['pairs'] = dataOut.pairsList
401 401
402 402
403 403 return data, meta
404 404
405 405 def plot(self):
406 406
407 407 self.x = self.data.times
408 408 self.y = self.data.yrange
409 409 self.z = self.data[self.CODE]
410 410
411 411 self.z = numpy.ma.masked_invalid(self.z)
412 412
413 413 if self.decimation is None:
414 414 x, y, z = self.fill_gaps(self.x, self.y, self.z)
415 415 else:
416 416 x, y, z = self.fill_gaps(*self.decimate())
417 417
418 418 for n, ax in enumerate(self.axes):
419 419 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
420 420 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
421 421 if ax.firsttime:
422 422 ax.plt = ax.pcolormesh(x, y, z[n].T,
423 423 vmin=self.zmin,
424 424 vmax=self.zmax,
425 425 cmap=plt.get_cmap(self.colormap)
426 426 )
427 427 if self.showprofile:
428 428 ax.plot_profile = self.pf_axes[n].plot(
429 429 self.data[self.CODE][n][-1], self.y)[0]
430 430 # ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
431 431 # color="k", linestyle="dashed", lw=1)[0]
432 432 else:
433 433 ax.collections.remove(ax.collections[0])
434 434 ax.plt = ax.pcolormesh(x, y, z[n].T,
435 435 vmin=self.zmin,
436 436 vmax=self.zmax,
437 437 cmap=plt.get_cmap(self.colormap)
438 438 )
439 439 if self.showprofile:
440 440 ax.plot_profile.set_data(self.data[self.CODE][n][-1], self.y)
441 441 # ax.plot_noise.set_data(numpy.repeat(
442 442 # self.data['noise'][n][-1], len(self.y)), self.y)
443 443
444 444
445 445
446 446 class PhasePlot(CoherencePlot):
447 447 '''
448 448 Plot for Phase map data
449 449 '''
450 450
451 451 CODE = 'phase'
452 452 colormap = 'seismic'
453 453
454 454 def update(self, dataOut):
455 455
456 456 data = {}
457 457 meta = {}
458 458 data['phase'] = dataOut.getCoherence(phase=True)
459 459 meta['pairs'] = dataOut.pairsList
460 460
461 461 return data, meta
462 462
463 463 class NoisePlot(Plot):
464 464 '''
465 465 Plot for noise
466 466 '''
467 467
468 468 CODE = 'noise'
469 469 plot_type = 'scatterbuffer'
470 470
471 471 def setup(self):
472 472 self.xaxis = 'time'
473 473 self.ncols = 1
474 474 self.nrows = 1
475 475 self.nplots = 1
476 476 self.ylabel = 'Intensity [dB]'
477 477 self.xlabel = 'Time'
478 478 self.titles = ['Noise']
479 479 self.colorbar = False
480 480 self.plots_adjust.update({'right': 0.85 })
481 481 #if not self.titles:
482 482 self.titles = ['Noise Plot']
483 483
484 484 def update(self, dataOut):
485 485
486 486 data = {}
487 487 meta = {}
488 488 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
489 489 noise = 10*numpy.log10(dataOut.getNoise())
490 490 noise = noise.reshape(dataOut.nChannels, 1)
491 491 data['noise'] = noise
492 492 meta['yrange'] = numpy.array([])
493 493
494 494 return data, meta
495 495
496 496 def plot(self):
497 497
498 498 x = self.data.times
499 499 xmin = self.data.min_time
500 500 xmax = xmin + self.xrange * 60 * 60
501 501 Y = self.data['noise']
502 502
503 503 if self.axes[0].firsttime:
504 504 if self.ymin is None: self.ymin = numpy.nanmin(Y) - 5
505 505 if self.ymax is None: self.ymax = numpy.nanmax(Y) + 5
506 506 for ch in self.data.channels:
507 507 y = Y[ch]
508 508 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
509 509 plt.legend(bbox_to_anchor=(1.18, 1.0))
510 510 else:
511 511 for ch in self.data.channels:
512 512 y = Y[ch]
513 513 self.axes[0].lines[ch].set_data(x, y)
514 514
515 515
516 516 class PowerProfilePlot(Plot):
517 517
518 518 CODE = 'pow_profile'
519 519 plot_type = 'scatter'
520 520
521 521 def setup(self):
522 522
523 523 self.ncols = 1
524 524 self.nrows = 1
525 525 self.nplots = 1
526 526 self.height = 4
527 527 self.width = 3
528 528 self.ylabel = 'Range [km]'
529 529 self.xlabel = 'Intensity [dB]'
530 530 self.titles = ['Power Profile']
531 531 self.colorbar = False
532 532
533 533 def update(self, dataOut):
534 534
535 535 data = {}
536 536 meta = {}
537 537 data[self.CODE] = dataOut.getPower()
538 538
539 539 return data, meta
540 540
541 541 def plot(self):
542 542
543 543 y = self.data.yrange
544 544 self.y = y
545 545
546 546 x = self.data[-1][self.CODE]
547 547
548 548 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
549 549 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
550 550
551 551 if self.axes[0].firsttime:
552 552 for ch in self.data.channels:
553 553 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
554 554 plt.legend()
555 555 else:
556 556 for ch in self.data.channels:
557 557 self.axes[0].lines[ch].set_data(x[ch], y)
558 558
559 559
560 560 class SpectraCutPlot(Plot):
561 561
562 562 CODE = 'spc_cut'
563 563 plot_type = 'scatter'
564 564 buffering = False
565 565 heights = []
566 566 channelList = []
567 567 maintitle = "Spectra Cuts"
568 568 flag_setIndex = False
569 569
570 570 def setup(self):
571 571
572 572 self.nplots = len(self.data.channels)
573 573 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
574 574 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
575 575 self.width = 4.5 * self.ncols + 2.5
576 576 self.height = 4.8 * self.nrows
577 577 self.ylabel = 'Power [dB]'
578 578 self.colorbar = False
579 579 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
580 580
581 581 if len(self.selectedHeightsList) > 0:
582 582 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
583 583
584 584
585 585
586 586 def update(self, dataOut):
587 587 if len(self.channelList) == 0:
588 588 self.channelList = dataOut.channelList
589 589
590 590 self.heights = dataOut.heightList
591 591 #print("sels: ",self.selectedHeightsList)
592 592 if len(self.selectedHeightsList)>0 and not self.flag_setIndex:
593 593
594 594 for sel_height in self.selectedHeightsList:
595 595 index_list = numpy.where(self.heights >= sel_height)
596 596 index_list = index_list[0]
597 597 self.height_index.append(index_list[0])
598 598 #print("sels i:"", self.height_index)
599 599 self.flag_setIndex = True
600 600 #print(self.height_index)
601 601 data = {}
602 602 meta = {}
603 603
604 604 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
605 605 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
606 606 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
607 607
608 608
609 609 z = []
610 610 for ch in range(dataOut.nChannels):
611 611 if hasattr(dataOut.normFactor,'shape'):
612 612 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
613 613 else:
614 614 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
615 615
616 616 z = numpy.asarray(z)
617 617 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
618 618 spc = 10*numpy.log10(z)
619 619
620 620
621 621 data['spc'] = spc - noise
622 622 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
623 623
624 624 return data, meta
625 625
626 626 def plot(self):
627 627 if self.xaxis == "frequency":
628 628 x = self.data.xrange[0][1:]
629 629 self.xlabel = "Frequency (kHz)"
630 630 elif self.xaxis == "time":
631 631 x = self.data.xrange[1]
632 632 self.xlabel = "Time (ms)"
633 633 else:
634 634 x = self.data.xrange[2]
635 635 self.xlabel = "Velocity (m/s)"
636 636
637 637 self.titles = []
638 638
639 639 y = self.data.yrange
640 640 z = self.data[-1]['spc']
641 641 #print(z.shape)
642 642 if len(self.height_index) > 0:
643 643 index = self.height_index
644 644 else:
645 645 index = numpy.arange(0, len(y), int((len(y))/9))
646 646 #print("inde x ", index, self.axes)
647 647
648 648 for n, ax in enumerate(self.axes):
649 649
650 650 if ax.firsttime:
651 651
652 652
653 653 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
654 654 self.xmin = self.xmin if self.xmin else -self.xmax
655 655 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
656 656 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
657 657
658 658
659 659 ax.plt = ax.plot(x, z[n, :, index].T)
660 660 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
661 661 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
662 662 ax.minorticks_on()
663 663 ax.grid(which='major', axis='both')
664 664 ax.grid(which='minor', axis='x')
665 665 else:
666 666 for i, line in enumerate(ax.plt):
667 667 line.set_data(x, z[n, :, index[i]])
668 668
669 669
670 670 self.titles.append('CH {}'.format(self.channelList[n]))
671 671 plt.suptitle(self.maintitle, fontsize=10)
672 672
673 673
674 674 class BeaconPhase(Plot):
675 675
676 676 __isConfig = None
677 677 __nsubplots = None
678 678
679 679 PREFIX = 'beacon_phase'
680 680
681 681 def __init__(self):
682 682 Plot.__init__(self)
683 683 self.timerange = 24*60*60
684 684 self.isConfig = False
685 685 self.__nsubplots = 1
686 686 self.counter_imagwr = 0
687 687 self.WIDTH = 800
688 688 self.HEIGHT = 400
689 689 self.WIDTHPROF = 120
690 690 self.HEIGHTPROF = 0
691 691 self.xdata = None
692 692 self.ydata = None
693 693
694 694 self.PLOT_CODE = BEACON_CODE
695 695
696 696 self.FTP_WEI = None
697 697 self.EXP_CODE = None
698 698 self.SUB_EXP_CODE = None
699 699 self.PLOT_POS = None
700 700
701 701 self.filename_phase = None
702 702
703 703 self.figfile = None
704 704
705 705 self.xmin = None
706 706 self.xmax = None
707 707
708 708 def getSubplots(self):
709 709
710 710 ncol = 1
711 711 nrow = 1
712 712
713 713 return nrow, ncol
714 714
715 715 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
716 716
717 717 self.__showprofile = showprofile
718 718 self.nplots = nplots
719 719
720 720 ncolspan = 7
721 721 colspan = 6
722 722 self.__nsubplots = 2
723 723
724 724 self.createFigure(id = id,
725 725 wintitle = wintitle,
726 726 widthplot = self.WIDTH+self.WIDTHPROF,
727 727 heightplot = self.HEIGHT+self.HEIGHTPROF,
728 728 show=show)
729 729
730 730 nrow, ncol = self.getSubplots()
731 731
732 732 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
733 733
734 734 def save_phase(self, filename_phase):
735 735 f = open(filename_phase,'w+')
736 736 f.write('\n\n')
737 737 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
738 738 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
739 739 f.close()
740 740
741 741 def save_data(self, filename_phase, data, data_datetime):
742 742 f=open(filename_phase,'a')
743 743 timetuple_data = data_datetime.timetuple()
744 744 day = str(timetuple_data.tm_mday)
745 745 month = str(timetuple_data.tm_mon)
746 746 year = str(timetuple_data.tm_year)
747 747 hour = str(timetuple_data.tm_hour)
748 748 minute = str(timetuple_data.tm_min)
749 749 second = str(timetuple_data.tm_sec)
750 750 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
751 751 f.close()
752 752
753 753 def plot(self):
754 754 log.warning('TODO: Not yet implemented...')
755 755
756 756 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
757 757 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
758 758 timerange=None,
759 759 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
760 760 server=None, folder=None, username=None, password=None,
761 761 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
762 762
763 763 if dataOut.flagNoData:
764 764 return dataOut
765 765
766 766 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
767 767 return
768 768
769 769 if pairsList == None:
770 770 pairsIndexList = dataOut.pairsIndexList[:10]
771 771 else:
772 772 pairsIndexList = []
773 773 for pair in pairsList:
774 774 if pair not in dataOut.pairsList:
775 775 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
776 776 pairsIndexList.append(dataOut.pairsList.index(pair))
777 777
778 778 if pairsIndexList == []:
779 779 return
780 780
781 781 # if len(pairsIndexList) > 4:
782 782 # pairsIndexList = pairsIndexList[0:4]
783 783
784 784 hmin_index = None
785 785 hmax_index = None
786 786
787 787 if hmin != None and hmax != None:
788 788 indexes = numpy.arange(dataOut.nHeights)
789 789 hmin_list = indexes[dataOut.heightList >= hmin]
790 790 hmax_list = indexes[dataOut.heightList <= hmax]
791 791
792 792 if hmin_list.any():
793 793 hmin_index = hmin_list[0]
794 794
795 795 if hmax_list.any():
796 796 hmax_index = hmax_list[-1]+1
797 797
798 798 x = dataOut.getTimeRange()
799 799
800 800 thisDatetime = dataOut.datatime
801 801
802 802 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
803 803 xlabel = "Local Time"
804 804 ylabel = "Phase (degrees)"
805 805
806 806 update_figfile = False
807 807
808 808 nplots = len(pairsIndexList)
809 809 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
810 810 phase_beacon = numpy.zeros(len(pairsIndexList))
811 811 for i in range(nplots):
812 812 pair = dataOut.pairsList[pairsIndexList[i]]
813 813 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
814 814 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
815 815 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
816 816 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
817 817 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
818 818
819 819 if dataOut.beacon_heiIndexList:
820 820 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
821 821 else:
822 822 phase_beacon[i] = numpy.average(phase)
823 823
824 824 if not self.isConfig:
825 825
826 826 nplots = len(pairsIndexList)
827 827
828 828 self.setup(id=id,
829 829 nplots=nplots,
830 830 wintitle=wintitle,
831 831 showprofile=showprofile,
832 832 show=show)
833 833
834 834 if timerange != None:
835 835 self.timerange = timerange
836 836
837 837 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
838 838
839 839 if ymin == None: ymin = 0
840 840 if ymax == None: ymax = 360
841 841
842 842 self.FTP_WEI = ftp_wei
843 843 self.EXP_CODE = exp_code
844 844 self.SUB_EXP_CODE = sub_exp_code
845 845 self.PLOT_POS = plot_pos
846 846
847 847 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
848 848 self.isConfig = True
849 849 self.figfile = figfile
850 850 self.xdata = numpy.array([])
851 851 self.ydata = numpy.array([])
852 852
853 853 update_figfile = True
854 854
855 855 #open file beacon phase
856 856 path = '%s%03d' %(self.PREFIX, self.id)
857 857 beacon_file = os.path.join(path,'%s.txt'%self.name)
858 858 self.filename_phase = os.path.join(figpath,beacon_file)
859 859 #self.save_phase(self.filename_phase)
860 860
861 861
862 862 #store data beacon phase
863 863 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
864 864
865 865 self.setWinTitle(title)
866 866
867 867
868 868 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
869 869
870 870 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
871 871
872 872 axes = self.axesList[0]
873 873
874 874 self.xdata = numpy.hstack((self.xdata, x[0:1]))
875 875
876 876 if len(self.ydata)==0:
877 877 self.ydata = phase_beacon.reshape(-1,1)
878 878 else:
879 879 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
880 880
881 881
882 882 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
883 883 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
884 884 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
885 885 XAxisAsTime=True, grid='both'
886 886 )
887 887
888 888 self.draw()
889 889
890 890 if dataOut.ltctime >= self.xmax:
891 891 self.counter_imagwr = wr_period
892 892 self.isConfig = False
893 893 update_figfile = True
894 894
895 895 self.save(figpath=figpath,
896 896 figfile=figfile,
897 897 save=save,
898 898 ftp=ftp,
899 899 wr_period=wr_period,
900 900 thisDatetime=thisDatetime,
901 901 update_figfile=update_figfile)
902 902
903 903 return dataOut
904 904
905 905 class NoiselessSpectraPlot(Plot):
906 906 '''
907 907 Plot for Spectra data, subtracting
908 908 the noise in all channels, using for
909 909 amisr-14 data
910 910 '''
911 911
912 912 CODE = 'noiseless_spc'
913 913 colormap = 'jet'
914 914 plot_type = 'pcolor'
915 915 buffering = False
916 916 channelList = []
917 917 last_noise = None
918 918
919 919 def setup(self):
920 920
921 921 self.nplots = len(self.data.channels)
922 922 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
923 923 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
924 924 self.height = 3.5 * self.nrows
925 925
926 926 self.cb_label = 'dB'
927 927 if self.showprofile:
928 928 self.width = 5.8 * self.ncols
929 929 else:
930 930 self.width = 4.8* self.ncols
931 931 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12})
932 932
933 933 self.ylabel = 'Range [km]'
934 934
935 935
936 936 def update_list(self,dataOut):
937 937 if len(self.channelList) == 0:
938 938 self.channelList = dataOut.channelList
939 939
940 940 def update(self, dataOut):
941 941
942 942 self.update_list(dataOut)
943 943 data = {}
944 944 meta = {}
945 945
946 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
947 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
946 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
947 n0 = (dataOut.getNoise()/norm)
948 948 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
949 noise = 10*numpy.log10(noise)
949 950
950 z = []
951 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
951 952 for ch in range(dataOut.nChannels):
952 if hasattr(dataOut.normFactor,'shape'):
953 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
953 if hasattr(dataOut.normFactor,'ndim'):
954 if dataOut.normFactor.ndim > 1:
955 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
954 956 else:
955 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
957 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
958 else:
959 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
956 960
957 z = numpy.asarray(z)
958 961 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
959 962 spc = 10*numpy.log10(z)
960 963
961 964
962 965 data['spc'] = spc - noise
963 966 #print(spc.shape)
964 967 data['rti'] = spc.mean(axis=1)
965 968 data['noise'] = noise
966 969
967 970
968 971
969 972 # data['noise'] = noise
970 973 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
971 974
972 975 return data, meta
973 976
974 977 def plot(self):
975 978 if self.xaxis == "frequency":
976 979 x = self.data.xrange[0]
977 980 self.xlabel = "Frequency (kHz)"
978 981 elif self.xaxis == "time":
979 982 x = self.data.xrange[1]
980 983 self.xlabel = "Time (ms)"
981 984 else:
982 985 x = self.data.xrange[2]
983 986 self.xlabel = "Velocity (m/s)"
984 987
985 988 self.titles = []
986 989 y = self.data.yrange
987 990 self.y = y
988 991
989 992 data = self.data[-1]
990 993 z = data['spc']
991 994
992 995 for n, ax in enumerate(self.axes):
993 996 #noise = data['noise'][n]
994 997
995 998 if ax.firsttime:
996 999 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
997 1000 self.xmin = self.xmin if self.xmin else -self.xmax
998 1001 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
999 1002 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1000 1003 ax.plt = ax.pcolormesh(x, y, z[n].T,
1001 1004 vmin=self.zmin,
1002 1005 vmax=self.zmax,
1003 1006 cmap=plt.get_cmap(self.colormap)
1004 1007 )
1005 1008
1006 1009 if self.showprofile:
1007 1010 ax.plt_profile = self.pf_axes[n].plot(
1008 1011 data['rti'][n], y)[0]
1009 1012
1010 1013
1011 1014 else:
1012 1015 ax.plt.set_array(z[n].T.ravel())
1013 1016 if self.showprofile:
1014 1017 ax.plt_profile.set_data(data['rti'][n], y)
1015 1018
1016 1019
1017 1020 self.titles.append('CH {}'.format(self.channelList[n]))
1018 1021
1019 1022
1020 1023 class NoiselessRTIPlot(RTIPlot):
1021 1024 '''
1022 1025 Plot for RTI data
1023 1026 '''
1024 1027
1025 1028 CODE = 'noiseless_rti'
1026 1029 colormap = 'jet'
1027 1030 plot_type = 'pcolorbuffer'
1028 1031 titles = None
1029 1032 channelList = []
1030 1033 elevationList = []
1031 1034 azimuthList = []
1032 1035 last_noise = None
1033 1036
1034 1037 def setup(self):
1035 1038 self.xaxis = 'time'
1036 1039 self.ncols = 1
1037 1040 #print("dataChannels ",self.data.channels)
1038 1041 self.nrows = len(self.data.channels)
1039 1042 self.nplots = len(self.data.channels)
1040 1043 self.ylabel = 'Range [km]'
1041 1044 self.xlabel = 'Time'
1042 1045 self.cb_label = 'dB'
1043 1046 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1044 1047 self.titles = ['{} Channel {}'.format(
1045 1048 self.CODE.upper(), x) for x in range(self.nplots)]
1046 1049
1047 1050 def update_list(self,dataOut):
1048 1051 if len(self.channelList) == 0:
1049 1052 self.channelList = dataOut.channelList
1050 1053 if len(self.elevationList) == 0:
1051 1054 self.elevationList = dataOut.elevationList
1052 1055 if len(self.azimuthList) == 0:
1053 1056 self.azimuthList = dataOut.azimuthList
1054 1057
1055 1058 def update(self, dataOut):
1056 1059 if len(self.channelList) == 0:
1057 1060 self.update_list(dataOut)
1058 1061
1059 1062 data = {}
1060 1063 meta = {}
1061 1064 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1062 1065 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt)
1063 1066 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1064
1065
1066 1067 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1067 1068 data['noise'] = n0
1068 1069 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1069 1070 noiseless_data = dataOut.getPower() - noise
1070 1071
1071 1072 #print("power, noise:", dataOut.getPower(), n0)
1072 1073 #print(noise)
1073 1074 #print(noiseless_data)
1074 1075
1075 1076 data['noiseless_rti'] = noiseless_data
1076 1077
1077 1078 return data, meta
1078 1079
1079 1080 def plot(self):
1080 1081 from matplotlib import pyplot as plt
1081 1082 self.x = self.data.times
1082 1083 self.y = self.data.yrange
1083 1084 self.z = self.data['noiseless_rti']
1084 1085 self.z = numpy.array(self.z, dtype=float)
1085 1086 self.z = numpy.ma.masked_invalid(self.z)
1086 1087
1087 1088
1088 1089 try:
1089 1090 if self.channelList != None:
1090 1091 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1091 1092 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1092 1093 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1093 1094 else:
1094 1095 self.titles = ['{} Channel {}'.format(
1095 1096 self.CODE.upper(), x) for x in self.channelList]
1096 1097 except:
1097 1098 if self.channelList.any() != None:
1098 1099 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1099 1100 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1100 1101 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1101 1102 else:
1102 1103 self.titles = ['{} Channel {}'.format(
1103 1104 self.CODE.upper(), x) for x in self.channelList]
1104 1105
1105 1106
1106 1107 if self.decimation is None:
1107 1108 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1108 1109 else:
1109 1110 x, y, z = self.fill_gaps(*self.decimate())
1110 1111
1111 1112 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
1112 1113 #print("plot shapes ", z.shape, x.shape, y.shape)
1113 1114 #print(self.axes)
1114 1115 for n, ax in enumerate(self.axes):
1115 1116
1116 1117
1117 1118 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
1118 1119 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1119 1120 data = self.data[-1]
1120 1121 if ax.firsttime:
1121 1122 ax.plt = ax.pcolormesh(x, y, z[n].T,
1122 1123 vmin=self.zmin,
1123 1124 vmax=self.zmax,
1124 1125 cmap=plt.get_cmap(self.colormap)
1125 1126 )
1126 1127 if self.showprofile:
1127 1128 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1128 1129
1129 1130 else:
1130 1131 ax.collections.remove(ax.collections[0])
1131 1132 ax.plt = ax.pcolormesh(x, y, z[n].T,
1132 1133 vmin=self.zmin,
1133 1134 vmax=self.zmax,
1134 1135 cmap=plt.get_cmap(self.colormap)
1135 1136 )
1136 1137 if self.showprofile:
1137 1138 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1138 1139 # if "noise" in self.data:
1139 1140 # #ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
1140 1141 # ax.plot_noise.set_data(data['noise'][n], self.y)
1141 1142
1142 1143
1143 1144 class OutliersRTIPlot(Plot):
1144 1145 '''
1145 1146 Plot for data_xxxx object
1146 1147 '''
1147 1148
1148 1149 CODE = 'outlier_rtc' # Range Time Counts
1149 1150 colormap = 'cool'
1150 1151 plot_type = 'pcolorbuffer'
1151 1152
1152 1153 def setup(self):
1153 1154 self.xaxis = 'time'
1154 1155 self.ncols = 1
1155 1156 self.nrows = self.data.shape('outlier_rtc')[0]
1156 1157 self.nplots = self.nrows
1157 1158 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1158 1159
1159 1160
1160 1161 if not self.xlabel:
1161 1162 self.xlabel = 'Time'
1162 1163
1163 1164 self.ylabel = 'Height [km]'
1164 1165 if not self.titles:
1165 1166 self.titles = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)]
1166 1167
1167 1168 def update(self, dataOut):
1168 1169
1169 1170 data = {}
1170 1171 data['outlier_rtc'] = dataOut.data_outlier
1171 1172
1172 1173 meta = {}
1173 1174
1174 1175 return data, meta
1175 1176
1176 1177 def plot(self):
1177 1178 # self.data.normalize_heights()
1178 1179 self.x = self.data.times
1179 1180 self.y = self.data.yrange
1180 1181 self.z = self.data['outlier_rtc']
1181 1182
1182 1183 #self.z = numpy.ma.masked_invalid(self.z)
1183 1184
1184 1185 if self.decimation is None:
1185 1186 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1186 1187 else:
1187 1188 x, y, z = self.fill_gaps(*self.decimate())
1188 1189
1189 1190 for n, ax in enumerate(self.axes):
1190 1191
1191 1192 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1192 1193 self.z[n])
1193 1194 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1194 1195 self.z[n])
1195 1196 data = self.data[-1]
1196 1197 if ax.firsttime:
1197 1198 if self.zlimits is not None:
1198 1199 self.zmin, self.zmax = self.zlimits[n]
1199 1200
1200 1201 ax.plt = ax.pcolormesh(x, y, z[n].T,
1201 1202 vmin=self.zmin,
1202 1203 vmax=self.zmax,
1203 1204 cmap=self.cmaps[n]
1204 1205 )
1205 1206 if self.showprofile:
1206 1207 ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0]
1207 1208 self.pf_axes[n].set_xlabel('')
1208 1209 else:
1209 1210 if self.zlimits is not None:
1210 1211 self.zmin, self.zmax = self.zlimits[n]
1211 1212 ax.collections.remove(ax.collections[0])
1212 1213 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1213 1214 vmin=self.zmin,
1214 1215 vmax=self.zmax,
1215 1216 cmap=self.cmaps[n]
1216 1217 )
1217 1218 if self.showprofile:
1218 1219 ax.plot_profile.set_data(data['outlier_rtc'][n], self.y)
1219 1220 self.pf_axes[n].set_xlabel('')
1220 1221
1221 1222 class NIncohIntRTIPlot(Plot):
1222 1223 '''
1223 1224 Plot for data_xxxx object
1224 1225 '''
1225 1226
1226 1227 CODE = 'integrations_rtc' # Range Time Counts
1227 1228 colormap = 'BuGn'
1228 1229 plot_type = 'pcolorbuffer'
1229 1230
1230 1231 def setup(self):
1231 1232 self.xaxis = 'time'
1232 1233 self.ncols = 1
1233 1234 self.nrows = self.data.shape('integrations_rtc')[0]
1234 1235 self.nplots = self.nrows
1235 1236 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1236 1237
1237 1238
1238 1239 if not self.xlabel:
1239 1240 self.xlabel = 'Time'
1240 1241
1241 1242 self.ylabel = 'Height [km]'
1242 1243 if not self.titles:
1243 1244 self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)]
1244 1245
1245 1246 def update(self, dataOut):
1246 1247
1247 1248 data = {}
1248 1249 data['integrations_rtc'] = dataOut.nIncohInt
1249 1250
1250 1251 meta = {}
1251 1252
1252 1253 return data, meta
1253 1254
1254 1255 def plot(self):
1255 1256 # self.data.normalize_heights()
1256 1257 self.x = self.data.times
1257 1258 self.y = self.data.yrange
1258 1259 self.z = self.data['integrations_rtc']
1259 1260
1260 1261 #self.z = numpy.ma.masked_invalid(self.z)
1261 1262
1262 1263 if self.decimation is None:
1263 1264 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1264 1265 else:
1265 1266 x, y, z = self.fill_gaps(*self.decimate())
1266 1267
1267 1268 for n, ax in enumerate(self.axes):
1268 1269
1269 1270 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1270 1271 self.z[n])
1271 1272 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1272 1273 self.z[n])
1273 1274 data = self.data[-1]
1274 1275 if ax.firsttime:
1275 1276 if self.zlimits is not None:
1276 1277 self.zmin, self.zmax = self.zlimits[n]
1277 1278
1278 1279 ax.plt = ax.pcolormesh(x, y, z[n].T,
1279 1280 vmin=self.zmin,
1280 1281 vmax=self.zmax,
1281 1282 cmap=self.cmaps[n]
1282 1283 )
1283 1284 if self.showprofile:
1284 1285 ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0]
1285 1286 self.pf_axes[n].set_xlabel('')
1286 1287 else:
1287 1288 if self.zlimits is not None:
1288 1289 self.zmin, self.zmax = self.zlimits[n]
1289 1290 ax.collections.remove(ax.collections[0])
1290 1291 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1291 1292 vmin=self.zmin,
1292 1293 vmax=self.zmax,
1293 1294 cmap=self.cmaps[n]
1294 1295 )
1295 1296 if self.showprofile:
1296 1297 ax.plot_profile.set_data(data['integrations_rtc'][n], self.y)
1297 1298 self.pf_axes[n].set_xlabel('')
1298 1299
1299 1300
1300 1301 import datetime
1301 1302 class NoiselessRTILinePlot(Plot):
1302 1303 '''
1303 1304 Plot for RTI data
1304 1305 '''
1305 1306
1306 1307 CODE = 'noiseless_rtiLine'
1307 1308
1308 1309 plot_type = 'scatter'
1309 1310 titles = None
1310 1311 channelList = []
1311 1312 elevationList = []
1312 1313 azimuthList = []
1313 1314 last_noise = None
1314 1315
1315 1316 def setup(self):
1316 1317 self.xaxis = 'Range (Km)'
1317 1318 self.nplots = len(self.data.channels)
1318 1319 self.nrows = int(numpy.ceil(self.nplots/2))
1319 1320 self.ncols = int(numpy.ceil(self.nplots/self.nrows))
1320 1321 self.ylabel = 'Intensity [dB]'
1321 1322 self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels]
1322 1323 self.colorbar = False
1323 1324 self.width = 6
1324 1325 self.height = 4
1325 1326
1326 1327 def update_list(self,dataOut):
1327 1328 if len(self.channelList) == 0:
1328 1329 self.channelList = dataOut.channelList
1329 1330 if len(self.elevationList) == 0:
1330 1331 self.elevationList = dataOut.elevationList
1331 1332 if len(self.azimuthList) == 0:
1332 1333 self.azimuthList = dataOut.azimuthList
1333 1334
1334 1335 def update(self, dataOut):
1335 1336 if len(self.channelList) == 0:
1336 1337 self.update_list(dataOut)
1337 1338
1338 1339 data = {}
1339 1340 meta = {}
1340 1341 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1341 1342 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt)
1342 1343 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1343 1344
1344 1345
1345 1346 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1346 1347 data['noise'] = n0
1347 1348
1348 1349 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1349 1350 noiseless_data = dataOut.getPower() - noise
1350 1351
1351 1352 #print("power, noise:", dataOut.getPower(), n0)
1352 1353 #print(noise)
1353 1354 #print(noiseless_data)
1354 1355
1355 1356 data['noiseless_rtiLine'] = noiseless_data
1356 1357
1357 1358 print(noiseless_data.shape, self.name)
1358 1359 data['time'] = dataOut.utctime
1359 1360
1360 1361 return data, meta
1361 1362
1362 1363 def plot(self):
1363 1364
1364 1365 self.x = self.data.times
1365 1366 self.y = self.data.yrange
1366 1367 print(self.data['noiseless_rtiLine'].shape, self.y.shape, self.name)
1367 1368 #ts = self.data['time'][0].squeeze()
1368 1369 if len(self.data['noiseless_rtiLine'])>2 :
1369 1370 self.z = self.data['noiseless_rtiLine'][:, -1,:]
1370 1371 else:
1371 1372 self.z = self.data['noiseless_rtiLine']
1372 1373 #print(self.z.shape, self.y.shape, ts)
1373 1374 #thisDatetime = datetime.datetime.utcfromtimestamp(ts)
1374 1375
1375 1376 for i,ax in enumerate(self.axes):
1376 1377 #self.titles[i] = "Channel {} {}".format(i, thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1377 1378
1378 1379
1379 1380 if ax.firsttime:
1380 1381 #self.xmin = min(self.z)
1381 1382 #self.xmax = max(self.z)
1382 1383 ax.plt_r = ax.plot(self.z[i], self.y)[0]
1383 1384 else:
1384 1385 ax.plt_r.set_data(self.z[i], self.y)
1385 1386
1386 1387
1387 1388
1388 1389 class GeneralProfilePlot(Plot):
1389 1390 '''
1390 1391 Plot for RTI data
1391 1392 '''
1392 1393
1393 1394 CODE = 'general_profilePlot'
1394 1395
1395 1396 plot_type = 'scatter'
1396 1397 titles = None
1397 1398 channelList = []
1398 1399 elevationList = []
1399 1400 azimuthList = []
1400 1401 last_noise = None
1401 1402
1402 1403 def setup(self):
1403 1404 self.xaxis = 'Range (Km)'
1404 1405 self.nplots = len(self.data.channels)
1405 1406 self.nrows = int(numpy.ceil(self.nplots/2))
1406 1407 self.ncols = int(numpy.ceil(self.nplots/self.nrows))
1407 1408 self.ylabel = 'Intensity [dB]'
1408 1409 self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels]
1409 1410 self.colorbar = False
1410 1411 self.width = 6
1411 1412 self.height = 4
1412 1413
1413 1414 def update_list(self,dataOut):
1414 1415 if len(self.channelList) == 0:
1415 1416 self.channelList = dataOut.channelList
1416 1417 if len(self.elevationList) == 0:
1417 1418 self.elevationList = dataOut.elevationList
1418 1419 if len(self.azimuthList) == 0:
1419 1420 self.azimuthList = dataOut.azimuthList
1420 1421
1421 1422 def update(self, dataOut):
1422 1423 if len(self.channelList) == 0:
1423 1424 self.update_list(dataOut)
1424 1425
1425 1426 data = {}
1426 1427 meta = {}
1427 1428
1428 1429 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1429 1430 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1430 1431 data['noise'] = n0
1431 1432
1432 1433 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1433 1434 noiseless_data = dataOut.getPower() - noise
1434 1435
1435 1436 data['noiseless_rtiLine'] = noiseless_data
1436 1437
1437 1438 print(noiseless_data.shape, self.name)
1438 1439 data['time'] = dataOut.utctime
1439 1440
1440 1441 return data, meta
1441 1442
1442 1443 def plot(self):
1443 1444
1444 1445 self.x = self.data.times
1445 1446 self.y = self.data.yrange
1446 1447 print(self.data['noiseless_rtiLine'].shape, self.y.shape, self.name)
1447 1448 #ts = self.data['time'][0].squeeze()
1448 1449 if len(self.data['noiseless_rtiLine'])>2 :
1449 1450 self.z = self.data['noiseless_rtiLine'][:, -1,:]
1450 1451 else:
1451 1452 self.z = self.data['noiseless_rtiLine']
1452 1453 #print(self.z.shape, self.y.shape, ts)
1453 1454 #thisDatetime = datetime.datetime.utcfromtimestamp(ts)
1454 1455
1455 1456 for i,ax in enumerate(self.axes):
1456 1457 #self.titles[i] = "Channel {} {}".format(i, thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1457 1458
1458 1459
1459 1460 if ax.firsttime:
1460 1461 #self.xmin = min(self.z)
1461 1462 #self.xmax = max(self.z)
1462 1463 ax.plt_r = ax.plot(self.z[i], self.y)[0]
1463 1464 else:
1464 1465 ax.plt_r.set_data(self.z[i], self.y) No newline at end of file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,2256 +1,2287
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.model.data import _noise
21 21
22 22 from schainpy.utils import log
23 23 import matplotlib.pyplot as plt
24 24 #from scipy.optimize import curve_fit
25 25 from schainpy.model.io.utilsIO import getHei_index
26 26 import datetime
27 27
28 28 class SpectraProc(ProcessingUnit):
29 29
30 30 def __init__(self):
31 31
32 32 ProcessingUnit.__init__(self)
33 33
34 34 self.buffer = None
35 35 self.firstdatatime = None
36 36 self.profIndex = 0
37 37 self.dataOut = Spectra()
38 38 self.id_min = None
39 39 self.id_max = None
40 40 self.setupReq = False #Agregar a todas las unidades de proc
41 41 self.nsamplesFFT = 0
42 42
43 43 def __updateSpecFromVoltage(self):
44 44
45 45
46 46
47 47 self.dataOut.timeZone = self.dataIn.timeZone
48 48 self.dataOut.dstFlag = self.dataIn.dstFlag
49 49 self.dataOut.errorCount = self.dataIn.errorCount
50 50 self.dataOut.useLocalTime = self.dataIn.useLocalTime
51 51
52 52 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
53 53 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
54 54 self.dataOut.ippSeconds = self.dataIn.ippSeconds
55 55 self.dataOut.ipp = self.dataIn.ipp
56 56 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
57 57 self.dataOut.channelList = self.dataIn.channelList
58 58 self.dataOut.heightList = self.dataIn.heightList
59 59 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
60 60 self.dataOut.nProfiles = self.dataOut.nFFTPoints
61 61 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
62 62 self.dataOut.utctime = self.firstdatatime
63 63 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
64 64 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
65 65 self.dataOut.flagShiftFFT = False
66 66 self.dataOut.nCohInt = self.dataIn.nCohInt
67 67 self.dataOut.nIncohInt = 1
68 68 self.dataOut.deltaHeight = self.dataIn.deltaHeight
69 69 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
70 70 self.dataOut.frequency = self.dataIn.frequency
71 71 self.dataOut.realtime = self.dataIn.realtime
72 72 self.dataOut.azimuth = self.dataIn.azimuth
73 73 self.dataOut.zenith = self.dataIn.zenith
74 74 self.dataOut.codeList = self.dataIn.codeList
75 75 self.dataOut.azimuthList = self.dataIn.azimuthList
76 76 self.dataOut.elevationList = self.dataIn.elevationList
77 77 self.dataOut.code = self.dataIn.code
78 78 self.dataOut.nCode = self.dataIn.nCode
79 self.dataOut.flagProfilesByRange = self.dataIn.flagProfilesByRange
80 self.dataOut.nProfilesByRange = self.dataIn.nProfilesByRange
79 81
80 82
81 83 def __getFft(self):
82 84 # print("fft donw")
83 85 """
84 86 Convierte valores de Voltaje a Spectra
85 87
86 88 Affected:
87 89 self.dataOut.data_spc
88 90 self.dataOut.data_cspc
89 91 self.dataOut.data_dc
90 92 self.dataOut.heightList
91 93 self.profIndex
92 94 self.buffer
93 95 self.dataOut.flagNoData
94 96 """
95 97 fft_volt = numpy.fft.fft(
96 98 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
97 99 fft_volt = fft_volt.astype(numpy.dtype('complex'))
98 100 dc = fft_volt[:, 0, :]
99 101
100 102 # calculo de self-spectra
101 103 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
102 104 spc = fft_volt * numpy.conjugate(fft_volt)
103 105 spc = spc.real
104 106
105 107 blocksize = 0
106 108 blocksize += dc.size
107 109 blocksize += spc.size
108 110
109 111 cspc = None
110 112 pairIndex = 0
111 113 if self.dataOut.pairsList != None:
112 114 # calculo de cross-spectra
113 115 cspc = numpy.zeros(
114 116 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
115 117 for pair in self.dataOut.pairsList:
116 118 if pair[0] not in self.dataOut.channelList:
117 119 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
118 120 str(pair), str(self.dataOut.channelList)))
119 121 if pair[1] not in self.dataOut.channelList:
120 122 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
121 123 str(pair), str(self.dataOut.channelList)))
122 124
123 125 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
124 126 numpy.conjugate(fft_volt[pair[1], :, :])
125 127 pairIndex += 1
126 128 blocksize += cspc.size
127 129
128 130 self.dataOut.data_spc = spc
129 131 self.dataOut.data_cspc = cspc
130 132 self.dataOut.data_dc = dc
131 133 self.dataOut.blockSize = blocksize
132 134 self.dataOut.flagShiftFFT = False
133 135
134 136 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, zeroPad=False):
135 137
136 138
137 139 try:
138 140 type = self.dataIn.type.decode("utf-8")
139 141 self.dataIn.type = type
140 142 except:
141 143 pass
142 144 if self.dataIn.type == "Spectra":
143 145
144 146 try:
145 147 self.dataOut.copy(self.dataIn)
146 148 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
147 149 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
148 150 self.dataOut.nProfiles = self.dataOut.nFFTPoints
149 151 #self.dataOut.nHeights = len(self.dataOut.heightList)
150 152 except Exception as e:
151 153 print("Error dataIn ",e)
152 154
153 155
154 156
155 157 if shift_fft:
156 158 #desplaza a la derecha en el eje 2 determinadas posiciones
157 159 shift = int(self.dataOut.nFFTPoints/2)
158 160 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
159 161
160 162 if self.dataOut.data_cspc is not None:
161 163 #desplaza a la derecha en el eje 2 determinadas posiciones
162 164 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
163 165 if pairsList:
164 166 self.__selectPairs(pairsList)
165 167
166 168
167 169 elif self.dataIn.type == "Voltage":
168 170
169 171 self.dataOut.flagNoData = True
170 172 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
171 173 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
172 174 if nFFTPoints == None:
173 175 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
174 176
175 177 if nProfiles == None:
176 178 nProfiles = nFFTPoints
177 179
178 180 if ippFactor == None:
179 181 self.dataOut.ippFactor = 1
180 182
181 183 self.dataOut.nFFTPoints = nFFTPoints
182 184 #print(" volts ch,prof, h: ", self.dataIn.data.shape)
183 185 if self.buffer is None:
184 186 if not zeroPad:
185 187 self.buffer = numpy.zeros((self.dataIn.nChannels,
186 188 nProfiles,
187 189 self.dataIn.nHeights),
188 190 dtype='complex')
189 191 else:
190 192 self.buffer = numpy.zeros((self.dataIn.nChannels,
191 193 nFFTPoints,
192 194 self.dataIn.nHeights),
193 195 dtype='complex')
194 196
195 197 if self.dataIn.flagDataAsBlock:
196 198 nVoltProfiles = self.dataIn.data.shape[1]
197 199
198 200 if nVoltProfiles == nProfiles or zeroPad:
199 201 self.buffer = self.dataIn.data.copy()
200 202 self.profIndex = nVoltProfiles
201 203
202 204 elif nVoltProfiles < nProfiles:
203 205
204 206 if self.profIndex == 0:
205 207 self.id_min = 0
206 208 self.id_max = nVoltProfiles
207 209
208 210 self.buffer[:, self.id_min:self.id_max,
209 211 :] = self.dataIn.data
210 212 self.profIndex += nVoltProfiles
211 213 self.id_min += nVoltProfiles
212 214 self.id_max += nVoltProfiles
213 215 else:
214 216 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
215 217 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
216 218 self.dataOut.flagNoData = True
217 219 else:
218 220 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
219 221 self.profIndex += 1
220 222
221 223 if self.firstdatatime == None:
222 224 self.firstdatatime = self.dataIn.utctime
223 225
224 226 if self.profIndex == nProfiles or zeroPad:
225 227
226 228 self.__updateSpecFromVoltage()
227 229
228 230 if pairsList == None:
229 231 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
230 232 else:
231 233 self.dataOut.pairsList = pairsList
232 234 self.__getFft()
233 235 self.dataOut.flagNoData = False
234 236 self.firstdatatime = None
235 237 self.nsamplesFFT = self.profIndex
236 238 self.profIndex = 0
237 239
238 240 #update Processing Header:
239 241 self.dataOut.processingHeaderObj.dtype = "Spectra"
240 242 self.dataOut.processingHeaderObj.nFFTPoints = self.dataOut.nFFTPoints
241 243 self.dataOut.processingHeaderObj.nSamplesFFT = self.nsamplesFFT
242 244 self.dataOut.processingHeaderObj.nIncohInt = 1
243 245
244 246
245 247 elif self.dataIn.type == "Parameters":
246 248
247 249 self.dataOut.data_spc = self.dataIn.data_spc
248 250 self.dataOut.data_cspc = self.dataIn.data_cspc
249 251 self.dataOut.data_outlier = self.dataIn.data_outlier
250 252 self.dataOut.nProfiles = self.dataIn.nProfiles
251 253 self.dataOut.nIncohInt = self.dataIn.nIncohInt
252 254 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
253 255 self.dataOut.ippFactor = self.dataIn.ippFactor
254 256 self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
255 257 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
256 258 self.dataOut.ProcessingHeader = self.dataIn.ProcessingHeader.copy()
257 259 self.dataOut.ippSeconds = self.dataIn.ippSeconds
258 260 self.dataOut.ipp = self.dataIn.ipp
259 261 #self.dataOut.abscissaList = self.dataIn.getVelRange(1)
260 262 #self.dataOut.spc_noise = self.dataIn.getNoise()
261 263 #self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
262 264 # self.dataOut.normFactor = self.dataIn.normFactor
263 265 if hasattr(self.dataIn, 'channelList'):
264 266 self.dataOut.channelList = self.dataIn.channelList
265 267 if hasattr(self.dataIn, 'pairsList'):
266 268 self.dataOut.pairsList = self.dataIn.pairsList
267 269 self.dataOut.groupList = self.dataIn.pairsList
268 270
269 271 self.dataOut.flagNoData = False
270 272
271 273 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
272 274 self.dataOut.ChanDist = self.dataIn.ChanDist
273 275 else: self.dataOut.ChanDist = None
274 276
275 277 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
276 278 # self.dataOut.VelRange = self.dataIn.VelRange
277 279 #else: self.dataOut.VelRange = None
278 280
279 281
280 282
281 283 else:
282 284 raise ValueError("The type of input object {} is not valid".format(
283 285 self.dataIn.type))
284 286
285 287
286 288
287 289
288 290 #print("spc proc Done", self.dataOut.data_spc.shape)
289 291 #print(self.dataOut.data_spc)
290 292 return
291 293
292 294 def __selectPairs(self, pairsList):
293 295
294 296 if not pairsList:
295 297 return
296 298
297 299 pairs = []
298 300 pairsIndex = []
299 301
300 302 for pair in pairsList:
301 303 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
302 304 continue
303 305 pairs.append(pair)
304 306 pairsIndex.append(pairs.index(pair))
305 307
306 308 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
307 309 self.dataOut.pairsList = pairs
308 310
309 311 return
310 312
311 313 def selectFFTs(self, minFFT, maxFFT ):
312 314 """
313 315 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
314 316 minFFT<= FFT <= maxFFT
315 317 """
316 318
317 319 if (minFFT > maxFFT):
318 320 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
319 321
320 322 if (minFFT < self.dataOut.getFreqRange()[0]):
321 323 minFFT = self.dataOut.getFreqRange()[0]
322 324
323 325 if (maxFFT > self.dataOut.getFreqRange()[-1]):
324 326 maxFFT = self.dataOut.getFreqRange()[-1]
325 327
326 328 minIndex = 0
327 329 maxIndex = 0
328 330 FFTs = self.dataOut.getFreqRange()
329 331
330 332 inda = numpy.where(FFTs >= minFFT)
331 333 indb = numpy.where(FFTs <= maxFFT)
332 334
333 335 try:
334 336 minIndex = inda[0][0]
335 337 except:
336 338 minIndex = 0
337 339
338 340 try:
339 341 maxIndex = indb[0][-1]
340 342 except:
341 343 maxIndex = len(FFTs)
342 344
343 345 self.selectFFTsByIndex(minIndex, maxIndex)
344 346
345 347 return 1
346 348
347 349 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
348 350 newheis = numpy.where(
349 351 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
350 352
351 353 if hei_ref != None:
352 354 newheis = numpy.where(self.dataOut.heightList > hei_ref)
353 355
354 356 minIndex = min(newheis[0])
355 357 maxIndex = max(newheis[0])
356 358 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
357 359 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
358 360
359 361 # determina indices
360 362 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
361 363 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
362 364 avg_dB = 10 * \
363 365 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
364 366 beacon_dB = numpy.sort(avg_dB)[-nheis:]
365 367 beacon_heiIndexList = []
366 368 for val in avg_dB.tolist():
367 369 if val >= beacon_dB[0]:
368 370 beacon_heiIndexList.append(avg_dB.tolist().index(val))
369 371
370 372 #data_spc = data_spc[:,:,beacon_heiIndexList]
371 373 data_cspc = None
372 374 if self.dataOut.data_cspc is not None:
373 375 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
374 376 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
375 377
376 378 data_dc = None
377 379 if self.dataOut.data_dc is not None:
378 380 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
379 381 #data_dc = data_dc[:,beacon_heiIndexList]
380 382
381 383 self.dataOut.data_spc = data_spc
382 384 self.dataOut.data_cspc = data_cspc
383 385 self.dataOut.data_dc = data_dc
384 386 self.dataOut.heightList = heightList
385 387 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
386 388
387 389 return 1
388 390
389 391 def selectFFTsByIndex(self, minIndex, maxIndex):
390 392 """
391 393
392 394 """
393 395
394 396 if (minIndex < 0) or (minIndex > maxIndex):
395 397 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
396 398
397 399 if (maxIndex >= self.dataOut.nProfiles):
398 400 maxIndex = self.dataOut.nProfiles-1
399 401
400 402 #Spectra
401 403 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
402 404
403 405 data_cspc = None
404 406 if self.dataOut.data_cspc is not None:
405 407 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
406 408
407 409 data_dc = None
408 410 if self.dataOut.data_dc is not None:
409 411 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
410 412
411 413 self.dataOut.data_spc = data_spc
412 414 self.dataOut.data_cspc = data_cspc
413 415 self.dataOut.data_dc = data_dc
414 416
415 417 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
416 418 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
417 419 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
418 420
419 421 return 1
420 422
421 423 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
422 424 # validacion de rango
423 425 if minHei == None:
424 426 minHei = self.dataOut.heightList[0]
425 427
426 428 if maxHei == None:
427 429 maxHei = self.dataOut.heightList[-1]
428 430
429 431 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
430 432 print('minHei: %.2f is out of the heights range' % (minHei))
431 433 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
432 434 minHei = self.dataOut.heightList[0]
433 435
434 436 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
435 437 print('maxHei: %.2f is out of the heights range' % (maxHei))
436 438 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
437 439 maxHei = self.dataOut.heightList[-1]
438 440
439 441 # validacion de velocidades
440 442 velrange = self.dataOut.getVelRange(1)
441 443
442 444 if minVel == None:
443 445 minVel = velrange[0]
444 446
445 447 if maxVel == None:
446 448 maxVel = velrange[-1]
447 449
448 450 if (minVel < velrange[0]) or (minVel > maxVel):
449 451 print('minVel: %.2f is out of the velocity range' % (minVel))
450 452 print('minVel is setting to %.2f' % (velrange[0]))
451 453 minVel = velrange[0]
452 454
453 455 if (maxVel > velrange[-1]) or (maxVel < minVel):
454 456 print('maxVel: %.2f is out of the velocity range' % (maxVel))
455 457 print('maxVel is setting to %.2f' % (velrange[-1]))
456 458 maxVel = velrange[-1]
457 459
458 460 # seleccion de indices para rango
459 461 minIndex = 0
460 462 maxIndex = 0
461 463 heights = self.dataOut.heightList
462 464
463 465 inda = numpy.where(heights >= minHei)
464 466 indb = numpy.where(heights <= maxHei)
465 467
466 468 try:
467 469 minIndex = inda[0][0]
468 470 except:
469 471 minIndex = 0
470 472
471 473 try:
472 474 maxIndex = indb[0][-1]
473 475 except:
474 476 maxIndex = len(heights)
475 477
476 478 if (minIndex < 0) or (minIndex > maxIndex):
477 479 raise ValueError("some value in (%d,%d) is not valid" % (
478 480 minIndex, maxIndex))
479 481
480 482 if (maxIndex >= self.dataOut.nHeights):
481 483 maxIndex = self.dataOut.nHeights - 1
482 484
483 485 # seleccion de indices para velocidades
484 486 indminvel = numpy.where(velrange >= minVel)
485 487 indmaxvel = numpy.where(velrange <= maxVel)
486 488 try:
487 489 minIndexVel = indminvel[0][0]
488 490 except:
489 491 minIndexVel = 0
490 492
491 493 try:
492 494 maxIndexVel = indmaxvel[0][-1]
493 495 except:
494 496 maxIndexVel = len(velrange)
495 497
496 498 # seleccion del espectro
497 499 data_spc = self.dataOut.data_spc[:,
498 500 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
499 501 # estimacion de ruido
500 502 noise = numpy.zeros(self.dataOut.nChannels)
501 503
502 504 for channel in range(self.dataOut.nChannels):
503 505 daux = data_spc[channel, :, :]
504 506 sortdata = numpy.sort(daux, axis=None)
505 507 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
506 508
507 509 self.dataOut.noise_estimation = noise.copy()
508 510
509 511 return 1
510 512
511 513 class removeDC(Operation):
512 514
513 515 def run(self, dataOut, mode=2):
514 516 self.dataOut = dataOut
515 517 jspectra = self.dataOut.data_spc
516 518 jcspectra = self.dataOut.data_cspc
517 519
518 520 num_chan = jspectra.shape[0]
519 521 num_hei = jspectra.shape[2]
520 522
521 523 if jcspectra is not None:
522 524 jcspectraExist = True
523 525 num_pairs = jcspectra.shape[0]
524 526 else:
525 527 jcspectraExist = False
526 528
527 529 freq_dc = int(jspectra.shape[1] / 2)
528 530 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
529 531 ind_vel = ind_vel.astype(int)
530 532
531 533 if ind_vel[0] < 0:
532 534 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
533 535
534 536 if mode == 1:
535 537 jspectra[:, freq_dc, :] = (
536 538 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
537 539
538 540 if jcspectraExist:
539 541 jcspectra[:, freq_dc, :] = (
540 542 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
541 543
542 544 if mode == 2:
543 545
544 546 vel = numpy.array([-2, -1, 1, 2])
545 547 xx = numpy.zeros([4, 4])
546 548
547 549 for fil in range(4):
548 550 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
549 551
550 552 xx_inv = numpy.linalg.inv(xx)
551 553 xx_aux = xx_inv[0, :]
552 554
553 555 for ich in range(num_chan):
554 556 yy = jspectra[ich, ind_vel, :]
555 557 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
556 558
557 559 junkid = jspectra[ich, freq_dc, :] <= 0
558 560 cjunkid = sum(junkid)
559 561
560 562 if cjunkid.any():
561 563 jspectra[ich, freq_dc, junkid.nonzero()] = (
562 564 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
563 565
564 566 if jcspectraExist:
565 567 for ip in range(num_pairs):
566 568 yy = jcspectra[ip, ind_vel, :]
567 569 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
568 570
569 571 self.dataOut.data_spc = jspectra
570 572 self.dataOut.data_cspc = jcspectra
571 573
572 574 return self.dataOut
573 575
574 576 class getNoiseB(Operation):
575 577
576 578 __slots__ =('offset','warnings', 'isConfig', 'minIndex','maxIndex','minIndexFFT','maxIndexFFT')
577 579 def __init__(self):
578 580
579 581 Operation.__init__(self)
580 582 self.isConfig = False
581 583
582 584 def setup(self, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
583 585
584 586 self.warnings = warnings
585 587 if minHei == None:
586 588 minHei = self.dataOut.heightList[0]
587 589
588 590 if maxHei == None:
589 591 maxHei = self.dataOut.heightList[-1]
590 592
591 593 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
592 594 if self.warnings:
593 595 print('minHei: %.2f is out of the heights range' % (minHei))
594 596 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
595 597 minHei = self.dataOut.heightList[0]
596 598
597 599 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
598 600 if self.warnings:
599 601 print('maxHei: %.2f is out of the heights range' % (maxHei))
600 602 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
601 603 maxHei = self.dataOut.heightList[-1]
602 604
603 605
604 606 #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
605 607 minIndexFFT = 0
606 608 maxIndexFFT = 0
607 609 # validacion de velocidades
608 610 indminPoint = None
609 611 indmaxPoint = None
610 612 if self.dataOut.type == 'Spectra':
611 613 if minVel == None and maxVel == None :
612 614
613 615 freqrange = self.dataOut.getFreqRange(1)
614 616
615 617 if minFreq == None:
616 618 minFreq = freqrange[0]
617 619
618 620 if maxFreq == None:
619 621 maxFreq = freqrange[-1]
620 622
621 623 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
622 624 if self.warnings:
623 625 print('minFreq: %.2f is out of the frequency range' % (minFreq))
624 626 print('minFreq is setting to %.2f' % (freqrange[0]))
625 627 minFreq = freqrange[0]
626 628
627 629 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
628 630 if self.warnings:
629 631 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
630 632 print('maxFreq is setting to %.2f' % (freqrange[-1]))
631 633 maxFreq = freqrange[-1]
632 634
633 635 indminPoint = numpy.where(freqrange >= minFreq)
634 636 indmaxPoint = numpy.where(freqrange <= maxFreq)
635 637
636 638 else:
637 639
638 640 velrange = self.dataOut.getVelRange(1)
639 641
640 642 if minVel == None:
641 643 minVel = velrange[0]
642 644
643 645 if maxVel == None:
644 646 maxVel = velrange[-1]
645 647
646 648 if (minVel < velrange[0]) or (minVel > maxVel):
647 649 if self.warnings:
648 650 print('minVel: %.2f is out of the velocity range' % (minVel))
649 651 print('minVel is setting to %.2f' % (velrange[0]))
650 652 minVel = velrange[0]
651 653
652 654 if (maxVel > velrange[-1]) or (maxVel < minVel):
653 655 if self.warnings:
654 656 print('maxVel: %.2f is out of the velocity range' % (maxVel))
655 657 print('maxVel is setting to %.2f' % (velrange[-1]))
656 658 maxVel = velrange[-1]
657 659
658 660 indminPoint = numpy.where(velrange >= minVel)
659 661 indmaxPoint = numpy.where(velrange <= maxVel)
660 662
661 663
662 664 # seleccion de indices para rango
663 665 minIndex = 0
664 666 maxIndex = 0
665 667 heights = self.dataOut.heightList
666 668
667 669 inda = numpy.where(heights >= minHei)
668 670 indb = numpy.where(heights <= maxHei)
669 671
670 672 try:
671 673 minIndex = inda[0][0]
672 674 except:
673 675 minIndex = 0
674 676
675 677 try:
676 678 maxIndex = indb[0][-1]
677 679 except:
678 680 maxIndex = len(heights)
679 681
680 682 if (minIndex < 0) or (minIndex > maxIndex):
681 683 raise ValueError("some value in (%d,%d) is not valid" % (
682 684 minIndex, maxIndex))
683 685
684 686 if (maxIndex >= self.dataOut.nHeights):
685 687 maxIndex = self.dataOut.nHeights - 1
686 688 #############################################################3
687 689 # seleccion de indices para velocidades
688 690 if self.dataOut.type == 'Spectra':
689 691 try:
690 692 minIndexFFT = indminPoint[0][0]
691 693 except:
692 694 minIndexFFT = 0
693 695
694 696 try:
695 697 maxIndexFFT = indmaxPoint[0][-1]
696 698 except:
697 699 maxIndexFFT = len( self.dataOut.getFreqRange(1))
698 700
699 701 self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT
700 702 self.isConfig = True
701 703 self.offset = 1
702 704 if offset!=None:
703 705 self.offset = 10**(offset/10)
704 706 #print("config getNoiseB Done")
705 707
706 708 def run(self, dataOut, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
707 709 self.dataOut = dataOut
708 710
709 711 if not self.isConfig:
710 712 self.setup(offset, minHei, maxHei,minVel, maxVel, minFreq, maxFreq, warnings)
711 713
712 714 self.dataOut.noise_estimation = None
713 715 noise = None
714 716 #print("data type: ",self.dataOut.type, self.dataOut.nIncohInt, self.dataOut.max_nIncohInt)
715 717 if self.dataOut.type == 'Voltage':
716 718 noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
717 719 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
718 720 elif self.dataOut.type == 'Spectra':
719 721 #print(self.dataOut.nChannels, self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.max_nIncohInt, self.dataOut.nIncohInt)
720 722 noise = numpy.zeros( self.dataOut.nChannels)
721 723 norm = 1
722 724
723 725 for channel in range( self.dataOut.nChannels):
724 726 if not hasattr(self.dataOut.nIncohInt,'__len__'):
725 727 norm = 1
726 728 else:
727 729 norm = self.dataOut.max_nIncohInt[channel]/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
728 730 #print("norm nIncoh: ", norm ,self.dataOut.data_spc.shape, self.dataOut.max_nIncohInt)
729 731 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
730 732 daux = numpy.multiply(daux, norm)
731 733 #print("offset: ", self.offset, 10*numpy.log10(self.offset))
732 734 # noise[channel] = self.getNoiseByMean(daux)/self.offset
733 735 #print(daux.shape, daux)
734 736 #noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
735 737 sortdata = numpy.sort(daux, axis=None)
736 738
737 739 noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt[channel])/self.offset
738 740 #print("noise shape", noise[channel], self.name)
739 741
740 742 #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
741 743 else:
742 744 noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
743 745 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
744 746 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
745 747 #print("2: ",self.dataOut.noise_estimation)
746 748 #print(self.dataOut.flagNoData)
747 749 #print("getNoise Done", noise, self.dataOut.nProfiles ,self.dataOut.ippFactor)
748 750 return self.dataOut
749 751
750 752 def getNoiseByMean(self,data):
751 753 #data debe estar ordenado
752 754 data = numpy.mean(data,axis=1)
753 755 sortdata = numpy.sort(data, axis=None)
754 756 #sortID=data.argsort()
755 757 #print(data.shape)
756 758
757 759 pnoise = None
758 760 j = 0
759 761
760 762 mean = numpy.mean(sortdata)
761 763 min = numpy.min(sortdata)
762 764 delta = mean - min
763 765 indexes = numpy.where(sortdata > (mean+delta))[0] #only array of indexes
764 766 #print(len(indexes))
765 767 if len(indexes)==0:
766 768 pnoise = numpy.mean(sortdata)
767 769 else:
768 770 j = indexes[0]
769 771 pnoise = numpy.mean(sortdata[0:j])
770 772
771 773 # from matplotlib import pyplot as plt
772 774 # plt.plot(sortdata)
773 775 # plt.vlines(j,(pnoise-delta),(pnoise+delta), color='r')
774 776 # plt.show()
775 777 #print("noise: ", 10*numpy.log10(pnoise))
776 778 return pnoise
777 779
778 780 def getNoiseByHS(self,data, navg):
779 781 #data debe estar ordenado
780 782 #data = numpy.mean(data,axis=1)
781 783 sortdata = numpy.sort(data, axis=None)
782 784
783 785 lenOfData = len(sortdata)
784 786 nums_min = lenOfData*0.2
785 787
786 788 if nums_min <= 5:
787 789
788 790 nums_min = 5
789 791
790 792 sump = 0.
791 793 sumq = 0.
792 794
793 795 j = 0
794 796 cont = 1
795 797
796 798 while((cont == 1)and(j < lenOfData)):
797 799
798 800 sump += sortdata[j]
799 801 sumq += sortdata[j]**2
800 802 #sumq -= sump**2
801 803 if j > nums_min:
802 804 rtest = float(j)/(j-1) + 1.0/navg
803 805 #if ((sumq*j) > (sump**2)):
804 806 if ((sumq*j) > (rtest*sump**2)):
805 807 j = j - 1
806 808 sump = sump - sortdata[j]
807 809 sumq = sumq - sortdata[j]**2
808 810 cont = 0
809 811
810 812 j += 1
811 813
812 814 lnoise = sump / j
813 815
814 816 return lnoise
815 817
816 818
817 819
818 820 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
819 821 z = (x - a1) / a2
820 822 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
821 823 return y
822 824
823 825
824 826 # class CleanRayleigh(Operation):
825 827 #
826 828 # def __init__(self):
827 829 #
828 830 # Operation.__init__(self)
829 831 # self.i=0
830 832 # self.isConfig = False
831 833 # self.__dataReady = False
832 834 # self.__profIndex = 0
833 835 # self.byTime = False
834 836 # self.byProfiles = False
835 837 #
836 838 # self.bloques = None
837 839 # self.bloque0 = None
838 840 #
839 841 # self.index = 0
840 842 #
841 843 # self.buffer = 0
842 844 # self.buffer2 = 0
843 845 # self.buffer3 = 0
844 846 #
845 847 #
846 848 # def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
847 849 #
848 850 # self.nChannels = dataOut.nChannels
849 851 # self.nProf = dataOut.nProfiles
850 852 # self.nPairs = dataOut.data_cspc.shape[0]
851 853 # self.pairsArray = numpy.array(dataOut.pairsList)
852 854 # self.spectra = dataOut.data_spc
853 855 # self.cspectra = dataOut.data_cspc
854 856 # self.heights = dataOut.heightList #alturas totales
855 857 # self.nHeights = len(self.heights)
856 858 # self.min_hei = min_hei
857 859 # self.max_hei = max_hei
858 860 # if (self.min_hei == None):
859 861 # self.min_hei = 0
860 862 # if (self.max_hei == None):
861 863 # self.max_hei = dataOut.heightList[-1]
862 864 # self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
863 865 # self.heightsClean = self.heights[self.hval] #alturas filtradas
864 866 # self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
865 867 # self.nHeightsClean = len(self.heightsClean)
866 868 # self.channels = dataOut.channelList
867 869 # self.nChan = len(self.channels)
868 870 # self.nIncohInt = dataOut.nIncohInt
869 871 # self.__initime = dataOut.utctime
870 872 # self.maxAltInd = self.hval[-1]+1
871 873 # self.minAltInd = self.hval[0]
872 874 #
873 875 # self.crosspairs = dataOut.pairsList
874 876 # self.nPairs = len(self.crosspairs)
875 877 # self.normFactor = dataOut.normFactor
876 878 # self.nFFTPoints = dataOut.nFFTPoints
877 879 # self.ippSeconds = dataOut.ippSeconds
878 880 # self.currentTime = self.__initime
879 881 # self.pairsArray = numpy.array(dataOut.pairsList)
880 882 # self.factor_stdv = factor_stdv
881 883 #
882 884 # if n != None :
883 885 # self.byProfiles = True
884 886 # self.nIntProfiles = n
885 887 # else:
886 888 # self.__integrationtime = timeInterval
887 889 #
888 890 # self.__dataReady = False
889 891 # self.isConfig = True
890 892 #
891 893 #
892 894 #
893 895 # def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
894 896 # #print("runing cleanRayleigh")
895 897 # if not self.isConfig :
896 898 #
897 899 # self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
898 900 #
899 901 # tini=dataOut.utctime
900 902 #
901 903 # if self.byProfiles:
902 904 # if self.__profIndex == self.nIntProfiles:
903 905 # self.__dataReady = True
904 906 # else:
905 907 # if (tini - self.__initime) >= self.__integrationtime:
906 908 #
907 909 # self.__dataReady = True
908 910 # self.__initime = tini
909 911 #
910 912 # #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
911 913 #
912 914 # if self.__dataReady:
913 915 #
914 916 # self.__profIndex = 0
915 917 # jspc = self.buffer
916 918 # jcspc = self.buffer2
917 919 # #jnoise = self.buffer3
918 920 # self.buffer = dataOut.data_spc
919 921 # self.buffer2 = dataOut.data_cspc
920 922 # #self.buffer3 = dataOut.noise
921 923 # self.currentTime = dataOut.utctime
922 924 # if numpy.any(jspc) :
923 925 # #print( jspc.shape, jcspc.shape)
924 926 # jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
925 927 # try:
926 928 # jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
927 929 # except:
928 930 # print("no cspc")
929 931 # self.__dataReady = False
930 932 # #print( jspc.shape, jcspc.shape)
931 933 # dataOut.flagNoData = False
932 934 # else:
933 935 # dataOut.flagNoData = True
934 936 # self.__dataReady = False
935 937 # return dataOut
936 938 # else:
937 939 # #print( len(self.buffer))
938 940 # if numpy.any(self.buffer):
939 941 # self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
940 942 # try:
941 943 # self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
942 944 # self.buffer3 += dataOut.data_dc
943 945 # except:
944 946 # pass
945 947 # else:
946 948 # self.buffer = dataOut.data_spc
947 949 # self.buffer2 = dataOut.data_cspc
948 950 # self.buffer3 = dataOut.data_dc
949 951 # #print self.index, self.fint
950 952 # #print self.buffer2.shape
951 953 # dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
952 954 # self.__profIndex += 1
953 955 # return dataOut ## NOTE: REV
954 956 #
955 957 #
956 958 # #index = tini.tm_hour*12+tini.tm_min/5
957 959 # '''
958 960 # #REVISAR
959 961 # '''
960 962 # # jspc = jspc/self.nFFTPoints/self.normFactor
961 963 # # jcspc = jcspc/self.nFFTPoints/self.normFactor
962 964 #
963 965 #
964 966 #
965 967 # tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
966 968 # dataOut.data_spc = tmp_spectra
967 969 # dataOut.data_cspc = tmp_cspectra
968 970 #
969 971 # #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
970 972 #
971 973 # dataOut.data_dc = self.buffer3
972 974 # dataOut.nIncohInt *= self.nIntProfiles
973 975 # dataOut.max_nIncohInt = self.nIntProfiles
974 976 # dataOut.utctime = self.currentTime #tiempo promediado
975 977 # #print("Time: ",time.localtime(dataOut.utctime))
976 978 # # dataOut.data_spc = sat_spectra
977 979 # # dataOut.data_cspc = sat_cspectra
978 980 # self.buffer = 0
979 981 # self.buffer2 = 0
980 982 # self.buffer3 = 0
981 983 #
982 984 # return dataOut
983 985 #
984 986 # def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
985 987 # print("OP cleanRayleigh")
986 988 # #import matplotlib.pyplot as plt
987 989 # #for k in range(149):
988 990 # #channelsProcssd = []
989 991 # #channelA_ok = False
990 992 # #rfunc = cspectra.copy() #self.bloques
991 993 # rfunc = spectra.copy()
992 994 # #rfunc = cspectra
993 995 # #val_spc = spectra*0.0 #self.bloque0*0.0
994 996 # #val_cspc = cspectra*0.0 #self.bloques*0.0
995 997 # #in_sat_spectra = spectra.copy() #self.bloque0
996 998 # #in_sat_cspectra = cspectra.copy() #self.bloques
997 999 #
998 1000 #
999 1001 # ###ONLY FOR TEST:
1000 1002 # raxs = math.ceil(math.sqrt(self.nPairs))
1001 1003 # if raxs == 0:
1002 1004 # raxs = 1
1003 1005 # caxs = math.ceil(self.nPairs/raxs)
1004 1006 # if self.nPairs <4:
1005 1007 # raxs = 2
1006 1008 # caxs = 2
1007 1009 # #print(raxs, caxs)
1008 1010 # fft_rev = 14 #nFFT to plot
1009 1011 # hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
1010 1012 # hei_rev = hei_rev[0]
1011 1013 # #print(hei_rev)
1012 1014 #
1013 1015 # #print numpy.absolute(rfunc[:,0,0,14])
1014 1016 #
1015 1017 # gauss_fit, covariance = None, None
1016 1018 # for ih in range(self.minAltInd,self.maxAltInd):
1017 1019 # for ifreq in range(self.nFFTPoints):
1018 1020 # '''
1019 1021 # ###ONLY FOR TEST:
1020 1022 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1021 1023 # fig, axs = plt.subplots(raxs, caxs)
1022 1024 # fig2, axs2 = plt.subplots(raxs, caxs)
1023 1025 # col_ax = 0
1024 1026 # row_ax = 0
1025 1027 # '''
1026 1028 # #print(self.nPairs)
1027 1029 # for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
1028 1030 # # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
1029 1031 # # continue
1030 1032 # # if not self.crosspairs[ii][0] in channelsProcssd:
1031 1033 # # channelA_ok = True
1032 1034 # #print("pair: ",self.crosspairs[ii])
1033 1035 # '''
1034 1036 # ###ONLY FOR TEST:
1035 1037 # if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
1036 1038 # col_ax = 0
1037 1039 # row_ax += 1
1038 1040 # '''
1039 1041 # func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
1040 1042 # #print(func2clean.shape)
1041 1043 # val = (numpy.isfinite(func2clean)==True).nonzero()
1042 1044 #
1043 1045 # if len(val)>0: #limitador
1044 1046 # min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
1045 1047 # if min_val <= -40 :
1046 1048 # min_val = -40
1047 1049 # max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
1048 1050 # if max_val >= 200 :
1049 1051 # max_val = 200
1050 1052 # #print min_val, max_val
1051 1053 # step = 1
1052 1054 # #print("Getting bins and the histogram")
1053 1055 # x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
1054 1056 # y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1055 1057 # #print(len(y_dist),len(binstep[:-1]))
1056 1058 # #print(row_ax,col_ax, " ..")
1057 1059 # #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
1058 1060 # mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
1059 1061 # sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
1060 1062 # parg = [numpy.amax(y_dist),mean,sigma]
1061 1063 #
1062 1064 # newY = None
1063 1065 #
1064 1066 # try :
1065 1067 # gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
1066 1068 # mode = gauss_fit[1]
1067 1069 # stdv = gauss_fit[2]
1068 1070 # #print(" FIT OK",gauss_fit)
1069 1071 # '''
1070 1072 # ###ONLY FOR TEST:
1071 1073 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1072 1074 # newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
1073 1075 # axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1074 1076 # axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1075 1077 # axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1076 1078 # '''
1077 1079 # except:
1078 1080 # mode = mean
1079 1081 # stdv = sigma
1080 1082 # #print("FIT FAIL")
1081 1083 # #continue
1082 1084 #
1083 1085 #
1084 1086 # #print(mode,stdv)
1085 1087 # #Removing echoes greater than mode + std_factor*stdv
1086 1088 # noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
1087 1089 # #noval tiene los indices que se van a remover
1088 1090 # #print("Chan ",ii," novals: ",len(noval[0]))
1089 1091 # if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
1090 1092 # novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
1091 1093 # #print(novall)
1092 1094 # #print(" ",self.pairsArray[ii])
1093 1095 # #cross_pairs = self.pairsArray[ii]
1094 1096 # #Getting coherent echoes which are removed.
1095 1097 # # if len(novall[0]) > 0:
1096 1098 # #
1097 1099 # # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
1098 1100 # # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
1099 1101 # # val_cspc[novall[0],ii,ifreq,ih] = 1
1100 1102 # #print("OUT NOVALL 1")
1101 1103 # try:
1102 1104 # pair = (self.channels[ii],self.channels[ii + 1])
1103 1105 # except:
1104 1106 # pair = (99,99)
1105 1107 # #print("par ", pair)
1106 1108 # if ( pair in self.crosspairs):
1107 1109 # q = self.crosspairs.index(pair)
1108 1110 # #print("está aqui: ", q, (ii,ii + 1))
1109 1111 # new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
1110 1112 # cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
1111 1113 #
1112 1114 # #if channelA_ok:
1113 1115 # #chA = self.channels.index(cross_pairs[0])
1114 1116 # new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
1115 1117 # spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
1116 1118 # #channelA_ok = False
1117 1119 #
1118 1120 # # chB = self.channels.index(cross_pairs[1])
1119 1121 # # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
1120 1122 # # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
1121 1123 # #
1122 1124 # # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
1123 1125 # # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
1124 1126 # '''
1125 1127 # ###ONLY FOR TEST:
1126 1128 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1127 1129 # func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
1128 1130 # y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1129 1131 # axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1130 1132 # axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1131 1133 # axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1132 1134 # '''
1133 1135 # '''
1134 1136 # ###ONLY FOR TEST:
1135 1137 # col_ax += 1 #contador de ploteo columnas
1136 1138 # ##print(col_ax)
1137 1139 # ###ONLY FOR TEST:
1138 1140 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1139 1141 # title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
1140 1142 # title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
1141 1143 # fig.suptitle(title)
1142 1144 # fig2.suptitle(title2)
1143 1145 # plt.show()
1144 1146 # '''
1145 1147 # ##################################################################################################
1146 1148 #
1147 1149 # #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
1148 1150 # out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
1149 1151 # out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
1150 1152 # for ih in range(self.nHeights):
1151 1153 # for ifreq in range(self.nFFTPoints):
1152 1154 # for ich in range(self.nChan):
1153 1155 # tmp = spectra[:,ich,ifreq,ih]
1154 1156 # valid = (numpy.isfinite(tmp[:])==True).nonzero()
1155 1157 #
1156 1158 # if len(valid[0]) >0 :
1157 1159 # out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1158 1160 #
1159 1161 # for icr in range(self.nPairs):
1160 1162 # tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
1161 1163 # valid = (numpy.isfinite(tmp)==True).nonzero()
1162 1164 # if len(valid[0]) > 0:
1163 1165 # out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1164 1166 #
1165 1167 # return out_spectra, out_cspectra
1166 1168 #
1167 1169 # def REM_ISOLATED_POINTS(self,array,rth):
1168 1170 # # import matplotlib.pyplot as plt
1169 1171 # if rth == None :
1170 1172 # rth = 4
1171 1173 # #print("REM ISO")
1172 1174 # num_prof = len(array[0,:,0])
1173 1175 # num_hei = len(array[0,0,:])
1174 1176 # n2d = len(array[:,0,0])
1175 1177 #
1176 1178 # for ii in range(n2d) :
1177 1179 # #print ii,n2d
1178 1180 # tmp = array[ii,:,:]
1179 1181 # #print tmp.shape, array[ii,101,:],array[ii,102,:]
1180 1182 #
1181 1183 # # fig = plt.figure(figsize=(6,5))
1182 1184 # # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1183 1185 # # ax = fig.add_axes([left, bottom, width, height])
1184 1186 # # x = range(num_prof)
1185 1187 # # y = range(num_hei)
1186 1188 # # cp = ax.contour(y,x,tmp)
1187 1189 # # ax.clabel(cp, inline=True,fontsize=10)
1188 1190 # # plt.show()
1189 1191 #
1190 1192 # #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
1191 1193 # tmp = numpy.reshape(tmp,num_prof*num_hei)
1192 1194 # indxs1 = (numpy.isfinite(tmp)==True).nonzero()
1193 1195 # indxs2 = (tmp > 0).nonzero()
1194 1196 #
1195 1197 # indxs1 = (indxs1[0])
1196 1198 # indxs2 = indxs2[0]
1197 1199 # #indxs1 = numpy.array(indxs1[0])
1198 1200 # #indxs2 = numpy.array(indxs2[0])
1199 1201 # indxs = None
1200 1202 # #print indxs1 , indxs2
1201 1203 # for iv in range(len(indxs2)):
1202 1204 # indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
1203 1205 # #print len(indxs2), indv
1204 1206 # if len(indv[0]) > 0 :
1205 1207 # indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
1206 1208 # # print indxs
1207 1209 # indxs = indxs[1:]
1208 1210 # #print(indxs, len(indxs))
1209 1211 # if len(indxs) < 4 :
1210 1212 # array[ii,:,:] = 0.
1211 1213 # return
1212 1214 #
1213 1215 # xpos = numpy.mod(indxs ,num_hei)
1214 1216 # ypos = (indxs / num_hei)
1215 1217 # sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
1216 1218 # #print sx
1217 1219 # xpos = xpos[sx]
1218 1220 # ypos = ypos[sx]
1219 1221 #
1220 1222 # # *********************************** Cleaning isolated points **********************************
1221 1223 # ic = 0
1222 1224 # while True :
1223 1225 # r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
1224 1226 # #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
1225 1227 # #plt.plot(r)
1226 1228 # #plt.show()
1227 1229 # no_coh1 = (numpy.isfinite(r)==True).nonzero()
1228 1230 # no_coh2 = (r <= rth).nonzero()
1229 1231 # #print r, no_coh1, no_coh2
1230 1232 # no_coh1 = numpy.array(no_coh1[0])
1231 1233 # no_coh2 = numpy.array(no_coh2[0])
1232 1234 # no_coh = None
1233 1235 # #print valid1 , valid2
1234 1236 # for iv in range(len(no_coh2)):
1235 1237 # indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
1236 1238 # if len(indv[0]) > 0 :
1237 1239 # no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
1238 1240 # no_coh = no_coh[1:]
1239 1241 # #print len(no_coh), no_coh
1240 1242 # if len(no_coh) < 4 :
1241 1243 # #print xpos[ic], ypos[ic], ic
1242 1244 # # plt.plot(r)
1243 1245 # # plt.show()
1244 1246 # xpos[ic] = numpy.nan
1245 1247 # ypos[ic] = numpy.nan
1246 1248 #
1247 1249 # ic = ic + 1
1248 1250 # if (ic == len(indxs)) :
1249 1251 # break
1250 1252 # #print( xpos, ypos)
1251 1253 #
1252 1254 # indxs = (numpy.isfinite(list(xpos))==True).nonzero()
1253 1255 # #print indxs[0]
1254 1256 # if len(indxs[0]) < 4 :
1255 1257 # array[ii,:,:] = 0.
1256 1258 # return
1257 1259 #
1258 1260 # xpos = xpos[indxs[0]]
1259 1261 # ypos = ypos[indxs[0]]
1260 1262 # for i in range(0,len(ypos)):
1261 1263 # ypos[i]=int(ypos[i])
1262 1264 # junk = tmp
1263 1265 # tmp = junk*0.0
1264 1266 #
1265 1267 # tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
1266 1268 # array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
1267 1269 #
1268 1270 # #print array.shape
1269 1271 # #tmp = numpy.reshape(tmp,(num_prof,num_hei))
1270 1272 # #print tmp.shape
1271 1273 #
1272 1274 # # fig = plt.figure(figsize=(6,5))
1273 1275 # # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1274 1276 # # ax = fig.add_axes([left, bottom, width, height])
1275 1277 # # x = range(num_prof)
1276 1278 # # y = range(num_hei)
1277 1279 # # cp = ax.contour(y,x,array[ii,:,:])
1278 1280 # # ax.clabel(cp, inline=True,fontsize=10)
1279 1281 # # plt.show()
1280 1282 # return array
1281 1283 #
1282 1284
1283 1285 class IntegrationFaradaySpectra(Operation):
1284 1286
1285 1287 __profIndex = 0
1286 1288 __withOverapping = False
1287 1289
1288 1290 __byTime = False
1289 1291 __initime = None
1290 1292 __lastdatatime = None
1291 1293 __integrationtime = None
1292 1294
1293 1295 __buffer_spc = None
1294 1296 __buffer_cspc = None
1295 1297 __buffer_dc = None
1296 1298
1297 1299 __dataReady = False
1298 1300
1299 1301 __timeInterval = None
1300 1302 n_ints = None #matriz de numero de integracions (CH,HEI)
1301 1303 n = None
1302 1304 minHei_ind = None
1303 1305 maxHei_ind = None
1304 1306 navg = 1.0
1305 1307 factor = 0.0
1306 1308 dataoutliers = None # (CHANNELS, HEIGHTS)
1307 1309
1310 _flagProfilesByRange = False
1311 _nProfilesByRange = 0
1312
1308 1313 def __init__(self):
1309 1314
1310 1315 Operation.__init__(self)
1311 1316
1312 1317 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, avg=1,factor=0.75):
1313 1318 """
1314 1319 Set the parameters of the integration class.
1315 1320
1316 1321 Inputs:
1317 1322
1318 1323 n : Number of coherent integrations
1319 1324 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1320 1325 overlapping :
1321 1326
1322 1327 """
1323 1328
1324 1329 self.__initime = None
1325 1330 self.__lastdatatime = 0
1326 1331
1327 1332 self.__buffer_spc = []
1328 1333 self.__buffer_cspc = []
1329 1334 self.__buffer_dc = 0
1330 1335
1331 1336 self.__profIndex = 0
1332 1337 self.__dataReady = False
1333 1338 self.__byTime = False
1334 1339
1335 1340 self.factor = factor
1336 1341 self.navg = avg
1337 1342 #self.ByLags = dataOut.ByLags ###REDEFINIR
1338 1343 self.ByLags = False
1339 1344 self.maxProfilesInt = 0
1340 1345 self.__nChannels = dataOut.nChannels
1341 1346 if DPL != None:
1342 1347 self.DPL=DPL
1343 1348 else:
1344 1349 #self.DPL=dataOut.DPL ###REDEFINIR
1345 1350 self.DPL=0
1346 1351
1347 1352 if n is None and timeInterval is None:
1348 1353 raise ValueError("n or timeInterval should be specified ...")
1349 1354
1350 1355 if n is not None:
1351 1356 self.n = int(n)
1352 1357 else:
1353 1358 self.__integrationtime = int(timeInterval)
1354 1359 self.n = None
1355 1360 self.__byTime = True
1356 1361
1357 1362
1358 1363 if minHei == None:
1359 1364 minHei = self.dataOut.heightList[0]
1360 1365
1361 1366 if maxHei == None:
1362 1367 maxHei = self.dataOut.heightList[-1]
1363 1368
1364 1369 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1365 1370 print('minHei: %.2f is out of the heights range' % (minHei))
1366 1371 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
1367 1372 minHei = self.dataOut.heightList[0]
1368 1373
1369 1374 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
1370 1375 print('maxHei: %.2f is out of the heights range' % (maxHei))
1371 1376 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
1372 1377 maxHei = self.dataOut.heightList[-1]
1373 1378
1374 1379 ind_list1 = numpy.where(self.dataOut.heightList >= minHei)
1375 1380 ind_list2 = numpy.where(self.dataOut.heightList <= maxHei)
1376 1381 self.minHei_ind = ind_list1[0][0]
1377 1382 self.maxHei_ind = ind_list2[0][-1]
1378 1383
1379 1384
1380 1385 def putData(self, data_spc, data_cspc, data_dc):
1381 1386 """
1382 1387 Add a profile to the __buffer_spc and increase in one the __profileIndex
1383 1388
1384 1389 """
1385 1390
1386 1391 self.__buffer_spc.append(data_spc)
1387 1392
1388 1393 if self.__nChannels < 2:
1389 1394 self.__buffer_cspc = None
1390 1395 else:
1391 1396 self.__buffer_cspc.append(data_cspc)
1392 1397
1393 1398 if data_dc is None:
1394 1399 self.__buffer_dc = None
1395 1400 else:
1396 1401 self.__buffer_dc += data_dc
1397 1402
1398 1403 self.__profIndex += 1
1399 1404
1400 1405 return
1401 1406
1402 1407 def hildebrand_sekhon_Integration(self,sortdata,navg, factor):
1403 1408 #data debe estar ordenado
1404 1409 #sortdata = numpy.sort(data, axis=None)
1405 1410 #sortID=data.argsort()
1406 1411 lenOfData = len(sortdata)
1407 1412 nums_min = lenOfData*factor
1408 1413 if nums_min <= 5:
1409 1414 nums_min = 5
1410 1415 sump = 0.
1411 1416 sumq = 0.
1412 1417 j = 0
1413 1418 cont = 1
1414 1419 while((cont == 1)and(j < lenOfData)):
1415 1420 sump += sortdata[j]
1416 1421 sumq += sortdata[j]**2
1417 1422 if j > nums_min:
1418 1423 rtest = float(j)/(j-1) + 1.0/navg
1419 1424 if ((sumq*j) > (rtest*sump**2)):
1420 1425 j = j - 1
1421 1426 sump = sump - sortdata[j]
1422 1427 sumq = sumq - sortdata[j]**2
1423 1428 cont = 0
1424 1429 j += 1
1425 1430 #lnoise = sump / j
1426 1431 #print("H S done")
1427 1432 #return j,sortID
1428 1433 return j
1429 1434
1430 1435
1431 1436 def pushData(self):
1432 1437 """
1433 1438 Return the sum of the last profiles and the profiles used in the sum.
1434 1439
1435 1440 Affected:
1436 1441
1437 1442 self.__profileIndex
1438 1443
1439 1444 """
1440 1445 bufferH=None
1441 1446 buffer=None
1442 1447 buffer1=None
1443 1448 buffer_cspc=None
1444 1449 #print("aes: ", self.__buffer_cspc)
1445 1450 self.__buffer_spc=numpy.array(self.__buffer_spc)
1446 1451 if self.__nChannels > 1 :
1447 1452 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1448 1453
1449 1454 #print("FREQ_DC",self.__buffer_spc.shape,self.__buffer_cspc.shape)
1450 1455
1451 1456 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1452 1457 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1453 1458
1454 1459 self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
1455 1460
1456 1461 for k in range(self.minHei_ind,self.maxHei_ind):
1457 1462 if self.__nChannels > 1:
1458 1463 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1459 1464
1460 1465 outliers_IDs_cspc=[]
1461 1466 cspc_outliers_exist=False
1462 1467 for i in range(self.nChannels):#dataOut.nChannels):
1463 1468
1464 1469 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1465 1470 indexes=[]
1466 1471 #sortIDs=[]
1467 1472 outliers_IDs=[]
1468 1473
1469 1474 for j in range(self.nProfiles): #frecuencias en el tiempo
1470 1475 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1471 1476 # continue
1472 1477 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1473 1478 # continue
1474 1479 buffer=buffer1[:,j]
1475 1480 sortdata = numpy.sort(buffer, axis=None)
1476 1481
1477 1482 sortID=buffer.argsort()
1478 1483 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1479 1484
1480 1485 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1481 1486
1482 1487 # fig,ax = plt.subplots()
1483 1488 # ax.set_title(str(k)+" "+str(j))
1484 1489 # x=range(len(sortdata))
1485 1490 # ax.scatter(x,sortdata)
1486 1491 # ax.axvline(index)
1487 1492 # plt.show()
1488 1493
1489 1494 indexes.append(index)
1490 1495 #sortIDs.append(sortID)
1491 1496 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1492 1497
1493 1498 #print("Outliers: ",outliers_IDs)
1494 1499 outliers_IDs=numpy.array(outliers_IDs)
1495 1500 outliers_IDs=outliers_IDs.ravel()
1496 1501 outliers_IDs=numpy.unique(outliers_IDs)
1497 1502 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1498 1503 indexes=numpy.array(indexes)
1499 1504 indexmin=numpy.min(indexes)
1500 1505
1501 1506
1502 1507 #print(indexmin,buffer1.shape[0], k)
1503 1508
1504 1509 # fig,ax = plt.subplots()
1505 1510 # ax.plot(sortdata)
1506 1511 # ax2 = ax.twinx()
1507 1512 # x=range(len(indexes))
1508 1513 # #plt.scatter(x,indexes)
1509 1514 # ax2.scatter(x,indexes)
1510 1515 # plt.show()
1511 1516
1512 1517 if indexmin != buffer1.shape[0]:
1513 1518 if self.__nChannels > 1:
1514 1519 cspc_outliers_exist= True
1515 1520
1516 1521 lt=outliers_IDs
1517 1522 #avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1518 1523
1519 1524 for p in list(outliers_IDs):
1520 1525 #buffer1[p,:]=avg
1521 1526 buffer1[p,:] = numpy.NaN
1522 1527
1523 1528 self.dataOutliers[i,k] = len(outliers_IDs)
1524 1529
1525 1530
1526 1531 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1527 1532
1528 1533
1529 1534 if self.__nChannels > 1:
1530 1535 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1531 1536
1532 1537
1533 1538 if self.__nChannels > 1:
1534 1539 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1535 1540 if cspc_outliers_exist:
1536 1541
1537 1542 lt=outliers_IDs_cspc
1538 1543
1539 1544 #avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1540 1545 for p in list(outliers_IDs_cspc):
1541 1546 #buffer_cspc[p,:]=avg
1542 1547 buffer_cspc[p,:] = numpy.NaN
1543 1548
1544 1549 if self.__nChannels > 1:
1545 1550 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1546 1551
1547 1552
1548 1553
1549 1554
1550 1555 nOutliers = len(outliers_IDs)
1551 1556 #print("Outliers n: ",self.dataOutliers,nOutliers)
1552 1557 buffer=None
1553 1558 bufferH=None
1554 1559 buffer1=None
1555 1560 buffer_cspc=None
1556 1561
1557 1562
1558 1563 buffer=None
1559 1564
1560 1565 #data_spc = numpy.sum(self.__buffer_spc,axis=0)
1561 1566 data_spc = numpy.nansum(self.__buffer_spc,axis=0)
1562 1567 if self.__nChannels > 1:
1563 1568 #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1564 1569 data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
1565 1570 else:
1566 1571 data_cspc = None
1567 1572 data_dc = self.__buffer_dc
1568 1573 #(CH, HEIGH)
1569 1574 self.maxProfilesInt = self.__profIndex - 1
1570 1575 n = self.__profIndex - self.dataOutliers # n becomes a matrix
1571 1576
1572 1577 self.__buffer_spc = []
1573 1578 self.__buffer_cspc = []
1574 1579 self.__buffer_dc = 0
1575 1580 self.__profIndex = 0
1576 1581 #print("cleaned ",data_cspc)
1577 1582 return data_spc, data_cspc, data_dc, n
1578 1583
1579 1584 def byProfiles(self, *args):
1580 1585
1581 1586 self.__dataReady = False
1582 1587 avgdata_spc = None
1583 1588 avgdata_cspc = None
1584 1589 avgdata_dc = None
1585 1590
1586 1591 self.putData(*args)
1587 1592
1588 1593 if self.__profIndex == self.n:
1589 1594
1590 1595 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1591 1596 self.n_ints = n
1592 1597 self.__dataReady = True
1593 1598
1594 1599 return avgdata_spc, avgdata_cspc, avgdata_dc
1595 1600
1596 1601 def byTime(self, datatime, *args):
1597 1602
1598 1603 self.__dataReady = False
1599 1604 avgdata_spc = None
1600 1605 avgdata_cspc = None
1601 1606 avgdata_dc = None
1602 1607
1603 1608 self.putData(*args)
1604 1609
1605 1610 if (datatime - self.__initime) >= self.__integrationtime:
1606 1611 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1607 1612 self.n_ints = n
1608 1613 self.__dataReady = True
1609 1614
1610 1615 return avgdata_spc, avgdata_cspc, avgdata_dc
1611 1616
1612 1617 def integrate(self, datatime, *args):
1613 1618
1614 1619 if self.__profIndex == 0:
1615 1620 self.__initime = datatime
1616 1621
1617 1622 if self.__byTime:
1618 1623 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1619 1624 datatime, *args)
1620 1625 else:
1621 1626 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1622 1627
1623 1628 if not self.__dataReady:
1624 1629 return None, None, None, None
1625 1630
1626 1631 #print("integrate", avgdata_cspc)
1627 1632 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1628 1633
1629 1634 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
1630 1635 self.dataOut = dataOut
1631 1636 if n == 1:
1632 1637 return self.dataOut
1633 1638 self.dataOut.processingHeaderObj.timeIncohInt = timeInterval
1634 #print("nchannels", self.dataOut.nChannels)
1639
1640 if dataOut.flagProfilesByRange == True:
1641 self._flagProfilesByRange = True
1642
1635 1643 if self.dataOut.nChannels == 1:
1636 1644 self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
1637 1645 #print("IN spc:", self.dataOut.data_spc.shape, self.dataOut.data_cspc)
1638 1646 if not self.isConfig:
1639 1647 self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor)
1640 1648 self.isConfig = True
1641 1649
1642 1650 if not self.ByLags:
1643 1651 self.nProfiles=self.dataOut.nProfiles
1644 1652 self.nChannels=self.dataOut.nChannels
1645 1653 self.nHeights=self.dataOut.nHeights
1646 1654 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1647 1655 self.dataOut.data_spc,
1648 1656 self.dataOut.data_cspc,
1649 1657 self.dataOut.data_dc)
1650 1658 else:
1651 1659 self.nProfiles=self.dataOut.nProfiles
1652 1660 self.nChannels=self.dataOut.nChannels
1653 1661 self.nHeights=self.dataOut.nHeights
1654 1662 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1655 1663 self.dataOut.dataLag_spc,
1656 1664 self.dataOut.dataLag_cspc,
1657 1665 self.dataOut.dataLag_dc)
1658 1666 self.dataOut.flagNoData = True
1667
1668 if self._flagProfilesByRange:
1669 dataOut.flagProfilesByRange = True
1670 self._nProfilesByRange += dataOut.nProfilesByRange
1671
1659 1672 if self.__dataReady:
1660 1673
1661 1674 if not self.ByLags:
1662 1675 if self.nChannels == 1:
1663 1676 #print("f int", avgdata_spc.shape)
1664 1677 self.dataOut.data_spc = avgdata_spc
1665 1678 self.dataOut.data_cspc = None
1666 1679 else:
1667 1680 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
1668 1681 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1669 1682 self.dataOut.data_dc = avgdata_dc
1670 1683 self.dataOut.data_outlier = self.dataOutliers
1671 1684
1685
1672 1686 else:
1673 1687 self.dataOut.dataLag_spc = avgdata_spc
1674 1688 self.dataOut.dataLag_cspc = avgdata_cspc
1675 1689 self.dataOut.dataLag_dc = avgdata_dc
1676 1690
1677 1691 self.dataOut.data_spc=self.dataOut.dataLag_spc[:,:,:,self.dataOut.LagPlot]
1678 1692 self.dataOut.data_cspc=self.dataOut.dataLag_cspc[:,:,:,self.dataOut.LagPlot]
1679 1693 self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot]
1680 1694
1681 1695
1682 1696 self.dataOut.nIncohInt *= self.n_ints
1683 1697 #print("maxProfilesInt: ",self.maxProfilesInt)
1684 1698
1685 1699 self.dataOut.utctime = avgdatatime
1686 1700 self.dataOut.flagNoData = False
1687 1701
1702 dataOut.nProfilesByRange = self._nProfilesByRange
1703 self._nProfilesByRange = 0
1704 self._flagProfilesByRange = False
1705
1688 1706 # #update Processing Header:
1689 1707 # self.dataOut.processingHeaderObj.nIncohInt =
1690 1708 # self.dataOut.processingHeaderObj.nFFTPoints = self.dataOut.nFFTPoints
1691 1709
1692 1710 #print("Faraday Integration DONE...", self.dataOut.data_cspc)
1693 1711 #print(self.dataOut.flagNoData)
1694 1712 return self.dataOut
1695 1713
1696 1714
1697 1715
1698 1716 class removeInterference(Operation):
1699 1717
1700 1718 def removeInterference3(self, min_hei = None, max_hei = None):
1701 1719
1702 1720 jspectra = self.dataOut.data_spc
1703 1721 #jcspectra = self.dataOut.data_cspc
1704 1722 jnoise = self.dataOut.getNoise()
1705 1723 num_incoh = self.dataOut.max_nIncohInt
1706 1724 #print(jspectra.shape)
1707 1725 num_channel, num_prof, num_hei = jspectra.shape
1708 1726 minHei = min_hei
1709 1727 maxHei = max_hei
1710 1728 ########################################################################
1711 1729 if minHei == None or (minHei < self.dataOut.heightList[0]):
1712 1730 minHei = self.dataOut.heightList[0]
1713 1731
1714 1732 if maxHei == None or (maxHei > self.dataOut.heightList[-1]):
1715 1733 maxHei = self.dataOut.heightList[-1]
1716 1734 minIndex = 0
1717 1735 maxIndex = 0
1718 1736 heights = self.dataOut.heightList
1719 1737
1720 1738 inda = numpy.where(heights >= minHei)
1721 1739 indb = numpy.where(heights <= maxHei)
1722 1740
1723 1741 try:
1724 1742 minIndex = inda[0][0]
1725 1743 except:
1726 1744 minIndex = 0
1727 1745 try:
1728 1746 maxIndex = indb[0][-1]
1729 1747 except:
1730 1748 maxIndex = len(heights)
1731 1749
1732 1750 if (minIndex < 0) or (minIndex > maxIndex):
1733 1751 raise ValueError("some value in (%d,%d) is not valid" % (
1734 1752 minIndex, maxIndex))
1735 1753 if (maxIndex >= self.dataOut.nHeights):
1736 1754 maxIndex = self.dataOut.nHeights - 1
1737 1755
1738 1756 ########################################################################
1739 1757
1740 1758
1741 1759 #dataOut.max_nIncohInt * dataOut.nCohInt
1742 1760 norm = self.dataOut.nIncohInt /self.dataOut.max_nIncohInt
1743 1761 #print(norm.shape)
1744 1762 # Subrutina de Remocion de la Interferencia
1745 1763 for ich in range(num_channel):
1746 1764 # Se ordena los espectros segun su potencia (menor a mayor)
1747 1765 #power = jspectra[ich, mask_prof, :]
1748 1766 interf = jspectra[ich, :, minIndex:maxIndex]
1749 1767 #print(interf.shape)
1750 1768 inttef = interf.mean(axis=1)
1751 1769
1752 1770 for hei in range(num_hei):
1753 1771 temp = jspectra[ich,:, hei]
1754 1772 temp -= inttef
1755 1773 temp += jnoise[ich]*norm[ich,hei]
1756 1774 jspectra[ich,:, hei] = temp
1757 1775
1758 1776 # Guardar Resultados
1759 1777 self.dataOut.data_spc = jspectra
1760 1778 #self.dataOut.data_cspc = jcspectra
1761 1779
1762 1780 return 1
1763 1781
1764 1782 def removeInterference2(self):
1765 1783
1766 1784 cspc = self.dataOut.data_cspc
1767 1785 spc = self.dataOut.data_spc
1768 1786 Heights = numpy.arange(cspc.shape[2])
1769 1787 realCspc = numpy.abs(cspc)
1770 1788
1771 1789 for i in range(cspc.shape[0]):
1772 1790 LinePower= numpy.sum(realCspc[i], axis=0)
1773 1791 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1774 1792 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1775 1793 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1776 1794 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1777 1795 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1778 1796
1779 1797
1780 1798 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1781 1799 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1782 1800 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1783 1801 cspc[i,InterferenceRange,:] = numpy.NaN
1784 1802
1785 1803 self.dataOut.data_cspc = cspc
1786 1804
1787 1805 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1788 1806
1789 1807 jspectra = self.dataOut.data_spc
1790 1808 jcspectra = self.dataOut.data_cspc
1791 1809 jnoise = self.dataOut.getNoise()
1792 1810 #num_incoh = self.dataOut.nIncohInt
1793 1811 num_incoh = self.dataOut.max_nIncohInt
1794 1812 #print("spc: ", jspectra.shape, jcspectra)
1795 1813 num_channel = jspectra.shape[0]
1796 1814 num_prof = jspectra.shape[1]
1797 1815 num_hei = jspectra.shape[2]
1798 1816
1799 1817 # hei_interf
1800 1818 if hei_interf is None:
1801 1819 count_hei = int(num_hei / 2) # a half of total ranges
1802 1820 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1803 1821 hei_interf = numpy.asarray(hei_interf)[0]
1804 1822 #print(hei_interf)
1805 1823 # nhei_interf
1806 1824 if (nhei_interf == None):
1807 1825 nhei_interf = 5
1808 1826 if (nhei_interf < 1):
1809 1827 nhei_interf = 1
1810 1828 if (nhei_interf > count_hei):
1811 1829 nhei_interf = count_hei
1812 1830 if (offhei_interf == None):
1813 1831 offhei_interf = 0
1814 1832
1815 1833 ind_hei = list(range(num_hei))
1816 1834 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1817 1835 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1818 1836 mask_prof = numpy.asarray(list(range(num_prof)))
1819 1837 num_mask_prof = mask_prof.size
1820 1838 comp_mask_prof = [0, num_prof / 2]
1821 1839
1822 1840 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1823 1841 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1824 1842 jnoise = numpy.nan
1825 1843 noise_exist = jnoise[0] < numpy.Inf
1826 1844
1827 1845 # Subrutina de Remocion de la Interferencia
1828 1846 for ich in range(num_channel):
1829 1847 # Se ordena los espectros segun su potencia (menor a mayor)
1830 1848 power = jspectra[ich, mask_prof, :]
1831 1849 power = power[:, hei_interf]
1832 1850 power = power.sum(axis=0)
1833 1851 psort = power.ravel().argsort()
1834 1852 print(hei_interf[psort[list(range(offhei_interf, nhei_interf + offhei_interf))]])
1835 1853 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1836 1854 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1837 1855 offhei_interf, nhei_interf + offhei_interf))]]]
1838 1856
1839 1857 if noise_exist:
1840 1858 # tmp_noise = jnoise[ich] / num_prof
1841 1859 tmp_noise = jnoise[ich]
1842 1860 junkspc_interf = junkspc_interf - tmp_noise
1843 1861 #junkspc_interf[:,comp_mask_prof] = 0
1844 1862 print(junkspc_interf.shape)
1845 1863 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1846 1864 jspc_interf = jspc_interf.transpose()
1847 1865 # Calculando el espectro de interferencia promedio
1848 1866 noiseid = numpy.where(jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1849 1867 noiseid = noiseid[0]
1850 1868 cnoiseid = noiseid.size
1851 1869 interfid = numpy.where(jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1852 1870 interfid = interfid[0]
1853 1871 cinterfid = interfid.size
1854 1872
1855 1873 if (cnoiseid > 0):
1856 1874 jspc_interf[noiseid] = 0
1857 1875 # Expandiendo los perfiles a limpiar
1858 1876 if (cinterfid > 0):
1859 1877 new_interfid = (
1860 1878 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1861 1879 new_interfid = numpy.asarray(new_interfid)
1862 1880 new_interfid = {x for x in new_interfid}
1863 1881 new_interfid = numpy.array(list(new_interfid))
1864 1882 new_cinterfid = new_interfid.size
1865 1883 else:
1866 1884 new_cinterfid = 0
1867 1885
1868 1886 for ip in range(new_cinterfid):
1869 1887 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1870 1888 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1871 1889
1872 1890 jspectra[ich, :, ind_hei] = jspectra[ich, :,ind_hei] - jspc_interf # Corregir indices
1873 1891
1874 1892 # Removiendo la interferencia del punto de mayor interferencia
1875 1893 ListAux = jspc_interf[mask_prof].tolist()
1876 1894 maxid = ListAux.index(max(ListAux))
1877 1895 print(cinterfid)
1878 1896 if cinterfid > 0:
1879 1897 for ip in range(cinterfid * (interf == 2) - 1):
1880 1898 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1881 1899 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1882 1900 cind = len(ind)
1883 1901
1884 1902 if (cind > 0):
1885 1903 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1886 1904 (1 + (numpy.random.uniform(cind) - 0.5) /
1887 1905 numpy.sqrt(num_incoh))
1888 1906
1889 1907 ind = numpy.array([-2, -1, 1, 2])
1890 1908 xx = numpy.zeros([4, 4])
1891 1909
1892 1910 for id1 in range(4):
1893 1911 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1894 1912 xx_inv = numpy.linalg.inv(xx)
1895 1913 xx = xx_inv[:, 0]
1896 1914 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1897 1915 yy = jspectra[ich, mask_prof[ind], :]
1898 1916 jspectra[ich, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1899 1917
1900 1918 indAux = (jspectra[ich, :, :] < tmp_noise *
1901 1919 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1902 1920 print(indAux)
1903 1921 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1904 1922 (1 - 1 / numpy.sqrt(num_incoh))
1905 1923
1906 1924 # Remocion de Interferencia en el Cross Spectra
1907 1925 if jcspectra is None:
1908 1926 return jspectra, jcspectra
1909 1927 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1910 1928 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1911 1929
1912 1930 for ip in range(num_pairs):
1913 1931
1914 1932 #-------------------------------------------
1915 1933
1916 1934 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1917 1935 cspower = cspower[:, hei_interf]
1918 1936 cspower = cspower.sum(axis=0)
1919 1937
1920 1938 cspsort = cspower.ravel().argsort()
1921 1939 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1922 1940 offhei_interf, nhei_interf + offhei_interf))]]]
1923 1941 junkcspc_interf = junkcspc_interf.transpose()
1924 1942 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1925 1943
1926 1944 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1927 1945
1928 1946 median_real = int(numpy.median(numpy.real(
1929 1947 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1930 1948 median_imag = int(numpy.median(numpy.imag(
1931 1949 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1932 1950 comp_mask_prof = [int(e) for e in comp_mask_prof]
1933 1951 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1934 1952 median_real, median_imag)
1935 1953
1936 1954 for iprof in range(num_prof):
1937 1955 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1938 1956 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1939 1957
1940 1958 # Removiendo la Interferencia
1941 1959 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1942 1960 :, ind_hei] - jcspc_interf
1943 1961
1944 1962 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1945 1963 maxid = ListAux.index(max(ListAux))
1946 1964
1947 1965 ind = numpy.array([-2, -1, 1, 2])
1948 1966 xx = numpy.zeros([4, 4])
1949 1967
1950 1968 for id1 in range(4):
1951 1969 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1952 1970
1953 1971 xx_inv = numpy.linalg.inv(xx)
1954 1972 xx = xx_inv[:, 0]
1955 1973
1956 1974 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1957 1975 yy = jcspectra[ip, mask_prof[ind], :]
1958 1976 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1959 1977
1960 1978 # Guardar Resultados
1961 1979 self.dataOut.data_spc = jspectra
1962 1980 self.dataOut.data_cspc = jcspectra
1963 1981
1964 1982 return 1
1965 1983
1966 1984 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1, minHei=None, maxHei=None):
1967 1985
1968 1986 self.dataOut = dataOut
1969 1987
1970 1988 if mode == 1:
1971 1989 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1972 1990 elif mode == 2:
1973 1991 self.removeInterference2()
1974 1992 elif mode == 3:
1975 1993 self.removeInterference3(min_hei=minHei, max_hei=maxHei)
1976 1994 return self.dataOut
1977 1995
1978 1996
1979 1997 class IncohInt(Operation):
1980 1998
1981 1999 __profIndex = 0
1982 2000 __withOverapping = False
1983 2001
1984 2002 __byTime = False
1985 2003 __initime = None
1986 2004 __lastdatatime = None
1987 2005 __integrationtime = None
1988 2006
1989 2007 __buffer_spc = None
1990 2008 __buffer_cspc = None
1991 2009 __buffer_dc = None
1992 2010
1993 2011 __dataReady = False
1994 2012
1995 2013 __timeInterval = None
1996 2014 incohInt = 0
1997 2015 nOutliers = 0
1998 2016 n = None
1999 2017
2018 _flagProfilesByRange = False
2019 _nProfilesByRange = 0
2000 2020 def __init__(self):
2001 2021
2002 2022 Operation.__init__(self)
2003 2023
2004 2024 def setup(self, n=None, timeInterval=None, overlapping=False):
2005 2025 """
2006 2026 Set the parameters of the integration class.
2007 2027
2008 2028 Inputs:
2009 2029
2010 2030 n : Number of coherent integrations
2011 2031 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
2012 2032 overlapping :
2013 2033
2014 2034 """
2015 2035
2016 2036 self.__initime = None
2017 2037 self.__lastdatatime = 0
2018 2038
2019 2039 self.__buffer_spc = 0
2020 2040 self.__buffer_cspc = 0
2021 2041 self.__buffer_dc = 0
2022 2042
2023 2043 self.__profIndex = 0
2024 2044 self.__dataReady = False
2025 2045 self.__byTime = False
2026 2046 self.incohInt = 0
2027 2047 self.nOutliers = 0
2028 2048 if n is None and timeInterval is None:
2029 2049 raise ValueError("n or timeInterval should be specified ...")
2030 2050
2031 2051 if n is not None:
2032 2052 self.n = int(n)
2033 2053 else:
2034 2054
2035 2055 self.__integrationtime = int(timeInterval)
2036 2056 self.n = None
2037 2057 self.__byTime = True
2038 2058
2039 2059
2040 2060 def putData(self, data_spc, data_cspc, data_dc):
2041 2061 """
2042 2062 Add a profile to the __buffer_spc and increase in one the __profileIndex
2043 2063
2044 2064 """
2045 2065 if data_spc.all() == numpy.nan :
2046 2066 print("nan ")
2047 2067 return
2048 2068 self.__buffer_spc += data_spc
2049 2069
2050 2070 if data_cspc is None:
2051 2071 self.__buffer_cspc = None
2052 2072 else:
2053 2073 self.__buffer_cspc += data_cspc
2054 2074
2055 2075 if data_dc is None:
2056 2076 self.__buffer_dc = None
2057 2077 else:
2058 2078 self.__buffer_dc += data_dc
2059 2079
2060 2080 self.__profIndex += 1
2061 2081
2062 2082 return
2063 2083
2064 2084 def pushData(self):
2065 2085 """
2066 2086 Return the sum of the last profiles and the profiles used in the sum.
2067 2087
2068 2088 Affected:
2069 2089
2070 2090 self.__profileIndex
2071 2091
2072 2092 """
2073 2093
2074 2094 data_spc = self.__buffer_spc
2075 2095 data_cspc = self.__buffer_cspc
2076 2096 data_dc = self.__buffer_dc
2077 2097 n = self.__profIndex
2078 2098
2079 2099 self.__buffer_spc = 0
2080 2100 self.__buffer_cspc = 0
2081 2101 self.__buffer_dc = 0
2082 2102
2083 2103
2084 2104 return data_spc, data_cspc, data_dc, n
2085 2105
2086 2106 def byProfiles(self, *args):
2087 2107
2088 2108 self.__dataReady = False
2089 2109 avgdata_spc = None
2090 2110 avgdata_cspc = None
2091 2111 avgdata_dc = None
2092 2112
2093 2113 self.putData(*args)
2094 2114
2095 2115 if self.__profIndex == self.n:
2096 2116
2097 2117 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
2098 2118 self.n = n
2099 2119 self.__dataReady = True
2100 2120
2101 2121 return avgdata_spc, avgdata_cspc, avgdata_dc
2102 2122
2103 2123 def byTime(self, datatime, *args):
2104 2124
2105 2125 self.__dataReady = False
2106 2126 avgdata_spc = None
2107 2127 avgdata_cspc = None
2108 2128 avgdata_dc = None
2109 2129
2110 2130 self.putData(*args)
2111 2131
2112 2132 if (datatime - self.__initime) >= self.__integrationtime:
2113 2133 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
2114 2134 self.n = n
2115 2135 self.__dataReady = True
2116 2136
2117 2137 return avgdata_spc, avgdata_cspc, avgdata_dc
2118 2138
2119 2139 def integrate(self, datatime, *args):
2120 2140
2121 2141 if self.__profIndex == 0:
2122 2142 self.__initime = datatime
2123 2143
2124 2144 if self.__byTime:
2125 2145 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
2126 2146 datatime, *args)
2127 2147 else:
2128 2148 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
2129 2149
2130 2150 if not self.__dataReady:
2131 2151 return None, None, None, None
2132 2152
2133 2153 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
2134 2154
2135 2155 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
2136 2156 if n == 1:
2137 2157 return dataOut
2138 2158
2139 2159 if dataOut.flagNoData == True:
2140 2160 return dataOut
2141 2161
2162 if dataOut.flagProfilesByRange == True:
2163 self._flagProfilesByRange = True
2164
2142 2165 dataOut.flagNoData = True
2143 2166 dataOut.processingHeaderObj.timeIncohInt = timeInterval
2144 2167 if not self.isConfig:
2145 2168 self.setup(n, timeInterval, overlapping)
2146 2169 self.isConfig = True
2147 2170
2148 2171
2149 2172 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
2150 2173 dataOut.data_spc,
2151 2174 dataOut.data_cspc,
2152 2175 dataOut.data_dc)
2153 2176
2154 2177 self.incohInt += dataOut.nIncohInt
2155 2178
2179
2156 2180 if isinstance(dataOut.data_outlier,numpy.ndarray) or isinstance(dataOut.data_outlier,int) or isinstance(dataOut.data_outlier, float):
2157 2181 self.nOutliers += dataOut.data_outlier
2158 2182
2183 if self._flagProfilesByRange:
2184 dataOut.flagProfilesByRange = True
2185 self._nProfilesByRange += dataOut.nProfilesByRange
2186
2159 2187 if self.__dataReady:
2160 2188 #print("prof: ",dataOut.max_nIncohInt,self.__profIndex)
2161 2189 dataOut.data_spc = avgdata_spc
2162 2190 dataOut.data_cspc = avgdata_cspc
2163 2191 dataOut.data_dc = avgdata_dc
2164 2192 dataOut.nIncohInt = self.incohInt
2165 2193 dataOut.data_outlier = self.nOutliers
2166 2194 dataOut.utctime = avgdatatime
2167 2195 dataOut.flagNoData = False
2168 2196 self.incohInt = 0
2169 2197 self.nOutliers = 0
2170 2198 self.__profIndex = 0
2199 dataOut.nProfilesByRange = self._nProfilesByRange
2200 self._nProfilesByRange = 0
2201 self._flagProfilesByRange = False
2171 2202 #print("IncohInt Done")
2172 2203 return dataOut
2173 2204
2174 2205 class dopplerFlip(Operation):
2175 2206
2176 2207 def run(self, dataOut):
2177 2208 # arreglo 1: (num_chan, num_profiles, num_heights)
2178 2209 self.dataOut = dataOut
2179 2210 # JULIA-oblicua, indice 2
2180 2211 # arreglo 2: (num_profiles, num_heights)
2181 2212 jspectra = self.dataOut.data_spc[2]
2182 2213 jspectra_tmp = numpy.zeros(jspectra.shape)
2183 2214 num_profiles = jspectra.shape[0]
2184 2215 freq_dc = int(num_profiles / 2)
2185 2216 # Flip con for
2186 2217 for j in range(num_profiles):
2187 2218 jspectra_tmp[num_profiles-j-1]= jspectra[j]
2188 2219 # Intercambio perfil de DC con perfil inmediato anterior
2189 2220 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
2190 2221 jspectra_tmp[freq_dc]= jspectra[freq_dc]
2191 2222 # canal modificado es re-escrito en el arreglo de canales
2192 2223 self.dataOut.data_spc[2] = jspectra_tmp
2193 2224
2194 2225 return self.dataOut
2195 2226
2196 2227
2197 2228
2198 2229
2199 2230
2200 2231
2201 2232 class cleanJULIAInterf(Operation):
2202 2233 """
2203 2234 Operación de prueba
2204 2235 """
2205 2236 __slots__ =('heights_indx', 'repeats','span' ,'step', 'factor', 'idate', 'idxs','isConfig','minHrefN', 'maxHrefN')
2206 2237 def __init__(self):
2207 2238 self.repeats = 0
2208 2239 self.factor=1
2209 2240 self.isConfig = False
2210 2241 self.idxs = None
2211 2242 self.heights_indx = None
2212 2243
2213 2244 def setup(self, dataOutHeightsList, heightsList, span=10, repeats=0, step=0, idate=None, startH=None, endH=None, minHref=None, maxHref=None):
2214 2245 totalHeihtList = dataOutHeightsList
2215 2246 heights = [float(hei) for hei in heightsList]
2216 2247 for r in range(repeats):
2217 2248 heights += [ (h+(step*(r+1))) for h in heights]
2218 2249 #print(heights)
2219 2250 self.heights_indx = [getHei_index(h,h,totalHeihtList)[0] for h in heights]
2220 2251
2221 2252 self.minHrefN, self.maxHrefN = getHei_index(minHref,maxHref,totalHeihtList)
2222 2253
2223 2254
2224 2255 self.config = True
2225 2256 self.span = span
2226 2257
2227 2258 def run(self, dataOut, heightsList, span=10, repeats=0, step=0, idate=None, startH=None, endH=None, minHref=None, maxHref=None):
2228 2259
2229 2260
2230 2261 self.dataOut = dataOut
2231 2262 startTime = datetime.datetime.combine(idate,startH)
2232 2263 endTime = datetime.datetime.combine(idate,endH)
2233 2264 currentTime = datetime.datetime.fromtimestamp(self.dataOut.utctime)
2234 2265
2235 2266 if currentTime < startTime or currentTime > endTime:
2236 2267 return self.dataOut
2237 2268
2238 2269 if not self.isConfig:
2239 2270 self.setup(self.dataOut.heightList,heightsList, span=span, repeats=repeats, step=step, idate=idate, startH=startH, endH=endH, minHref=minHref, maxHref=maxHref )
2240 2271
2241 2272 for ch in range(self.dataOut.data_spc.shape[0]):
2242 2273 i = 0
2243 2274 N_ref = self.dataOut.data_spc[ch, :, self.minHrefN: self.maxHrefN].mean()
2244 2275 mn = self.heights_indx[-1] - self.span/2
2245 2276 mx = self.heights_indx[-1] + self.span/2
2246 2277 J_lev = self.dataOut.data_spc[ch, :, mn: mx].mean() - N_ref
2247 2278
2248 2279 for hei in self.heights_indx:
2249 2280 h = hei - 1
2250 2281 mn_i = hei - self.span/2
2251 2282 mx_i = hei + self.span/2
2252 2283 self.dataOut.data_spc[ch, :,mn_i:mx_i ] -= J_lev
2253 2284 i += 1
2254 2285
2255 2286
2256 2287 return self.dataOut No newline at end of file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now