##// END OF EJS Templates
schain
RSS
test 2
joabAM -
r1751:2d2f96339643
parent child
Show More
Side by Side Unified
Context file:
r1751:2d2f96339643 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1163 +1,1167
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
25 25 def getNumpyDtype(dataTypeCode):
26 26
27 27 if dataTypeCode == 0:
28 28 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
29 29 elif dataTypeCode == 1:
30 30 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
31 31 elif dataTypeCode == 2:
32 32 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
33 33 elif dataTypeCode == 3:
34 34 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
35 35 elif dataTypeCode == 4:
36 36 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
37 37 elif dataTypeCode == 5:
38 38 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
39 39 else:
40 40 raise ValueError('dataTypeCode was not defined')
41 41
42 42 return numpyDtype
43 43
44 44
45 45 def getDataTypeCode(numpyDtype):
46 46
47 47 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
48 48 datatype = 0
49 49 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
50 50 datatype = 1
51 51 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
52 52 datatype = 2
53 53 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
54 54 datatype = 3
55 55 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
56 56 datatype = 4
57 57 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
58 58 datatype = 5
59 59 else:
60 60 datatype = None
61 61
62 62 return datatype
63 63
64 64
65 65 def hildebrand_sekhon(data, navg):
66 66 """
67 67 This method is for the objective determination of the noise level in Doppler spectra. This
68 68 implementation technique is based on the fact that the standard deviation of the spectral
69 69 densities is equal to the mean spectral density for white Gaussian noise
70 70
71 71 Inputs:
72 72 Data : heights
73 73 navg : numbers of averages
74 74
75 75 Return:
76 76 mean : noise's level
77 77 """
78 78
79 79 sortdata = numpy.sort(data, axis=None)
80 80 '''
81 81 lenOfData = len(sortdata)
82 82 nums_min = lenOfData*0.2
83 83
84 84 if nums_min <= 5:
85 85
86 86 nums_min = 5
87 87
88 88 sump = 0.
89 89 sumq = 0.
90 90
91 91 j = 0
92 92 cont = 1
93 93
94 94 while((cont == 1)and(j < lenOfData)):
95 95
96 96 sump += sortdata[j]
97 97 sumq += sortdata[j]**2
98 98
99 99 if j > nums_min:
100 100 rtest = float(j)/(j-1) + 1.0/navg
101 101 if ((sumq*j) > (rtest*sump**2)):
102 102 j = j - 1
103 103 sump = sump - sortdata[j]
104 104 sumq = sumq - sortdata[j]**2
105 105 cont = 0
106 106
107 107 j += 1
108 108
109 109 lnoise = sump / j
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 class GenericData(object):
123 123
124 124 flagNoData = True
125 125
126 126 def copy(self, inputObj=None):
127 127
128 128 if inputObj == None:
129 129 return copy.deepcopy(self)
130 130
131 131 for key in list(inputObj.__dict__.keys()):
132 132
133 133 attribute = inputObj.__dict__[key]
134 134
135 135 # If this attribute is a tuple or list
136 136 if type(inputObj.__dict__[key]) in (tuple, list):
137 137 self.__dict__[key] = attribute[:]
138 138 continue
139 139
140 140 # If this attribute is another object or instance
141 141 if hasattr(attribute, '__dict__'):
142 142 self.__dict__[key] = attribute.copy()
143 143 continue
144 144
145 145 self.__dict__[key] = inputObj.__dict__[key]
146 146
147 147 def deepcopy(self):
148 148
149 149 return copy.deepcopy(self)
150 150
151 151 def isEmpty(self):
152 152
153 153 return self.flagNoData
154 154
155 155 def isReady(self):
156 156
157 157 return not self.flagNoData
158 158
159 159
160 160 class JROData(GenericData):
161 161
162 162 systemHeaderObj = SystemHeader()
163 163 radarControllerHeaderObj = RadarControllerHeader()
164 164 type = None
165 165 datatype = None # dtype but in string
166 166 nProfiles = None
167 167 heightList = None
168 168 channelList = None
169 169 flagDiscontinuousBlock = False
170 170 useLocalTime = False
171 171 utctime = None
172 172 timeZone = None
173 173 dstFlag = None
174 174 errorCount = None
175 175 blocksize = None
176 176 flagDecodeData = False # asumo q la data no esta decodificada
177 177 flagDeflipData = False # asumo q la data no esta sin flip
178 178 flagShiftFFT = False
179 179 nCohInt = None
180 180 windowOfFilter = 1
181 181 C = 3e8
182 182 frequency = 49.92e6
183 183 realtime = False
184 184 beacon_heiIndexList = None
185 185 last_block = None
186 186 blocknow = None
187 187 azimuth = None
188 188 zenith = None
189 189 beam = Beam()
190 190 profileIndex = None
191 191 error = None
192 192 data = None
193 193 nmodes = None
194 194 metadata_list = ['heightList', 'timeZone', 'type']
195 195
196 196 ippFactor = 1 #Added to correct the freq and vel range for AMISR data
197 197 useInputBuffer = False
198 198 buffer_empty = True
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
216 216 def __str__(self):
217 217
218 218 return '{} - {}'.format(self.type, self.datatime())
219 219
220 220 def getNoise(self):
221 221
222 222 raise NotImplementedError
223 223
224 224 @property
225 225 def nChannels(self):
226 226
227 227 return len(self.channelList)
228 228
229 229 @property
230 230 def channelIndexList(self):
231 231
232 232 return list(range(self.nChannels))
233 233
234 234 @property
235 235 def nHeights(self):
236 236
237 237 return len(self.heightList)
238 238
239 239 def getDeltaH(self):
240 240
241 241 return self.heightList[1] - self.heightList[0]
242 242
243 243 @property
244 244 def ltctime(self):
245
245 try:
246 self.timeZone = self.timeZone.decode("utf-8")
247 except Exception as e:
248 pass
249
246 250 if self.useLocalTime:
247 251 if self.timeZone =='lt':
248 252 return self.utctime - 300 * 60
249 253 elif self.timeZone =='ut':
250 254 return self.utctime
251 255 else:
252 log.error("No valid timeZone detected")
256 log.error("No valid timeZone detected:{}".format(self.timeZone))
253 257 return self.utctime
254 258
255 259 @property
256 260 def datatime(self):
257 261
258 262 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
259 263 return datatimeValue
260 264
261 265 def getTimeRange(self):
262 266
263 267 datatime = []
264 268
265 269 datatime.append(self.ltctime)
266 270 datatime.append(self.ltctime + self.timeInterval + 1)
267 271
268 272 datatime = numpy.array(datatime)
269 273
270 274 return datatime
271 275
272 276 def getFmaxTimeResponse(self):
273 277
274 278 period = (10**-6) * self.getDeltaH() / (0.15)
275 279
276 280 PRF = 1. / (period * self.nCohInt)
277 281
278 282 fmax = PRF
279 283
280 284 return fmax
281 285
282 286 def getFmax(self):
283 287 PRF = 1. / (self.__ippSeconds * self.nCohInt)
284 288
285 289 fmax = PRF
286 290 return fmax
287 291
288 292 def getVmax(self):
289 293
290 294 _lambda = self.C / self.frequency
291 295
292 296 vmax = self.getFmax() * _lambda / 2
293 297
294 298 return vmax
295 299
296 300 ## Radar Controller Header must be immutable
297 301 @property
298 302 def ippSeconds(self):
299 303 '''
300 304 '''
301 305 #return self.radarControllerHeaderObj.ippSeconds
302 306 return self.__ippSeconds
303 307
304 308 @ippSeconds.setter
305 309 def ippSeconds(self, ippSeconds):
306 310 '''
307 311 '''
308 312 #self.radarControllerHeaderObj.ippSeconds = ippSeconds
309 313 self.__ippSeconds = ippSeconds
310 314 self.__ipp = ippSeconds*SPEED_OF_LIGHT/2000.0
311 315
312 316 @property
313 317 def code(self):
314 318 '''
315 319 '''
316 320 # return self.radarControllerHeaderObj.code
317 321 return self.__code
318 322
319 323 @code.setter
320 324 def code(self, code):
321 325 '''
322 326 '''
323 327 # self.radarControllerHeaderObj.code = code
324 328 self.__code = code
325 329
326 330 @property
327 331 def nCode(self):
328 332 '''
329 333 '''
330 334 # return self.radarControllerHeaderObj.nCode
331 335 return self.__nCode
332 336
333 337 @nCode.setter
334 338 def nCode(self, ncode):
335 339 '''
336 340 '''
337 341 # self.radarControllerHeaderObj.nCode = ncode
338 342 self.__nCode = ncode
339 343
340 344 @property
341 345 def nBaud(self):
342 346 '''
343 347 '''
344 348 # return self.radarControllerHeaderObj.nBaud
345 349 return self.__nBaud
346 350
347 351 @nBaud.setter
348 352 def nBaud(self, nbaud):
349 353 '''
350 354 '''
351 355 # self.radarControllerHeaderObj.nBaud = nbaud
352 356 self.__nBaud = nbaud
353 357
354 358 @property
355 359 def ipp(self):
356 360 '''
357 361 '''
358 362 # return self.radarControllerHeaderObj.ipp
359 363 return self.__ipp
360 364
361 365 @ipp.setter
362 366 def ipp(self, ipp):
363 367 '''
364 368 '''
365 369 # self.radarControllerHeaderObj.ipp = ipp
366 370 self.__ipp = ipp
367 371
368 372 @property
369 373 def metadata(self):
370 374 '''
371 375 '''
372 376
373 377 return {attr: getattr(self, attr) for attr in self.metadata_list}
374 378
375 379
376 380 class Voltage(JROData):
377 381
378 382 dataPP_POW = None
379 383 dataPP_DOP = None
380 384 dataPP_WIDTH = None
381 385 dataPP_SNR = None
382 386
383 387 # To use oper
384 388 flagProfilesByRange = False
385 389 nProfilesByRange = None
386 390 max_nIncohInt = 1
387 391
388 392 def __init__(self):
389 393 '''
390 394 Constructor
391 395 '''
392 396
393 397 self.useLocalTime = True
394 398 self.radarControllerHeaderObj = RadarControllerHeader()
395 399 self.systemHeaderObj = SystemHeader()
396 400 self.processingHeaderObj = ProcessingHeader()
397 401 self.type = "Voltage"
398 402 self.data = None
399 403 self.nProfiles = None
400 404 self.heightList = None
401 405 self.channelList = None
402 406 self.flagNoData = True
403 407 self.flagDiscontinuousBlock = False
404 408 self.utctime = None
405 409 self.timeZone = 0
406 410 self.dstFlag = None
407 411 self.errorCount = None
408 412 self.nCohInt = None
409 413 self.blocksize = None
410 414 self.flagCohInt = False
411 415 self.flagDecodeData = False # asumo q la data no esta decodificada
412 416 self.flagDeflipData = False # asumo q la data no esta sin flip
413 417 self.flagShiftFFT = False
414 418 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
415 419 self.profileIndex = 0
416 420 self.ippFactor=1
417 421 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
418 422 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
419 423
420 424 def getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None):
421 425 """
422 426 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
423 427
424 428 Return:
425 429 noiselevel
426 430 """
427 431
428 432 if channel != None:
429 433 data = self.data[channel,ymin_index:ymax_index]
430 434 nChannels = 1
431 435 else:
432 436 data = self.data[:,ymin_index:ymax_index]
433 437 nChannels = self.nChannels
434 438
435 439 noise = numpy.zeros(nChannels)
436 440 power = data * numpy.conjugate(data)
437 441
438 442 for thisChannel in range(nChannels):
439 443 if nChannels == 1:
440 444 daux = power[:].real
441 445 else:
442 446 daux = power[thisChannel, :].real
443 447 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
444 448
445 449 return noise
446 450
447 451 def getNoise(self, type=1, channel=None,ymin_index=None, ymax_index=None):
448 452
449 453 if type == 1:
450 454 noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index)
451 455
452 456 return noise
453 457
454 458 def getPower(self, channel=None):
455 459
456 460 if channel != None:
457 461 data = self.data[channel]
458 462 else:
459 463 data = self.data
460 464
461 465 power = data * numpy.conjugate(data)
462 466 powerdB = 10 * numpy.log10(power.real)
463 467 powerdB = numpy.squeeze(powerdB)
464 468
465 469 return powerdB
466 470 @property
467 471 def data_pow(self):
468 472 return self.getPower()
469 473
470 474 @property
471 475 def timeInterval(self):
472 476
473 477 return self.ippSeconds * self.nCohInt
474 478
475 479 noise = property(getNoise, "I'm the 'nHeights' property.")
476 480
477 481
478 482 class Spectra(JROData):
479 483
480 484 data_outlier = None
481 485 flagProfilesByRange = False
482 486 nProfilesByRange = None
483 487
484 488 def __init__(self):
485 489 '''
486 490 Constructor
487 491 '''
488 492
489 493 self.data_dc = None
490 494 self.data_spc = None
491 495 self.data_cspc = None
492 496 self.useLocalTime = True
493 497 self.radarControllerHeaderObj = RadarControllerHeader()
494 498 self.systemHeaderObj = SystemHeader()
495 499 self.processingHeaderObj = ProcessingHeader()
496 500 self.type = "Spectra"
497 501 self.timeZone = 0
498 502 self.nProfiles = None
499 503 self.heightList = None
500 504 self.channelList = None
501 505 self.pairsList = None
502 506 self.flagNoData = True
503 507 self.flagDiscontinuousBlock = False
504 508 self.utctime = None
505 509 self.nCohInt = None
506 510 self.nIncohInt = None
507 511 self.blocksize = None
508 512 self.nFFTPoints = None
509 513 self.wavelength = None
510 514 self.flagDecodeData = False # asumo q la data no esta decodificada
511 515 self.flagDeflipData = False # asumo q la data no esta sin flip
512 516 self.flagShiftFFT = False
513 517 self.ippFactor = 1
514 518 self.beacon_heiIndexList = []
515 519 self.noise_estimation = None
516 520 self.codeList = []
517 521 self.azimuthList = []
518 522 self.elevationList = []
519 523 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
520 524 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
521 525
522 526 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
523 527 """
524 528 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
525 529
526 530 Return:
527 531 noiselevel
528 532 """
529 533
530 534 noise = numpy.zeros(self.nChannels)
531 535
532 536 for channel in range(self.nChannels):
533 537 daux = self.data_spc[channel,
534 538 xmin_index:xmax_index, ymin_index:ymax_index]
535 539 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
536 540 noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt[channel])
537 541
538 542 return noise
539 543
540 544 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
541 545
542 546 if self.noise_estimation is not None:
543 547 # this was estimated by getNoise Operation defined in jroproc_spectra.py
544 548 return self.noise_estimation
545 549 else:
546 550 noise = self.getNoisebyHildebrand(
547 551 xmin_index, xmax_index, ymin_index, ymax_index)
548 552 return noise
549 553
550 554 def getFreqRangeTimeResponse(self, extrapoints=0):
551 555
552 556 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
553 557 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
554 558
555 559 return freqrange
556 560
557 561 def getAcfRange(self, extrapoints=0):
558 562
559 563 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
560 564 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
561 565
562 566 return freqrange
563 567
564 568 def getFreqRange(self, extrapoints=0):
565 569
566 570 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
567 571 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
568 572
569 573 return freqrange
570 574
571 575 def getVelRange(self, extrapoints=0):
572 576
573 577 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
574 578 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
575 579
576 580 if self.nmodes:
577 581 return velrange/self.nmodes
578 582 else:
579 583 return velrange
580 584
581 585 @property
582 586 def nPairs(self):
583 587
584 588 return len(self.pairsList)
585 589
586 590 @property
587 591 def pairsIndexList(self):
588 592
589 593 return list(range(self.nPairs))
590 594
591 595 @property
592 596 def normFactor(self):
593 597
594 598 pwcode = 1
595 599 if self.flagDecodeData:
596 600 try:
597 601 pwcode = numpy.sum(self.code[0]**2)
598 602 except Exception as e:
599 603 log.warning("Failed pwcode read, setting to 1")
600 604 pwcode = 1
601 605 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
602 606 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
603 607 if self.flagProfilesByRange:
604 608 normFactor *= (self.nProfilesByRange/self.nProfilesByRange.max())
605 609 return normFactor
606 610
607 611 @property
608 612 def flag_cspc(self):
609 613
610 614 if self.data_cspc is None:
611 615 return True
612 616
613 617 return False
614 618
615 619 @property
616 620 def flag_dc(self):
617 621
618 622 if self.data_dc is None:
619 623 return True
620 624
621 625 return False
622 626
623 627 @property
624 628 def timeInterval(self):
625 629
626 630 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
627 631 if self.nmodes:
628 632 return self.nmodes*timeInterval
629 633 else:
630 634 return timeInterval
631 635
632 636 def getPower(self):
633 637
634 638 factor = self.normFactor
635 639 power = numpy.zeros( (self.nChannels,self.nHeights) )
636 640 for ch in range(self.nChannels):
637 641 z = None
638 642 if hasattr(factor,'shape'):
639 643 if factor.ndim > 1:
640 644 z = self.data_spc[ch]/factor[ch]
641 645 else:
642 646 z = self.data_spc[ch]/factor
643 647 else:
644 648 z = self.data_spc[ch]/factor
645 649 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
646 650 avg = numpy.average(z, axis=0)
647 651 power[ch] = 10 * numpy.log10(avg)
648 652 return power
649 653
650 654 @property
651 655 def max_nIncohInt(self):
652 656
653 657 ints = numpy.zeros(self.nChannels)
654 658 for ch in range(self.nChannels):
655 659 if hasattr(self.nIncohInt,'shape'):
656 660 if self.nIncohInt.ndim > 1:
657 661 ints[ch,] = self.nIncohInt[ch].max()
658 662 else:
659 663 ints[ch,] = self.nIncohInt
660 664 self.nIncohInt = int(self.nIncohInt)
661 665 else:
662 666 ints[ch,] = self.nIncohInt
663 667
664 668 return ints
665 669
666 670 def getCoherence(self, pairsList=None, phase=False):
667 671
668 672 z = []
669 673 if pairsList is None:
670 674 pairsIndexList = self.pairsIndexList
671 675 else:
672 676 pairsIndexList = []
673 677 for pair in pairsList:
674 678 if pair not in self.pairsList:
675 679 raise ValueError("Pair %s is not in dataOut.pairsList" % (
676 680 pair))
677 681 pairsIndexList.append(self.pairsList.index(pair))
678 682 for i in range(len(pairsIndexList)):
679 683 pair = self.pairsList[pairsIndexList[i]]
680 684 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
681 685 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
682 686 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
683 687 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
684 688 if phase:
685 689 data = numpy.arctan2(avgcoherenceComplex.imag,
686 690 avgcoherenceComplex.real) * 180 / numpy.pi
687 691 else:
688 692 data = numpy.abs(avgcoherenceComplex)
689 693
690 694 z.append(data)
691 695
692 696 return numpy.array(z)
693 697
694 698 def setValue(self, value):
695 699
696 700 print("This property should not be initialized", value)
697 701
698 702 return
699 703
700 704 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
701 705
702 706
703 707 class SpectraHeis(Spectra):
704 708
705 709 def __init__(self):
706 710
707 711 self.radarControllerHeaderObj = RadarControllerHeader()
708 712 self.systemHeaderObj = SystemHeader()
709 713 self.type = "SpectraHeis"
710 714 self.nProfiles = None
711 715 self.heightList = None
712 716 self.channelList = None
713 717 self.flagNoData = True
714 718 self.flagDiscontinuousBlock = False
715 719 self.utctime = None
716 720 self.blocksize = None
717 721 self.profileIndex = 0
718 722 self.nCohInt = 1
719 723 self.nIncohInt = 1
720 724
721 725 @property
722 726 def normFactor(self):
723 727 pwcode = 1
724 728 if self.flagDecodeData:
725 729 pwcode = numpy.sum(self.code[0]**2)
726 730
727 731 normFactor = self.nIncohInt * self.nCohInt * pwcode
728 732
729 733 return normFactor
730 734
731 735 @property
732 736 def timeInterval(self):
733 737
734 738 return self.ippSeconds * self.nCohInt * self.nIncohInt
735 739
736 740
737 741 class Fits(JROData):
738 742
739 743 def __init__(self):
740 744
741 745 self.type = "Fits"
742 746 self.nProfiles = None
743 747 self.heightList = None
744 748 self.channelList = None
745 749 self.flagNoData = True
746 750 self.utctime = None
747 751 self.nCohInt = 1
748 752 self.nIncohInt = 1
749 753 self.useLocalTime = True
750 754 self.profileIndex = 0
751 755 self.timeZone = 0
752 756
753 757 def getTimeRange(self):
754 758
755 759 datatime = []
756 760
757 761 datatime.append(self.ltctime)
758 762 datatime.append(self.ltctime + self.timeInterval)
759 763
760 764 datatime = numpy.array(datatime)
761 765
762 766 return datatime
763 767
764 768 def getChannelIndexList(self):
765 769
766 770 return list(range(self.nChannels))
767 771
768 772 def getNoise(self, type=1):
769 773
770 774
771 775 if type == 1:
772 776 noise = self.getNoisebyHildebrand()
773 777
774 778 if type == 2:
775 779 noise = self.getNoisebySort()
776 780
777 781 if type == 3:
778 782 noise = self.getNoisebyWindow()
779 783
780 784 return noise
781 785
782 786 @property
783 787 def timeInterval(self):
784 788
785 789 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
786 790
787 791 return timeInterval
788 792
789 793 @property
790 794 def ippSeconds(self):
791 795 '''
792 796 '''
793 797 return self.ipp_sec
794 798
795 799 noise = property(getNoise, "I'm the 'nHeights' property.")
796 800
797 801
798 802 class Correlation(JROData):
799 803
800 804 def __init__(self):
801 805 '''
802 806 Constructor
803 807 '''
804 808 self.radarControllerHeaderObj = RadarControllerHeader()
805 809 self.systemHeaderObj = SystemHeader()
806 810 self.type = "Correlation"
807 811 self.data = None
808 812 self.dtype = None
809 813 self.nProfiles = None
810 814 self.heightList = None
811 815 self.channelList = None
812 816 self.flagNoData = True
813 817 self.flagDiscontinuousBlock = False
814 818 self.utctime = None
815 819 self.timeZone = 0
816 820 self.dstFlag = None
817 821 self.errorCount = None
818 822 self.blocksize = None
819 823 self.flagDecodeData = False # asumo q la data no esta decodificada
820 824 self.flagDeflipData = False # asumo q la data no esta sin flip
821 825 self.pairsList = None
822 826 self.nPoints = None
823 827
824 828 def getPairsList(self):
825 829
826 830 return self.pairsList
827 831
828 832 def getNoise(self, mode=2):
829 833
830 834 indR = numpy.where(self.lagR == 0)[0][0]
831 835 indT = numpy.where(self.lagT == 0)[0][0]
832 836
833 837 jspectra0 = self.data_corr[:, :, indR, :]
834 838 jspectra = copy.copy(jspectra0)
835 839
836 840 num_chan = jspectra.shape[0]
837 841 num_hei = jspectra.shape[2]
838 842
839 843 freq_dc = jspectra.shape[1] / 2
840 844 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
841 845
842 846 if ind_vel[0] < 0:
843 847 ind_vel[list(range(0, 1))] = ind_vel[list(
844 848 range(0, 1))] + self.num_prof
845 849
846 850 if mode == 1:
847 851 jspectra[:, freq_dc, :] = (
848 852 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
849 853
850 854 if mode == 2:
851 855
852 856 vel = numpy.array([-2, -1, 1, 2])
853 857 xx = numpy.zeros([4, 4])
854 858
855 859 for fil in range(4):
856 860 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
857 861
858 862 xx_inv = numpy.linalg.inv(xx)
859 863 xx_aux = xx_inv[0, :]
860 864
861 865 for ich in range(num_chan):
862 866 yy = jspectra[ich, ind_vel, :]
863 867 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
864 868
865 869 junkid = jspectra[ich, freq_dc, :] <= 0
866 870 cjunkid = sum(junkid)
867 871
868 872 if cjunkid.any():
869 873 jspectra[ich, freq_dc, junkid.nonzero()] = (
870 874 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
871 875
872 876 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
873 877
874 878 return noise
875 879
876 880 @property
877 881 def timeInterval(self):
878 882
879 883 return self.ippSeconds * self.nCohInt * self.nProfiles
880 884
881 885 def splitFunctions(self):
882 886
883 887 pairsList = self.pairsList
884 888 ccf_pairs = []
885 889 acf_pairs = []
886 890 ccf_ind = []
887 891 acf_ind = []
888 892 for l in range(len(pairsList)):
889 893 chan0 = pairsList[l][0]
890 894 chan1 = pairsList[l][1]
891 895
892 896 # Obteniendo pares de Autocorrelacion
893 897 if chan0 == chan1:
894 898 acf_pairs.append(chan0)
895 899 acf_ind.append(l)
896 900 else:
897 901 ccf_pairs.append(pairsList[l])
898 902 ccf_ind.append(l)
899 903
900 904 data_acf = self.data_cf[acf_ind]
901 905 data_ccf = self.data_cf[ccf_ind]
902 906
903 907 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
904 908
905 909 @property
906 910 def normFactor(self):
907 911 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
908 912 acf_pairs = numpy.array(acf_pairs)
909 913 normFactor = numpy.zeros((self.nPairs, self.nHeights))
910 914
911 915 for p in range(self.nPairs):
912 916 pair = self.pairsList[p]
913 917
914 918 ch0 = pair[0]
915 919 ch1 = pair[1]
916 920
917 921 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
918 922 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
919 923 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
920 924
921 925 return normFactor
922 926
923 927
924 928 class Parameters(Spectra):
925 929
926 930 groupList = None # List of Pairs, Groups, etc
927 931 data_param = None # Parameters obtained
928 932 data_pre = None # Data Pre Parametrization
929 933 data_SNR = None # Signal to Noise Ratio
930 934 abscissaList = None # Abscissa, can be velocities, lags or time
931 935 utctimeInit = None # Initial UTC time
932 936 paramInterval = None # Time interval to calculate Parameters in seconds
933 937 useLocalTime = True
934 938 # Fitting
935 939 data_error = None # Error of the estimation
936 940 constants = None
937 941 library = None
938 942 # Output signal
939 943 outputInterval = None # Time interval to calculate output signal in seconds
940 944 data_output = None # Out signal
941 945 nAvg = None
942 946 noise_estimation = None
943 947 GauSPC = None # Fit gaussian SPC
944 948
945 949 data_outlier = None
946 950 data_vdrift = None
947 951 radarControllerHeaderTxt=None #header Controller like text
948 952 txPower = None
949 953 flagProfilesByRange = False
950 954 nProfilesByRange = None
951 955
952 956
953 957 def __init__(self):
954 958 '''
955 959 Constructor
956 960 '''
957 961 self.radarControllerHeaderObj = RadarControllerHeader()
958 962 self.systemHeaderObj = SystemHeader()
959 963 self.processingHeaderObj = ProcessingHeader()
960 964 self.type = "Parameters"
961 965 self.timeZone = 0
962 966
963 967 def getTimeRange1(self, interval):
964 968
965 969 datatime = []
966 970
967 971 if self.useLocalTime:
968 972 time1 = self.utctimeInit - self.timeZone * 60
969 973 else:
970 974 time1 = self.utctimeInit
971 975
972 976 datatime.append(time1)
973 977 datatime.append(time1 + interval)
974 978 datatime = numpy.array(datatime)
975 979
976 980 return datatime
977 981
978 982 @property
979 983 def timeInterval(self):
980 984
981 985 if hasattr(self, 'timeInterval1'):
982 986 return self.timeInterval1
983 987 else:
984 988 return self.paramInterval
985 989
986 990 def setValue(self, value):
987 991
988 992 print("This property should not be initialized")
989 993
990 994 return
991 995
992 996 def getNoise(self):
993 997
994 998 return self.spc_noise
995 999
996 1000 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
997 1001
998 1002
999 1003 class PlotterData(object):
1000 1004 '''
1001 1005 Object to hold data to be plotted
1002 1006 '''
1003 1007
1004 1008 MAXNUMX = 200
1005 1009 MAXNUMY = 200
1006 1010
1007 1011 def __init__(self, code, exp_code, localtime=True):
1008 1012
1009 1013 self.key = code
1010 1014 self.exp_code = exp_code
1011 1015 self.ready = False
1012 1016 self.flagNoData = False
1013 1017 self.localtime = localtime
1014 1018 self.data = {}
1015 1019 self.meta = {}
1016 1020 self.__heights = []
1017 1021
1018 1022 def __str__(self):
1019 1023 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1020 1024 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1021 1025
1022 1026 def __len__(self):
1023 1027 return len(self.data)
1024 1028
1025 1029 def __getitem__(self, key):
1026 1030 if isinstance(key, int):
1027 1031 return self.data[self.times[key]]
1028 1032 elif isinstance(key, str):
1029 1033 ret = numpy.array([self.data[x][key] for x in self.times])
1030 1034 if ret.ndim > 1:
1031 1035 ret = numpy.swapaxes(ret, 0, 1)
1032 1036 return ret
1033 1037
1034 1038 def __contains__(self, key):
1035 1039 return key in self.data[self.min_time]
1036 1040
1037 1041 def setup(self):
1038 1042 '''
1039 1043 Configure object
1040 1044 '''
1041 1045 self.type = ''
1042 1046 self.ready = False
1043 1047 del self.data
1044 1048 self.data = {}
1045 1049 self.__heights = []
1046 1050 self.__all_heights = set()
1047 1051
1048 1052 def shape(self, key):
1049 1053 '''
1050 1054 Get the shape of the one-element data for the given key
1051 1055 '''
1052 1056
1053 1057 if len(self.data[self.min_time][key]):
1054 1058 return self.data[self.min_time][key].shape
1055 1059 return (0,)
1056 1060
1057 1061 def update(self, data, tm, meta={}):
1058 1062 '''
1059 1063 Update data object with new dataOut
1060 1064 '''
1061 1065
1062 1066 self.data[tm] = data
1063 1067
1064 1068 for key, value in meta.items():
1065 1069 setattr(self, key, value)
1066 1070
1067 1071 def normalize_heights(self):
1068 1072 '''
1069 1073 Ensure same-dimension of the data for different heighList
1070 1074 '''
1071 1075
1072 1076 H = numpy.array(list(self.__all_heights))
1073 1077 H.sort()
1074 1078 for key in self.data:
1075 1079 shape = self.shape(key)[:-1] + H.shape
1076 1080 for tm, obj in list(self.data[key].items()):
1077 1081 h = self.__heights[self.times.tolist().index(tm)]
1078 1082 if H.size == h.size:
1079 1083 continue
1080 1084 index = numpy.where(numpy.in1d(H, h))[0]
1081 1085 dummy = numpy.zeros(shape) + numpy.nan
1082 1086 if len(shape) == 2:
1083 1087 dummy[:, index] = obj
1084 1088 else:
1085 1089 dummy[index] = obj
1086 1090 self.data[key][tm] = dummy
1087 1091
1088 1092 self.__heights = [H for tm in self.times]
1089 1093
1090 1094 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1091 1095 '''
1092 1096 Convert data to json
1093 1097 '''
1094 1098
1095 1099 meta = {}
1096 1100 meta['xrange'] = []
1097 1101 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1098 1102 tmp = self.data[tm][self.key]
1099 1103 shape = tmp.shape
1100 1104 if len(shape) == 2:
1101 1105 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1102 1106 elif len(shape) == 3:
1103 1107 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1104 1108 data = self.roundFloats(
1105 1109 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1106 1110 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1107 1111 else:
1108 1112 data = self.roundFloats(self.data[tm][self.key].tolist())
1109 1113
1110 1114 ret = {
1111 1115 'plot': plot_name,
1112 1116 'code': self.exp_code,
1113 1117 'time': float(tm),
1114 1118 'data': data,
1115 1119 }
1116 1120 meta['type'] = plot_type
1117 1121 meta['interval'] = float(self.interval)
1118 1122 meta['localtime'] = self.localtime
1119 1123 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1120 1124 meta.update(self.meta)
1121 1125 ret['metadata'] = meta
1122 1126 return json.dumps(ret)
1123 1127
1124 1128 @property
1125 1129 def times(self):
1126 1130 '''
1127 1131 Return the list of times of the current data
1128 1132 '''
1129 1133
1130 1134 ret = [t for t in self.data]
1131 1135 ret.sort()
1132 1136 return numpy.array(ret)
1133 1137
1134 1138 @property
1135 1139 def min_time(self):
1136 1140 '''
1137 1141 Return the minimun time value
1138 1142 '''
1139 1143
1140 1144 return self.times[0]
1141 1145
1142 1146 @property
1143 1147 def max_time(self):
1144 1148 '''
1145 1149 Return the maximun time value
1146 1150 '''
1147 1151
1148 1152 return self.times[-1]
1149 1153
1150 1154 # @property
1151 1155 # def heights(self):
1152 1156 # '''
1153 1157 # Return the list of heights of the current data
1154 1158 # '''
1155 1159
1156 1160 # return numpy.array(self.__heights[-1])
1157 1161
1158 1162 @staticmethod
1159 1163 def roundFloats(obj):
1160 1164 if isinstance(obj, list):
1161 1165 return list(map(PlotterData.roundFloats, obj))
1162 1166 elif isinstance(obj, float):
1163 1167 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,437 +1,437
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
7 7 from schainpy.utils import log
8 8
9 9 EARTH_RADIUS = 6.3710e3
10 10
11 11
12 12 def ll2xy(lat1, lon1, lat2, lon2):
13 13
14 14 p = 0.017453292519943295
15 15 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
16 16 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
17 17 r = 12742 * numpy.arcsin(numpy.sqrt(a))
18 18 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
19 19 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
20 20 theta = -theta + numpy.pi/2
21 21 return r*numpy.cos(theta), r*numpy.sin(theta)
22 22
23 23
24 24 def km2deg(km):
25 25 '''
26 26 Convert distance in km to degrees
27 27 '''
28 28
29 29 return numpy.rad2deg(km/EARTH_RADIUS)
30 30
31 31
32 32
33 33 class SpectralMomentsPlot(SpectraPlot):
34 34 '''
35 35 Plot for Spectral Moments
36 36 '''
37 37 CODE = 'spc_moments'
38 38 # colormap = 'jet'
39 39 # plot_type = 'pcolor'
40 40
41 41 class DobleGaussianPlot(SpectraPlot):
42 42 '''
43 43 Plot for Double Gaussian Plot
44 44 '''
45 45 CODE = 'gaussian_fit'
46 46 # colormap = 'jet'
47 47 # plot_type = 'pcolor'
48 48
49 49 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
50 50 '''
51 51 Plot SpectraCut with Double Gaussian Fit
52 52 '''
53 53 CODE = 'cut_gaussian_fit'
54 54
55 55
56 56 class SpectralFitObliquePlot(SpectraPlot):
57 57 '''
58 58 Plot for Spectral Oblique
59 59 '''
60 60 CODE = 'spc_moments'
61 61 colormap = 'jet'
62 62 plot_type = 'pcolor'
63 63
64 64
65 65 class SnrPlot(RTIPlot):
66 66 '''
67 67 Plot for SNR Data
68 68 '''
69 69
70 70 CODE = 'snr'
71 71 colormap = 'jet'
72 72
73 73 def update(self, dataOut):
74 74 if len(self.channelList) == 0:
75 75 self.update_list(dataOut)
76 76
77 77 meta = {}
78 78 data = {
79 79 'snr': 10 * numpy.log10(dataOut.data_snr)
80 80 }
81 81 return data, meta
82 82
83 83 class DopplerPlot(RTIPlot):
84 84 '''
85 85 Plot for DOPPLER Data (1st moment)
86 86 '''
87 87
88 88 CODE = 'dop'
89 89 colormap = 'RdBu_r'
90 90
91 91 def update(self, dataOut):
92 92 self.update_list(dataOut)
93 93 data = {
94 94 'dop': dataOut.data_dop
95 95 }
96 96
97 97 return data, {}
98 98
99 99 class PowerPlot(RTIPlot):
100 100 '''
101 101 Plot for Power Data (0 moment)
102 102 '''
103 103
104 104 CODE = 'pow'
105 105 colormap = 'jet'
106 106
107 107 def update(self, dataOut):
108 108 self.update_list(dataOut)
109 109 data = {
110 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
110 'pow': 10*numpy.log10(dataOut.data_pow)
111 111 }
112 112 try:
113 113 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
114 114 except:
115 115 pass
116 116 return data, {}
117 117
118 118 class SpectralWidthPlot(RTIPlot):
119 119 '''
120 120 Plot for Spectral Width Data (2nd moment)
121 121 '''
122 122
123 123 CODE = 'width'
124 124 colormap = 'jet'
125 125
126 126 def update(self, dataOut):
127 127 self.update_list(dataOut)
128 128 data = {
129 129 'width': dataOut.data_width
130 130 }
131 131 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
132 132 return data, {}
133 133
134 134 class SkyMapPlot(Plot):
135 135 '''
136 136 Plot for meteors detection data
137 137 '''
138 138
139 139 CODE = 'param'
140 140
141 141 def setup(self):
142 142
143 143 self.ncols = 1
144 144 self.nrows = 1
145 145 self.width = 7.2
146 146 self.height = 7.2
147 147 self.nplots = 1
148 148 self.xlabel = 'Zonal Zenith Angle (deg)'
149 149 self.ylabel = 'Meridional Zenith Angle (deg)'
150 150 self.polar = True
151 151 self.ymin = -180
152 152 self.ymax = 180
153 153 self.colorbar = False
154 154
155 155 def plot(self):
156 156
157 157 arrayParameters = numpy.concatenate(self.data['param'])
158 158 error = arrayParameters[:, -1]
159 159 indValid = numpy.where(error == 0)[0]
160 160 finalMeteor = arrayParameters[indValid, :]
161 161 finalAzimuth = finalMeteor[:, 3]
162 162 finalZenith = finalMeteor[:, 4]
163 163
164 164 x = finalAzimuth * numpy.pi / 180
165 165 y = finalZenith
166 166
167 167 ax = self.axes[0]
168 168
169 169 if ax.firsttime:
170 170 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
171 171 else:
172 172 ax.plot.set_data(x, y)
173 173
174 174 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
175 175 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
176 176 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
177 177 dt2,
178 178 len(x))
179 179 self.titles[0] = title
180 180
181 181 class GenericRTIPlot(Plot):
182 182 '''
183 183 Plot for data_xxxx object
184 184 '''
185 185
186 186 CODE = 'param'
187 187 colormap = 'viridis'
188 188 plot_type = 'pcolorbuffer'
189 189
190 190 def setup(self):
191 191 self.xaxis = 'time'
192 192 self.ncols = 1
193 193 self.nrows = self.data.shape('param')[0]
194 194 self.nplots = self.nrows
195 195 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
196 196
197 197 if not self.xlabel:
198 198 self.xlabel = 'Time'
199 199
200 200 self.ylabel = 'Range [km]'
201 201 if not self.titles:
202 202 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
203 203
204 204 def update(self, dataOut):
205 205
206 206 data = {
207 207 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
208 208 }
209 209
210 210 meta = {}
211 211
212 212 return data, meta
213 213
214 214 def plot(self):
215 215 # self.data.normalize_heights()
216 216 self.x = self.data.times
217 217 self.y = self.data.yrange
218 218 self.z = self.data['param']
219 219
220 220 self.z = numpy.ma.masked_invalid(self.z)
221 221
222 222 if self.decimation is None:
223 223 x, y, z = self.fill_gaps(self.x, self.y, self.z)
224 224 else:
225 225 x, y, z = self.fill_gaps(*self.decimate())
226 226
227 227 for n, ax in enumerate(self.axes):
228 228
229 229 self.zmax = self.zmax if self.zmax is not None else numpy.max(
230 230 self.z[n])
231 231 self.zmin = self.zmin if self.zmin is not None else numpy.min(
232 232 self.z[n])
233 233
234 234 if ax.firsttime:
235 235 if self.zlimits is not None:
236 236 self.zmin, self.zmax = self.zlimits[n]
237 237
238 238 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
239 239 vmin=self.zmin,
240 240 vmax=self.zmax,
241 241 cmap=self.cmaps[n]
242 242 )
243 243 else:
244 244 if self.zlimits is not None:
245 245 self.zmin, self.zmax = self.zlimits[n]
246 246 try:
247 247 ax.collections.remove(ax.collections[0])
248 248 except:
249 249 pass
250 250 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
251 251 vmin=self.zmin,
252 252 vmax=self.zmax,
253 253 cmap=self.cmaps[n]
254 254 )
255 255
256 256
257 257 class PolarMapPlot(Plot):
258 258 '''
259 259 Plot for weather radar
260 260 '''
261 261
262 262 CODE = 'param'
263 263 colormap = 'seismic'
264 264
265 265 def setup(self):
266 266 self.ncols = 1
267 267 self.nrows = 1
268 268 self.width = 9
269 269 self.height = 8
270 270 self.mode = self.data.meta['mode']
271 271 if self.channels is not None:
272 272 self.nplots = len(self.channels)
273 273 self.nrows = len(self.channels)
274 274 else:
275 275 self.nplots = self.data.shape(self.CODE)[0]
276 276 self.nrows = self.nplots
277 277 self.channels = list(range(self.nplots))
278 278 if self.mode == 'E':
279 279 self.xlabel = 'Longitude'
280 280 self.ylabel = 'Latitude'
281 281 else:
282 282 self.xlabel = 'Range (km)'
283 283 self.ylabel = 'Height (km)'
284 284 self.bgcolor = 'white'
285 285 self.cb_labels = self.data.meta['units']
286 286 self.lat = self.data.meta['latitude']
287 287 self.lon = self.data.meta['longitude']
288 288 self.xmin, self.xmax = float(
289 289 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
290 290 self.ymin, self.ymax = float(
291 291 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
292 292 # self.polar = True
293 293
294 294 def plot(self):
295 295
296 296 for n, ax in enumerate(self.axes):
297 297 data = self.data['param'][self.channels[n]]
298 298
299 299 zeniths = numpy.linspace(
300 300 0, self.data.meta['max_range'], data.shape[1])
301 301 if self.mode == 'E':
302 302 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
303 303 r, theta = numpy.meshgrid(zeniths, azimuths)
304 304 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
305 305 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
306 306 x = km2deg(x) + self.lon
307 307 y = km2deg(y) + self.lat
308 308 else:
309 309 azimuths = numpy.radians(self.data.yrange)
310 310 r, theta = numpy.meshgrid(zeniths, azimuths)
311 311 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
312 312 self.y = zeniths
313 313
314 314 if ax.firsttime:
315 315 if self.zlimits is not None:
316 316 self.zmin, self.zmax = self.zlimits[n]
317 317 ax.plt = ax.pcolormesh(# r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
318 318 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
319 319 vmin=self.zmin,
320 320 vmax=self.zmax,
321 321 cmap=self.cmaps[n])
322 322 else:
323 323 if self.zlimits is not None:
324 324 self.zmin, self.zmax = self.zlimits[n]
325 325 ax.collections.remove(ax.collections[0])
326 326 ax.plt = ax.pcolormesh(# r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
327 327 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
328 328 vmin=self.zmin,
329 329 vmax=self.zmax,
330 330 cmap=self.cmaps[n])
331 331
332 332 if self.mode == 'A':
333 333 continue
334 334
335 335 # plot district names
336 336 f = open('/data/workspace/schain_scripts/distrito.csv')
337 337 for line in f:
338 338 label, lon, lat = [s.strip() for s in line.split(',') if s]
339 339 lat = float(lat)
340 340 lon = float(lon)
341 341 # ax.plot(lon, lat, '.b', ms=2)
342 342 ax.text(lon, lat, label.decode('utf8'), ha='center',
343 343 va='bottom', size='8', color='black')
344 344
345 345 # plot limites
346 346 limites = []
347 347 tmp = []
348 348 for line in open('/data/workspace/schain_scripts/lima.csv'):
349 349 if '#' in line:
350 350 if tmp:
351 351 limites.append(tmp)
352 352 tmp = []
353 353 continue
354 354 values = line.strip().split(',')
355 355 tmp.append((float(values[0]), float(values[1])))
356 356 for points in limites:
357 357 ax.add_patch(
358 358 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
359 359
360 360 # plot Cuencas
361 361 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
362 362 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
363 363 values = [line.strip().split(',') for line in f]
364 364 points = [(float(s[0]), float(s[1])) for s in values]
365 365 ax.add_patch(Polygon(points, ec='b', fc='none'))
366 366
367 367 # plot grid
368 368 for r in (15, 30, 45, 60):
369 369 ax.add_artist(plt.Circle((self.lon, self.lat),
370 370 km2deg(r), color='0.6', fill=False, lw=0.2))
371 371 ax.text(
372 372 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
373 373 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
374 374 '{}km'.format(r),
375 375 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
376 376
377 377 if self.mode == 'E':
378 378 title = 'El={}\N{DEGREE SIGN}'.format(self.data.meta['elevation'])
379 379 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
380 380 else:
381 381 title = 'Az={}\N{DEGREE SIGN}'.format(self.data.meta['azimuth'])
382 382 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
383 383
384 384 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
385 385 self.titles = ['{} {}'.format(
386 386 self.data.parameters[x], title) for x in self.channels]
387 387
388 388
389 389
390 390 class TxPowerPlot(Plot):
391 391 '''
392 392 Plot for TX Power from external file
393 393 '''
394 394
395 395 CODE = 'tx_power'
396 396 plot_type = 'scatterbuffer'
397 397
398 398 def setup(self):
399 399 self.xaxis = 'time'
400 400 self.ncols = 1
401 401 self.nrows = 1
402 402 self.nplots = 1
403 403 self.ylabel = 'Power [kW]'
404 404 self.xlabel = 'Time'
405 405 self.titles = ['TX power']
406 406 self.colorbar = False
407 407 self.plots_adjust.update({'right': 0.85 })
408 408 #if not self.titles:
409 409 self.titles = ['TX Power Plot']
410 410
411 411 def update(self, dataOut):
412 412
413 413 data = {}
414 414 meta = {}
415 415
416 416 data['tx_power'] = dataOut.txPower/1000
417 417 meta['yrange'] = numpy.array([])
418 418 #print(dataOut.txPower/1000)
419 419 return data, meta
420 420
421 421 def plot(self):
422 422
423 423 x = self.data.times
424 424 xmin = self.data.min_time
425 425 xmax = xmin + self.xrange * 60 * 60
426 426 Y = self.data['tx_power']
427 427
428 428 if self.axes[0].firsttime:
429 429 if self.ymin is None: self.ymin = 0
430 430 if self.ymax is None: self.ymax = numpy.nanmax(Y) + 5
431 431 if self.ymax == 5:
432 432 self.ymax = 250
433 433 self.ymin = 100
434 434 self.axes[0].plot(x, Y, lw=1, label='Power')
435 435 plt.legend(bbox_to_anchor=(1.18, 1.0))
436 436 else:
437 437 self.axes[0].lines[0].set_data(x, Y) No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1935 +1,1934
1 1 # Copyright (c) 2012-2021 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 import datetime
12 12
13 13 from schainpy.model.graphics.jroplot_base import Plot, plt, log
14 14 from itertools import combinations
15 15 from matplotlib.ticker import LinearLocator
16 16
17 17 from schainpy.model.utils.BField import BField
18 18 from scipy.interpolate import splrep
19 19 from scipy.interpolate import splev
20 20
21 21 from matplotlib import __version__ as plt_version
22 22
23 23 if plt_version >='3.3.4':
24 24 EXTRA_POINTS = 0
25 25 else:
26 26 EXTRA_POINTS = 1
27 27 class SpectraPlot(Plot):
28 28 '''
29 29 Plot for Spectra data
30 30 '''
31 31
32 32 CODE = 'spc'
33 33 colormap = 'jet'
34 34 plot_type = 'pcolor'
35 35 buffering = False
36 36 channelList = []
37 37 elevationList = []
38 38 azimuthList = []
39 39
40 40 def setup(self):
41 41
42 42 self.nplots = len(self.data.channels)
43 43 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
44 44 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
45 45 self.height = 3.4 * self.nrows
46 46 self.cb_label = 'dB'
47 47 if self.showprofile:
48 48 self.width = 5.2 * self.ncols
49 49 else:
50 50 self.width = 4.2* self.ncols
51 51 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.12})
52 52 self.ylabel = 'Range [km]'
53 53
54 54 def update_list(self,dataOut):
55 55
56 56 if len(self.channelList) == 0:
57 57 self.channelList = dataOut.channelList
58 58 if len(self.elevationList) == 0:
59 59 self.elevationList = dataOut.elevationList
60 60 if len(self.azimuthList) == 0:
61 61 self.azimuthList = dataOut.azimuthList
62 62
63 63 def update(self, dataOut):
64 64
65 65 self.update_list(dataOut)
66 66 data = {}
67 67 meta = {}
68 68 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
69 69 if dataOut.type == "Parameters":
70 70 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
71 71 spc = 10*numpy.log10(dataOut.data_spc/(dataOut.nProfiles))
72 72 else:
73 73 noise = 10*numpy.log10(dataOut.getNoise()/norm)
74 74
75 75 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
76 76 for ch in range(dataOut.nChannels):
77 77 if hasattr(dataOut.normFactor,'ndim'):
78 78 if dataOut.normFactor.ndim > 1:
79 79 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
80 80
81 81 else:
82 82 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
83 83 else:
84 84 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
85 85
86 86 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
87 87 spc = 10*numpy.log10(z)
88 88
89 89 data['spc'] = spc
90 90 data['rti'] = spc.mean(axis=1)
91 91 data['noise'] = noise
92 92 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
93 93 if self.CODE == 'spc_moments':
94 94 data['moments'] = dataOut.moments
95 95
96 96 return data, meta
97 97
98 98 def plot(self):
99 99
100 100 if self.xaxis == "frequency":
101 101 x = self.data.xrange[0]
102 102 self.xlabel = "Frequency (kHz)"
103 103 elif self.xaxis == "time":
104 104 x = self.data.xrange[1]
105 105 self.xlabel = "Time (ms)"
106 106 else:
107 107 x = self.data.xrange[2]
108 108 self.xlabel = "Velocity (m/s)"
109 109
110 110 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
111 111 x = self.data.xrange[2]
112 112 self.xlabel = "Velocity (m/s)"
113 113
114 114 self.titles = []
115 115
116 116 y = self.data.yrange
117 117 self.y = y
118 118
119 119 data = self.data[-1]
120 120 z = data['spc']
121 121
122 122 for n, ax in enumerate(self.axes):
123 123 noise = self.data['noise'][n][0]
124 124 # noise = data['noise'][n]
125 125
126 126 if self.CODE == 'spc_moments':
127 127 mean = data['moments'][n, 1]
128 128 if self.CODE == 'gaussian_fit':
129 129 gau0 = data['gaussfit'][n][2,:,0]
130 130 gau1 = data['gaussfit'][n][2,:,1]
131 131 if ax.firsttime:
132 132 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
133 133 self.xmin = self.xmin if self.xmin else -self.xmax
134 134 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
135 135 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
136 136 ax.plt = ax.pcolormesh(x, y, z[n].T,
137 137 vmin=self.zmin,
138 138 vmax=self.zmax,
139 139 cmap=plt.get_cmap(self.colormap)
140 140 )
141 141
142 142 if self.showprofile:
143 143 ax.plt_profile = self.pf_axes[n].plot(
144 144 data['rti'][n], y)[0]
145 145 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
146 146 color="k", linestyle="dashed", lw=1)[0]
147 147 if self.CODE == 'spc_moments':
148 148 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
149 149 if self.CODE == 'gaussian_fit':
150 150 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
151 151 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
152 152 else:
153 153 ax.plt.set_array(z[n].T.ravel())
154 154 if self.showprofile:
155 155 ax.plt_profile.set_data(data['rti'][n], y)
156 156 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
157 157 if self.CODE == 'spc_moments':
158 158 ax.plt_mean.set_data(mean, y)
159 159 if self.CODE == 'gaussian_fit':
160 160 ax.plt_gau0.set_data(gau0, y)
161 161 ax.plt_gau1.set_data(gau1, y)
162 162 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
163 163 self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
164 164 else:
165 165 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
166 166
167 167 class SpectraObliquePlot(Plot):
168 168 '''
169 169 Plot for Spectra data
170 170 '''
171 171
172 172 CODE = 'spc_oblique'
173 173 colormap = 'jet'
174 174 plot_type = 'pcolor'
175 175
176 176 def setup(self):
177 177 self.xaxis = "oblique"
178 178 self.nplots = len(self.data.channels)
179 179 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
180 180 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
181 181 self.height = 2.6 * self.nrows
182 182 self.cb_label = 'dB'
183 183 if self.showprofile:
184 184 self.width = 4 * self.ncols
185 185 else:
186 186 self.width = 3.5 * self.ncols
187 187 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
188 188 self.ylabel = 'Range [km]'
189 189
190 190 def update(self, dataOut):
191 191
192 192 data = {}
193 193 meta = {}
194 194
195 195 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
196 196 data['spc'] = spc
197 197 data['rti'] = dataOut.getPower()
198 198 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
199 199 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
200 200
201 201 data['shift1'] = dataOut.Dop_EEJ_T1[0]
202 202 data['shift2'] = dataOut.Dop_EEJ_T2[0]
203 203 data['max_val_2'] = dataOut.Oblique_params[0,-1,:]
204 204 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
205 205 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
206 206
207 207 return data, meta
208 208
209 209 def plot(self):
210 210
211 211 if self.xaxis == "frequency":
212 212 x = self.data.xrange[0]
213 213 self.xlabel = "Frequency (kHz)"
214 214 elif self.xaxis == "time":
215 215 x = self.data.xrange[1]
216 216 self.xlabel = "Time (ms)"
217 217 else:
218 218 x = self.data.xrange[2]
219 219 self.xlabel = "Velocity (m/s)"
220 220
221 221 self.titles = []
222 222
223 223 y = self.data.yrange
224 224 self.y = y
225 225
226 226 data = self.data[-1]
227 227 z = data['spc']
228 228
229 229 for n, ax in enumerate(self.axes):
230 230 noise = self.data['noise'][n][-1]
231 231 shift1 = data['shift1']
232 232 shift2 = data['shift2']
233 233 max_val_2 = data['max_val_2']
234 234 err1 = data['shift1_error']
235 235 err2 = data['shift2_error']
236 236 if ax.firsttime:
237 237
238 238 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
239 239 self.xmin = self.xmin if self.xmin else -self.xmax
240 240 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
241 241 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
242 242 ax.plt = ax.pcolormesh(x, y, z[n].T,
243 243 vmin=self.zmin,
244 244 vmax=self.zmax,
245 245 cmap=plt.get_cmap(self.colormap)
246 246 )
247 247
248 248 if self.showprofile:
249 249 ax.plt_profile = self.pf_axes[n].plot(
250 250 self.data['rti'][n][-1], y)[0]
251 251 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
252 252 color="k", linestyle="dashed", lw=1)[0]
253 253
254 254 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
255 255 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
256 256 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
257 257
258 258 else:
259 259 self.ploterr1.remove()
260 260 self.ploterr2.remove()
261 261 self.ploterr3.remove()
262 262 ax.plt.set_array(z[n].T.ravel())
263 263 if self.showprofile:
264 264 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
265 265 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
266 266 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
267 267 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
268 268 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
269 269
270 270 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
271 271
272 272
273 273 class CrossSpectraPlot(Plot):
274 274
275 275 CODE = 'cspc'
276 276 colormap = 'jet'
277 277 plot_type = 'pcolor'
278 278 zmin_coh = None
279 279 zmax_coh = None
280 280 zmin_phase = None
281 281 zmax_phase = None
282 282 realChannels = None
283 283 crossPairs = None
284 284
285 285 def setup(self):
286 286
287 287 self.ncols = 4
288 288 self.nplots = len(self.data.pairs) * 2
289 289 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
290 290 self.width = 3.1 * self.ncols
291 291 self.height = 2.6 * self.nrows
292 292 self.ylabel = 'Range [km]'
293 293 self.showprofile = False
294 294 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
295 295
296 296 def update(self, dataOut):
297 297
298 298 data = {}
299 299 meta = {}
300 300
301 301 spc = dataOut.data_spc
302 302 cspc = dataOut.data_cspc
303 303 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
304 304 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
305 305 meta['pairs'] = rawPairs
306 306 if self.crossPairs == None:
307 307 self.crossPairs = dataOut.pairsList
308 308 tmp = []
309 309
310 310 for n, pair in enumerate(meta['pairs']):
311 311 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
312 312 coh = numpy.abs(out)
313 313 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
314 314 tmp.append(coh)
315 315 tmp.append(phase)
316 316
317 317 data['cspc'] = numpy.array(tmp)
318 318
319 319 return data, meta
320 320
321 321 def plot(self):
322 322
323 323 if self.xaxis == "frequency":
324 324 x = self.data.xrange[0]
325 325 self.xlabel = "Frequency (kHz)"
326 326 elif self.xaxis == "time":
327 327 x = self.data.xrange[1]
328 328 self.xlabel = "Time (ms)"
329 329 else:
330 330 x = self.data.xrange[2]
331 331 self.xlabel = "Velocity (m/s)"
332 332
333 333 self.titles = []
334 334
335 335 y = self.data.yrange
336 336 self.y = y
337 337
338 338 data = self.data[-1]
339 339 cspc = data['cspc']
340 340
341 341 for n in range(len(self.data.pairs)):
342 342 pair = self.crossPairs[n]
343 343 coh = cspc[n*2]
344 344 phase = cspc[n*2+1]
345 345 ax = self.axes[2 * n]
346 346 if ax.firsttime:
347 347 ax.plt = ax.pcolormesh(x, y, coh.T,
348 348 vmin=self.zmin_coh,
349 349 vmax=self.zmax_coh,
350 350 cmap=plt.get_cmap(self.colormap_coh)
351 351 )
352 352 else:
353 353 ax.plt.set_array(coh.T.ravel())
354 354 self.titles.append(
355 355 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
356 356
357 357 ax = self.axes[2 * n + 1]
358 358 if ax.firsttime:
359 359 ax.plt = ax.pcolormesh(x, y, phase.T,
360 360 vmin=-180,
361 361 vmax=180,
362 362 cmap=plt.get_cmap(self.colormap_phase)
363 363 )
364 364 else:
365 365 ax.plt.set_array(phase.T.ravel())
366 366 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
367 367
368 368
369 369 class CrossSpectra4Plot(Plot):
370 370
371 371 CODE = 'cspc'
372 372 colormap = 'jet'
373 373 plot_type = 'pcolor'
374 374 zmin_coh = None
375 375 zmax_coh = None
376 376 zmin_phase = None
377 377 zmax_phase = None
378 378
379 379 def setup(self):
380 380
381 381 self.ncols = 4
382 382 self.nrows = len(self.data.pairs)
383 383 self.nplots = self.nrows * 4
384 384 self.width = 3.1 * self.ncols
385 385 self.height = 5 * self.nrows
386 386 self.ylabel = 'Range [km]'
387 387 self.showprofile = False
388 388 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
389 389
390 390 def plot(self):
391 391
392 392 if self.xaxis == "frequency":
393 393 x = self.data.xrange[0]
394 394 self.xlabel = "Frequency (kHz)"
395 395 elif self.xaxis == "time":
396 396 x = self.data.xrange[1]
397 397 self.xlabel = "Time (ms)"
398 398 else:
399 399 x = self.data.xrange[2]
400 400 self.xlabel = "Velocity (m/s)"
401 401
402 402 self.titles = []
403 403
404 404
405 405 y = self.data.heights
406 406 self.y = y
407 407 nspc = self.data['spc']
408 408 spc = self.data['cspc'][0]
409 409 cspc = self.data['cspc'][1]
410 410
411 411 for n in range(self.nrows):
412 412 noise = self.data['noise'][:,-1]
413 413 pair = self.data.pairs[n]
414 414
415 415 ax = self.axes[4 * n]
416 416 if ax.firsttime:
417 417 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
418 418 self.xmin = self.xmin if self.xmin else -self.xmax
419 419 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
420 420 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
421 421 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
422 422 vmin=self.zmin,
423 423 vmax=self.zmax,
424 424 cmap=plt.get_cmap(self.colormap)
425 425 )
426 426 else:
427 427
428 428 ax.plt.set_array(nspc[pair[0]].T.ravel())
429 429 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
430 430
431 431 ax = self.axes[4 * n + 1]
432 432
433 433 if ax.firsttime:
434 434 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
435 435 vmin=self.zmin,
436 436 vmax=self.zmax,
437 437 cmap=plt.get_cmap(self.colormap)
438 438 )
439 439 else:
440 440
441 441 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
442 442 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
443 443
444 444 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
445 445 coh = numpy.abs(out)
446 446 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
447 447
448 448 ax = self.axes[4 * n + 2]
449 449 if ax.firsttime:
450 450 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
451 451 vmin=0,
452 452 vmax=1,
453 453 cmap=plt.get_cmap(self.colormap_coh)
454 454 )
455 455 else:
456 456 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
457 457 self.titles.append(
458 458 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
459 459
460 460 ax = self.axes[4 * n + 3]
461 461 if ax.firsttime:
462 462 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
463 463 vmin=-180,
464 464 vmax=180,
465 465 cmap=plt.get_cmap(self.colormap_phase)
466 466 )
467 467 else:
468 468 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
469 469 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
470 470
471 471
472 472 class CrossSpectra2Plot(Plot):
473 473
474 474 CODE = 'cspc'
475 475 colormap = 'jet'
476 476 plot_type = 'pcolor'
477 477 zmin_coh = None
478 478 zmax_coh = None
479 479 zmin_phase = None
480 480 zmax_phase = None
481 481
482 482 def setup(self):
483 483
484 484 self.ncols = 1
485 485 self.nrows = len(self.data.pairs)
486 486 self.nplots = self.nrows * 1
487 487 self.width = 3.1 * self.ncols
488 488 self.height = 5 * self.nrows
489 489 self.ylabel = 'Range [km]'
490 490 self.showprofile = False
491 491 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
492 492
493 493 def plot(self):
494 494
495 495 if self.xaxis == "frequency":
496 496 x = self.data.xrange[0]
497 497 self.xlabel = "Frequency (kHz)"
498 498 elif self.xaxis == "time":
499 499 x = self.data.xrange[1]
500 500 self.xlabel = "Time (ms)"
501 501 else:
502 502 x = self.data.xrange[2]
503 503 self.xlabel = "Velocity (m/s)"
504 504
505 505 self.titles = []
506 506
507 507
508 508 y = self.data.heights
509 509 self.y = y
510 510 cspc = self.data['cspc'][1]
511 511
512 512 for n in range(self.nrows):
513 513 noise = self.data['noise'][:,-1]
514 514 pair = self.data.pairs[n]
515 515 out = cspc[n]
516 516 cross = numpy.abs(out)
517 517 z = cross/self.data.nFactor
518 518 cross = 10*numpy.log10(z)
519 519
520 520 ax = self.axes[1 * n]
521 521 if ax.firsttime:
522 522 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
523 523 self.xmin = self.xmin if self.xmin else -self.xmax
524 524 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
525 525 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
526 526 ax.plt = ax.pcolormesh(x, y, cross.T,
527 527 vmin=self.zmin,
528 528 vmax=self.zmax,
529 529 cmap=plt.get_cmap(self.colormap)
530 530 )
531 531 else:
532 532 ax.plt.set_array(cross.T.ravel())
533 533 self.titles.append(
534 534 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
535 535
536 536
537 537 class CrossSpectra3Plot(Plot):
538 538
539 539 CODE = 'cspc'
540 540 colormap = 'jet'
541 541 plot_type = 'pcolor'
542 542 zmin_coh = None
543 543 zmax_coh = None
544 544 zmin_phase = None
545 545 zmax_phase = None
546 546
547 547 def setup(self):
548 548
549 549 self.ncols = 3
550 550 self.nrows = len(self.data.pairs)
551 551 self.nplots = self.nrows * 3
552 552 self.width = 3.1 * self.ncols
553 553 self.height = 5 * self.nrows
554 554 self.ylabel = 'Range [km]'
555 555 self.showprofile = False
556 556 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
557 557
558 558 def plot(self):
559 559
560 560 if self.xaxis == "frequency":
561 561 x = self.data.xrange[0]
562 562 self.xlabel = "Frequency (kHz)"
563 563 elif self.xaxis == "time":
564 564 x = self.data.xrange[1]
565 565 self.xlabel = "Time (ms)"
566 566 else:
567 567 x = self.data.xrange[2]
568 568 self.xlabel = "Velocity (m/s)"
569 569
570 570 self.titles = []
571 571
572 572
573 573 y = self.data.heights
574 574 self.y = y
575 575
576 576 cspc = self.data['cspc'][1]
577 577
578 578 for n in range(self.nrows):
579 579 noise = self.data['noise'][:,-1]
580 580 pair = self.data.pairs[n]
581 581 out = cspc[n]
582 582
583 583 cross = numpy.abs(out)
584 584 z = cross/self.data.nFactor
585 585 cross = 10*numpy.log10(z)
586 586
587 587 out_r= out.real/self.data.nFactor
588 588
589 589 out_i= out.imag/self.data.nFactor
590 590
591 591 ax = self.axes[3 * n]
592 592 if ax.firsttime:
593 593 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
594 594 self.xmin = self.xmin if self.xmin else -self.xmax
595 595 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
596 596 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
597 597 ax.plt = ax.pcolormesh(x, y, cross.T,
598 598 vmin=self.zmin,
599 599 vmax=self.zmax,
600 600 cmap=plt.get_cmap(self.colormap)
601 601 )
602 602 else:
603 603 ax.plt.set_array(cross.T.ravel())
604 604 self.titles.append(
605 605 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
606 606
607 607 ax = self.axes[3 * n + 1]
608 608 if ax.firsttime:
609 609 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
610 610 self.xmin = self.xmin if self.xmin else -self.xmax
611 611 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
612 612 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
613 613 ax.plt = ax.pcolormesh(x, y, out_r.T,
614 614 vmin=-1.e6,
615 615 vmax=0,
616 616 cmap=plt.get_cmap(self.colormap)
617 617 )
618 618 else:
619 619 ax.plt.set_array(out_r.T.ravel())
620 620 self.titles.append(
621 621 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
622 622
623 623 ax = self.axes[3 * n + 2]
624 624
625 625
626 626 if ax.firsttime:
627 627 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
628 628 self.xmin = self.xmin if self.xmin else -self.xmax
629 629 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
630 630 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
631 631 ax.plt = ax.pcolormesh(x, y, out_i.T,
632 632 vmin=-1.e6,
633 633 vmax=1.e6,
634 634 cmap=plt.get_cmap(self.colormap)
635 635 )
636 636 else:
637 637 ax.plt.set_array(out_i.T.ravel())
638 638 self.titles.append(
639 639 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
640 640
641 641 class RTIPlot(Plot):
642 642 '''
643 643 Plot for RTI data
644 644 '''
645 645
646 646 CODE = 'rti'
647 647 colormap = 'jet'
648 648 plot_type = 'pcolorbuffer'
649 649 titles = None
650 650 channelList = []
651 651 elevationList = []
652 652 azimuthList = []
653 653
654 654 def setup(self):
655 655 self.xaxis = 'time'
656 656 self.ncols = 1
657 657 self.nrows = len(self.data.channels)
658 658 self.nplots = len(self.data.channels)
659 659 self.ylabel = 'Range [km]'
660 660 #self.xlabel = 'Time'
661 661 self.cb_label = 'dB'
662 662 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
663 663 self.titles = ['{} Channel {}'.format(
664 664 self.CODE.upper(), x) for x in range(self.nplots)]
665 665
666 666 def update_list(self,dataOut):
667 667
668 668 if len(self.channelList) == 0:
669 669 self.channelList = dataOut.channelList
670 670 if len(self.elevationList) == 0:
671 671 self.elevationList = dataOut.elevationList
672 672 if len(self.azimuthList) == 0:
673 673 self.azimuthList = dataOut.azimuthList
674 674
675 675
676 676 def update(self, dataOut):
677 677
678 678 if len(self.channelList) == 0:
679 679 self.update_list(dataOut)
680 680 data = {}
681 681 meta = {}
682 682 data['rti'] = dataOut.getPower()
683 683 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
684 684 noise = 10*numpy.log10(dataOut.getNoise()/norm)
685 685 data['noise'] = noise
686 686
687 687 return data, meta
688 688
689 689 def plot(self):
690 690
691 691 self.x = self.data.times
692 692 self.y = self.data.yrange
693 693 self.z = self.data[self.CODE]
694 694 self.z = numpy.array(self.z, dtype=float)
695 695 self.z = numpy.ma.masked_invalid(self.z)
696 696
697 697 try:
698 698 if self.channelList != None:
699 699 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
700 700 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
701 701 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
702 702 else:
703 703 self.titles = ['{} Channel {}'.format(
704 704 self.CODE.upper(), x) for x in self.channelList]
705 705 except:
706 706 if self.channelList.any() != None:
707 707 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
708 708 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
709 709 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
710 710 else:
711 711 self.titles = ['{} Channel {}'.format(
712 712 self.CODE.upper(), x) for x in self.channelList]
713 713
714 714 if self.decimation is None:
715 715 x, y, z = self.fill_gaps(self.x, self.y, self.z)
716 716 else:
717 717 x, y, z = self.fill_gaps(*self.decimate())
718 718
719 719 for n, ax in enumerate(self.axes):
720 720
721 721 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
722 722 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
723 723 data = self.data[-1]
724 724 if ax.firsttime:
725 725 ax.plt = ax.pcolormesh(x, y, z[n].T,
726 726 vmin=self.zmin,
727 727 vmax=self.zmax,
728 728 cmap=plt.get_cmap(self.colormap)
729 729 )
730 730 if self.showprofile:
731 ax.plot_profile = self.pf_axes[n].plot(
732 data[self.CODE][n], self.y)[0]
731 ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0]
733 732 if "noise" in self.data:
733
734 734 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
735 735 color="k", linestyle="dashed", lw=1)[0]
736 736 else:
737 ax.collections.remove(ax.collections[0]) # error while running in 3.12
737 ax.collections.remove(ax.collections[0])
738 738 ax.plt = ax.pcolormesh(x, y, z[n].T,
739 739 vmin=self.zmin,
740 740 vmax=self.zmax,
741 741 cmap=plt.get_cmap(self.colormap)
742 742 )
743 743 if self.showprofile:
744 744 ax.plot_profile.set_data(data[self.CODE][n], self.y)
745 745 if "noise" in self.data:
746 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
747 color="k", linestyle="dashed", lw=1)[0]
746 ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
748 747
749 748 class SpectrogramPlot(Plot):
750 749 '''
751 750 Plot for Spectrogram data
752 751 '''
753 752
754 753 CODE = 'Spectrogram_Profile'
755 754 colormap = 'binary'
756 755 plot_type = 'pcolorbuffer'
757 756
758 757 def setup(self):
759 758 self.xaxis = 'time'
760 759 self.ncols = 1
761 760 self.nrows = len(self.data.channels)
762 761 self.nplots = len(self.data.channels)
763 762 self.xlabel = 'Time'
764 763 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
765 764 self.titles = []
766 765
767 766 self.titles = ['{} Channel {}'.format(
768 767 self.CODE.upper(), x) for x in range(self.nrows)]
769 768
770 769
771 770 def update(self, dataOut):
772 771 data = {}
773 772 meta = {}
774 773
775 774 maxHei = 1620#+12000
776 775 indb = numpy.where(dataOut.heightList <= maxHei)
777 776 hei = indb[0][-1]
778 777
779 778 factor = dataOut.nIncohInt
780 779 z = dataOut.data_spc[:,:,hei] / factor
781 780 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
782 781
783 782 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
784 783 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
785 784
786 785 data['hei'] = hei
787 786 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
788 787 data['nProfiles'] = dataOut.nProfiles
789 788
790 789 return data, meta
791 790
792 791 def plot(self):
793 792
794 793 self.x = self.data.times
795 794 self.z = self.data[self.CODE]
796 795 self.y = self.data.xrange[0]
797 796
798 797 hei = self.data['hei'][-1]
799 798 DH = self.data['DH'][-1]
800 799 nProfiles = self.data['nProfiles'][-1]
801 800
802 801 self.ylabel = "Frequency (kHz)"
803 802
804 803 self.z = numpy.ma.masked_invalid(self.z)
805 804
806 805 if self.decimation is None:
807 806 x, y, z = self.fill_gaps(self.x, self.y, self.z)
808 807 else:
809 808 x, y, z = self.fill_gaps(*self.decimate())
810 809
811 810 for n, ax in enumerate(self.axes):
812 811 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
813 812 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
814 813 data = self.data[-1]
815 814 if ax.firsttime:
816 815 ax.plt = ax.pcolormesh(x, y, z[n].T,
817 816 vmin=self.zmin,
818 817 vmax=self.zmax,
819 818 cmap=plt.get_cmap(self.colormap)
820 819 )
821 820 else:
822 # ax.collections.remove(ax.collections[0]) # error while running
821 ax.collections.remove(ax.collections[0]) # error while running
823 822 ax.plt = ax.pcolormesh(x, y, z[n].T,
824 823 vmin=self.zmin,
825 824 vmax=self.zmax,
826 825 cmap=plt.get_cmap(self.colormap)
827 826 )
828 827
829 828
830 829
831 830 class CoherencePlot(RTIPlot):
832 831 '''
833 832 Plot for Coherence data
834 833 '''
835 834
836 835 CODE = 'coh'
837 836 titles = None
838 837
839 838 def setup(self):
840 839 self.xaxis = 'time'
841 840 self.ncols = 1
842 841 self.nrows = len(self.data.pairs)
843 842 self.nplots = len(self.data.pairs)
844 843 self.ylabel = 'Range [km]'
845 844 self.xlabel = 'Time'
846 845 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
847 846 if self.CODE == 'coh':
848 847 self.cb_label = ''
849 848 self.titles = [
850 849 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
851 850 else:
852 851 self.cb_label = 'Degrees'
853 852 self.titles = [
854 853 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
855 854
856 855 def update(self, dataOut):
857 856
858 857 data = {}
859 858 meta = {}
860 859 data['coh'] = dataOut.getCoherence()
861 860 meta['pairs'] = dataOut.pairsList
862 861
863 862 return data, meta
864 863
865 864 class PhasePlot(CoherencePlot):
866 865 '''
867 866 Plot for Phase map data
868 867 '''
869 868
870 869 CODE = 'phase'
871 870 colormap = 'seismic'
872 871
873 872 def update(self, dataOut):
874 873
875 874 data = {}
876 875 meta = {}
877 876 data['phase'] = dataOut.getCoherence(phase=True)
878 877 meta['pairs'] = dataOut.pairsList
879 878
880 879 return data, meta
881 880
882 881 class NoisePlot(Plot):
883 882 '''
884 883 Plot for noise
885 884 '''
886 885
887 886 CODE = 'noise'
888 887 plot_type = 'scatterbuffer'
889 888
890 889 def setup(self):
891 890 self.xaxis = 'time'
892 891 self.ncols = 1
893 892 self.nrows = 1
894 893 self.nplots = 1
895 894 self.ylabel = 'Intensity [dB]'
896 895 self.xlabel = 'Time'
897 896 self.titles = ['Noise']
898 897 self.colorbar = False
899 898 self.plots_adjust.update({'right': 0.85 })
900 899 self.titles = ['Noise Plot']
901 900
902 901 def update(self, dataOut):
903 902
904 903 data = {}
905 904 meta = {}
906 905 noise = 10*numpy.log10(dataOut.getNoise())
907 906 noise = noise.reshape(dataOut.nChannels, 1)
908 907 data['noise'] = noise
909 908 meta['yrange'] = numpy.array([])
910 909
911 910 return data, meta
912 911
913 912 def plot(self):
914 913
915 914 x = self.data.times
916 915 xmin = self.data.min_time
917 916 xmax = xmin + self.xrange * 60 * 60
918 917 Y = self.data['noise']
919 918
920 919 if self.axes[0].firsttime:
921 920 self.ymin = numpy.nanmin(Y) - 5
922 921 self.ymax = numpy.nanmax(Y) + 5
923 922 for ch in self.data.channels:
924 923 y = Y[ch]
925 924 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
926 925 plt.legend(bbox_to_anchor=(1.18, 1.0))
927 926 else:
928 927 for ch in self.data.channels:
929 928 y = Y[ch]
930 929 self.axes[0].lines[ch].set_data(x, y)
931 930
932 931 class PowerProfilePlot(Plot):
933 932
934 933 CODE = 'pow_profile'
935 934 plot_type = 'scatter'
936 935
937 936 def setup(self):
938 937
939 938 self.ncols = 1
940 939 self.nrows = 1
941 940 self.nplots = 1
942 941 self.height = 4
943 942 self.width = 3
944 943 self.ylabel = 'Range [km]'
945 944 self.xlabel = 'Intensity [dB]'
946 945 self.titles = ['Power Profile']
947 946 self.colorbar = False
948 947
949 948 def update(self, dataOut):
950 949
951 950 data = {}
952 951 meta = {}
953 952 data[self.CODE] = dataOut.getPower()
954 953
955 954 return data, meta
956 955
957 956 def plot(self):
958 957
959 958 y = self.data.yrange
960 959 self.y = y
961 960
962 961 x = self.data[-1][self.CODE]
963 962
964 963 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
965 964 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
966 965
967 966 if self.axes[0].firsttime:
968 967 for ch in self.data.channels:
969 968 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
970 969 plt.legend()
971 970 else:
972 971 for ch in self.data.channels:
973 972 self.axes[0].lines[ch].set_data(x[ch], y)
974 973
975 974
976 975 class SpectraCutPlot(Plot):
977 976
978 977 CODE = 'spc_cut'
979 978 plot_type = 'scatter'
980 979 buffering = False
981 980 heights = []
982 981 channelList = []
983 982 maintitle = "Spectra Cuts"
984 983 flag_setIndex = False
985 984
986 985 def setup(self):
987 986
988 987 self.nplots = len(self.data.channels)
989 988 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
990 989 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
991 990 self.width = 4.5 * self.ncols + 2.5
992 991 self.height = 4.8 * self.nrows
993 992 self.ylabel = 'Power [dB]'
994 993 self.colorbar = False
995 994 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
996 995
997 996 if len(self.selectedHeightsList) > 0:
998 997 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
999 998
1000 999
1001 1000
1002 1001 def update(self, dataOut):
1003 1002 if len(self.channelList) == 0:
1004 1003 self.channelList = dataOut.channelList
1005 1004
1006 1005 self.heights = dataOut.heightList
1007 1006 #print("sels: ",self.selectedHeightsList)
1008 1007 if len(self.selectedHeightsList)>0 and not self.flag_setIndex:
1009 1008
1010 1009 for sel_height in self.selectedHeightsList:
1011 1010 index_list = numpy.where(self.heights >= sel_height)
1012 1011 index_list = index_list[0]
1013 1012 self.height_index.append(index_list[0])
1014 1013 #print("sels i:"", self.height_index)
1015 1014 self.flag_setIndex = True
1016 1015 #print(self.height_index)
1017 1016 data = {}
1018 1017 meta = {}
1019 1018
1020 1019 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
1021 1020 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1022 1021 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1023 1022
1024 1023
1025 1024 z = []
1026 1025 for ch in range(dataOut.nChannels):
1027 1026 if hasattr(dataOut.normFactor,'shape'):
1028 1027 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1029 1028 else:
1030 1029 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1031 1030
1032 1031 z = numpy.asarray(z)
1033 1032 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1034 1033 spc = 10*numpy.log10(z)
1035 1034
1036 1035
1037 1036 data['spc'] = spc - noise
1038 1037 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1039 1038
1040 1039 return data, meta
1041 1040
1042 1041 def plot(self):
1043 1042 if self.xaxis == "frequency":
1044 1043 x = self.data.xrange[0][0:]
1045 1044 self.xlabel = "Frequency (kHz)"
1046 1045 elif self.xaxis == "time":
1047 1046 x = self.data.xrange[1]
1048 1047 self.xlabel = "Time (ms)"
1049 1048 else:
1050 1049 x = self.data.xrange[2]
1051 1050 self.xlabel = "Velocity (m/s)"
1052 1051
1053 1052 self.titles = []
1054 1053
1055 1054 y = self.data.yrange
1056 1055 z = self.data[-1]['spc']
1057 1056 #print(z.shape)
1058 1057 if len(self.height_index) > 0:
1059 1058 index = self.height_index
1060 1059 else:
1061 1060 index = numpy.arange(0, len(y), int((len(y))/9))
1062 1061 #print("inde x ", index, self.axes)
1063 1062
1064 1063 for n, ax in enumerate(self.axes):
1065 1064
1066 1065 if ax.firsttime:
1067 1066
1068 1067
1069 1068 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1070 1069 self.xmin = self.xmin if self.xmin else -self.xmax
1071 1070 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
1072 1071 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
1073 1072
1074 1073
1075 1074 ax.plt = ax.plot(x, z[n, :, index].T)
1076 1075 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1077 1076 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
1078 1077 ax.minorticks_on()
1079 1078 ax.grid(which='major', axis='both')
1080 1079 ax.grid(which='minor', axis='x')
1081 1080 else:
1082 1081 for i, line in enumerate(ax.plt):
1083 1082 line.set_data(x, z[n, :, index[i]])
1084 1083
1085 1084
1086 1085 self.titles.append('CH {}'.format(self.channelList[n]))
1087 1086 plt.suptitle(self.maintitle, fontsize=10)
1088 1087
1089 1088
1090 1089 class BeaconPhase(Plot):
1091 1090
1092 1091 __isConfig = None
1093 1092 __nsubplots = None
1094 1093
1095 1094 PREFIX = 'beacon_phase'
1096 1095
1097 1096 def __init__(self):
1098 1097 Plot.__init__(self)
1099 1098 self.timerange = 24*60*60
1100 1099 self.isConfig = False
1101 1100 self.__nsubplots = 1
1102 1101 self.counter_imagwr = 0
1103 1102 self.WIDTH = 800
1104 1103 self.HEIGHT = 400
1105 1104 self.WIDTHPROF = 120
1106 1105 self.HEIGHTPROF = 0
1107 1106 self.xdata = None
1108 1107 self.ydata = None
1109 1108
1110 1109 self.PLOT_CODE = BEACON_CODE
1111 1110
1112 1111 self.FTP_WEI = None
1113 1112 self.EXP_CODE = None
1114 1113 self.SUB_EXP_CODE = None
1115 1114 self.PLOT_POS = None
1116 1115
1117 1116 self.filename_phase = None
1118 1117
1119 1118 self.figfile = None
1120 1119
1121 1120 self.xmin = None
1122 1121 self.xmax = None
1123 1122
1124 1123 def getSubplots(self):
1125 1124
1126 1125 ncol = 1
1127 1126 nrow = 1
1128 1127
1129 1128 return nrow, ncol
1130 1129
1131 1130 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1132 1131
1133 1132 self.__showprofile = showprofile
1134 1133 self.nplots = nplots
1135 1134
1136 1135 ncolspan = 7
1137 1136 colspan = 6
1138 1137 self.__nsubplots = 2
1139 1138
1140 1139 self.createFigure(id = id,
1141 1140 wintitle = wintitle,
1142 1141 widthplot = self.WIDTH+self.WIDTHPROF,
1143 1142 heightplot = self.HEIGHT+self.HEIGHTPROF,
1144 1143 show=show)
1145 1144
1146 1145 nrow, ncol = self.getSubplots()
1147 1146
1148 1147 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1149 1148
1150 1149 def save_phase(self, filename_phase):
1151 1150 f = open(filename_phase,'w+')
1152 1151 f.write('\n\n')
1153 1152 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1154 1153 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1155 1154 f.close()
1156 1155
1157 1156 def save_data(self, filename_phase, data, data_datetime):
1158 1157 f=open(filename_phase,'a')
1159 1158 timetuple_data = data_datetime.timetuple()
1160 1159 day = str(timetuple_data.tm_mday)
1161 1160 month = str(timetuple_data.tm_mon)
1162 1161 year = str(timetuple_data.tm_year)
1163 1162 hour = str(timetuple_data.tm_hour)
1164 1163 minute = str(timetuple_data.tm_min)
1165 1164 second = str(timetuple_data.tm_sec)
1166 1165 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1167 1166 f.close()
1168 1167
1169 1168 def plot(self):
1170 1169 log.warning('TODO: Not yet implemented...')
1171 1170
1172 1171 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1173 1172 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1174 1173 timerange=None,
1175 1174 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1176 1175 server=None, folder=None, username=None, password=None,
1177 1176 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1178 1177
1179 1178 if dataOut.flagNoData:
1180 1179 return dataOut
1181 1180
1182 1181 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1183 1182 return
1184 1183
1185 1184 if pairsList == None:
1186 1185 pairsIndexList = dataOut.pairsIndexList[:10]
1187 1186 else:
1188 1187 pairsIndexList = []
1189 1188 for pair in pairsList:
1190 1189 if pair not in dataOut.pairsList:
1191 1190 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1192 1191 pairsIndexList.append(dataOut.pairsList.index(pair))
1193 1192
1194 1193 if pairsIndexList == []:
1195 1194 return
1196 1195
1197 1196 # if len(pairsIndexList) > 4:
1198 1197 # pairsIndexList = pairsIndexList[0:4]
1199 1198
1200 1199 hmin_index = None
1201 1200 hmax_index = None
1202 1201
1203 1202 if hmin != None and hmax != None:
1204 1203 indexes = numpy.arange(dataOut.nHeights)
1205 1204 hmin_list = indexes[dataOut.heightList >= hmin]
1206 1205 hmax_list = indexes[dataOut.heightList <= hmax]
1207 1206
1208 1207 if hmin_list.any():
1209 1208 hmin_index = hmin_list[0]
1210 1209
1211 1210 if hmax_list.any():
1212 1211 hmax_index = hmax_list[-1]+1
1213 1212
1214 1213 x = dataOut.getTimeRange()
1215 1214
1216 1215 thisDatetime = dataOut.datatime
1217 1216
1218 1217 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1219 1218 xlabel = "Local Time"
1220 1219 ylabel = "Phase (degrees)"
1221 1220
1222 1221 update_figfile = False
1223 1222
1224 1223 nplots = len(pairsIndexList)
1225 1224 phase_beacon = numpy.zeros(len(pairsIndexList))
1226 1225 for i in range(nplots):
1227 1226 pair = dataOut.pairsList[pairsIndexList[i]]
1228 1227 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1229 1228 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1230 1229 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1231 1230 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1232 1231 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1233 1232
1234 1233 if dataOut.beacon_heiIndexList:
1235 1234 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1236 1235 else:
1237 1236 phase_beacon[i] = numpy.average(phase)
1238 1237
1239 1238 if not self.isConfig:
1240 1239
1241 1240 nplots = len(pairsIndexList)
1242 1241
1243 1242 self.setup(id=id,
1244 1243 nplots=nplots,
1245 1244 wintitle=wintitle,
1246 1245 showprofile=showprofile,
1247 1246 show=show)
1248 1247
1249 1248 if timerange != None:
1250 1249 self.timerange = timerange
1251 1250
1252 1251 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1253 1252
1254 1253 if ymin == None: ymin = 0
1255 1254 if ymax == None: ymax = 360
1256 1255
1257 1256 self.FTP_WEI = ftp_wei
1258 1257 self.EXP_CODE = exp_code
1259 1258 self.SUB_EXP_CODE = sub_exp_code
1260 1259 self.PLOT_POS = plot_pos
1261 1260
1262 1261 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1263 1262 self.isConfig = True
1264 1263 self.figfile = figfile
1265 1264 self.xdata = numpy.array([])
1266 1265 self.ydata = numpy.array([])
1267 1266
1268 1267 update_figfile = True
1269 1268
1270 1269 #open file beacon phase
1271 1270 path = '%s%03d' %(self.PREFIX, self.id)
1272 1271 beacon_file = os.path.join(path,'%s.txt'%self.name)
1273 1272 self.filename_phase = os.path.join(figpath,beacon_file)
1274 1273
1275 1274 self.setWinTitle(title)
1276 1275
1277 1276
1278 1277 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1279 1278
1280 1279 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1281 1280
1282 1281 axes = self.axesList[0]
1283 1282
1284 1283 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1285 1284
1286 1285 if len(self.ydata)==0:
1287 1286 self.ydata = phase_beacon.reshape(-1,1)
1288 1287 else:
1289 1288 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1290 1289
1291 1290
1292 1291 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1293 1292 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1294 1293 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1295 1294 XAxisAsTime=True, grid='both'
1296 1295 )
1297 1296
1298 1297 self.draw()
1299 1298
1300 1299 if dataOut.ltctime >= self.xmax:
1301 1300 self.counter_imagwr = wr_period
1302 1301 self.isConfig = False
1303 1302 update_figfile = True
1304 1303
1305 1304 self.save(figpath=figpath,
1306 1305 figfile=figfile,
1307 1306 save=save,
1308 1307 ftp=ftp,
1309 1308 wr_period=wr_period,
1310 1309 thisDatetime=thisDatetime,
1311 1310 update_figfile=update_figfile)
1312 1311
1313 1312 return dataOut
1314 1313
1315 1314 #####################################
1316 1315 class NoiselessSpectraPlot(Plot):
1317 1316 '''
1318 1317 Plot for Spectra data, subtracting
1319 1318 the noise in all channels, using for
1320 1319 amisr-14 data
1321 1320 '''
1322 1321
1323 1322 CODE = 'noiseless_spc'
1324 1323 colormap = 'jet'
1325 1324 plot_type = 'pcolor'
1326 1325 buffering = False
1327 1326 channelList = []
1328 1327 last_noise = None
1329 1328
1330 1329 def setup(self):
1331 1330
1332 1331 self.nplots = len(self.data.channels)
1333 1332 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
1334 1333 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
1335 1334 self.height = 3.5 * self.nrows
1336 1335
1337 1336 self.cb_label = 'dB'
1338 1337 if self.showprofile:
1339 1338 self.width = 5.8 * self.ncols
1340 1339 else:
1341 1340 self.width = 4.8* self.ncols
1342 1341 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12})
1343 1342
1344 1343 self.ylabel = 'Range [km]'
1345 1344
1346 1345
1347 1346 def update_list(self,dataOut):
1348 1347 if len(self.channelList) == 0:
1349 1348 self.channelList = dataOut.channelList
1350 1349
1351 1350 def update(self, dataOut):
1352 1351
1353 1352 self.update_list(dataOut)
1354 1353 data = {}
1355 1354 meta = {}
1356 1355
1357 1356 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1358 1357 n0 = (dataOut.getNoise()/norm)
1359 1358 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1360 1359 noise = 10*numpy.log10(noise)
1361 1360
1362 1361 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
1363 1362 for ch in range(dataOut.nChannels):
1364 1363 if hasattr(dataOut.normFactor,'ndim'):
1365 1364 if dataOut.normFactor.ndim > 1:
1366 1365 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1367 1366 else:
1368 1367 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1369 1368 else:
1370 1369 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1371 1370
1372 1371 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1373 1372 spc = 10*numpy.log10(z)
1374 1373
1375 1374
1376 1375 data['spc'] = spc - noise
1377 1376 #print(spc.shape)
1378 1377 data['rti'] = spc.mean(axis=1)
1379 1378 data['noise'] = noise
1380 1379
1381 1380
1382 1381
1383 1382 # data['noise'] = noise
1384 1383 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1385 1384
1386 1385 return data, meta
1387 1386
1388 1387 def plot(self):
1389 1388 if self.xaxis == "frequency":
1390 1389 x = self.data.xrange[0]
1391 1390 self.xlabel = "Frequency (kHz)"
1392 1391 elif self.xaxis == "time":
1393 1392 x = self.data.xrange[1]
1394 1393 self.xlabel = "Time (ms)"
1395 1394 else:
1396 1395 x = self.data.xrange[2]
1397 1396 self.xlabel = "Velocity (m/s)"
1398 1397
1399 1398 self.titles = []
1400 1399 y = self.data.yrange
1401 1400 self.y = y
1402 1401
1403 1402 data = self.data[-1]
1404 1403 z = data['spc']
1405 1404
1406 1405 for n, ax in enumerate(self.axes):
1407 1406 #noise = data['noise'][n]
1408 1407
1409 1408 if ax.firsttime:
1410 1409 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1411 1410 self.xmin = self.xmin if self.xmin else -self.xmax
1412 1411 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1413 1412 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1414 1413 ax.plt = ax.pcolormesh(x, y, z[n].T,
1415 1414 vmin=self.zmin,
1416 1415 vmax=self.zmax,
1417 1416 cmap=plt.get_cmap(self.colormap)
1418 1417 )
1419 1418
1420 1419 if self.showprofile:
1421 1420 ax.plt_profile = self.pf_axes[n].plot(
1422 1421 data['rti'][n], y)[0]
1423 1422
1424 1423
1425 1424 else:
1426 1425 ax.plt.set_array(z[n].T.ravel())
1427 1426 if self.showprofile:
1428 1427 ax.plt_profile.set_data(data['rti'][n], y)
1429 1428
1430 1429
1431 1430 self.titles.append('CH {}'.format(self.channelList[n]))
1432 1431
1433 1432
1434 1433 class NoiselessRTIPlot(RTIPlot):
1435 1434 '''
1436 1435 Plot for RTI data
1437 1436 '''
1438 1437
1439 1438 CODE = 'noiseless_rti'
1440 1439 colormap = 'jet'
1441 1440 plot_type = 'pcolorbuffer'
1442 1441 titles = None
1443 1442 channelList = []
1444 1443 elevationList = []
1445 1444 azimuthList = []
1446 1445 last_noise = None
1447 1446
1448 1447 def setup(self):
1449 1448 self.xaxis = 'time'
1450 1449 self.ncols = 1
1451 1450 #print("dataChannels ",self.data.channels)
1452 1451 self.nrows = len(self.data.channels)
1453 1452 self.nplots = len(self.data.channels)
1454 1453 self.ylabel = 'Range [km]'
1455 1454 #self.xlabel = 'Time'
1456 1455 self.cb_label = 'dB'
1457 1456 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1458 1457 self.titles = ['{} Channel {}'.format(
1459 1458 self.CODE.upper(), x) for x in range(self.nplots)]
1460 1459
1461 1460 def update_list(self,dataOut):
1462 1461 if len(self.channelList) == 0:
1463 1462 self.channelList = dataOut.channelList
1464 1463 if len(self.elevationList) == 0:
1465 1464 self.elevationList = dataOut.elevationList
1466 1465 if len(self.azimuthList) == 0:
1467 1466 self.azimuthList = dataOut.azimuthList
1468 1467
1469 1468 def update(self, dataOut):
1470 1469 if len(self.channelList) == 0:
1471 1470 self.update_list(dataOut)
1472 1471
1473 1472 data = {}
1474 1473 meta = {}
1475 1474 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1476 1475 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt
1477 1476 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1478 1477 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1479 1478 data['noise'] = n0
1480 1479 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1481 1480 noiseless_data = dataOut.getPower() - noise
1482 1481
1483 1482 #print("power, noise:", dataOut.getPower(), n0)
1484 1483 #print(noise)
1485 1484 #print(noiseless_data)
1486 1485
1487 1486 data['noiseless_rti'] = noiseless_data
1488 1487
1489 1488 return data, meta
1490 1489
1491 1490 def plot(self):
1492 1491 from matplotlib import pyplot as plt
1493 1492 self.x = self.data.times
1494 1493 self.y = self.data.yrange
1495 1494 self.z = self.data['noiseless_rti']
1496 1495 self.z = numpy.array(self.z, dtype=float)
1497 1496 self.z = numpy.ma.masked_invalid(self.z)
1498 1497
1499 1498
1500 1499 try:
1501 1500 if self.channelList != None:
1502 1501 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1503 1502 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1504 1503 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1505 1504 else:
1506 1505 self.titles = ['{} Channel {}'.format(
1507 1506 self.CODE.upper(), x) for x in self.channelList]
1508 1507 except:
1509 1508 if self.channelList.any() != None:
1510 1509 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1511 1510 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1512 1511 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1513 1512 else:
1514 1513 self.titles = ['{} Channel {}'.format(
1515 1514 self.CODE.upper(), x) for x in self.channelList]
1516 1515
1517 1516
1518 1517 if self.decimation is None:
1519 1518 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1520 1519 else:
1521 1520 x, y, z = self.fill_gaps(*self.decimate())
1522 1521
1523 1522 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
1524 1523 #print("plot shapes ", z.shape, x.shape, y.shape)
1525 1524 #print(self.axes)
1526 1525 for n, ax in enumerate(self.axes):
1527 1526
1528 1527
1529 1528 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
1530 1529 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1531 1530 data = self.data[-1]
1532 1531 if ax.firsttime:
1533 1532 if (n+1) == len(self.channelList):
1534 1533 ax.set_xlabel('Time')
1535 1534 ax.plt = ax.pcolormesh(x, y, z[n].T,
1536 1535 vmin=self.zmin,
1537 1536 vmax=self.zmax,
1538 1537 cmap=plt.get_cmap(self.colormap)
1539 1538 )
1540 1539 if self.showprofile:
1541 1540 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1542 1541
1543 1542 else:
1544 # ax.collections.remove(ax.collections[0]) # error while running
1543 ax.collections.remove(ax.collections[0]) # error while running
1545 1544 ax.plt = ax.pcolormesh(x, y, z[n].T,
1546 1545 vmin=self.zmin,
1547 1546 vmax=self.zmax,
1548 1547 cmap=plt.get_cmap(self.colormap)
1549 1548 )
1550 1549 if self.showprofile:
1551 1550 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1552 1551 # if "noise" in self.data:
1553 1552 # #ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
1554 1553 # ax.plot_noise.set_data(data['noise'][n], self.y)
1555 1554
1556 1555
1557 1556 class OutliersRTIPlot(Plot):
1558 1557 '''
1559 1558 Plot for data_xxxx object
1560 1559 '''
1561 1560
1562 1561 CODE = 'outlier_rtc' # Range Time Counts
1563 1562 colormap = 'cool'
1564 1563 plot_type = 'pcolorbuffer'
1565 1564
1566 1565 def setup(self):
1567 1566 self.xaxis = 'time'
1568 1567 self.ncols = 1
1569 1568 self.nrows = self.data.shape('outlier_rtc')[0]
1570 1569 self.nplots = self.nrows
1571 1570 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1572 1571
1573 1572
1574 1573 if not self.xlabel:
1575 1574 self.xlabel = 'Time'
1576 1575
1577 1576 self.ylabel = 'Height [km]'
1578 1577 if not self.titles:
1579 1578 self.titles = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)]
1580 1579
1581 1580 def update(self, dataOut):
1582 1581
1583 1582 data = {}
1584 1583 data['outlier_rtc'] = dataOut.data_outlier
1585 1584
1586 1585 meta = {}
1587 1586
1588 1587 return data, meta
1589 1588
1590 1589 def plot(self):
1591 1590 # self.data.normalize_heights()
1592 1591 self.x = self.data.times
1593 1592 self.y = self.data.yrange
1594 1593 self.z = self.data['outlier_rtc']
1595 1594
1596 1595 #self.z = numpy.ma.masked_invalid(self.z)
1597 1596
1598 1597 if self.decimation is None:
1599 1598 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1600 1599 else:
1601 1600 x, y, z = self.fill_gaps(*self.decimate())
1602 1601
1603 1602 for n, ax in enumerate(self.axes):
1604 1603
1605 1604 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1606 1605 self.z[n])
1607 1606 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1608 1607 self.z[n])
1609 1608 data = self.data[-1]
1610 1609 if ax.firsttime:
1611 1610 if self.zlimits is not None:
1612 1611 self.zmin, self.zmax = self.zlimits[n]
1613 1612
1614 1613 ax.plt = ax.pcolormesh(x, y, z[n].T,
1615 1614 vmin=self.zmin,
1616 1615 vmax=self.zmax,
1617 1616 cmap=self.cmaps[n]
1618 1617 )
1619 1618 if self.showprofile:
1620 1619 ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0]
1621 1620 self.pf_axes[n].set_xlabel('')
1622 1621 else:
1623 1622 if self.zlimits is not None:
1624 1623 self.zmin, self.zmax = self.zlimits[n]
1625 # ax.collections.remove(ax.collections[0]) # error while running
1624 ax.collections.remove(ax.collections[0]) # error while running
1626 1625 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1627 1626 vmin=self.zmin,
1628 1627 vmax=self.zmax,
1629 1628 cmap=self.cmaps[n]
1630 1629 )
1631 1630 if self.showprofile:
1632 1631 ax.plot_profile.set_data(data['outlier_rtc'][n], self.y)
1633 1632 self.pf_axes[n].set_xlabel('')
1634 1633
1635 1634 class NIncohIntRTIPlot(Plot):
1636 1635 '''
1637 1636 Plot for data_xxxx object
1638 1637 '''
1639 1638
1640 1639 CODE = 'integrations_rtc' # Range Time Counts
1641 1640 colormap = 'BuGn'
1642 1641 plot_type = 'pcolorbuffer'
1643 1642
1644 1643 def setup(self):
1645 1644 self.xaxis = 'time'
1646 1645 self.ncols = 1
1647 1646 self.nrows = self.data.shape('integrations_rtc')[0]
1648 1647 self.nplots = self.nrows
1649 1648 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1650 1649
1651 1650
1652 1651 if not self.xlabel:
1653 1652 self.xlabel = 'Time'
1654 1653
1655 1654 self.ylabel = 'Height [km]'
1656 1655 if not self.titles:
1657 1656 self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)]
1658 1657
1659 1658 def update(self, dataOut):
1660 1659
1661 1660 data = {}
1662 1661 data['integrations_rtc'] = dataOut.nIncohInt
1663 1662
1664 1663 meta = {}
1665 1664
1666 1665 return data, meta
1667 1666
1668 1667 def plot(self):
1669 1668 # self.data.normalize_heights()
1670 1669 self.x = self.data.times
1671 1670 self.y = self.data.yrange
1672 1671 self.z = self.data['integrations_rtc']
1673 1672
1674 1673 #self.z = numpy.ma.masked_invalid(self.z)
1675 1674
1676 1675 if self.decimation is None:
1677 1676 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1678 1677 else:
1679 1678 x, y, z = self.fill_gaps(*self.decimate())
1680 1679
1681 1680 for n, ax in enumerate(self.axes):
1682 1681
1683 1682 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1684 1683 self.z[n])
1685 1684 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1686 1685 self.z[n])
1687 1686 data = self.data[-1]
1688 1687 if ax.firsttime:
1689 1688 if self.zlimits is not None:
1690 1689 self.zmin, self.zmax = self.zlimits[n]
1691 1690
1692 1691 ax.plt = ax.pcolormesh(x, y, z[n].T,
1693 1692 vmin=self.zmin,
1694 1693 vmax=self.zmax,
1695 1694 cmap=self.cmaps[n]
1696 1695 )
1697 1696 if self.showprofile:
1698 1697 ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0]
1699 1698 self.pf_axes[n].set_xlabel('')
1700 1699 else:
1701 1700 if self.zlimits is not None:
1702 1701 self.zmin, self.zmax = self.zlimits[n]
1703 # ax.collections.remove(ax.collections[0]) # error while running
1702 ax.collections.remove(ax.collections[0]) # error while running
1704 1703 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1705 1704 vmin=self.zmin,
1706 1705 vmax=self.zmax,
1707 1706 cmap=self.cmaps[n]
1708 1707 )
1709 1708 if self.showprofile:
1710 1709 ax.plot_profile.set_data(data['integrations_rtc'][n], self.y)
1711 1710 self.pf_axes[n].set_xlabel('')
1712 1711
1713 1712
1714 1713
1715 1714 class RTIMapPlot(Plot):
1716 1715 '''
1717 1716 Plot for RTI data
1718 1717
1719 1718 Example:
1720 1719
1721 1720 controllerObj = Project()
1722 1721 controllerObj.setup(id = '11', name='eej_proc', description=desc)
1723 1722 ##.......................................................................................
1724 1723 ##.......................................................................................
1725 1724 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', path=inPath, startDate='2023/05/24',endDate='2023/05/24',
1726 1725 startTime='12:00:00',endTime='12:45:59',walk=1,timezone='lt',margin_days=1,code = code,nCode = nCode,
1727 1726 nBaud = nBaud,nOsamp = nosamp,nChannels=nChannels,nFFT=NFFT,
1728 1727 syncronization=False,shiftChannels=0)
1729 1728
1730 1729 volts_proc = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
1731 1730
1732 1731 opObj01 = volts_proc.addOperation(name='Decoder', optype='other')
1733 1732 opObj01.addParameter(name='code', value=code, format='floatlist')
1734 1733 opObj01.addParameter(name='nCode', value=1, format='int')
1735 1734 opObj01.addParameter(name='nBaud', value=nBaud, format='int')
1736 1735 opObj01.addParameter(name='osamp', value=nosamp, format='int')
1737 1736
1738 1737 opObj12 = volts_proc.addOperation(name='selectHeights', optype='self')
1739 1738 opObj12.addParameter(name='minHei', value='90', format='float')
1740 1739 opObj12.addParameter(name='maxHei', value='150', format='float')
1741 1740
1742 1741 proc_spc = controllerObj.addProcUnit(datatype='SpectraProc', inputId=volts_proc.getId())
1743 1742 proc_spc.addParameter(name='nFFTPoints', value='8', format='int')
1744 1743
1745 1744 opObj11 = proc_spc.addOperation(name='IncohInt', optype='other')
1746 1745 opObj11.addParameter(name='n', value='1', format='int')
1747 1746
1748 1747 beamMapFile = "/home/japaza/Documents/AMISR_sky_mapper/UMET_beamcodes.csv"
1749 1748
1750 1749 opObj12 = proc_spc.addOperation(name='RTIMapPlot', optype='external')
1751 1750 opObj12.addParameter(name='selectedHeightsList', value='95, 100, 105, 110 ', format='int')
1752 1751 opObj12.addParameter(name='bField', value='100', format='int')
1753 1752 opObj12.addParameter(name='filename', value=beamMapFile, format='str')
1754 1753
1755 1754 '''
1756 1755
1757 1756 CODE = 'rti_skymap'
1758 1757
1759 1758 plot_type = 'scatter'
1760 1759 titles = None
1761 1760 colormap = 'jet'
1762 1761 channelList = []
1763 1762 elevationList = []
1764 1763 azimuthList = []
1765 1764 last_noise = None
1766 1765 flag_setIndex = False
1767 1766 heights = []
1768 1767 dcosx = []
1769 1768 dcosy = []
1770 1769 fullDcosy = None
1771 1770 fullDcosy = None
1772 1771 hindex = []
1773 1772 mapFile = False
1774 1773 ##### BField ####
1775 1774 flagBField = False
1776 1775 dcosxB = []
1777 1776 dcosyB = []
1778 1777 Bmarker = ['+','*','D','x','s','>','o','^']
1779 1778
1780 1779
1781 1780 def setup(self):
1782 1781
1783 1782 self.xaxis = 'Range (Km)'
1784 1783 if len(self.selectedHeightsList) > 0:
1785 1784 self.nplots = len(self.selectedHeightsList)
1786 1785 else:
1787 1786 self.nplots = 4
1788 1787 self.ncols = int(numpy.ceil(self.nplots/2))
1789 1788 self.nrows = int(numpy.ceil(self.nplots/self.ncols))
1790 1789 self.ylabel = 'dcosy'
1791 1790 self.xlabel = 'dcosx'
1792 1791 self.colorbar = True
1793 1792 self.width = 6 + 4.1*self.nrows
1794 1793 self.height = 3 + 3.5*self.ncols
1795 1794
1796 1795
1797 1796 if self.extFile!=None:
1798 1797 try:
1799 1798 pointings = numpy.genfromtxt(self.extFile, delimiter=',')
1800 1799 full_azi = pointings[:,1]
1801 1800 full_elev = pointings[:,2]
1802 1801 self.fullDcosx = numpy.cos(numpy.radians(full_elev))*numpy.sin(numpy.radians(full_azi))
1803 1802 self.fullDcosy = numpy.cos(numpy.radians(full_elev))*numpy.cos(numpy.radians(full_azi))
1804 1803 mapFile = True
1805 1804 except Exception as e:
1806 1805 self.extFile = None
1807 1806 print(e)
1808 1807
1809 1808
1810 1809 def update_list(self,dataOut):
1811 1810 if len(self.channelList) == 0:
1812 1811 self.channelList = dataOut.channelList
1813 1812 if len(self.elevationList) == 0:
1814 1813 self.elevationList = dataOut.elevationList
1815 1814 if len(self.azimuthList) == 0:
1816 1815 self.azimuthList = dataOut.azimuthList
1817 1816 a = numpy.radians(numpy.asarray(self.azimuthList))
1818 1817 e = numpy.radians(numpy.asarray(self.elevationList))
1819 1818 self.heights = dataOut.heightList
1820 1819 self.dcosx = numpy.cos(e)*numpy.sin(a)
1821 1820 self.dcosy = numpy.cos(e)*numpy.cos(a)
1822 1821
1823 1822 if len(self.bFieldList)>0:
1824 1823 datetObj = datetime.datetime.fromtimestamp(dataOut.utctime)
1825 1824 doy = datetObj.timetuple().tm_yday
1826 1825 year = datetObj.year
1827 1826 # self.dcosxB, self.dcosyB
1828 1827 ObjB = BField(year=year,doy=doy,site=2,heights=self.bFieldList)
1829 1828 [dcos, alpha, nlon, nlat] = ObjB.getBField()
1830 1829
1831 1830 alpha_location = numpy.zeros((nlon,2,len(self.bFieldList)))
1832 1831 for ih in range(len(self.bFieldList)):
1833 1832 alpha_location[:,0,ih] = dcos[:,0,ih,0]
1834 1833 for ilon in numpy.arange(nlon):
1835 1834 myx = (alpha[ilon,:,ih])[::-1]
1836 1835 myy = (dcos[ilon,:,ih,0])[::-1]
1837 1836 tck = splrep(myx,myy,s=0)
1838 1837 mydcosx = splev(ObjB.alpha_i,tck,der=0)
1839 1838
1840 1839 myx = (alpha[ilon,:,ih])[::-1]
1841 1840 myy = (dcos[ilon,:,ih,1])[::-1]
1842 1841 tck = splrep(myx,myy,s=0)
1843 1842 mydcosy = splev(ObjB.alpha_i,tck,der=0)
1844 1843 alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy])
1845 1844 self.dcosxB.append(alpha_location[:,0,ih])
1846 1845 self.dcosyB.append(alpha_location[:,1,ih])
1847 1846 self.flagBField = True
1848 1847
1849 1848 if len(self.celestialList)>0:
1850 1849 #getBField(self.bFieldList, date)
1851 1850 #pass = kwargs.get('celestial', [])
1852 1851 pass
1853 1852
1854 1853
1855 1854 def update(self, dataOut):
1856 1855
1857 1856 if len(self.channelList) == 0:
1858 1857 self.update_list(dataOut)
1859 1858
1860 1859 if not self.flag_setIndex:
1861 1860 if len(self.selectedHeightsList)>0:
1862 1861 for sel_height in self.selectedHeightsList:
1863 1862 index_list = numpy.where(self.heights >= sel_height)
1864 1863 index_list = index_list[0]
1865 1864 self.hindex.append(index_list[0])
1866 1865 self.flag_setIndex = True
1867 1866
1868 1867 data = {}
1869 1868 meta = {}
1870 1869
1871 1870 data['rti_skymap'] = dataOut.getPower()
1872 1871 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1873 1872 noise = 10*numpy.log10(dataOut.getNoise()/norm)
1874 1873 data['noise'] = noise
1875 1874
1876 1875 return data, meta
1877 1876
1878 1877 def plot(self):
1879 1878
1880 1879 self.x = self.dcosx
1881 1880 self.y = self.dcosy
1882 1881 self.z = self.data[-1]['rti_skymap']
1883 1882 self.z = numpy.array(self.z, dtype=float)
1884 1883
1885 1884 if len(self.hindex) > 0:
1886 1885 index = self.hindex
1887 1886 else:
1888 1887 index = numpy.arange(0, len(self.heights), int((len(self.heights))/4.2))
1889 1888
1890 1889 self.titles = ['Height {:.2f} km '.format(self.heights[i])+" " for i in index]
1891 1890 for n, ax in enumerate(self.axes):
1892 1891
1893 1892 if ax.firsttime:
1894 1893
1895 1894 self.xmax = self.xmax if self.xmax else numpy.nanmax(self.x)
1896 1895 self.xmin = self.xmin if self.xmin else numpy.nanmin(self.x)
1897 1896 self.ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
1898 1897 self.ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
1899 1898 self.zmax = self.zmax if self.zmax else numpy.nanmax(self.z)
1900 1899 self.zmin = self.zmin if self.zmin else numpy.nanmin(self.z)
1901 1900
1902 1901 if self.extFile!=None:
1903 1902 ax.scatter(self.fullDcosx, self.fullDcosy, marker="+", s=20)
1904 1903
1905 1904 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1906 1905 s=60, marker="s", vmax = self.zmax)
1907 1906
1908 1907
1909 1908 ax.minorticks_on()
1910 1909 ax.grid(which='major', axis='both')
1911 1910 ax.grid(which='minor', axis='x')
1912 1911
1913 1912 if self.flagBField :
1914 1913
1915 1914 for ih in range(len(self.bFieldList)):
1916 1915 label = str(self.bFieldList[ih]) + ' km'
1917 1916 ax.plot(self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1918 1917 label=label, linestyle='--', ms=4.0,lw=0.5)
1919 1918 handles, labels = ax.get_legend_handles_labels()
1920 1919 a = -0.05
1921 1920 b = 1.15 - 1.19*(self.nrows)
1922 1921 self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥')
1923 1922
1924 1923 else:
1925 1924
1926 1925 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1927 1926 s=80, marker="s", vmax = self.zmax)
1928 1927
1929 1928 if self.flagBField :
1930 1929 for ih in range(len(self.bFieldList)):
1931 1930 ax.plot (self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1932 1931 linestyle='--', ms=4.0,lw=0.5)
1933 1932
1934 1933
1935 1934
Show file before | Show file after | Hide comments
@@ -1,819 +1,820
1 1 import os
2 2 import time
3 3 import datetime
4 4
5 5 import numpy
6 6 import h5py
7 7
8 8 import schainpy.admin
9 9 from schainpy.model.data.jrodata import *
10 10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 11 from schainpy.model.io.jroIO_base import *
12 12 from schainpy.utils import log
13 13
14 14
15 15 class HDFReader(Reader, ProcessingUnit):
16 16 """Processing unit to read HDF5 format files
17 17
18 18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 20 attributes.
21 21 It is possible to read any HDF5 file by given the structure in the `description`
22 22 parameter, also you can add extra values to metadata with the parameter `extras`.
23 23
24 24 Parameters:
25 25 -----------
26 26 path : str
27 27 Path where files are located.
28 28 startDate : date
29 29 Start date of the files
30 30 endDate : list
31 31 End date of the files
32 32 startTime : time
33 33 Start time of the files
34 34 endTime : time
35 35 End time of the files
36 36 description : dict, optional
37 37 Dictionary with the description of the HDF5 file
38 38 extras : dict, optional
39 39 Dictionary with extra metadata to be be added to `dataOut`
40 40
41 41 Attention: Be carefull, add attribute utcoffset, in the last part of reader in order to work in Local Time without time problems.
42 42
43 43 -----------
44 44 utcoffset='-18000'
45 45
46 46
47 47 Examples
48 48 --------
49 49
50 50 desc = {
51 51 'Data': {
52 52 'data_output': ['u', 'v', 'w'],
53 53 'utctime': 'timestamps',
54 54 } ,
55 55 'Metadata': {
56 56 'heightList': 'heights'
57 57 }
58 58 }
59 59
60 60 desc = {
61 61 'Data': {
62 62 'data_output': 'winds',
63 63 'utctime': 'timestamps'
64 64 },
65 65 'Metadata': {
66 66 'heightList': 'heights'
67 67 }
68 68 }
69 69
70 70 extras = {
71 71 'timeZone': 300
72 72 }
73 73
74 74 reader = project.addReadUnit(
75 75 name='HDFReader',
76 76 path='/path/to/files',
77 77 startDate='2019/01/01',
78 78 endDate='2019/01/31',
79 79 startTime='00:00:00',
80 80 endTime='23:59:59',
81 81 utcoffset='-18000'
82 82 # description=json.dumps(desc),
83 83 # extras=json.dumps(extras),
84 84 )
85 85
86 86 """
87 87
88 88 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
89 89
90 90 def __init__(self):
91 91
92 92 ProcessingUnit.__init__(self)
93 93 self.ext = ".hdf5"
94 94 self.optchar = "D"
95 95 self.meta = {}
96 96 self.data = {}
97 97 self.open_file = h5py.File
98 98 self.open_mode = 'r'
99 99 self.description = {}
100 100 self.extras = {}
101 101 self.filefmt = "*%Y%j***"
102 102 self.folderfmt = "*%Y%j"
103 103 self.utcoffset = 0
104 104 self.flagUpdateDataOut = False
105 105 self.dataOut = Parameters()
106 106 self.dataOut.error=False ## NOTE: Importante definir esto antes inicio
107 107 self.dataOut.flagNoData = True
108 108
109 109 def setup(self, **kwargs):
110 110
111 111 self.set_kwargs(**kwargs)
112 112 if not self.ext.startswith('.'):
113 113 self.ext = '.{}'.format(self.ext)
114 114
115 115 if self.online:
116 116 log.log("Searching files in online mode...", self.name)
117 117
118 118 for nTries in range(self.nTries):
119 119 fullpath = self.searchFilesOnLine(self.path, self.startDate,
120 120 self.endDate, self.expLabel, self.ext, self.walk,
121 121 self.filefmt, self.folderfmt)
122 122 pathname, filename = os.path.split(fullpath)
123 123 try:
124 124 fullpath = next(fullpath)
125 125 except:
126 126 fullpath = None
127 127
128 128 if fullpath:
129 129 break
130 130
131 131 log.warning(
132 132 'Waiting {} sec for a valid file in {}: try {} ...'.format(
133 133 self.delay, self.path, nTries + 1),
134 134 self.name)
135 135 time.sleep(self.delay)
136 136
137 137 if not(fullpath):
138 138 raise schainpy.admin.SchainError(
139 139 'There isn\'t any valid file in {}'.format(self.path))
140 140
141 141 pathname, filename = os.path.split(fullpath)
142 142 self.year = int(filename[1:5])
143 143 self.doy = int(filename[5:8])
144 144 self.set = int(filename[8:11]) - 1
145 145 else:
146 146 log.log("Searching files in {}".format(self.path), self.name)
147 147 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
148 148 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
149 149
150 150 self.setNextFile()
151 151
152 152 return
153 153
154 154 # def readFirstHeader(self):
155 155 # '''Read metadata and data'''
156 156
157 157 # self.__readMetadata()
158 158 # self.__readData()
159 159 # self.__setBlockList()
160 160
161 161 # if 'type' in self.meta:
162 162 # self.dataOut = eval(self.meta['type'])()
163 163
164 164 # for attr in self.meta:
165 165 # setattr(self.dataOut, attr, self.meta[attr])
166 166
167 167 # self.blockIndex = 0
168 168
169 169 # return
170 170
171 171 def readFirstHeader(self):
172 172 '''Read metadata and data'''
173 173
174 174 self.__readMetadata2()
175 175 self.__readData()
176 176 self.__setBlockList()
177 if 'type' in self.meta:
178 self.dataOut = eval(self.meta['type'])()
177 # if 'type' in self.meta:
178 # self.dataOut = eval(self.meta['type'])()
179 179
180 180 for attr in self.meta:
181 181 if "processingHeaderObj" in attr:
182 182 self.flagUpdateDataOut=True
183 183 at = attr.split('.')
184 184 if len(at) > 1:
185 185 setattr(eval("self.dataOut."+at[0]),at[1], self.meta[attr])
186 186 else:
187 187 setattr(self.dataOut, attr, self.meta[attr])
188 188 self.blockIndex = 0
189 189
190 190 if self.flagUpdateDataOut:
191 191 self.updateDataOut()
192 192
193 193 return
194 194
195 195 def updateDataOut(self):
196 196
197 197 self.dataOut.azimuthList = self.dataOut.processingHeaderObj.azimuthList
198 198 self.dataOut.elevationList = self.dataOut.processingHeaderObj.elevationList
199 199 self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
200 200 self.dataOut.ippSeconds = self.dataOut.processingHeaderObj.ipp
201 201 self.dataOut.elevationList = self.dataOut.processingHeaderObj.elevationList
202 202 self.dataOut.channelList = self.dataOut.processingHeaderObj.channelList
203 203 self.dataOut.nCohInt = self.dataOut.processingHeaderObj.nCohInt
204 204 self.dataOut.nFFTPoints = self.dataOut.processingHeaderObj.nFFTPoints
205 205 self.flagUpdateDataOut = False
206 206 self.dataOut.frequency = self.dataOut.radarControllerHeaderObj.frequency
207 207 #self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
208 208
209 209 def __setBlockList(self):
210 210 '''
211 211 Selects the data within the times defined
212 212
213 213 self.fp
214 214 self.startTime
215 215 self.endTime
216 216 self.blockList
217 217 self.blocksPerFile
218 218
219 219 '''
220 220
221 221 startTime = self.startTime
222 222 endTime = self.endTime
223 223 thisUtcTime = self.data['utctime'] + self.utcoffset
224 224 # self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
225 225 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
226 226 self.startFileDatetime = thisDatetime
227 227 thisDate = thisDatetime.date()
228 228 thisTime = thisDatetime.time()
229 229 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
230 230 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
231 231 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
232 232
233 233 self.blockList = ind
234 234 self.blocksPerFile = len(ind)
235 235 # self.blocksPerFile = len(thisUtcTime)
236 236 if len(ind)==0:
237 237 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.blockIndex,
238 238 self.blocksPerFile,
239 239 thisDatetime))
240 240 self.setNextFile()
241 241
242 242 return
243 243
244 244 def __readMetadata(self):
245 245 '''
246 246 Reads Metadata
247 247 '''
248 248
249 249 meta = {}
250 250
251 251 if self.description:
252 252 for key, value in self.description['Metadata'].items():
253 253 meta[key] = self.fp[value][()]
254 254 else:
255 255 grp = self.fp['Metadata']
256 256 for name in grp:
257 257 meta[name] = grp[name][()]
258 258
259 259 if self.extras:
260 260 for key, value in self.extras.items():
261 261 meta[key] = value
262 262 self.meta = meta
263 263
264 264 return
265 265
266 266 def __readMetadata2(self):
267 267 '''
268 268 Reads Metadata
269 269 '''
270 270 meta = {}
271
271 272 if self.description:
272 273 for key, value in self.description['Metadata'].items():
273 274 meta[key] = self.fp[value][()]
274 275 else:
275 276 grp = self.fp['Metadata']
276 277 for item in grp.values():
277 278 name = item.name
278 279 if isinstance(item, h5py.Dataset):
279 280 name = name.split("/")[-1]
280 281 meta[name] = item[()]
281 282 else:
282 283 grp2 = self.fp[name]
283 284 Obj = name.split("/")[-1]
284 285
285 286 for item2 in grp2.values():
286 287 name2 = Obj+"."+item2.name.split("/")[-1]
287 288 meta[name2] = item2[()]
288 289
289 290 if self.extras:
290 291 for key, value in self.extras.items():
291 292 meta[key] = value
292 293 self.meta = meta
293 294
294 295 return
295 296
296 297 def __readData(self):
297 298
298 299 data = {}
299 300
300 301 if self.description:
301 302 for key, value in self.description['Data'].items():
302 303 if isinstance(value, str):
303 304 if isinstance(self.fp[value], h5py.Dataset):
304 305 data[key] = self.fp[value][()]
305 306 elif isinstance(self.fp[value], h5py.Group):
306 307 array = []
307 308 for ch in self.fp[value]:
308 309 array.append(self.fp[value][ch][()])
309 310 data[key] = numpy.array(array)
310 311 elif isinstance(value, list):
311 312 array = []
312 313 for ch in value:
313 314 array.append(self.fp[ch][()])
314 315 data[key] = numpy.array(array)
315 316 else:
316 317 grp = self.fp['Data']
317 318 for name in grp:
318 319 if isinstance(grp[name], h5py.Dataset):
319 320 array = grp[name][()]
320 321 elif isinstance(grp[name], h5py.Group):
321 322 array = []
322 323 for ch in grp[name]:
323 324 array.append(grp[name][ch][()])
324 325 array = numpy.array(array)
325 326 else:
326 327 log.warning('Unknown type: {}'.format(name))
327 328
328 329 if name in self.description:
329 330 key = self.description[name]
330 331 else:
331 332 key = name
332 333 data[key] = array
333 334
334 335 self.data = data
335 336 return
336 337
337 338 def getData(self):
338 339
339 340 if not self.isDateTimeInRange(self.startFileDatetime, self.startDate, self.endDate, self.startTime, self.endTime):
340 341 self.dataOut.flagNoData = True
341 342 self.blockIndex = self.blocksPerFile
342 343 self.dataOut.error = True # TERMINA EL PROGRAMA
343 344 return
344 345 for attr in self.data:
345 346
346 347 if self.data[attr].ndim == 1:
347 348 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
348 349 else:
349 350 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
350 351
351 352
352 353 self.blockIndex += 1
353 354
354 355 if self.blockIndex == 1:
355 356 log.log("Block No. {}/{} -> {}".format(
356 357 self.blockIndex,
357 358 self.blocksPerFile,
358 359 self.dataOut.datatime.ctime()), self.name)
359 360 else:
360 361 log.log("Block No. {}/{} ".format(
361 362 self.blockIndex,
362 363 self.blocksPerFile),self.name)
363 364
364 365 if self.blockIndex == self.blocksPerFile:
365 366 self.setNextFile()
366 367
367 368 self.dataOut.flagNoData = False
368 369
369 370 return
370 371
371 372 def run(self, **kwargs):
372 373
373 374 if not(self.isConfig):
374 375 self.setup(**kwargs)
375 376 self.isConfig = True
376 377
377 378 if self.blockIndex == self.blocksPerFile:
378 379 self.setNextFile()
379 380
380 381 self.getData()
381 382
382 383 return
383 384
384 385 @MPDecorator
385 386 class HDFWriter(Operation):
386 387 """Operation to write HDF5 files.
387 388
388 389 The HDF5 file contains by default two groups Data and Metadata where
389 390 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
390 391 parameters, data attributes are normaly time dependent where the metadata
391 392 are not.
392 393 It is possible to customize the structure of the HDF5 file with the
393 394 optional description parameter see the examples.
394 395
395 396 Parameters:
396 397 -----------
397 398 path : str
398 399 Path where files will be saved.
399 400 blocksPerFile : int
400 401 Number of blocks per file
401 402 metadataList : list
402 403 List of the dataOut attributes that will be saved as metadata
403 404 dataList : int
404 405 List of the dataOut attributes that will be saved as data
405 406 setType : bool
406 407 If True the name of the files corresponds to the timestamp of the data
407 408 description : dict, optional
408 409 Dictionary with the desired description of the HDF5 file
409 410
410 411 Examples
411 412 --------
412 413
413 414 desc = {
414 415 'data_output': {'winds': ['z', 'w', 'v']},
415 416 'utctime': 'timestamps',
416 417 'heightList': 'heights'
417 418 }
418 419 desc = {
419 420 'data_output': ['z', 'w', 'v'],
420 421 'utctime': 'timestamps',
421 422 'heightList': 'heights'
422 423 }
423 424 desc = {
424 425 'Data': {
425 426 'data_output': 'winds',
426 427 'utctime': 'timestamps'
427 428 },
428 429 'Metadata': {
429 430 'heightList': 'heights'
430 431 }
431 432 }
432 433
433 434 writer = proc_unit.addOperation(name='HDFWriter')
434 435 writer.addParameter(name='path', value='/path/to/file')
435 436 writer.addParameter(name='blocksPerFile', value='32')
436 437 writer.addParameter(name='metadataList', value='heightList,timeZone')
437 438 writer.addParameter(name='dataList',value='data_output,utctime')
438 439 # writer.addParameter(name='description',value=json.dumps(desc))
439 440
440 441 """
441 442
442 443 ext = ".hdf5"
443 444 optchar = "D"
444 445 filename = None
445 446 path = None
446 447 setFile = None
447 448 fp = None
448 449 ds = None
449 450 firsttime = True
450 451 #Configurations
451 452 blocksPerFile = None
452 453 blockIndex = None
453 454 dataOut = None #eval ??????
454 455 #Data Arrays
455 456 dataList = None
456 457 metadataList = None
457 458 currentDay = None
458 459 lastTime = None
459 460 timeZone = "ut"
460 461 hourLimit = 3
461 462 breakDays = True
462 463
463 464 def __init__(self):
464 465
465 466 Operation.__init__(self)
466 467 return
467 468
468 469 def set_kwargs(self, **kwargs):
469 470
470 471 for key, value in kwargs.items():
471 472 setattr(self, key, value)
472 473
473 474 def set_kwargs_obj(self, obj, **kwargs):
474 475
475 476 for key, value in kwargs.items():
476 477 setattr(obj, key, value)
477 478
478 479 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None,
479 480 description={},timeZone = "ut",hourLimit = 3, breakDays=True, **kwargs):
480 481 self.path = path
481 482 self.blocksPerFile = blocksPerFile
482 483 self.metadataList = metadataList
483 484 self.dataList = [s.strip() for s in dataList]
484 485 self.setType = setType
485 486 self.description = description
486 487 self.timeZone = timeZone
487 488 self.hourLimit = hourLimit
488 489 self.breakDays = breakDays
489 490 self.set_kwargs(**kwargs)
490 491
491 492 if self.metadataList is None:
492 493 self.metadataList = self.dataOut.metadata_list
493 494
494 495 self.metadataList = list(set(self.metadataList))
495 496
496 497 tableList = []
497 498 dsList = []
498 499
499 500 for i in range(len(self.dataList)):
500 501 dsDict = {}
501 502 if hasattr(self.dataOut, self.dataList[i]):
502 503 dataAux = getattr(self.dataOut, self.dataList[i])
503 504 dsDict['variable'] = self.dataList[i]
504 505 else:
505 506 log.warning('Attribute {} not found in dataOut'.format(self.dataList[i]),self.name)
506 507 continue
507 508
508 509 if dataAux is None:
509 510 continue
510 511 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float_)):
511 512 dsDict['nDim'] = 0
512 513 else:
513 514 dsDict['nDim'] = len(dataAux.shape)
514 515 dsDict['shape'] = dataAux.shape
515 516 dsDict['dsNumber'] = dataAux.shape[0]
516 517 dsDict['dtype'] = dataAux.dtype
517 518
518 519 dsList.append(dsDict)
519 520
520 521 self.blockIndex = 0
521 522 self.dsList = dsList
522 523 self.currentDay = self.dataOut.datatime.date()
523 524
524 525 def timeFlag(self):
525 526 currentTime = self.dataOut.utctime
526 527 timeTuple = None
527 528 if self.timeZone == "lt":
528 529 timeTuple = time.localtime(currentTime)
529 530 else :
530 531 timeTuple = time.gmtime(currentTime)
531 532 dataDay = timeTuple.tm_yday
532 533
533 534 if self.lastTime is None:
534 535 self.lastTime = currentTime
535 536 self.currentDay = dataDay
536 537 return False
537 538
538 539 timeDiff = currentTime - self.lastTime
539 540
540 541 # Si el dia es diferente o si la diferencia entre un
541 542 # dato y otro supera self.hourLimit
542 543 if (dataDay != self.currentDay) and self.breakDays:
543 544 self.currentDay = dataDay
544 545 return True
545 546 elif timeDiff > self.hourLimit*60*60:
546 547 self.lastTime = currentTime
547 548 return True
548 549 else:
549 550 self.lastTime = currentTime
550 551 return False
551 552
552 553 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
553 554 dataList=[], setType=None, description={}, **kwargs):
554 555
555 556 self.dataOut = dataOut
556 557 self.set_kwargs_obj(self.dataOut, **kwargs)
557 558 if not(self.isConfig):
558 559 self.setup(path=path, blocksPerFile=blocksPerFile,
559 560 metadataList=metadataList, dataList=dataList,
560 561 setType=setType, description=description, **kwargs)
561 562
562 563 self.isConfig = True
563 564 self.setNextFile()
564 565
565 566 self.putData()
566 567 return
567 568
568 569 def setNextFile(self):
569 570
570 571 ext = self.ext
571 572 path = self.path
572 573 setFile = self.setFile
573 574 timeTuple = None
574 575 if self.timeZone == "lt":
575 576 timeTuple = time.localtime(self.dataOut.utctime)
576 577 elif self.timeZone == "ut":
577 578 timeTuple = time.gmtime(self.dataOut.utctime)
578 579 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
579 580 fullpath = os.path.join(path, subfolder)
580 581
581 582 if os.path.exists(fullpath):
582 583 filesList = os.listdir(fullpath)
583 584 filesList = [k for k in filesList if k.startswith(self.optchar)]
584 585 if len(filesList) > 0:
585 586 filesList = sorted(filesList, key=str.lower)
586 587 filen = filesList[-1]
587 588 # el filename debera tener el siguiente formato
588 589 # 0 1234 567 89A BCDE (hex)
589 590 # x YYYY DDD SSS .ext
590 591 if isNumber(filen[8:11]):
591 592 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
592 593 else:
593 594 setFile = -1
594 595 else:
595 596 setFile = -1 #inicializo mi contador de seteo
596 597 else:
597 598 os.makedirs(fullpath)
598 599 setFile = -1 #inicializo mi contador de seteo
599 600
600 601 if self.setType is None:
601 602 setFile += 1
602 603 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
603 604 timeTuple.tm_year,
604 605 timeTuple.tm_yday,
605 606 setFile,
606 607 ext)
607 608 else:
608 609 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
609 610 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
610 611 timeTuple.tm_year,
611 612 timeTuple.tm_yday,
612 613 setFile,
613 614 ext)
614 615
615 616 self.filename = os.path.join(path, subfolder, file)
616 617
617 618
618 619
619 620 def getLabel(self, name, x=None):
620 621
621 622 if x is None:
622 623 if 'Data' in self.description:
623 624 data = self.description['Data']
624 625 if 'Metadata' in self.description:
625 626 data.update(self.description['Metadata'])
626 627 else:
627 628 data = self.description
628 629 if name in data:
629 630 if isinstance(data[name], str):
630 631 return data[name]
631 632 elif isinstance(data[name], list):
632 633 return None
633 634 elif isinstance(data[name], dict):
634 635 for key, value in data[name].items():
635 636 return key
636 637 return name
637 638 else:
638 639 if 'Metadata' in self.description:
639 640 meta = self.description['Metadata']
640 641 else:
641 642 meta = self.description
642 643 if name in meta:
643 644 if isinstance(meta[name], list):
644 645 return meta[name][x]
645 646 elif isinstance(meta[name], dict):
646 647 for key, value in meta[name].items():
647 648 return value[x]
648 649 if 'cspc' in name:
649 650 return 'pair{:02d}'.format(x)
650 651 else:
651 652 return 'channel{:02d}'.format(x)
652 653
653 654 def writeMetadata(self, fp):
654 655
655 656 if self.description:
656 657 if 'Metadata' in self.description:
657 658 grp = fp.create_group('Metadata')
658 659 else:
659 660 grp = fp
660 661 else:
661 662 grp = fp.create_group('Metadata')
662 663
663 664 for i in range(len(self.metadataList)):
664 665 if not hasattr(self.dataOut, self.metadataList[i]):
665 666 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
666 667 continue
667 668 value = getattr(self.dataOut, self.metadataList[i])
668 669 if isinstance(value, bool):
669 670 if value is True:
670 671 value = 1
671 672 else:
672 673 value = 0
673 674 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
674 675 return
675 676
676 677 def writeMetadata2(self, fp):
677 678
678 679 if self.description:
679 680 if 'Metadata' in self.description:
680 681 grp = fp.create_group('Metadata')
681 682 else:
682 683 grp = fp
683 684 else:
684 685 grp = fp.create_group('Metadata')
685 686
686 687 for i in range(len(self.metadataList)):
687 688
688 689 attribute = self.metadataList[i]
689 690 attr = attribute.split('.')
690 691 if len(attr) > 1:
691 692 if not hasattr(eval("self.dataOut."+attr[0]),attr[1]):
692 693 log.warning('Metadata: {}.{} not found'.format(attr[0],attr[1]), self.name)
693 694 continue
694 695 value = getattr(eval("self.dataOut."+attr[0]),attr[1])
695 696 if isinstance(value, bool):
696 697 if value is True:
697 698 value = 1
698 699 else:
699 700 value = 0
700 701 if isinstance(value,type(None)):
701 702 log.warning("Invalid value detected, {} is None".format(attribute), self.name)
702 703 value = 0
703 704 grp2 = None
704 705 if not 'Metadata/'+attr[0] in fp:
705 706 grp2 = fp.create_group('Metadata/'+attr[0])
706 707 else:
707 708 grp2 = fp['Metadata/'+attr[0]]
708 709 grp2.create_dataset(attr[1], data=value)
709 710
710 711 else:
711 712 if not hasattr(self.dataOut, attr[0] ):
712 713 log.warning('Metadata: `{}` not found'.format(attribute), self.name)
713 714 continue
714 715 value = getattr(self.dataOut, attr[0])
715 716 if isinstance(value, bool):
716 717 if value is True:
717 718 value = 1
718 719 else:
719 720 value = 0
720 721 if isinstance(value, type(None)):
721 722 log.error("Value {} is None".format(attribute),self.name)
722 723
723 724 grp.create_dataset(self.getLabel(attribute), data=value)
724 725
725 726 return
726 727
727 728 def writeData(self, fp):
728 729
729 730 if self.description:
730 731 if 'Data' in self.description:
731 732 grp = fp.create_group('Data')
732 733 else:
733 734 grp = fp
734 735 else:
735 736 grp = fp.create_group('Data')
736 737
737 738 dtsets = []
738 739 data = []
739 740
740 741 for dsInfo in self.dsList:
741 742 if dsInfo['nDim'] == 0:
742 743 ds = grp.create_dataset(
743 744 self.getLabel(dsInfo['variable']),
744 745 (self.blocksPerFile,),
745 746 chunks=True,
746 747 dtype=numpy.float64)
747 748 dtsets.append(ds)
748 749 data.append((dsInfo['variable'], -1))
749 750 else:
750 751 label = self.getLabel(dsInfo['variable'])
751 752 if label is not None:
752 753 sgrp = grp.create_group(label)
753 754 else:
754 755 sgrp = grp
755 756 for i in range(dsInfo['dsNumber']):
756 757 ds = sgrp.create_dataset(
757 758 self.getLabel(dsInfo['variable'], i),
758 759 (self.blocksPerFile,) + dsInfo['shape'][1:],
759 760 chunks=True,
760 761 dtype=dsInfo['dtype'])
761 762 dtsets.append(ds)
762 763 data.append((dsInfo['variable'], i))
763 764 fp.flush()
764 765
765 766 log.log('Creating file: {}'.format(fp.filename), self.name)
766 767
767 768 self.ds = dtsets
768 769 self.data = data
769 770 self.firsttime = True
770 771
771 772 return
772 773
773 774 def putData(self):
774 775
775 776 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
776 777 self.closeFile()
777 778 self.setNextFile()
778 779 self.dataOut.flagNoData = False
779 780 self.blockIndex = 0
780 781
781 782 if self.blockIndex == 0:
782 783 #Setting HDF5 File
783 784 self.fp = h5py.File(self.filename, 'w')
784 785 #write metadata
785 786 self.writeMetadata2(self.fp)
786 787 #Write data
787 788 self.writeData(self.fp)
788 789 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
789 790 elif (self.blockIndex % 10 ==0):
790 791 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
791 792 else:
792 793
793 794 log.log('Block No. {}/{}'.format(self.blockIndex+1, self.blocksPerFile), self.name)
794 795
795 796 for i, ds in enumerate(self.ds):
796 797 attr, ch = self.data[i]
797 798 if ch == -1:
798 799 ds[self.blockIndex] = getattr(self.dataOut, attr)
799 800 else:
800 801 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
801 802
802 803 self.blockIndex += 1
803 804
804 805 self.fp.flush()
805 806 self.dataOut.flagNoData = True
806 807
807 808 def closeFile(self):
808 809
809 810 if self.blockIndex != self.blocksPerFile:
810 811 for ds in self.ds:
811 812 ds.resize(self.blockIndex, axis=0)
812 813
813 814 if self.fp:
814 815 self.fp.flush()
815 816 self.fp.close()
816 817
817 818 def close(self):
818 819
819 820 self.closeFile()
Show file before | Show file after | Hide comments
@@ -1,565 +1,576
1 1 """
2 2 Utilities for IO modules
3 3 @modified: Joab Apaza
4 4 @email: roj-op01@igp.gob.pe, joab.apaza32@gmail.com
5 5 """
6 6 ################################################################################
7 7 ################################################################################
8 8 import os
9 9 from datetime import datetime
10 10 import numpy
11 11 from schainpy.model.data.jrodata import Parameters
12 12 import itertools
13 13 import numpy
14 14 import h5py
15 15 import re
16 16 import time
17 17 from schainpy.utils import log
18 18 ################################################################################
19 19 ################################################################################
20 20 ################################################################################
21 21 def folder_in_range(folder, start_date, end_date, pattern):
22 22 """
23 23 Check whether folder is bettwen start_date and end_date
24 24 Args:
25 25 folder (str): Folder to check
26 26 start_date (date): Initial date
27 27 end_date (date): Final date
28 28 pattern (str): Datetime format of the folder
29 29 Returns:
30 30 bool: True for success, False otherwise
31 31 """
32 32 try:
33 33 dt = datetime.strptime(folder, pattern)
34 34 except:
35 35 raise ValueError('Folder {} does not match {} format'.format(folder, pattern))
36 36 return start_date <= dt.date() <= end_date
37 37 ################################################################################
38 38 ################################################################################
39 39 ################################################################################
40 40 def getHei_index( minHei, maxHei, heightList):
41 41 try:
42 42 if (minHei < heightList[0]):
43 43 minHei = heightList[0]
44 44 if (maxHei > heightList[-1]):
45 45 maxHei = heightList[-1]
46 46 minIndex = 0
47 47 maxIndex = 0
48 48 heights = numpy.asarray(heightList)
49 49 inda = numpy.where(heights >= minHei)
50 50 indb = numpy.where(heights <= maxHei)
51 51 try:
52 52 minIndex = inda[0][0]
53 53 except:
54 54 minIndex = 0
55 55 try:
56 56 maxIndex = indb[0][-1]
57 57 except:
58 58 maxIndex = len(heightList)
59 59 return minIndex,maxIndex
60 60 except Exception as e:
61 61 log.error("In getHei_index: ", __name__)
62 62 log.error(e , __name__)
63 63 ################################################################################
64 64 ################################################################################
65 65 ################################################################################
66 66 class MergeH5(object):
67 67 """Processing unit to read HDF5 format files
68 68 This unit reads HDF5 files created with `HDFWriter` operation when channels area
69 69 processed by separated. Then merge all channels in a single files.
70 70 "example"
71 71 nChannels = 4
72 72 pathOut = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/merged"
73 73 p0 = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/d2022240_Ch0"
74 74 p1 = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/d2022240_Ch1"
75 75 p2 = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/d2022240_Ch2"
76 76 p3 = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/d2022240_Ch3"
77 77 list = ['data_spc','data_cspc','nIncohInt','utctime']
78 78 merger = MergeH5(nChannels,pathOut,list, p0, p1,p2,p3)
79 79 merger.run()
80 80 The file example_FULLmultiprocessing_merge.txt show an application for AMISR data
81 81 """
82 82 # #__attrs__ = ['paths', 'nChannels']
83 83 isConfig = False
84 84 inPaths = None
85 85 nChannels = None
86 86 ch_dataIn = []
87 87 channelList = []
88 88 def __init__(self,nChannels, pOut, dataList, *args):
89 89 self.inPaths = [p for p in args]
90 90 #print(self.inPaths)
91 91 if len(self.inPaths) != nChannels:
92 92 print("ERROR, number of channels different from iput paths {} != {}".format(nChannels, len(args)))
93 93 return
94 94 self.pathOut = pOut
95 95 self.dataList = dataList
96 96 self.nChannels = len(self.inPaths)
97 97 self.ch_dataIn = [Parameters() for p in args]
98 98 self.dataOut = Parameters()
99 99 self.channelList = [n for n in range(nChannels)]
100 100 self.blocksPerFile = None
101 101 self.date = None
102 102 self.ext = ".hdf5$"
103 103 self.dataList = dataList
104 104 self.optchar = "D"
105 105 self.meta = {}
106 106 self.data = {}
107 107 self.open_file = h5py.File
108 108 self.open_mode = 'r'
109 109 self.description = {}
110 110 self.extras = {}
111 111 self.filefmt = "*%Y%j***"
112 112 self.folderfmt = "*%Y%j"
113 113 self.flag_spc = False
114 114 self.flag_pow = False
115 115 self.flag_snr = False
116 116 self.flag_nIcoh = False
117 117 self.flagProcessingHeader = False
118 118 self.flagControllerHeader = False
119 119 def setup(self):
120 120 # if not self.ext.startswith('.'):
121 121 # self.ext = '.{}'.format(self.ext)
122 122 self.filenameList = self.searchFiles(self.inPaths, None)
123 123 self.nfiles = len(self.filenameList[0])
124 124 def searchFiles(self, paths, date, walk=True):
125 125 # self.paths = path
126 126 #self.date = startDate
127 127 #self.walk = walk
128 128 filenameList = [[] for n in range(self.nChannels)]
129 129 ch = 0
130 130 for path in paths:
131 131 if os.path.exists(path):
132 132 print("Searching files in {}".format(path))
133 133 filenameList[ch] = self.getH5files(path, walk)
134 134 print("Found: ")
135 135 for f in filenameList[ch]:
136 136 print(f)
137 137 else:
138 138 self.status = 0
139 139 print('Path:%s does not exists'%path)
140 140 return 0
141 141 ch+=1
142 142 return filenameList
143 143 def getH5files(self, path, walk):
144 144 dirnameList = []
145 145 pat = '(\d)+.'+self.ext
146 146 if walk:
147 147 for root, dirs, files in os.walk(path):
148 148 for dir in dirs:
149 149 #print(os.path.join(root,dir))
150 150 files = [re.search(pat,x) for x in os.listdir(os.path.join(root,dir))]
151 151 #print(files)
152 152 files = [x for x in files if x!=None]
153 153 files = [x.string for x in files]
154 154 files = [os.path.join(root,dir,x) for x in files if x!=None]
155 155 files.sort()
156 156 dirnameList += files
157 157 return dirnameList
158 158 else:
159 159 dirnameList = [re.search(pat,x) for x in os.listdir(path)]
160 160 dirnameList = [x for x in dirnameList if x!=None]
161 161 dirnameList = [x.string for x in dirnameList]
162 162 dirnameList = [x for x in dirnameList if x!=None]
163 163 dirnameList.sort()
164 164 return dirnameList
165 165 def readFile(self,fp,ch):
166 166 '''Read metadata and data'''
167 167 self.readMetadata(fp,ch)
168 #print(self.metadataList)
168 # print(self.metadataList)
169 169 data = self.readData(fp)
170 170 for attr in self.meta:
171 171 if "processingHeaderObj" in attr:
172 172 self.flagProcessingHeader=True
173 173 if "radarControllerHeaderObj" in attr:
174 174 self.flagControllerHeader=True
175 175 at = attr.split('.')
176 #print("AT ", at)
176 # print("AT ", at)
177 177 if len(at) > 1:
178 178 setattr(eval("self.ch_dataIn[ch]."+at[0]),at[1], self.meta[attr])
179 179 else:
180 180 setattr(self.ch_dataIn[ch], attr, self.meta[attr])
181 181 self.fill_dataIn(data, self.ch_dataIn[ch])
182 182 return
183 183 def readMetadata(self, fp, ch):
184 184 '''
185 185 Reads Metadata
186 186 '''
187 187 meta = {}
188 188 self.metadataList = []
189 189 grp = fp['Metadata']
190 190 for item in grp.values():
191 191 name = item.name
192 192 if isinstance(item, h5py.Dataset):
193 193 name = name.split("/")[-1]
194 194 if 'List' in name:
195 195 meta[name] = item[()].tolist()
196 196 else:
197 197 meta[name] = item[()]
198 198 self.metadataList.append(name)
199 199 else:
200 200 grp2 = fp[name]
201 201 Obj = name.split("/")[-1]
202 202 #print(Obj)
203 203 for item2 in grp2.values():
204 204 name2 = Obj+"."+item2.name.split("/")[-1]
205 205 if 'List' in name2:
206 206 meta[name2] = item2[()].tolist()
207 207 else:
208 208 meta[name2] = item2[()]
209 209 self.metadataList.append(name2)
210 210 if not self.meta:
211 211 self.meta = meta.copy()
212 212 for key in list(self.meta.keys()):
213 213 if "channelList" in key:
214 214 self.meta["channelList"] =[n for n in range(self.nChannels)]
215 215 if "processingHeaderObj" in key:
216 216 self.meta["processingHeaderObj.channelList"] =[n for n in range(self.nChannels)]
217 217 if "radarControllerHeaderObj" in key:
218 218 self.meta["radarControllerHeaderObj.channelList"] =[n for n in range(self.nChannels)]
219 219 return 1
220 220 else:
221 221 for k in list(self.meta.keys()):
222 222 if 'List' in k and 'channel' not in k and "height" not in k and "radarControllerHeaderObj" not in k:
223 223 self.meta[k] += meta[k]
224 224 #print("Metadata: ",self.meta)
225 225 return 1
226 226 def fill_dataIn(self,data, dataIn):
227 227 for attr in data:
228 228 if data[attr].ndim == 1:
229 229 setattr(dataIn, attr, data[attr][:])
230 230 else:
231 231 setattr(dataIn, attr, numpy.squeeze(data[attr][:,:]))
232 #print("shape in", dataIn.data_spc.shape, len(dataIn.data_spc))
232 # print("shape in", dataIn.data_spc.shape, len(dataIn.data_spc))
233 233 if self.flag_spc:
234 234 if dataIn.data_spc.ndim > 3:
235 235 dataIn.data_spc = dataIn.data_spc[0]
236 236 #print("shape in", dataIn.data_spc.shape)
237 237 def getBlocksPerFile(self):
238 238 b = numpy.zeros(self.nChannels)
239 239 for i in range(self.nChannels):
240 240 if self.flag_spc:
241 241 b[i] = self.ch_dataIn[i].data_spc.shape[0] #number of blocks
242 242 elif self.flag_pow:
243 243 b[i] = self.ch_dataIn[i].data_pow.shape[0] #number of blocks
244 244 elif self.flag_snr:
245 245 b[i] = self.ch_dataIn[i].data_snr.shape[0] #number of blocks
246 246 self.blocksPerFile = int(b.min())
247 247 iresh_ch = numpy.where(b > self.blocksPerFile)[0]
248 248 if len(iresh_ch) > 0:
249 249 for ich in iresh_ch:
250 250 for i in range(len(self.dataList)):
251 251 if hasattr(self.ch_dataIn[ich], self.dataList[i]):
252 252 # print("reshaping ", self.dataList[i])
253 253 # print(getattr(self.ch_dataIn[ich], self.dataList[i]).shape)
254 254 dataAux = getattr(self.ch_dataIn[ich], self.dataList[i])
255 255 setattr(self.ch_dataIn[ich], self.dataList[i], None)
256 256 setattr(self.ch_dataIn[ich], self.dataList[i], dataAux[0:self.blocksPerFile])
257 257 # print(getattr(self.ch_dataIn[ich], self.dataList[i]).shape)
258 258 else:
259 # log.error("Channels number error,iresh_ch=", iresh_ch)
259 260 return
260 261 def getLabel(self, name, x=None):
261 262 if x is None:
262 263 if 'Data' in self.description:
263 264 data = self.description['Data']
264 265 if 'Metadata' in self.description:
265 266 data.update(self.description['Metadata'])
266 267 else:
267 268 data = self.description
268 269 if name in data:
269 270 if isinstance(data[name], str):
270 271 return data[name]
271 272 elif isinstance(data[name], list):
272 273 return None
273 274 elif isinstance(data[name], dict):
274 275 for key, value in data[name].items():
275 276 return key
276 277 return name
277 278 else:
278 279 if 'Metadata' in self.description:
279 280 meta = self.description['Metadata']
280 281 else:
281 282 meta = self.description
282 283 if name in meta:
283 284 if isinstance(meta[name], list):
284 285 return meta[name][x]
285 286 elif isinstance(meta[name], dict):
286 287 for key, value in meta[name].items():
287 288 return value[x]
288 289 if 'cspc' in name:
289 290 return 'pair{:02d}'.format(x)
290 291 else:
291 292 return 'channel{:02d}'.format(x)
293
292 294 def readData(self, fp):
293 #print("read fp: ", fp)
295 # print("read fp: ", fp)
294 296 data = {}
295 297 grp = fp['Data']
296 298 for name in grp:
297 299 if "spc" in name:
298 300 self.flag_spc = True
299 301 if "pow" in name:
300 302 self.flag_pow = True
301 303 if "snr" in name:
302 304 self.flag_snr = True
303 305 if "nIncohInt" in name:
304 306 self.flag_nIcoh = True
305
307 # print("spc:",self.flag_spc," pow:",self.flag_pow," snr:", self.flag_snr)
306 308 if isinstance(grp[name], h5py.Dataset):
307 309 array = grp[name][()]
308 310 elif isinstance(grp[name], h5py.Group):
309 311 array = []
310 312 for ch in grp[name]:
311 313 array.append(grp[name][ch][()])
312 314 array = numpy.array(array)
313 315 else:
314 316 print('Unknown type: {}'.format(name))
315 317 data[name] = array
316 318 return data
319
317 320 def getDataOut(self):
321 # print("Getting DataOut")
318 322 self.dataOut = self.ch_dataIn[0].copy() #dataIn #blocks, fft, hei for metadata
319 323 if self.flagProcessingHeader:
320 324 self.dataOut.processingHeaderObj = self.ch_dataIn[0].processingHeaderObj.copy()
321 325 self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
322 326 self.dataOut.ippSeconds = self.dataOut.processingHeaderObj.ipp
323 327 self.dataOut.channelList = self.dataOut.processingHeaderObj.channelList
324 328 self.dataOut.nCohInt = self.dataOut.processingHeaderObj.nCohInt
325 329 self.dataOut.nFFTPoints = self.dataOut.processingHeaderObj.nFFTPoints
326 330 if self.flagControllerHeader:
327 331 self.dataOut.radarControllerHeaderObj = self.ch_dataIn[0].radarControllerHeaderObj.copy()
328 332 self.dataOut.frequency = self.dataOut.radarControllerHeaderObj.frequency
329 333 #--------------------------------------------------------------------
330 334 #--------------------------------------------------------------------
331 335 if self.flag_spc:
332 336 if self.dataOut.data_spc.ndim < 3:
333 337 print("shape spc in: ",self.dataOut.data_spc.shape )
334 338 return 0
335 339 if self.flag_pow:
336 340 if self.dataOut.data_pow.ndim < 2:
337 341 print("shape pow in: ",self.dataOut.data_pow.shape )
338 342 return 0
339 343 if self.flag_snr:
340 344 if self.dataOut.data_snr.ndim < 2:
341 345 print("shape snr in: ",self.dataOut.data_snr.shape )
342 346 return 0
343 347 self.dataOut.data_spc = None
344 348 self.dataOut.data_cspc = None
345 349 self.dataOut.data_pow = None
346 350 self.dataOut.data_snr = None
347 351 self.dataOut.utctime = None
348 352 self.dataOut.nIncohInt = None
349 353 #--------------------------------------------------------------------
350 354 if self.flag_spc:
351 355 spc = [data.data_spc for data in self.ch_dataIn]
352 356 self.dataOut.data_spc = numpy.stack(spc, axis=1) #blocks, ch, fft, hei
353 357 #--------------------------------------------------------------------
354 358 if self.flag_pow:
355 359 pow = [data.data_pow for data in self.ch_dataIn]
356 360 self.dataOut.data_pow = numpy.stack(pow, axis=1) #blocks, ch, fft, hei
357 361 #--------------------------------------------------------------------
358 362 if self.flag_snr:
359 363 snr = [data.data_snr for data in self.ch_dataIn]
360 364 self.dataOut.data_snr = numpy.stack(snr, axis=1) #blocks, ch, fft, hei
361 365 #--------------------------------------------------------------------
362 366 time = [data.utctime for data in self.ch_dataIn]
363 367 time = numpy.asarray(time).mean(axis=0)
364 368 time = numpy.squeeze(time)
365 369 self.dataOut.utctime = time
366 370 #--------------------------------------------------------------------
367 371 if self.flag_nIcoh:
368 372 ints = [data.nIncohInt for data in self.ch_dataIn]
369 373 self.dataOut.nIncohInt = numpy.stack(ints, axis=1)
370 374 if self.dataOut.nIncohInt.ndim > 3:
371 375 aux = self.dataOut.nIncohInt
372 376 self.dataOut.nIncohInt = None
373 377 self.dataOut.nIncohInt = aux[0]
374 378 if self.dataOut.nIncohInt.ndim < 3:
375 379 nIncohInt = numpy.repeat(self.dataOut.nIncohInt, self.dataOut.nHeights).reshape(self.blocksPerFile,self.nChannels, self.dataOut.nHeights)
376 380 #nIncohInt = numpy.reshape(nIncohInt, (self.blocksPerFile,self.nChannels, self.dataOut.nHeights))
377 381 self.dataOut.nIncohInt = None
378 382 self.dataOut.nIncohInt = nIncohInt
379 383 if (self.dataOut.nIncohInt.shape)[0]==self.nChannels: ## ch,blocks, hei
380 384 self.dataOut.nIncohInt = numpy.swapaxes(self.dataOut.nIncohInt, 0, 1) ## blocks,ch, hei
381 385 else:
382 386 self.dataOut.nIncohInt = self.ch_dataIn[0].nIncohInt
383 387 #--------------------------------------------------------------------
384 #print("utcTime: ", time.shape)
385 #print("data_spc ",self.dataOut.data_spc.shape)
388 # print("utcTime: ", time.shape)
389 # print("data_spc ",self.dataOut.data_spc.shape)
386 390 if "data_cspc" in self.dataList:
387 391 pairsList = [pair for pair in itertools.combinations(self.channelList, 2)]
388 392 #print("PairsList: ", pairsList)
389 393 self.dataOut.pairsList = pairsList
390 394 cspc = []
391 395 for i, j in pairsList:
392 396 cspc.append(self.ch_dataIn[i].data_spc*numpy.conjugate(self.ch_dataIn[j].data_spc)) #blocks, fft, hei
393 397 cspc = numpy.asarray(cspc) # # pairs, blocks, fft, hei
394 398 #print("cspc: ",cspc.shape)
395 399 self.dataOut.data_cspc = numpy.swapaxes(cspc, 0, 1) ## blocks, pairs, fft, hei
396 400 #print("dataOut.data_cspc: ",self.dataOut.data_cspc.shape)
397 401 #if "data_pow" in self.dataList:
398 402 return 1
399 403 # def writeMetadata(self, fp):
400 404 #
401 405 #
402 406 # grp = fp.create_group('Metadata')
403 407 #
404 408 # for i in range(len(self.metadataList)):
405 409 # if not hasattr(self.dataOut, self.metadataList[i]):
406 410 # print('Metadata: `{}` not found'.format(self.metadataList[i]))
407 411 # continue
408 412 # value = getattr(self.dataOut, self.metadataList[i])
409 413 # if isinstance(value, bool):
410 414 # if value is True:
411 415 # value = 1
412 416 # else:
413 417 # value = 0
414 418 # grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
415 419 # return
416 420 def writeMetadata(self, fp):
417 421 grp = fp.create_group('Metadata')
418 422 for i in range(len(self.metadataList)):
419 423 attribute = self.metadataList[i]
420 424 attr = attribute.split('.')
421 425 if '' in attr:
422 426 attr.remove('')
423 427 #print(attr)
424 428 if len(attr) > 1:
425 429 if not hasattr(eval("self.dataOut."+attr[0]),attr[1]):
426 430 print('Metadata: {}.{} not found'.format(attr[0],attr[1]))
427 431 continue
428 432 value = getattr(eval("self.dataOut."+attr[0]),attr[1])
429 433 if isinstance(value, bool):
430 434 if value is True:
431 435 value = 1
432 436 else:
433 437 value = 0
434 438 grp2 = None
435 439 if not 'Metadata/'+attr[0] in fp:
436 440 grp2 = fp.create_group('Metadata/'+attr[0])
437 441 else:
438 442 grp2 = fp['Metadata/'+attr[0]]
439 443 grp2.create_dataset(attr[1], data=value)
440 444 else:
441 445 if not hasattr(self.dataOut, attr[0] ):
442 446 print('Metadata: `{}` not found'.format(attribute))
443 447 continue
444 448 value = getattr(self.dataOut, attr[0])
445 449 if isinstance(value, bool):
446 450 if value is True:
447 451 value = 1
448 452 else:
449 453 value = 0
450 454 if isinstance(value, type(None)):
451 455 print("------ERROR, value {} is None".format(attribute))
452 456
453 457 grp.create_dataset(self.getLabel(attribute), data=value)
454 458 return
459
455 460 def getDsList(self):
461 # print("Getting DS List", self.dataList)
456 462 dsList =[]
463 dataAux = None
457 464 for i in range(len(self.dataList)):
458 465 dsDict = {}
459 466 if hasattr(self.dataOut, self.dataList[i]):
460 467 dataAux = getattr(self.dataOut, self.dataList[i])
461 468 dsDict['variable'] = self.dataList[i]
462 469 else:
463 470 print('Attribute {} not found in dataOut'.format(self.dataList[i]))
464 471 continue
465 472 if dataAux is None:
466 473 continue
467 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
474 elif isinstance(dataAux, (int, float, numpy.int_, numpy.float_)):
468 475 dsDict['nDim'] = 0
469 476 else:
470 477 dsDict['nDim'] = len(dataAux.shape) -1
471 478 dsDict['shape'] = dataAux.shape
472 479 if len(dsDict['shape'])>=2:
473 480 dsDict['dsNumber'] = dataAux.shape[1]
474 481 else:
475 482 dsDict['dsNumber'] = 1
476 483 dsDict['dtype'] = dataAux.dtype
477 484 # if len(dataAux.shape) == 4:
478 485 # dsDict['nDim'] = len(dataAux.shape) -1
479 486 # dsDict['shape'] = dataAux.shape
480 487 # dsDict['dsNumber'] = dataAux.shape[1]
481 488 # dsDict['dtype'] = dataAux.dtype
482 489 # else:
483 490 # dsDict['nDim'] = len(dataAux.shape)
484 491 # dsDict['shape'] = dataAux.shape
485 492 # dsDict['dsNumber'] = dataAux.shape[0]
486 493 # dsDict['dtype'] = dataAux.dtype
487 494 dsList.append(dsDict)
488 #print(dsList)
495 # print("dsList: ", dsList)
489 496 self.dsList = dsList
497
490 498 def clean_dataIn(self):
491 499 for ch in range(self.nChannels):
492 500 self.ch_dataIn[ch].data_spc = None
493 501 self.ch_dataIn[ch].utctime = None
494 502 self.ch_dataIn[ch].nIncohInt = None
495 503 self.meta ={}
496 504 self.blocksPerFile = None
505
497 506 def writeData(self, outFilename):
498 507 self.getDsList()
499 508 fp = h5py.File(outFilename, 'w')
509 # print("--> Merged file: ",fp)
500 510 self.writeMetadata(fp)
501 511 grp = fp.create_group('Data')
502 512 dtsets = []
503 513 data = []
504 514 for dsInfo in self.dsList:
505 515 if dsInfo['nDim'] == 0:
506 516 ds = grp.create_dataset(
507 517 self.getLabel(dsInfo['variable']),(self.blocksPerFile, ),chunks=True,dtype=numpy.float64)
508 518 dtsets.append(ds)
509 519 data.append((dsInfo['variable'], -1))
510 520 else:
511 521 label = self.getLabel(dsInfo['variable'])
512 522 if label is not None:
513 523 sgrp = grp.create_group(label)
514 524 else:
515 525 sgrp = grp
516 526 k = -1*(dsInfo['nDim'] - 1)
517 #print(k, dsInfo['shape'], dsInfo['shape'][k:])
527 # print(k, dsInfo['shape'], dsInfo['shape'][k:])
518 528 for i in range(dsInfo['dsNumber']):
519 529 ds = sgrp.create_dataset(
520 530 self.getLabel(dsInfo['variable'], i),(self.blocksPerFile, ) + dsInfo['shape'][k:],
521 531 chunks=True,
522 532 dtype=dsInfo['dtype'])
523 533 dtsets.append(ds)
524 534 data.append((dsInfo['variable'], i))
525 535 #print("\n",dtsets)
526 536 print('Creating merged file: {}'.format(fp.filename))
527 537 for i, ds in enumerate(dtsets):
528 538 attr, ch = data[i]
529 539 if ch == -1:
530 540 ds[:] = getattr(self.dataOut, attr)
531 541 else:
532 542 #print(ds, getattr(self.dataOut, attr)[ch].shape)
533 543 aux = getattr(self.dataOut, attr)# block, ch, ...
534 544 aux = numpy.swapaxes(aux,0,1) # ch, blocks, ...
535 545 #print(ds.shape, aux.shape)
536 546 #ds[:] = getattr(self.dataOut, attr)[ch]
537 547 ds[:] = aux[ch]
538 548 fp.flush()
539 549 fp.close()
540 550 self.clean_dataIn()
541 551 return
552
542 553 def run(self):
543 554 if not(self.isConfig):
544 555 self.setup()
545 556 self.isConfig = True
546 557 for nf in range(self.nfiles):
547 558 name = None
548 559 for ch in range(self.nChannels):
549 560 name = self.filenameList[ch][nf]
550 561 filename = os.path.join(self.inPaths[ch], name)
551 562 fp = h5py.File(filename, 'r')
552 #print("Opening file: ",filename)
563 print("Opening file: ",filename)
553 564 self.readFile(fp,ch)
554 565 fp.close()
555 566 if self.blocksPerFile == None:
556 567 self.getBlocksPerFile()
557 568 print("blocks per file: ", self.blocksPerFile)
558 569 if not self.getDataOut():
559 570 print("Error getting DataOut invalid number of blocks")
560 571 return
561 572 name = name[-16:]
562 #print("Final name out: ", name)
573 # print("Final name out: ", name)
563 574 outFile = os.path.join(self.pathOut, name)
564 #print("Outfile: ", outFile)
575 # print("Outfile: ", outFile)
565 576 self.writeData(outFile) No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1737 +1,1738
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
16 16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 17 from schainpy.model.data.jrodata import Spectra
18 18 from schainpy.model.data.jrodata import hildebrand_sekhon
19 19 from schainpy.model.data import _noise
20 20 from schainpy.utils import log
21 21 import matplotlib.pyplot as plt
22 22 from schainpy.model.io.utilsIO import getHei_index
23 23 import datetime
24 24
25 25 class SpectraProc(ProcessingUnit):
26 26
27 27 def __init__(self):
28 28
29 29 ProcessingUnit.__init__(self)
30 30
31 31 self.buffer = None
32 32 self.firstdatatime = None
33 33 self.profIndex = 0
34 34 self.dataOut = Spectra()
35 35 self.dataOut.error=False
36 36 self.id_min = None
37 37 self.id_max = None
38 38 self.setupReq = False #Agregar a todas las unidades de proc
39 39 self.nsamplesFFT = 0
40 40
41 41 def __updateSpecFromVoltage(self):
42 42
43 43 self.dataOut.timeZone = self.dataIn.timeZone
44 44 self.dataOut.dstFlag = self.dataIn.dstFlag
45 45 self.dataOut.errorCount = self.dataIn.errorCount
46 46 self.dataOut.useLocalTime = self.dataIn.useLocalTime
47 47 try:
48 48 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
49 49 except:
50 50 pass
51 51 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
52 52 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
53 53 self.dataOut.ippSeconds = self.dataIn.ippSeconds
54 54 self.dataOut.ipp = self.dataIn.ipp
55 55 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
56 56 self.dataOut.channelList = self.dataIn.channelList
57 57 self.dataOut.heightList = self.dataIn.heightList
58 58 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
59 59 self.dataOut.nProfiles = self.dataOut.nFFTPoints
60 60 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
61 61 self.dataOut.utctime = self.firstdatatime
62 62 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
63 63 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
64 64 self.dataOut.flagShiftFFT = False
65 65 self.dataOut.nCohInt = self.dataIn.nCohInt
66 66 self.dataOut.nIncohInt = 1
67 67 self.dataOut.deltaHeight = self.dataIn.deltaHeight
68 68 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
69 69 self.dataOut.frequency = self.dataIn.frequency
70 70 self.dataOut.realtime = self.dataIn.realtime
71 71 self.dataOut.azimuth = self.dataIn.azimuth
72 72 self.dataOut.zenith = self.dataIn.zenith
73 73 self.dataOut.codeList = self.dataIn.codeList
74 74 self.dataOut.azimuthList = self.dataIn.azimuthList
75 75 self.dataOut.elevationList = self.dataIn.elevationList
76 76 self.dataOut.code = self.dataIn.code
77 77 self.dataOut.nCode = self.dataIn.nCode
78 78 self.dataOut.flagProfilesByRange = self.dataIn.flagProfilesByRange
79 79 self.dataOut.nProfilesByRange = self.dataIn.nProfilesByRange
80 80 self.dataOut.runNextUnit = self.dataIn.runNextUnit
81 81 try:
82 82 self.dataOut.step = self.dataIn.step
83 83 except:
84 84 pass
85 85
86 86 def __getFft(self):
87 87 """
88 88 Convierte valores de Voltaje a Spectra
89 89
90 90 Affected:
91 91 self.dataOut.data_spc
92 92 self.dataOut.data_cspc
93 93 self.dataOut.data_dc
94 94 self.dataOut.heightList
95 95 self.profIndex
96 96 self.buffer
97 97 self.dataOut.flagNoData
98 98 """
99 99 fft_volt = numpy.fft.fft(
100 100 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
101 101 fft_volt = fft_volt.astype(numpy.dtype('complex'))
102 102 dc = fft_volt[:, 0, :]
103 103
104 104 # calculo de self-spectra
105 105 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
106 106 spc = fft_volt * numpy.conjugate(fft_volt)
107 107 spc = spc.real
108 108
109 109 blocksize = 0
110 110 blocksize += dc.size
111 111 blocksize += spc.size
112 112
113 113 cspc = None
114 114 pairIndex = 0
115 115 if self.dataOut.pairsList != None:
116 116 # calculo de cross-spectra
117 117 cspc = numpy.zeros(
118 118 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
119 119 for pair in self.dataOut.pairsList:
120 120 if pair[0] not in self.dataOut.channelList:
121 121 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
122 122 str(pair), str(self.dataOut.channelList)))
123 123 if pair[1] not in self.dataOut.channelList:
124 124 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
125 125 str(pair), str(self.dataOut.channelList)))
126 126
127 127 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
128 128 numpy.conjugate(fft_volt[pair[1], :, :])
129 129 pairIndex += 1
130 130 blocksize += cspc.size
131 131
132 132 self.dataOut.data_spc = spc
133 133 self.dataOut.data_cspc = cspc
134 134 self.dataOut.data_dc = dc
135 135 self.dataOut.blockSize = blocksize
136 136 self.dataOut.flagShiftFFT = False
137 137
138 138 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False,
139 139 zeroPad=False, zeroPoints=0, runNextUnit=0):
140
140 141 self.dataIn.runNextUnit = runNextUnit
141 142 try:
142 type = self.dataIn.type.decode("utf-8")
143 self.dataIn.type = type
143 _type = self.dataIn.type.decode("utf-8")
144 self.dataIn.type = _type
144 145 except Exception as e:
145 # print("spc -> ",e)
146 #print("spc -> ",self.dataIn.type, e)
146 147 pass
147 148
148 149 if self.dataIn.type == "Spectra":
149 150 #print("AQUI")
150 151 try:
151 152 self.dataOut.copy(self.dataIn)
152 153 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
153 154 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
154 155 self.dataOut.nProfiles = self.dataOut.nFFTPoints
155 156 #self.dataOut.nHeights = len(self.dataOut.heightList)
156 157 except Exception as e:
157 158 print("Error dataIn ",e)
158 159
159 160 if shift_fft:
160 161 #desplaza a la derecha en el eje 2 determinadas posiciones
161 162 shift = int(self.dataOut.nFFTPoints/2)
162 163 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
163 164
164 165 if self.dataOut.data_cspc is not None:
165 166 #desplaza a la derecha en el eje 2 determinadas posiciones
166 167 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
167 168 if pairsList:
168 169 self.__selectPairs(pairsList)
169 170
170 171 elif self.dataIn.type == "Voltage":
171 172
172 173 self.dataOut.flagNoData = True
173 174 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
174 175 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
175 176
176 177 if nFFTPoints == None:
177 178 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
178 179
179 180 if nProfiles == None:
180 181 nProfiles = nFFTPoints
181 182
182 183 if ippFactor == None:
183 184 self.dataOut.ippFactor = self.dataIn.ippFactor
184 185 else:
185 186 self.dataOut.ippFactor = ippFactor
186 187
187 188 if self.buffer is None:
188 189 if not zeroPad:
189 190 self.buffer = numpy.zeros((self.dataIn.nChannels,
190 191 nProfiles,
191 192 self.dataIn.nHeights),
192 193 dtype='complex')
193 194 zeroPoints = 0
194 195 else:
195 196 self.buffer = numpy.zeros((self.dataIn.nChannels,
196 197 nFFTPoints+int(zeroPoints),
197 198 self.dataIn.nHeights),
198 199 dtype='complex')
199 200
200 201 self.dataOut.nFFTPoints = nFFTPoints
201 202
202 203 if self.buffer is None:
203 204 self.buffer = numpy.zeros((self.dataIn.nChannels,
204 205 nProfiles,
205 206 self.dataIn.nHeights),
206 207 dtype='complex')
207 208
208 209 if self.dataIn.flagDataAsBlock:
209 210 nVoltProfiles = self.dataIn.data.shape[1]
210 211 zeroPoints = 0
211 212 if nVoltProfiles == nProfiles or zeroPad:
212 213 self.buffer = self.dataIn.data.copy()
213 214 self.profIndex = nVoltProfiles
214 215
215 216 elif nVoltProfiles < nProfiles:
216 217
217 218 if self.profIndex == 0:
218 219 self.id_min = 0
219 220 self.id_max = nVoltProfiles
220 221
221 222 self.buffer[:, self.id_min:self.id_max,
222 223 :] = self.dataIn.data
223 224 self.profIndex += nVoltProfiles
224 225 self.id_min += nVoltProfiles
225 226 self.id_max += nVoltProfiles
226 227 elif nVoltProfiles > nProfiles:
227 228 self.reader.bypass = True
228 229 if self.profIndex == 0:
229 230 self.id_min = 0
230 231 self.id_max = nProfiles
231 232
232 233 self.buffer = self.dataIn.data[:, self.id_min:self.id_max,:]
233 234 self.profIndex += nProfiles
234 235 self.id_min += nProfiles
235 236 self.id_max += nProfiles
236 237 if self.id_max == nVoltProfiles:
237 238 self.reader.bypass = False
238 239
239 240 else:
240 241 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
241 242 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
242 243 self.dataOut.flagNoData = True
243 244 else:
244 245 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
245 246 self.profIndex += 1
246 247
247 248 if self.firstdatatime == None:
248 249 self.firstdatatime = self.dataIn.utctime
249 250
250 251 if self.profIndex == nProfiles or (zeroPad and zeroPoints==0):
251 252
252 253 self.__updateSpecFromVoltage()
253 254 if pairsList == None:
254 255 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
255 256 else:
256 257 self.dataOut.pairsList = pairsList
257 258 self.__getFft()
258 259 self.dataOut.flagNoData = False
259 260 self.firstdatatime = None
260 261 self.nsamplesFFT = self.profIndex
261 262 #if not self.reader.bypass:
262 263 self.profIndex = 0
263 264 #update Processing Header:
264 265 self.dataOut.processingHeaderObj.dtype = "Spectra"
265 266 self.dataOut.processingHeaderObj.nFFTPoints = self.dataOut.nFFTPoints
266 267 self.dataOut.processingHeaderObj.nSamplesFFT = self.nsamplesFFT
267 268 self.dataOut.processingHeaderObj.nIncohInt = 1
268 269
269 270 elif self.dataIn.type == "Parameters": #when get data from h5 spc file
270 271
271 272 self.dataOut.data_spc = self.dataIn.data_spc
272 273 self.dataOut.data_cspc = self.dataIn.data_cspc
273 274 self.dataOut.data_outlier = self.dataIn.data_outlier
274 275 self.dataOut.nProfiles = self.dataIn.nProfiles
275 276 self.dataOut.nIncohInt = self.dataIn.nIncohInt
276 277 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
277 278 self.dataOut.ippFactor = self.dataIn.ippFactor
278 279 self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
279 280 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
280 281 self.dataOut.ProcessingHeader = self.dataIn.ProcessingHeader.copy()
281 282 self.dataOut.ippSeconds = self.dataIn.ippSeconds
282 283 self.dataOut.ipp = self.dataIn.ipp
283 284 #self.dataOut.abscissaList = self.dataIn.getVelRange(1)
284 285 #self.dataOut.spc_noise = self.dataIn.getNoise()
285 286 #self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
286 287 # self.dataOut.normFactor = self.dataIn.normFactor
287 288 if hasattr(self.dataIn, 'channelList'):
288 289 self.dataOut.channelList = self.dataIn.channelList
289 290 if hasattr(self.dataIn, 'pairsList'):
290 291 self.dataOut.pairsList = self.dataIn.pairsList
291 292 self.dataOut.groupList = self.dataIn.pairsList
292 293
293 294 self.dataOut.flagNoData = False
294 295
295 296 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
296 297 self.dataOut.ChanDist = self.dataIn.ChanDist
297 298 else: self.dataOut.ChanDist = None
298 299
299 300 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
300 301 # self.dataOut.VelRange = self.dataIn.VelRange
301 302 #else: self.dataOut.VelRange = None
302 303
303 304 else:
304 305 raise ValueError("The type of input object '%s' is not valid".format(
305 306 self.dataIn.type))
306 307 # print("SPC done")
307 308
308 309 def __selectPairs(self, pairsList):
309 310
310 311 if not pairsList:
311 312 return
312 313
313 314 pairs = []
314 315 pairsIndex = []
315 316
316 317 for pair in pairsList:
317 318 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
318 319 continue
319 320 pairs.append(pair)
320 321 pairsIndex.append(pairs.index(pair))
321 322
322 323 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
323 324 self.dataOut.pairsList = pairs
324 325
325 326 return
326 327
327 328 def selectFFTs(self, minFFT, maxFFT ):
328 329 """
329 330 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
330 331 minFFT<= FFT <= maxFFT
331 332 """
332 333
333 334 if (minFFT > maxFFT):
334 335 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
335 336
336 337 if (minFFT < self.dataOut.getFreqRange()[0]):
337 338 minFFT = self.dataOut.getFreqRange()[0]
338 339
339 340 if (maxFFT > self.dataOut.getFreqRange()[-1]):
340 341 maxFFT = self.dataOut.getFreqRange()[-1]
341 342
342 343 minIndex = 0
343 344 maxIndex = 0
344 345 FFTs = self.dataOut.getFreqRange()
345 346
346 347 inda = numpy.where(FFTs >= minFFT)
347 348 indb = numpy.where(FFTs <= maxFFT)
348 349
349 350 try:
350 351 minIndex = inda[0][0]
351 352 except:
352 353 minIndex = 0
353 354
354 355 try:
355 356 maxIndex = indb[0][-1]
356 357 except:
357 358 maxIndex = len(FFTs)
358 359
359 360 self.selectFFTsByIndex(minIndex, maxIndex)
360 361
361 362 return 1
362 363
363 364 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
364 365 newheis = numpy.where(
365 366 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
366 367
367 368 if hei_ref != None:
368 369 newheis = numpy.where(self.dataOut.heightList > hei_ref)
369 370
370 371 minIndex = min(newheis[0])
371 372 maxIndex = max(newheis[0])
372 373 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
373 374 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
374 375
375 376 # determina indices
376 377 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
377 378 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
378 379 avg_dB = 10 * \
379 380 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
380 381 beacon_dB = numpy.sort(avg_dB)[-nheis:]
381 382 beacon_heiIndexList = []
382 383 for val in avg_dB.tolist():
383 384 if val >= beacon_dB[0]:
384 385 beacon_heiIndexList.append(avg_dB.tolist().index(val))
385 386
386 387 data_cspc = None
387 388 if self.dataOut.data_cspc is not None:
388 389 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
389 390
390 391 data_dc = None
391 392 if self.dataOut.data_dc is not None:
392 393 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
393 394
394 395 self.dataOut.data_spc = data_spc
395 396 self.dataOut.data_cspc = data_cspc
396 397 self.dataOut.data_dc = data_dc
397 398 self.dataOut.heightList = heightList
398 399 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
399 400
400 401 return 1
401 402
402 403 def selectFFTsByIndex(self, minIndex, maxIndex):
403 404 """
404 405
405 406 """
406 407
407 408 if (minIndex < 0) or (minIndex > maxIndex):
408 409 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
409 410
410 411 if (maxIndex >= self.dataOut.nProfiles):
411 412 maxIndex = self.dataOut.nProfiles-1
412 413
413 414 #Spectra
414 415 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
415 416
416 417 data_cspc = None
417 418 if self.dataOut.data_cspc is not None:
418 419 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
419 420
420 421 data_dc = None
421 422 if self.dataOut.data_dc is not None:
422 423 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
423 424
424 425 self.dataOut.data_spc = data_spc
425 426 self.dataOut.data_cspc = data_cspc
426 427 self.dataOut.data_dc = data_dc
427 428
428 429 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
429 430 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
430 431 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
431 432
432 433 return 1
433 434
434 435 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
435 436 # validacion de rango
436 437 if minHei == None:
437 438 minHei = self.dataOut.heightList[0]
438 439
439 440 if maxHei == None:
440 441 maxHei = self.dataOut.heightList[-1]
441 442
442 443 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
443 444 print('minHei: %.2f is out of the heights range' % (minHei))
444 445 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
445 446 minHei = self.dataOut.heightList[0]
446 447
447 448 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
448 449 print('maxHei: %.2f is out of the heights range' % (maxHei))
449 450 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
450 451 maxHei = self.dataOut.heightList[-1]
451 452
452 453 # validacion de velocidades
453 454 velrange = self.dataOut.getVelRange(1)
454 455
455 456 if minVel == None:
456 457 minVel = velrange[0]
457 458
458 459 if maxVel == None:
459 460 maxVel = velrange[-1]
460 461
461 462 if (minVel < velrange[0]) or (minVel > maxVel):
462 463 print('minVel: %.2f is out of the velocity range' % (minVel))
463 464 print('minVel is setting to %.2f' % (velrange[0]))
464 465 minVel = velrange[0]
465 466
466 467 if (maxVel > velrange[-1]) or (maxVel < minVel):
467 468 print('maxVel: %.2f is out of the velocity range' % (maxVel))
468 469 print('maxVel is setting to %.2f' % (velrange[-1]))
469 470 maxVel = velrange[-1]
470 471
471 472 # seleccion de indices para rango
472 473 minIndex = 0
473 474 maxIndex = 0
474 475 heights = self.dataOut.heightList
475 476
476 477 inda = numpy.where(heights >= minHei)
477 478 indb = numpy.where(heights <= maxHei)
478 479
479 480 try:
480 481 minIndex = inda[0][0]
481 482 except:
482 483 minIndex = 0
483 484
484 485 try:
485 486 maxIndex = indb[0][-1]
486 487 except:
487 488 maxIndex = len(heights)
488 489
489 490 if (minIndex < 0) or (minIndex > maxIndex):
490 491 raise ValueError("some value in (%d,%d) is not valid" % (
491 492 minIndex, maxIndex))
492 493
493 494 if (maxIndex >= self.dataOut.nHeights):
494 495 maxIndex = self.dataOut.nHeights - 1
495 496
496 497 # seleccion de indices para velocidades
497 498 indminvel = numpy.where(velrange >= minVel)
498 499 indmaxvel = numpy.where(velrange <= maxVel)
499 500 try:
500 501 minIndexVel = indminvel[0][0]
501 502 except:
502 503 minIndexVel = 0
503 504
504 505 try:
505 506 maxIndexVel = indmaxvel[0][-1]
506 507 except:
507 508 maxIndexVel = len(velrange)
508 509
509 510 # seleccion del espectro
510 511 data_spc = self.dataOut.data_spc[:,
511 512 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
512 513 # estimacion de ruido
513 514 noise = numpy.zeros(self.dataOut.nChannels)
514 515
515 516 for channel in range(self.dataOut.nChannels):
516 517 daux = data_spc[channel, :, :]
517 518 sortdata = numpy.sort(daux, axis=None)
518 519 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
519 520
520 521 self.dataOut.noise_estimation = noise.copy()
521 522
522 523 return 1
523 524
524 525 class GetSNR(Operation):
525 526 '''
526 527 Written by R. Flores
527 528 '''
528 529 """Operation to get SNR.
529 530
530 531 Parameters:
531 532 -----------
532 533
533 534 Example
534 535 --------
535 536
536 537 op = proc_unit.addOperation(name='GetSNR', optype='other')
537 538
538 539 """
539 540
540 541 def __init__(self, **kwargs):
541 542
542 543 Operation.__init__(self, **kwargs)
543 544
544 545 def run(self,dataOut):
545 546
546 547 noise = dataOut.getNoise(ymin_index=-10) #Región superior donde solo debería de haber ruido
547 548 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
548 549 dataOut.snl = numpy.log10(dataOut.data_snr)
549 550 dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
550 551
551 552 return dataOut
552 553
553 554 class removeDC(Operation):
554 555
555 556 def run(self, dataOut, mode=2):
556 557 self.dataOut = dataOut
557 558 jspectra = self.dataOut.data_spc
558 559 jcspectra = self.dataOut.data_cspc
559 560
560 561 num_chan = jspectra.shape[0]
561 562 num_hei = jspectra.shape[2]
562 563
563 564 if jcspectra is not None:
564 565 jcspectraExist = True
565 566 num_pairs = jcspectra.shape[0]
566 567 else:
567 568 jcspectraExist = False
568 569
569 570 freq_dc = int(jspectra.shape[1] / 2)
570 571 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
571 572 ind_vel = ind_vel.astype(int)
572 573
573 574 if ind_vel[0] < 0:
574 575 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
575 576
576 577 if mode == 1:
577 578 jspectra[:, freq_dc, :] = (
578 579 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
579 580
580 581 if jcspectraExist:
581 582 jcspectra[:, freq_dc, :] = (
582 583 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
583 584
584 585 if mode == 2:
585 586
586 587 vel = numpy.array([-2, -1, 1, 2])
587 588 xx = numpy.zeros([4, 4])
588 589
589 590 for fil in range(4):
590 591 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
591 592
592 593 xx_inv = numpy.linalg.inv(xx)
593 594 xx_aux = xx_inv[0, :]
594 595
595 596 for ich in range(num_chan):
596 597 yy = jspectra[ich, ind_vel, :]
597 598 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
598 599
599 600 junkid = jspectra[ich, freq_dc, :] <= 0
600 601 cjunkid = sum(junkid)
601 602
602 603 if cjunkid.any():
603 604 jspectra[ich, freq_dc, junkid.nonzero()] = (
604 605 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
605 606
606 607 if jcspectraExist:
607 608 for ip in range(num_pairs):
608 609 yy = jcspectra[ip, ind_vel, :]
609 610 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
610 611
611 612 self.dataOut.data_spc = jspectra
612 613 self.dataOut.data_cspc = jcspectra
613 614
614 615 return self.dataOut
615 616 class getNoiseB(Operation):
616 617 """
617 618 Get noise from custom heights and frequency ranges,
618 619 offset for additional manual correction
619 620 J. Apaza -> developed to amisr isr spectra
620 621
621 622 """
622 623 __slots__ =('offset','warnings', 'isConfig', 'minIndex','maxIndex','minIndexFFT','maxIndexFFT')
623 624 def __init__(self):
624 625
625 626 Operation.__init__(self)
626 627 self.isConfig = False
627 628
628 629 def setup(self, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
629 630
630 631 self.warnings = warnings
631 632 if minHei == None:
632 633 minHei = self.dataOut.heightList[0]
633 634
634 635 if maxHei == None:
635 636 maxHei = self.dataOut.heightList[-1]
636 637
637 638 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
638 639 if self.warnings:
639 640 print('minHei: %.2f is out of the heights range' % (minHei))
640 641 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
641 642 minHei = self.dataOut.heightList[0]
642 643
643 644 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
644 645 if self.warnings:
645 646 print('maxHei: %.2f is out of the heights range' % (maxHei))
646 647 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
647 648 maxHei = self.dataOut.heightList[-1]
648 649
649 650
650 651 #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
651 652 minIndexFFT = 0
652 653 maxIndexFFT = 0
653 654 # validacion de velocidades
654 655 indminPoint = None
655 656 indmaxPoint = None
656 657 if self.dataOut.type == 'Spectra':
657 658 if minVel == None and maxVel == None :
658 659
659 660 freqrange = self.dataOut.getFreqRange(1)
660 661
661 662 if minFreq == None:
662 663 minFreq = freqrange[0]
663 664
664 665 if maxFreq == None:
665 666 maxFreq = freqrange[-1]
666 667
667 668 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
668 669 if self.warnings:
669 670 print('minFreq: %.2f is out of the frequency range' % (minFreq))
670 671 print('minFreq is setting to %.2f' % (freqrange[0]))
671 672 minFreq = freqrange[0]
672 673
673 674 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
674 675 if self.warnings:
675 676 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
676 677 print('maxFreq is setting to %.2f' % (freqrange[-1]))
677 678 maxFreq = freqrange[-1]
678 679
679 680 indminPoint = numpy.where(freqrange >= minFreq)
680 681 indmaxPoint = numpy.where(freqrange <= maxFreq)
681 682
682 683 else:
683 684
684 685 velrange = self.dataOut.getVelRange(1)
685 686
686 687 if minVel == None:
687 688 minVel = velrange[0]
688 689
689 690 if maxVel == None:
690 691 maxVel = velrange[-1]
691 692
692 693 if (minVel < velrange[0]) or (minVel > maxVel):
693 694 if self.warnings:
694 695 print('minVel: %.2f is out of the velocity range' % (minVel))
695 696 print('minVel is setting to %.2f' % (velrange[0]))
696 697 minVel = velrange[0]
697 698
698 699 if (maxVel > velrange[-1]) or (maxVel < minVel):
699 700 if self.warnings:
700 701 print('maxVel: %.2f is out of the velocity range' % (maxVel))
701 702 print('maxVel is setting to %.2f' % (velrange[-1]))
702 703 maxVel = velrange[-1]
703 704
704 705 indminPoint = numpy.where(velrange >= minVel)
705 706 indmaxPoint = numpy.where(velrange <= maxVel)
706 707
707 708
708 709 # seleccion de indices para rango REEMPLAZAR FOR FUNCION EXTERNA LUEGO
709 710 # minIndex = 0
710 711 # maxIndex = 0
711 712 # heights = self.dataOut.heightList
712 713 # inda = numpy.where(heights >= minHei)
713 714 # indb = numpy.where(heights <= maxHei)
714 715 # try:
715 716 # minIndex = inda[0][0]
716 717 # except:
717 718 # minIndex = 0
718 719 # try:
719 720 # maxIndex = indb[0][-1]
720 721 # except:
721 722 # maxIndex = len(heights)
722 723 # if (minIndex < 0) or (minIndex > maxIndex):
723 724 # raise ValueError("some value in (%d,%d) is not valid" % (
724 725 # minIndex, maxIndex))
725 726 # if (maxIndex >= self.dataOut.nHeights):
726 727 # maxIndex = self.dataOut.nHeights - 1
727 728
728 729 minIndex, maxIndex = getHei_index(minHei,maxHei,self.dataOut.heightList)
729 730
730 731
731 732 #############################################################3
732 733 # seleccion de indices para velocidades
733 734 if self.dataOut.type == 'Spectra':
734 735 try:
735 736 minIndexFFT = indminPoint[0][0]
736 737 except:
737 738 minIndexFFT = 0
738 739
739 740 try:
740 741 maxIndexFFT = indmaxPoint[0][-1]
741 742 except:
742 743 maxIndexFFT = len( self.dataOut.getFreqRange(1))
743 744
744 745 self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT
745 746 self.isConfig = True
746 747 self.offset = 1
747 748 if offset!=None:
748 749 self.offset = 10**(offset/10)
749 750
750 751
751 752 def run(self, dataOut, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
752 753 self.dataOut = dataOut
753 754
754 755 if not self.isConfig:
755 756 self.setup(offset, minHei, maxHei,minVel, maxVel, minFreq, maxFreq, warnings)
756 757
757 758 self.dataOut.noise_estimation = None
758 759 noise = None
759 760 if self.dataOut.type == 'Voltage':
760 761 noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
761 762 elif self.dataOut.type == 'Spectra':
762 763 noise = numpy.zeros( self.dataOut.nChannels)
763 764 norm = 1
764 765
765 766 for channel in range( self.dataOut.nChannels):
766 767 if not hasattr(self.dataOut.nIncohInt,'__len__'):
767 768 norm = 1
768 769 else:
769 770 norm = self.dataOut.max_nIncohInt[channel]/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
770 771
771 772 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
772 773 daux = numpy.multiply(daux, norm)
773 774 sortdata = numpy.sort(daux, axis=None)
774 775 noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt[channel])/self.offset
775 776
776 777 else:
777 778 noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
778 779
779 780 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
780 781
781 782 return self.dataOut
782 783
783 784 def getNoiseByMean(self,data):
784 785 #data debe estar ordenado
785 786 data = numpy.mean(data,axis=1)
786 787 sortdata = numpy.sort(data, axis=None)
787 788 pnoise = None
788 789 j = 0
789 790
790 791 mean = numpy.mean(sortdata)
791 792 min = numpy.min(sortdata)
792 793 delta = mean - min
793 794 indexes = numpy.where(sortdata > (mean+delta))[0] #only array of indexes
794 795 #print(len(indexes))
795 796 if len(indexes)==0:
796 797 pnoise = numpy.mean(sortdata)
797 798 else:
798 799 j = indexes[0]
799 800 pnoise = numpy.mean(sortdata[0:j])
800 801
801 802 return pnoise
802 803
803 804 def getNoiseByHS(self,data, navg):
804 805 #data debe estar ordenado
805 806 #data = numpy.mean(data,axis=1)
806 807 sortdata = numpy.sort(data, axis=None)
807 808
808 809 lenOfData = len(sortdata)
809 810 nums_min = lenOfData*0.2
810 811
811 812 if nums_min <= 5:
812 813
813 814 nums_min = 5
814 815
815 816 sump = 0.
816 817 sumq = 0.
817 818
818 819 j = 0
819 820 cont = 1
820 821
821 822 while((cont == 1)and(j < lenOfData)):
822 823
823 824 sump += sortdata[j]
824 825 sumq += sortdata[j]**2
825 826 #sumq -= sump**2
826 827 if j > nums_min:
827 828 rtest = float(j)/(j-1) + 1.0/navg
828 829 #if ((sumq*j) > (sump**2)):
829 830 if ((sumq*j) > (rtest*sump**2)):
830 831 j = j - 1
831 832 sump = sump - sortdata[j]
832 833 sumq = sumq - sortdata[j]**2
833 834 cont = 0
834 835
835 836 j += 1
836 837
837 838 lnoise = sump / j
838 839
839 840 return lnoise
840 841
841 842 class removeInterference(Operation):
842 843
843 844 def removeInterference2(self):
844 845
845 846 cspc = self.dataOut.data_cspc
846 847 spc = self.dataOut.data_spc
847 848 Heights = numpy.arange(cspc.shape[2])
848 849 realCspc = numpy.abs(cspc)
849 850
850 851 for i in range(cspc.shape[0]):
851 852 LinePower= numpy.sum(realCspc[i], axis=0)
852 853 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
853 854 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
854 855 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
855 856 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
856 857 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
857 858
858 859
859 860 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
860 861 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
861 862 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
862 863 cspc[i,InterferenceRange,:] = numpy.NaN
863 864
864 865 self.dataOut.data_cspc = cspc
865 866
866 867 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
867 868
868 869 jspectra = self.dataOut.data_spc
869 870 jcspectra = self.dataOut.data_cspc
870 871 jnoise = self.dataOut.getNoise()
871 872 num_incoh = self.dataOut.nIncohInt
872 873
873 874 num_channel = jspectra.shape[0]
874 875 num_prof = jspectra.shape[1]
875 876 num_hei = jspectra.shape[2]
876 877
877 878 # hei_interf
878 879 if hei_interf is None:
879 880 count_hei = int(num_hei / 2)
880 881 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
881 882 hei_interf = numpy.asarray(hei_interf)[0]
882 883 # nhei_interf
883 884 if (nhei_interf == None):
884 885 nhei_interf = 5
885 886 if (nhei_interf < 1):
886 887 nhei_interf = 1
887 888 if (nhei_interf > count_hei):
888 889 nhei_interf = count_hei
889 890 if (offhei_interf == None):
890 891 offhei_interf = 0
891 892
892 893 ind_hei = list(range(num_hei))
893 894 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
894 895 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
895 896 mask_prof = numpy.asarray(list(range(num_prof)))
896 897 num_mask_prof = mask_prof.size
897 898 comp_mask_prof = [0, num_prof / 2]
898 899
899 900 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
900 901 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
901 902 jnoise = numpy.nan
902 903 noise_exist = jnoise[0] < numpy.Inf
903 904
904 905 # Subrutina de Remocion de la Interferencia
905 906 for ich in range(num_channel):
906 907 # Se ordena los espectros segun su potencia (menor a mayor)
907 908 power = jspectra[ich, mask_prof, :]
908 909 power = power[:, hei_interf]
909 910 power = power.sum(axis=0)
910 911 psort = power.ravel().argsort()
911 912
912 913 # Se estima la interferencia promedio en los Espectros de Potencia empleando
913 914 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
914 915 offhei_interf, nhei_interf + offhei_interf))]]]
915 916
916 917 if noise_exist:
917 918 # tmp_noise = jnoise[ich] / num_prof
918 919 tmp_noise = jnoise[ich]
919 920 junkspc_interf = junkspc_interf - tmp_noise
920 921 #junkspc_interf[:,comp_mask_prof] = 0
921 922
922 923 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
923 924 jspc_interf = jspc_interf.transpose()
924 925 # Calculando el espectro de interferencia promedio
925 926 noiseid = numpy.where(
926 927 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
927 928 noiseid = noiseid[0]
928 929 cnoiseid = noiseid.size
929 930 interfid = numpy.where(
930 931 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
931 932 interfid = interfid[0]
932 933 cinterfid = interfid.size
933 934
934 935 if (cnoiseid > 0):
935 936 jspc_interf[noiseid] = 0
936 937
937 938 # Expandiendo los perfiles a limpiar
938 939 if (cinterfid > 0):
939 940 new_interfid = (
940 941 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
941 942 new_interfid = numpy.asarray(new_interfid)
942 943 new_interfid = {x for x in new_interfid}
943 944 new_interfid = numpy.array(list(new_interfid))
944 945 new_cinterfid = new_interfid.size
945 946 else:
946 947 new_cinterfid = 0
947 948
948 949 for ip in range(new_cinterfid):
949 950 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
950 951 jspc_interf[new_interfid[ip]
951 952 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
952 953
953 954 jspectra[ich, :, ind_hei] = jspectra[ich, :,
954 955 ind_hei] - jspc_interf # Corregir indices
955 956
956 957 # Removiendo la interferencia del punto de mayor interferencia
957 958 ListAux = jspc_interf[mask_prof].tolist()
958 959 maxid = ListAux.index(max(ListAux))
959 960
960 961 if cinterfid > 0:
961 962 for ip in range(cinterfid * (interf == 2) - 1):
962 963 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
963 964 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
964 965 cind = len(ind)
965 966
966 967 if (cind > 0):
967 968 jspectra[ich, interfid[ip], ind] = tmp_noise * \
968 969 (1 + (numpy.random.uniform(cind) - 0.5) /
969 970 numpy.sqrt(num_incoh))
970 971
971 972 ind = numpy.array([-2, -1, 1, 2])
972 973 xx = numpy.zeros([4, 4])
973 974
974 975 for id1 in range(4):
975 976 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
976 977
977 978 xx_inv = numpy.linalg.inv(xx)
978 979 xx = xx_inv[:, 0]
979 980 ind = (ind + maxid + num_mask_prof) % num_mask_prof
980 981 yy = jspectra[ich, mask_prof[ind], :]
981 982 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
982 983 yy.transpose(), xx)
983 984
984 985 indAux = (jspectra[ich, :, :] < tmp_noise *
985 986 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
986 987 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
987 988 (1 - 1 / numpy.sqrt(num_incoh))
988 989
989 990 # Remocion de Interferencia en el Cross Spectra
990 991 if jcspectra is None:
991 992 return jspectra, jcspectra
992 993 num_pairs = int(jcspectra.size / (num_prof * num_hei))
993 994 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
994 995
995 996 for ip in range(num_pairs):
996 997
997 998 #-------------------------------------------
998 999
999 1000 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1000 1001 cspower = cspower[:, hei_interf]
1001 1002 cspower = cspower.sum(axis=0)
1002 1003
1003 1004 cspsort = cspower.ravel().argsort()
1004 1005 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1005 1006 offhei_interf, nhei_interf + offhei_interf))]]]
1006 1007 junkcspc_interf = junkcspc_interf.transpose()
1007 1008 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1008 1009
1009 1010 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1010 1011
1011 1012 median_real = int(numpy.median(numpy.real(
1012 1013 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1013 1014 median_imag = int(numpy.median(numpy.imag(
1014 1015 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1015 1016 comp_mask_prof = [int(e) for e in comp_mask_prof]
1016 1017 junkcspc_interf[comp_mask_prof, :] = numpy.complex_(
1017 1018 median_real, median_imag)
1018 1019
1019 1020 for iprof in range(num_prof):
1020 1021 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1021 1022 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1022 1023
1023 1024 # Removiendo la Interferencia
1024 1025 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1025 1026 :, ind_hei] - jcspc_interf
1026 1027
1027 1028 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1028 1029 maxid = ListAux.index(max(ListAux))
1029 1030
1030 1031 ind = numpy.array([-2, -1, 1, 2])
1031 1032 xx = numpy.zeros([4, 4])
1032 1033
1033 1034 for id1 in range(4):
1034 1035 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1035 1036
1036 1037 xx_inv = numpy.linalg.inv(xx)
1037 1038 xx = xx_inv[:, 0]
1038 1039
1039 1040 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1040 1041 yy = jcspectra[ip, mask_prof[ind], :]
1041 1042 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1042 1043
1043 1044 # Guardar Resultados
1044 1045 self.dataOut.data_spc = jspectra
1045 1046 self.dataOut.data_cspc = jcspectra
1046 1047
1047 1048 return 1
1048 1049
1049 1050
1050 1051 def run(self, dataOut, interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None, mode=1):
1051 1052
1052 1053 self.dataOut = dataOut
1053 1054
1054 1055 if mode == 1:
1055 1056 self.removeInterference(interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None)
1056 1057 elif mode == 2:
1057 1058 self.removeInterference2()
1058 1059
1059 1060 return self.dataOut
1060 1061
1061 1062
1062 1063 class deflip(Operation):
1063 1064
1064 1065 def run(self, dataOut):
1065 1066 # arreglo 1: (num_chan, num_profiles, num_heights)
1066 1067 self.dataOut = dataOut
1067 1068
1068 1069 # JULIA-oblicua, indice 2
1069 1070 # arreglo 2: (num_profiles, num_heights)
1070 1071 jspectra = self.dataOut.data_spc[2]
1071 1072 jspectra_tmp=numpy.zeros(jspectra.shape)
1072 1073 num_profiles=jspectra.shape[0]
1073 1074 freq_dc = int(num_profiles / 2)
1074 1075 # Flip con for
1075 1076 for j in range(num_profiles):
1076 1077 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1077 1078 # Intercambio perfil de DC con perfil inmediato anterior
1078 1079 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1079 1080 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1080 1081 # canal modificado es re-escrito en el arreglo de canales
1081 1082 self.dataOut.data_spc[2] = jspectra_tmp
1082 1083
1083 1084 return self.dataOut
1084 1085
1085 1086
1086 1087 class IncohInt(Operation):
1087 1088
1088 1089 __profIndex = 0
1089 1090 __withOverapping = False
1090 1091
1091 1092 __byTime = False
1092 1093 __initime = None
1093 1094 __lastdatatime = None
1094 1095 __integrationtime = None
1095 1096
1096 1097 __buffer_spc = None
1097 1098 __buffer_cspc = None
1098 1099 __buffer_dc = None
1099 1100
1100 1101 __dataReady = False
1101 1102
1102 1103 __timeInterval = None
1103 1104 incohInt = 0
1104 1105 nOutliers = 0
1105 1106 n = None
1106 1107
1107 1108 _flagProfilesByRange = False
1108 1109 _nProfilesByRange = 0
1109 1110 def __init__(self):
1110 1111
1111 1112 Operation.__init__(self)
1112 1113
1113 1114 def setup(self, n=None, timeInterval=None, overlapping=False):
1114 1115 """
1115 1116 Set the parameters of the integration class.
1116 1117
1117 1118 Inputs:
1118 1119
1119 1120 n : Number of coherent integrations
1120 1121 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1121 1122 overlapping :
1122 1123
1123 1124 """
1124 1125
1125 1126 self.__initime = None
1126 1127 self.__lastdatatime = 0
1127 1128
1128 1129 self.__buffer_spc = 0
1129 1130 self.__buffer_cspc = 0
1130 1131 self.__buffer_dc = 0
1131 1132
1132 1133 self.__profIndex = 0
1133 1134 self.__dataReady = False
1134 1135 self.__byTime = False
1135 1136 self.incohInt = 0
1136 1137 self.nOutliers = 0
1137 1138 if n is None and timeInterval is None:
1138 1139 raise ValueError("n or timeInterval should be specified ...")
1139 1140
1140 1141 if n is not None:
1141 1142 self.n = int(n)
1142 1143 else:
1143 1144
1144 1145 self.__integrationtime = int(timeInterval)
1145 1146 self.n = None
1146 1147 self.__byTime = True
1147 1148
1148 1149
1149 1150
1150 1151 def putData(self, data_spc, data_cspc, data_dc):
1151 1152 """
1152 1153 Add a profile to the __buffer_spc and increase in one the __profileIndex
1153 1154
1154 1155 """
1155 1156 if data_spc.all() == numpy.nan :
1156 1157 print("nan ")
1157 1158 return
1158 1159 self.__buffer_spc += data_spc
1159 1160
1160 1161 if data_cspc is None:
1161 1162 self.__buffer_cspc = None
1162 1163 else:
1163 1164 self.__buffer_cspc += data_cspc
1164 1165
1165 1166 if data_dc is None:
1166 1167 self.__buffer_dc = None
1167 1168 else:
1168 1169 self.__buffer_dc += data_dc
1169 1170
1170 1171 self.__profIndex += 1
1171 1172
1172 1173 return
1173 1174
1174 1175 def pushData(self):
1175 1176 """
1176 1177 Return the sum of the last profiles and the profiles used in the sum.
1177 1178
1178 1179 Affected:
1179 1180
1180 1181 self.__profileIndex
1181 1182
1182 1183 """
1183 1184
1184 1185 data_spc = self.__buffer_spc
1185 1186 data_cspc = self.__buffer_cspc
1186 1187 data_dc = self.__buffer_dc
1187 1188 n = self.__profIndex
1188 1189
1189 1190 self.__buffer_spc = 0
1190 1191 self.__buffer_cspc = 0
1191 1192 self.__buffer_dc = 0
1192 1193
1193 1194
1194 1195 return data_spc, data_cspc, data_dc, n
1195 1196
1196 1197 def byProfiles(self, *args):
1197 1198
1198 1199 self.__dataReady = False
1199 1200 avgdata_spc = None
1200 1201 avgdata_cspc = None
1201 1202 avgdata_dc = None
1202 1203
1203 1204 self.putData(*args)
1204 1205
1205 1206 if self.__profIndex == self.n:
1206 1207
1207 1208 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1208 1209 self.n = n
1209 1210 self.__dataReady = True
1210 1211
1211 1212 return avgdata_spc, avgdata_cspc, avgdata_dc
1212 1213
1213 1214 def byTime(self, datatime, *args):
1214 1215
1215 1216 self.__dataReady = False
1216 1217 avgdata_spc = None
1217 1218 avgdata_cspc = None
1218 1219 avgdata_dc = None
1219 1220
1220 1221 self.putData(*args)
1221 1222
1222 1223 if (datatime - self.__initime) >= self.__integrationtime:
1223 1224 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1224 1225 self.n = n
1225 1226 self.__dataReady = True
1226 1227
1227 1228 return avgdata_spc, avgdata_cspc, avgdata_dc
1228 1229
1229 1230 def integrate(self, datatime, *args):
1230 1231
1231 1232 if self.__profIndex == 0:
1232 1233 self.__initime = datatime
1233 1234
1234 1235 if self.__byTime:
1235 1236 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1236 1237 datatime, *args)
1237 1238 else:
1238 1239 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1239 1240
1240 1241 if not self.__dataReady:
1241 1242 return None, None, None, None
1242 1243
1243 1244 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1244 1245
1245 1246 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1246 1247 if n == 1:
1247 1248 return dataOut
1248 1249
1249 1250 if dataOut.flagNoData == True:
1250 1251 return dataOut
1251 1252
1252 1253 if dataOut.flagProfilesByRange == True:
1253 1254 self._flagProfilesByRange = True
1254 1255
1255 1256 dataOut.flagNoData = True
1256 1257 dataOut.processingHeaderObj.timeIncohInt = timeInterval
1257 1258 if not self.isConfig:
1258 1259 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
1259 1260 self.setup(n, timeInterval, overlapping)
1260 1261 self.isConfig = True
1261 1262
1262 1263
1263 1264 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1264 1265 dataOut.data_spc,
1265 1266 dataOut.data_cspc,
1266 1267 dataOut.data_dc)
1267 1268
1268 1269 self.incohInt += dataOut.nIncohInt
1269 1270
1270 1271
1271 1272 if isinstance(dataOut.data_outlier,numpy.ndarray) or isinstance(dataOut.data_outlier,int) or isinstance(dataOut.data_outlier, float):
1272 1273 self.nOutliers += dataOut.data_outlier
1273 1274
1274 1275 if self._flagProfilesByRange:
1275 1276 dataOut.flagProfilesByRange = True
1276 1277 self._nProfilesByRange += dataOut.nProfilesByRange
1277 1278
1278 1279 if self.__dataReady:
1279 1280 #print("prof: ",dataOut.max_nIncohInt,self.__profIndex)
1280 1281 dataOut.data_spc = avgdata_spc
1281 1282 dataOut.data_cspc = avgdata_cspc
1282 1283 dataOut.data_dc = avgdata_dc
1283 1284 dataOut.nIncohInt = self.incohInt
1284 1285 dataOut.data_outlier = self.nOutliers
1285 1286 dataOut.utctime = avgdatatime
1286 1287 dataOut.flagNoData = False
1287 1288 self.incohInt = 0
1288 1289 self.nOutliers = 0
1289 1290 self.__profIndex = 0
1290 1291 dataOut.nProfilesByRange = self._nProfilesByRange
1291 1292 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
1292 1293 self._flagProfilesByRange = False
1293 1294 # print("IncohInt Done")
1294 1295 return dataOut
1295 1296
1296 1297
1297 1298 class IntegrationFaradaySpectra(Operation):
1298 1299
1299 1300 __profIndex = 0
1300 1301 __withOverapping = False
1301 1302
1302 1303 __byTime = False
1303 1304 __initime = None
1304 1305 __lastdatatime = None
1305 1306 __integrationtime = None
1306 1307
1307 1308 __buffer_spc = None
1308 1309 __buffer_cspc = None
1309 1310 __buffer_dc = None
1310 1311
1311 1312 __dataReady = False
1312 1313
1313 1314 __timeInterval = None
1314 1315 n_ints = None #matriz de numero de integracions (CH,HEI)
1315 1316 n = None
1316 1317 minHei_ind = None
1317 1318 maxHei_ind = None
1318 1319 navg = 1.0
1319 1320 factor = 0.0
1320 1321 dataoutliers = None # (CHANNELS, HEIGHTS)
1321 1322
1322 1323 _flagProfilesByRange = False
1323 1324 _nProfilesByRange = 0
1324 1325
1325 1326 def __init__(self):
1326 1327
1327 1328 Operation.__init__(self)
1328 1329
1329 1330 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, avg=1,factor=0.75):
1330 1331 """
1331 1332 Set the parameters of the integration class.
1332 1333
1333 1334 Inputs:
1334 1335
1335 1336 n : Number of coherent integrations
1336 1337 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1337 1338 overlapping :
1338 1339
1339 1340 """
1340 1341
1341 1342 self.__initime = None
1342 1343 self.__lastdatatime = 0
1343 1344
1344 1345 self.__buffer_spc = []
1345 1346 self.__buffer_cspc = []
1346 1347 self.__buffer_dc = 0
1347 1348
1348 1349 self.__profIndex = 0
1349 1350 self.__dataReady = False
1350 1351 self.__byTime = False
1351 1352
1352 1353 self.factor = factor
1353 1354 self.navg = avg
1354 1355 #self.ByLags = dataOut.ByLags ###REDEFINIR
1355 1356 self.ByLags = False
1356 1357 self.maxProfilesInt = 0
1357 1358 self.__nChannels = dataOut.nChannels
1358 1359 if DPL != None:
1359 1360 self.DPL=DPL
1360 1361 else:
1361 1362 #self.DPL=dataOut.DPL ###REDEFINIR
1362 1363 self.DPL=0
1363 1364
1364 1365 if n is None and timeInterval is None:
1365 1366 raise ValueError("n or timeInterval should be specified ...")
1366 1367
1367 1368 if n is not None:
1368 1369 self.n = int(n)
1369 1370 else:
1370 1371 self.__integrationtime = int(timeInterval)
1371 1372 self.n = None
1372 1373 self.__byTime = True
1373 1374
1374 1375
1375 1376 if minHei == None:
1376 1377 minHei = self.dataOut.heightList[0]
1377 1378
1378 1379 if maxHei == None:
1379 1380 maxHei = self.dataOut.heightList[-1]
1380 1381
1381 1382 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1382 1383 print('minHei: %.2f is out of the heights range' % (minHei))
1383 1384 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
1384 1385 minHei = self.dataOut.heightList[0]
1385 1386
1386 1387 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
1387 1388 print('maxHei: %.2f is out of the heights range' % (maxHei))
1388 1389 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
1389 1390 maxHei = self.dataOut.heightList[-1]
1390 1391
1391 1392 ind_list1 = numpy.where(self.dataOut.heightList >= minHei)
1392 1393 ind_list2 = numpy.where(self.dataOut.heightList <= maxHei)
1393 1394 self.minHei_ind = ind_list1[0][0]
1394 1395 self.maxHei_ind = ind_list2[0][-1]
1395 1396
1396 1397 def putData(self, data_spc, data_cspc, data_dc):
1397 1398 """
1398 1399 Add a profile to the __buffer_spc and increase in one the __profileIndex
1399 1400
1400 1401 """
1401 1402
1402 1403 self.__buffer_spc.append(data_spc)
1403 1404
1404 1405 if self.__nChannels < 2:
1405 1406 self.__buffer_cspc = None
1406 1407 else:
1407 1408 self.__buffer_cspc.append(data_cspc)
1408 1409
1409 1410 if data_dc is None:
1410 1411 self.__buffer_dc = None
1411 1412 else:
1412 1413 self.__buffer_dc += data_dc
1413 1414
1414 1415 self.__profIndex += 1
1415 1416
1416 1417 return
1417 1418
1418 1419 def hildebrand_sekhon_Integration(self,sortdata,navg, factor):
1419 1420 #data debe estar ordenado
1420 1421 #sortdata = numpy.sort(data, axis=None)
1421 1422 #sortID=data.argsort()
1422 1423 lenOfData = len(sortdata)
1423 1424 nums_min = lenOfData*factor
1424 1425 if nums_min <= 5:
1425 1426 nums_min = 5
1426 1427 sump = 0.
1427 1428 sumq = 0.
1428 1429 j = 0
1429 1430 cont = 1
1430 1431 while((cont == 1)and(j < lenOfData)):
1431 1432 sump += sortdata[j]
1432 1433 sumq += sortdata[j]**2
1433 1434 if j > nums_min:
1434 1435 rtest = float(j)/(j-1) + 1.0/navg
1435 1436 if ((sumq*j) > (rtest*sump**2)):
1436 1437 j = j - 1
1437 1438 sump = sump - sortdata[j]
1438 1439 sumq = sumq - sortdata[j]**2
1439 1440 cont = 0
1440 1441 j += 1
1441 1442 #lnoise = sump / j
1442 1443 #print("H S done")
1443 1444 #return j,sortID
1444 1445 return j
1445 1446
1446 1447
1447 1448 def pushData(self):
1448 1449 """
1449 1450 Return the sum of the last profiles and the profiles used in the sum.
1450 1451
1451 1452 Affected:
1452 1453
1453 1454 self.__profileIndex
1454 1455
1455 1456 """
1456 1457 bufferH=None
1457 1458 buffer=None
1458 1459 buffer1=None
1459 1460 buffer_cspc=None
1460 1461 #print("aes: ", self.__buffer_cspc)
1461 1462 self.__buffer_spc=numpy.array(self.__buffer_spc)
1462 1463 if self.__nChannels > 1 :
1463 1464 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1464 1465
1465 1466 #print("FREQ_DC",self.__buffer_spc.shape,self.__buffer_cspc.shape)
1466 1467
1467 1468 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1468 1469 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1469 1470
1470 1471 self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
1471 1472
1472 1473 for k in range(self.minHei_ind,self.maxHei_ind):
1473 1474 if self.__nChannels > 1:
1474 1475 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1475 1476
1476 1477 outliers_IDs_cspc=[]
1477 1478 cspc_outliers_exist=False
1478 1479 for i in range(self.nChannels):#dataOut.nChannels):
1479 1480
1480 1481 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1481 1482 indexes=[]
1482 1483 #sortIDs=[]
1483 1484 outliers_IDs=[]
1484 1485
1485 1486 for j in range(self.nProfiles): #frecuencias en el tiempo
1486 1487 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1487 1488 # continue
1488 1489 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1489 1490 # continue
1490 1491 buffer=buffer1[:,j]
1491 1492 sortdata = numpy.sort(buffer, axis=None)
1492 1493
1493 1494 sortID=buffer.argsort()
1494 1495 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1495 1496
1496 1497 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1497 1498
1498 1499 # fig,ax = plt.subplots()
1499 1500 # ax.set_title(str(k)+" "+str(j))
1500 1501 # x=range(len(sortdata))
1501 1502 # ax.scatter(x,sortdata)
1502 1503 # ax.axvline(index)
1503 1504 # plt.show()
1504 1505
1505 1506 indexes.append(index)
1506 1507 #sortIDs.append(sortID)
1507 1508 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1508 1509
1509 1510 #print("Outliers: ",outliers_IDs)
1510 1511 outliers_IDs=numpy.array(outliers_IDs)
1511 1512 outliers_IDs=outliers_IDs.ravel()
1512 1513 outliers_IDs=numpy.unique(outliers_IDs)
1513 1514 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1514 1515 indexes=numpy.array(indexes)
1515 1516 indexmin=numpy.min(indexes)
1516 1517
1517 1518
1518 1519 #print(indexmin,buffer1.shape[0], k)
1519 1520
1520 1521 # fig,ax = plt.subplots()
1521 1522 # ax.plot(sortdata)
1522 1523 # ax2 = ax.twinx()
1523 1524 # x=range(len(indexes))
1524 1525 # #plt.scatter(x,indexes)
1525 1526 # ax2.scatter(x,indexes)
1526 1527 # plt.show()
1527 1528
1528 1529 if indexmin != buffer1.shape[0]:
1529 1530 if self.__nChannels > 1:
1530 1531 cspc_outliers_exist= True
1531 1532
1532 1533 lt=outliers_IDs
1533 1534 #avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1534 1535
1535 1536 for p in list(outliers_IDs):
1536 1537 #buffer1[p,:]=avg
1537 1538 buffer1[p,:] = numpy.NaN
1538 1539
1539 1540 self.dataOutliers[i,k] = len(outliers_IDs)
1540 1541
1541 1542
1542 1543 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1543 1544
1544 1545
1545 1546 if self.__nChannels > 1:
1546 1547 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1547 1548
1548 1549
1549 1550 if self.__nChannels > 1:
1550 1551 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1551 1552 if cspc_outliers_exist:
1552 1553
1553 1554 lt=outliers_IDs_cspc
1554 1555
1555 1556 #avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1556 1557 for p in list(outliers_IDs_cspc):
1557 1558 #buffer_cspc[p,:]=avg
1558 1559 buffer_cspc[p,:] = numpy.NaN
1559 1560
1560 1561 if self.__nChannels > 1:
1561 1562 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1562 1563
1563 1564
1564 1565
1565 1566
1566 1567 nOutliers = len(outliers_IDs)
1567 1568 #print("Outliers n: ",self.dataOutliers,nOutliers)
1568 1569 buffer=None
1569 1570 bufferH=None
1570 1571 buffer1=None
1571 1572 buffer_cspc=None
1572 1573
1573 1574
1574 1575 buffer=None
1575 1576
1576 1577 #data_spc = numpy.sum(self.__buffer_spc,axis=0)
1577 1578 data_spc = numpy.nansum(self.__buffer_spc,axis=0)
1578 1579 if self.__nChannels > 1:
1579 1580 #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1580 1581 data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
1581 1582 else:
1582 1583 data_cspc = None
1583 1584 data_dc = self.__buffer_dc
1584 1585 #(CH, HEIGH)
1585 1586 self.maxProfilesInt = self.__profIndex - 1
1586 1587 n = self.__profIndex - self.dataOutliers # n becomes a matrix
1587 1588
1588 1589 self.__buffer_spc = []
1589 1590 self.__buffer_cspc = []
1590 1591 self.__buffer_dc = 0
1591 1592 self.__profIndex = 0
1592 1593 #print("cleaned ",data_cspc)
1593 1594 return data_spc, data_cspc, data_dc, n
1594 1595
1595 1596 def byProfiles(self, *args):
1596 1597
1597 1598 self.__dataReady = False
1598 1599 avgdata_spc = None
1599 1600 avgdata_cspc = None
1600 1601 avgdata_dc = None
1601 1602
1602 1603 self.putData(*args)
1603 1604
1604 1605 if self.__profIndex >= self.n:
1605 1606
1606 1607 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1607 1608 self.n_ints = n
1608 1609 self.__dataReady = True
1609 1610
1610 1611 return avgdata_spc, avgdata_cspc, avgdata_dc
1611 1612
1612 1613 def byTime(self, datatime, *args):
1613 1614
1614 1615 self.__dataReady = False
1615 1616 avgdata_spc = None
1616 1617 avgdata_cspc = None
1617 1618 avgdata_dc = None
1618 1619
1619 1620 self.putData(*args)
1620 1621
1621 1622 if (datatime - self.__initime) >= self.__integrationtime:
1622 1623 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1623 1624 self.n_ints = n
1624 1625 self.__dataReady = True
1625 1626
1626 1627 return avgdata_spc, avgdata_cspc, avgdata_dc
1627 1628
1628 1629 def integrate(self, datatime, *args):
1629 1630
1630 1631 if self.__profIndex == 0:
1631 1632 self.__initime = datatime
1632 1633
1633 1634 if self.__byTime:
1634 1635 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1635 1636 datatime, *args)
1636 1637 else:
1637 1638 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1638 1639
1639 1640 if not self.__dataReady:
1640 1641 return None, None, None, None
1641 1642
1642 1643 #print("integrate", avgdata_cspc)
1643 1644 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1644 1645
1645 1646 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
1646 1647 self.dataOut = dataOut
1647 1648 if n == 1:
1648 1649 return self.dataOut
1649 1650 self.dataOut.processingHeaderObj.timeIncohInt = timeInterval
1650 1651
1651 1652 if dataOut.flagProfilesByRange:
1652 1653 self._flagProfilesByRange = True
1653 1654
1654 1655 if self.dataOut.nChannels == 1:
1655 1656 self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
1656 1657 #print("IN spc:", self.dataOut.data_spc.shape, self.dataOut.data_cspc)
1657 1658 if not self.isConfig:
1658 1659 self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor)
1659 1660 self.isConfig = True
1660 1661
1661 1662 if not self.ByLags:
1662 1663 self.nProfiles=self.dataOut.nProfiles
1663 1664 self.nChannels=self.dataOut.nChannels
1664 1665 self.nHeights=self.dataOut.nHeights
1665 1666 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1666 1667 self.dataOut.data_spc,
1667 1668 self.dataOut.data_cspc,
1668 1669 self.dataOut.data_dc)
1669 1670 else:
1670 1671 self.nProfiles=self.dataOut.nProfiles
1671 1672 self.nChannels=self.dataOut.nChannels
1672 1673 self.nHeights=self.dataOut.nHeights
1673 1674 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1674 1675 self.dataOut.dataLag_spc,
1675 1676 self.dataOut.dataLag_cspc,
1676 1677 self.dataOut.dataLag_dc)
1677 1678 self.dataOut.flagNoData = True
1678 1679
1679 1680 if self._flagProfilesByRange:
1680 1681 dataOut.flagProfilesByRange = True
1681 1682 self._nProfilesByRange += dataOut.nProfilesByRange
1682 1683
1683 1684 if self.__dataReady:
1684 1685
1685 1686 if not self.ByLags:
1686 1687 if self.nChannels == 1:
1687 1688 #print("f int", avgdata_spc.shape)
1688 1689 self.dataOut.data_spc = avgdata_spc
1689 1690 self.dataOut.data_cspc = None
1690 1691 else:
1691 1692 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
1692 1693 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1693 1694 self.dataOut.data_dc = avgdata_dc
1694 1695 self.dataOut.data_outlier = self.dataOutliers
1695 1696
1696 1697
1697 1698 else:
1698 1699 self.dataOut.dataLag_spc = avgdata_spc
1699 1700 self.dataOut.dataLag_cspc = avgdata_cspc
1700 1701 self.dataOut.dataLag_dc = avgdata_dc
1701 1702
1702 1703 self.dataOut.data_spc=self.dataOut.dataLag_spc[:,:,:,self.dataOut.LagPlot]
1703 1704 self.dataOut.data_cspc=self.dataOut.dataLag_cspc[:,:,:,self.dataOut.LagPlot]
1704 1705 self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot]
1705 1706
1706 1707 self.dataOut.nIncohInt *= self.n_ints
1707 1708
1708 1709 self.dataOut.utctime = avgdatatime
1709 1710 self.dataOut.flagNoData = False
1710 1711
1711 1712 dataOut.nProfilesByRange = self._nProfilesByRange
1712 1713 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
1713 1714 self._flagProfilesByRange = False
1714 1715
1715 1716 return self.dataOut
1716 1717
1717 1718 class dopplerFlip(Operation):
1718 1719
1719 1720 def run(self, dataOut, chann = None):
1720 1721 # arreglo 1: (num_chan, num_profiles, num_heights)
1721 1722 self.dataOut = dataOut
1722 1723 # JULIA-oblicua, indice 2
1723 1724 # arreglo 2: (num_profiles, num_heights)
1724 1725 jspectra = self.dataOut.data_spc[chann]
1725 1726 jspectra_tmp = numpy.zeros(jspectra.shape)
1726 1727 num_profiles = jspectra.shape[0]
1727 1728 freq_dc = int(num_profiles / 2)
1728 1729 # Flip con for
1729 1730 for j in range(num_profiles):
1730 1731 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1731 1732 # Intercambio perfil de DC con perfil inmediato anterior
1732 1733 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1733 1734 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1734 1735 # canal modificado es re-escrito en el arreglo de canales
1735 1736 self.dataOut.data_spc[chann] = jspectra_tmp
1736 1737
1737 1738 return self.dataOut No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now