##// END OF EJS Templates
schain
RSS
Bug in jrodata ncode nbaud
jespinoza -
r1340:4f5903c4ea97
parent child
Show More
Side by Side Unified
Context file:
r1340:4f5903c4ea97 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1193 +1,1193
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Definition of diferent Data objects for different types of data
6 6
7 7 Here you will find the diferent data objects for the different types
8 8 of data, this data objects must be used as dataIn or dataOut objects in
9 9 processing units and operations. Currently the supported data objects are:
10 10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 11 """
12 12
13 13 import copy
14 14 import numpy
15 15 import datetime
16 16 import json
17 17
18 18 import schainpy.admin
19 19 from schainpy.utils import log
20 20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 21 from schainpy.model.data import _noise
22 22
23 23
24 24 def getNumpyDtype(dataTypeCode):
25 25
26 26 if dataTypeCode == 0:
27 27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 28 elif dataTypeCode == 1:
29 29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 30 elif dataTypeCode == 2:
31 31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 32 elif dataTypeCode == 3:
33 33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 34 elif dataTypeCode == 4:
35 35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 36 elif dataTypeCode == 5:
37 37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 38 else:
39 39 raise ValueError('dataTypeCode was not defined')
40 40
41 41 return numpyDtype
42 42
43 43
44 44 def getDataTypeCode(numpyDtype):
45 45
46 46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 47 datatype = 0
48 48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 49 datatype = 1
50 50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 51 datatype = 2
52 52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 53 datatype = 3
54 54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 55 datatype = 4
56 56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 57 datatype = 5
58 58 else:
59 59 datatype = None
60 60
61 61 return datatype
62 62
63 63
64 64 def hildebrand_sekhon(data, navg):
65 65 """
66 66 This method is for the objective determination of the noise level in Doppler spectra. This
67 67 implementation technique is based on the fact that the standard deviation of the spectral
68 68 densities is equal to the mean spectral density for white Gaussian noise
69 69
70 70 Inputs:
71 71 Data : heights
72 72 navg : numbers of averages
73 73
74 74 Return:
75 75 mean : noise's level
76 76 """
77 77
78 78 sortdata = numpy.sort(data, axis=None)
79 79 '''
80 80 lenOfData = len(sortdata)
81 81 nums_min = lenOfData*0.2
82 82
83 83 if nums_min <= 5:
84 84
85 85 nums_min = 5
86 86
87 87 sump = 0.
88 88 sumq = 0.
89 89
90 90 j = 0
91 91 cont = 1
92 92
93 93 while((cont == 1)and(j < lenOfData)):
94 94
95 95 sump += sortdata[j]
96 96 sumq += sortdata[j]**2
97 97
98 98 if j > nums_min:
99 99 rtest = float(j)/(j-1) + 1.0/navg
100 100 if ((sumq*j) > (rtest*sump**2)):
101 101 j = j - 1
102 102 sump = sump - sortdata[j]
103 103 sumq = sumq - sortdata[j]**2
104 104 cont = 0
105 105
106 106 j += 1
107 107
108 108 lnoise = sump / j
109 109 '''
110 110 return _noise.hildebrand_sekhon(sortdata, navg)
111 111
112 112
113 113 class Beam:
114 114
115 115 def __init__(self):
116 116 self.codeList = []
117 117 self.azimuthList = []
118 118 self.zenithList = []
119 119
120 120
121 121 class GenericData(object):
122 122
123 123 flagNoData = True
124 124
125 125 def copy(self, inputObj=None):
126 126
127 127 if inputObj == None:
128 128 return copy.deepcopy(self)
129 129
130 130 for key in list(inputObj.__dict__.keys()):
131 131
132 132 attribute = inputObj.__dict__[key]
133 133
134 134 # If this attribute is a tuple or list
135 135 if type(inputObj.__dict__[key]) in (tuple, list):
136 136 self.__dict__[key] = attribute[:]
137 137 continue
138 138
139 139 # If this attribute is another object or instance
140 140 if hasattr(attribute, '__dict__'):
141 141 self.__dict__[key] = attribute.copy()
142 142 continue
143 143
144 144 self.__dict__[key] = inputObj.__dict__[key]
145 145
146 146 def deepcopy(self):
147 147
148 148 return copy.deepcopy(self)
149 149
150 150 def isEmpty(self):
151 151
152 152 return self.flagNoData
153 153
154 154 def isReady(self):
155 155
156 156 return not self.flagNoData
157 157
158 158
159 159 class JROData(GenericData):
160 160
161 161 systemHeaderObj = SystemHeader()
162 162 radarControllerHeaderObj = RadarControllerHeader()
163 163 type = None
164 164 datatype = None # dtype but in string
165 165 nProfiles = None
166 166 heightList = None
167 167 channelList = None
168 168 flagDiscontinuousBlock = False
169 169 useLocalTime = False
170 170 utctime = None
171 171 timeZone = None
172 172 dstFlag = None
173 173 errorCount = None
174 174 blocksize = None
175 175 flagDecodeData = False # asumo q la data no esta decodificada
176 176 flagDeflipData = False # asumo q la data no esta sin flip
177 177 flagShiftFFT = False
178 178 nCohInt = None
179 179 windowOfFilter = 1
180 180 C = 3e8
181 181 frequency = 49.92e6
182 182 realtime = False
183 183 beacon_heiIndexList = None
184 184 last_block = None
185 185 blocknow = None
186 186 azimuth = None
187 187 zenith = None
188 188 beam = Beam()
189 189 profileIndex = None
190 190 error = None
191 191 data = None
192 192 nmodes = None
193 193 metadata_list = ['heightList', 'timeZone', 'type']
194 194
195 195 def __str__(self):
196 196
197 197 return '{} - {}'.format(self.type, self.datatime())
198 198
199 199 def getNoise(self):
200 200
201 201 raise NotImplementedError
202 202
203 203 @property
204 204 def nChannels(self):
205 205
206 206 return len(self.channelList)
207 207
208 208 @property
209 209 def channelIndexList(self):
210 210
211 211 return list(range(self.nChannels))
212 212
213 213 @property
214 214 def nHeights(self):
215 215
216 216 return len(self.heightList)
217 217
218 218 def getDeltaH(self):
219 219
220 220 return self.heightList[1] - self.heightList[0]
221 221
222 222 @property
223 223 def ltctime(self):
224 224
225 225 if self.useLocalTime:
226 226 return self.utctime - self.timeZone * 60
227 227
228 228 return self.utctime
229 229
230 230 @property
231 231 def datatime(self):
232 232
233 233 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
234 234 return datatimeValue
235 235
236 236 def getTimeRange(self):
237 237
238 238 datatime = []
239 239
240 240 datatime.append(self.ltctime)
241 241 datatime.append(self.ltctime + self.timeInterval + 1)
242 242
243 243 datatime = numpy.array(datatime)
244 244
245 245 return datatime
246 246
247 247 def getFmaxTimeResponse(self):
248 248
249 249 period = (10**-6) * self.getDeltaH() / (0.15)
250 250
251 251 PRF = 1. / (period * self.nCohInt)
252 252
253 253 fmax = PRF
254 254
255 255 return fmax
256 256
257 257 def getFmax(self):
258 258 PRF = 1. / (self.ippSeconds * self.nCohInt)
259 259
260 260 fmax = PRF
261 261 return fmax
262 262
263 263 def getVmax(self):
264 264
265 265 _lambda = self.C / self.frequency
266 266
267 267 vmax = self.getFmax() * _lambda / 2
268 268
269 269 return vmax
270 270
271 271 @property
272 272 def ippSeconds(self):
273 273 '''
274 274 '''
275 275 return self.radarControllerHeaderObj.ippSeconds
276 276
277 277 @ippSeconds.setter
278 278 def ippSeconds(self, ippSeconds):
279 279 '''
280 280 '''
281 281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
282 282
283 283 @property
284 284 def code(self):
285 285 '''
286 286 '''
287 287 return self.radarControllerHeaderObj.code
288 288
289 289 @code.setter
290 290 def code(self, code):
291 291 '''
292 292 '''
293 293 self.radarControllerHeaderObj.code = code
294 294
295 295 @property
296 def ncode(self):
296 def nCode(self):
297 297 '''
298 298 '''
299 299 return self.radarControllerHeaderObj.nCode
300 300
301 @ncode.setter
302 def ncode(self, ncode):
301 @nCode.setter
302 def nCode(self, ncode):
303 303 '''
304 304 '''
305 305 self.radarControllerHeaderObj.nCode = ncode
306 306
307 307 @property
308 def nbaud(self):
308 def nBaud(self):
309 309 '''
310 310 '''
311 311 return self.radarControllerHeaderObj.nBaud
312 312
313 @nbaud.setter
314 def nbaud(self, nbaud):
313 @nBaud.setter
314 def nBaud(self, nbaud):
315 315 '''
316 316 '''
317 317 self.radarControllerHeaderObj.nBaud = nbaud
318 318
319 319 @property
320 320 def ipp(self):
321 321 '''
322 322 '''
323 323 return self.radarControllerHeaderObj.ipp
324 324
325 325 @ipp.setter
326 326 def ipp(self, ipp):
327 327 '''
328 328 '''
329 329 self.radarControllerHeaderObj.ipp = ipp
330 330
331 331 @property
332 332 def metadata(self):
333 333 '''
334 334 '''
335 335
336 336 return {attr: getattr(self, attr) for attr in self.metadata_list}
337 337
338 338
339 339 class Voltage(JROData):
340 340
341 341 dataPP_POW = None
342 342 dataPP_DOP = None
343 343 dataPP_WIDTH = None
344 344 dataPP_SNR = None
345 345
346 346 def __init__(self):
347 347 '''
348 348 Constructor
349 349 '''
350 350
351 351 self.useLocalTime = True
352 352 self.radarControllerHeaderObj = RadarControllerHeader()
353 353 self.systemHeaderObj = SystemHeader()
354 354 self.type = "Voltage"
355 355 self.data = None
356 356 self.nProfiles = None
357 357 self.heightList = None
358 358 self.channelList = None
359 359 self.flagNoData = True
360 360 self.flagDiscontinuousBlock = False
361 361 self.utctime = None
362 362 self.timeZone = 0
363 363 self.dstFlag = None
364 364 self.errorCount = None
365 365 self.nCohInt = None
366 366 self.blocksize = None
367 367 self.flagCohInt = False
368 368 self.flagDecodeData = False # asumo q la data no esta decodificada
369 369 self.flagDeflipData = False # asumo q la data no esta sin flip
370 370 self.flagShiftFFT = False
371 371 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
372 372 self.profileIndex = 0
373 373 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
374 'code', 'ncode', 'nbaud', 'ippSeconds', 'ipp']
374 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
375 375
376 376 def getNoisebyHildebrand(self, channel=None):
377 377 """
378 378 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
379 379
380 380 Return:
381 381 noiselevel
382 382 """
383 383
384 384 if channel != None:
385 385 data = self.data[channel]
386 386 nChannels = 1
387 387 else:
388 388 data = self.data
389 389 nChannels = self.nChannels
390 390
391 391 noise = numpy.zeros(nChannels)
392 392 power = data * numpy.conjugate(data)
393 393
394 394 for thisChannel in range(nChannels):
395 395 if nChannels == 1:
396 396 daux = power[:].real
397 397 else:
398 398 daux = power[thisChannel, :].real
399 399 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
400 400
401 401 return noise
402 402
403 403 def getNoise(self, type=1, channel=None):
404 404
405 405 if type == 1:
406 406 noise = self.getNoisebyHildebrand(channel)
407 407
408 408 return noise
409 409
410 410 def getPower(self, channel=None):
411 411
412 412 if channel != None:
413 413 data = self.data[channel]
414 414 else:
415 415 data = self.data
416 416
417 417 power = data * numpy.conjugate(data)
418 418 powerdB = 10 * numpy.log10(power.real)
419 419 powerdB = numpy.squeeze(powerdB)
420 420
421 421 return powerdB
422 422
423 423 @property
424 424 def timeInterval(self):
425 425
426 426 return self.ippSeconds * self.nCohInt
427 427
428 428 noise = property(getNoise, "I'm the 'nHeights' property.")
429 429
430 430
431 431 class Spectra(JROData):
432 432
433 433 def __init__(self):
434 434 '''
435 435 Constructor
436 436 '''
437 437
438 438 self.useLocalTime = True
439 439 self.radarControllerHeaderObj = RadarControllerHeader()
440 440 self.systemHeaderObj = SystemHeader()
441 441 self.type = "Spectra"
442 442 self.timeZone = 0
443 443 self.nProfiles = None
444 444 self.heightList = None
445 445 self.channelList = None
446 446 self.pairsList = None
447 447 self.flagNoData = True
448 448 self.flagDiscontinuousBlock = False
449 449 self.utctime = None
450 450 self.nCohInt = None
451 451 self.nIncohInt = None
452 452 self.blocksize = None
453 453 self.nFFTPoints = None
454 454 self.wavelength = None
455 455 self.flagDecodeData = False # asumo q la data no esta decodificada
456 456 self.flagDeflipData = False # asumo q la data no esta sin flip
457 457 self.flagShiftFFT = False
458 458 self.ippFactor = 1
459 459 self.beacon_heiIndexList = []
460 460 self.noise_estimation = None
461 461 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
462 'code', 'ncode', 'nbaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
462 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
463 463
464 464 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
465 465 """
466 466 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
467 467
468 468 Return:
469 469 noiselevel
470 470 """
471 471
472 472 noise = numpy.zeros(self.nChannels)
473 473
474 474 for channel in range(self.nChannels):
475 475 daux = self.data_spc[channel,
476 476 xmin_index:xmax_index, ymin_index:ymax_index]
477 477 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
478 478
479 479 return noise
480 480
481 481 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
482 482
483 483 if self.noise_estimation is not None:
484 484 # this was estimated by getNoise Operation defined in jroproc_spectra.py
485 485 return self.noise_estimation
486 486 else:
487 487 noise = self.getNoisebyHildebrand(
488 488 xmin_index, xmax_index, ymin_index, ymax_index)
489 489 return noise
490 490
491 491 def getFreqRangeTimeResponse(self, extrapoints=0):
492 492
493 493 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
494 494 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
495 495
496 496 return freqrange
497 497
498 498 def getAcfRange(self, extrapoints=0):
499 499
500 500 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
501 501 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
502 502
503 503 return freqrange
504 504
505 505 def getFreqRange(self, extrapoints=0):
506 506
507 507 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
508 508 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
509 509
510 510 return freqrange
511 511
512 512 def getVelRange(self, extrapoints=0):
513 513
514 514 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
515 515 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
516 516
517 517 if self.nmodes:
518 518 return velrange/self.nmodes
519 519 else:
520 520 return velrange
521 521
522 522 @property
523 523 def nPairs(self):
524 524
525 525 return len(self.pairsList)
526 526
527 527 @property
528 528 def pairsIndexList(self):
529 529
530 530 return list(range(self.nPairs))
531 531
532 532 @property
533 533 def normFactor(self):
534 534
535 535 pwcode = 1
536 536
537 537 if self.flagDecodeData:
538 538 pwcode = numpy.sum(self.code[0]**2)
539 539 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
540 540 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
541 541
542 542 return normFactor
543 543
544 544 @property
545 545 def flag_cspc(self):
546 546
547 547 if self.data_cspc is None:
548 548 return True
549 549
550 550 return False
551 551
552 552 @property
553 553 def flag_dc(self):
554 554
555 555 if self.data_dc is None:
556 556 return True
557 557
558 558 return False
559 559
560 560 @property
561 561 def timeInterval(self):
562 562
563 563 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
564 564 if self.nmodes:
565 565 return self.nmodes*timeInterval
566 566 else:
567 567 return timeInterval
568 568
569 569 def getPower(self):
570 570
571 571 factor = self.normFactor
572 572 z = self.data_spc / factor
573 573 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
574 574 avg = numpy.average(z, axis=1)
575 575
576 576 return 10 * numpy.log10(avg)
577 577
578 578 def getCoherence(self, pairsList=None, phase=False):
579 579
580 580 z = []
581 581 if pairsList is None:
582 582 pairsIndexList = self.pairsIndexList
583 583 else:
584 584 pairsIndexList = []
585 585 for pair in pairsList:
586 586 if pair not in self.pairsList:
587 587 raise ValueError("Pair %s is not in dataOut.pairsList" % (
588 588 pair))
589 589 pairsIndexList.append(self.pairsList.index(pair))
590 590 for i in range(len(pairsIndexList)):
591 591 pair = self.pairsList[pairsIndexList[i]]
592 592 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
593 593 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
594 594 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
595 595 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
596 596 if phase:
597 597 data = numpy.arctan2(avgcoherenceComplex.imag,
598 598 avgcoherenceComplex.real) * 180 / numpy.pi
599 599 else:
600 600 data = numpy.abs(avgcoherenceComplex)
601 601
602 602 z.append(data)
603 603
604 604 return numpy.array(z)
605 605
606 606 def setValue(self, value):
607 607
608 608 print("This property should not be initialized")
609 609
610 610 return
611 611
612 612 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
613 613
614 614
615 615 class SpectraHeis(Spectra):
616 616
617 617 def __init__(self):
618 618
619 619 self.radarControllerHeaderObj = RadarControllerHeader()
620 620 self.systemHeaderObj = SystemHeader()
621 621 self.type = "SpectraHeis"
622 622 self.nProfiles = None
623 623 self.heightList = None
624 624 self.channelList = None
625 625 self.flagNoData = True
626 626 self.flagDiscontinuousBlock = False
627 627 self.utctime = None
628 628 self.blocksize = None
629 629 self.profileIndex = 0
630 630 self.nCohInt = 1
631 631 self.nIncohInt = 1
632 632
633 633 @property
634 634 def normFactor(self):
635 635 pwcode = 1
636 636 if self.flagDecodeData:
637 637 pwcode = numpy.sum(self.code[0]**2)
638 638
639 639 normFactor = self.nIncohInt * self.nCohInt * pwcode
640 640
641 641 return normFactor
642 642
643 643 @property
644 644 def timeInterval(self):
645 645
646 646 return self.ippSeconds * self.nCohInt * self.nIncohInt
647 647
648 648
649 649 class Fits(JROData):
650 650
651 651 def __init__(self):
652 652
653 653 self.type = "Fits"
654 654 self.nProfiles = None
655 655 self.heightList = None
656 656 self.channelList = None
657 657 self.flagNoData = True
658 658 self.utctime = None
659 659 self.nCohInt = 1
660 660 self.nIncohInt = 1
661 661 self.useLocalTime = True
662 662 self.profileIndex = 0
663 663 self.timeZone = 0
664 664
665 665 def getTimeRange(self):
666 666
667 667 datatime = []
668 668
669 669 datatime.append(self.ltctime)
670 670 datatime.append(self.ltctime + self.timeInterval)
671 671
672 672 datatime = numpy.array(datatime)
673 673
674 674 return datatime
675 675
676 676 def getChannelIndexList(self):
677 677
678 678 return list(range(self.nChannels))
679 679
680 680 def getNoise(self, type=1):
681 681
682 682
683 683 if type == 1:
684 684 noise = self.getNoisebyHildebrand()
685 685
686 686 if type == 2:
687 687 noise = self.getNoisebySort()
688 688
689 689 if type == 3:
690 690 noise = self.getNoisebyWindow()
691 691
692 692 return noise
693 693
694 694 @property
695 695 def timeInterval(self):
696 696
697 697 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
698 698
699 699 return timeInterval
700 700
701 701 @property
702 702 def ippSeconds(self):
703 703 '''
704 704 '''
705 705 return self.ipp_sec
706 706
707 707 noise = property(getNoise, "I'm the 'nHeights' property.")
708 708
709 709
710 710 class Correlation(JROData):
711 711
712 712 def __init__(self):
713 713 '''
714 714 Constructor
715 715 '''
716 716 self.radarControllerHeaderObj = RadarControllerHeader()
717 717 self.systemHeaderObj = SystemHeader()
718 718 self.type = "Correlation"
719 719 self.data = None
720 720 self.dtype = None
721 721 self.nProfiles = None
722 722 self.heightList = None
723 723 self.channelList = None
724 724 self.flagNoData = True
725 725 self.flagDiscontinuousBlock = False
726 726 self.utctime = None
727 727 self.timeZone = 0
728 728 self.dstFlag = None
729 729 self.errorCount = None
730 730 self.blocksize = None
731 731 self.flagDecodeData = False # asumo q la data no esta decodificada
732 732 self.flagDeflipData = False # asumo q la data no esta sin flip
733 733 self.pairsList = None
734 734 self.nPoints = None
735 735
736 736 def getPairsList(self):
737 737
738 738 return self.pairsList
739 739
740 740 def getNoise(self, mode=2):
741 741
742 742 indR = numpy.where(self.lagR == 0)[0][0]
743 743 indT = numpy.where(self.lagT == 0)[0][0]
744 744
745 745 jspectra0 = self.data_corr[:, :, indR, :]
746 746 jspectra = copy.copy(jspectra0)
747 747
748 748 num_chan = jspectra.shape[0]
749 749 num_hei = jspectra.shape[2]
750 750
751 751 freq_dc = jspectra.shape[1] / 2
752 752 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
753 753
754 754 if ind_vel[0] < 0:
755 755 ind_vel[list(range(0, 1))] = ind_vel[list(
756 756 range(0, 1))] + self.num_prof
757 757
758 758 if mode == 1:
759 759 jspectra[:, freq_dc, :] = (
760 760 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
761 761
762 762 if mode == 2:
763 763
764 764 vel = numpy.array([-2, -1, 1, 2])
765 765 xx = numpy.zeros([4, 4])
766 766
767 767 for fil in range(4):
768 768 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
769 769
770 770 xx_inv = numpy.linalg.inv(xx)
771 771 xx_aux = xx_inv[0, :]
772 772
773 773 for ich in range(num_chan):
774 774 yy = jspectra[ich, ind_vel, :]
775 775 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
776 776
777 777 junkid = jspectra[ich, freq_dc, :] <= 0
778 778 cjunkid = sum(junkid)
779 779
780 780 if cjunkid.any():
781 781 jspectra[ich, freq_dc, junkid.nonzero()] = (
782 782 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
783 783
784 784 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
785 785
786 786 return noise
787 787
788 788 @property
789 789 def timeInterval(self):
790 790
791 791 return self.ippSeconds * self.nCohInt * self.nProfiles
792 792
793 793 def splitFunctions(self):
794 794
795 795 pairsList = self.pairsList
796 796 ccf_pairs = []
797 797 acf_pairs = []
798 798 ccf_ind = []
799 799 acf_ind = []
800 800 for l in range(len(pairsList)):
801 801 chan0 = pairsList[l][0]
802 802 chan1 = pairsList[l][1]
803 803
804 804 # Obteniendo pares de Autocorrelacion
805 805 if chan0 == chan1:
806 806 acf_pairs.append(chan0)
807 807 acf_ind.append(l)
808 808 else:
809 809 ccf_pairs.append(pairsList[l])
810 810 ccf_ind.append(l)
811 811
812 812 data_acf = self.data_cf[acf_ind]
813 813 data_ccf = self.data_cf[ccf_ind]
814 814
815 815 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
816 816
817 817 @property
818 818 def normFactor(self):
819 819 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
820 820 acf_pairs = numpy.array(acf_pairs)
821 821 normFactor = numpy.zeros((self.nPairs, self.nHeights))
822 822
823 823 for p in range(self.nPairs):
824 824 pair = self.pairsList[p]
825 825
826 826 ch0 = pair[0]
827 827 ch1 = pair[1]
828 828
829 829 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
830 830 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
831 831 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
832 832
833 833 return normFactor
834 834
835 835
836 836 class Parameters(Spectra):
837 837
838 838 groupList = None # List of Pairs, Groups, etc
839 839 data_param = None # Parameters obtained
840 840 data_pre = None # Data Pre Parametrization
841 841 data_SNR = None # Signal to Noise Ratio
842 842 abscissaList = None # Abscissa, can be velocities, lags or time
843 843 utctimeInit = None # Initial UTC time
844 844 paramInterval = None # Time interval to calculate Parameters in seconds
845 845 useLocalTime = True
846 846 # Fitting
847 847 data_error = None # Error of the estimation
848 848 constants = None
849 849 library = None
850 850 # Output signal
851 851 outputInterval = None # Time interval to calculate output signal in seconds
852 852 data_output = None # Out signal
853 853 nAvg = None
854 854 noise_estimation = None
855 855 GauSPC = None # Fit gaussian SPC
856 856
857 857 def __init__(self):
858 858 '''
859 859 Constructor
860 860 '''
861 861 self.radarControllerHeaderObj = RadarControllerHeader()
862 862 self.systemHeaderObj = SystemHeader()
863 863 self.type = "Parameters"
864 864 self.timeZone = 0
865 865
866 866 def getTimeRange1(self, interval):
867 867
868 868 datatime = []
869 869
870 870 if self.useLocalTime:
871 871 time1 = self.utctimeInit - self.timeZone * 60
872 872 else:
873 873 time1 = self.utctimeInit
874 874
875 875 datatime.append(time1)
876 876 datatime.append(time1 + interval)
877 877 datatime = numpy.array(datatime)
878 878
879 879 return datatime
880 880
881 881 @property
882 882 def timeInterval(self):
883 883
884 884 if hasattr(self, 'timeInterval1'):
885 885 return self.timeInterval1
886 886 else:
887 887 return self.paramInterval
888 888
889 889 def setValue(self, value):
890 890
891 891 print("This property should not be initialized")
892 892
893 893 return
894 894
895 895 def getNoise(self):
896 896
897 897 return self.spc_noise
898 898
899 899 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
900 900
901 901
902 902 class PlotterData(object):
903 903 '''
904 904 Object to hold data to be plotted
905 905 '''
906 906
907 907 MAXNUMX = 200
908 908 MAXNUMY = 200
909 909
910 910 def __init__(self, code, throttle_value, exp_code, localtime=True, buffering=True, snr=False):
911 911
912 912 self.key = code
913 913 self.throttle = throttle_value
914 914 self.exp_code = exp_code
915 915 self.buffering = buffering
916 916 self.ready = False
917 917 self.flagNoData = False
918 918 self.localtime = localtime
919 919 self.data = {}
920 920 self.meta = {}
921 921 self.__heights = []
922 922
923 923 if 'snr' in code:
924 924 self.plottypes = ['snr']
925 925 elif code == 'spc':
926 926 self.plottypes = ['spc', 'noise', 'rti']
927 927 elif code == 'cspc':
928 928 self.plottypes = ['cspc', 'spc', 'noise', 'rti']
929 929 elif code == 'rti':
930 930 self.plottypes = ['noise', 'rti']
931 931 else:
932 932 self.plottypes = [code]
933 933
934 934 if 'snr' not in self.plottypes and snr:
935 935 self.plottypes.append('snr')
936 936
937 937 for plot in self.plottypes:
938 938 self.data[plot] = {}
939 939
940 940 def __str__(self):
941 941 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
942 942 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
943 943
944 944 def __len__(self):
945 945 return len(self.data[self.key])
946 946
947 947 def __getitem__(self, key):
948 948
949 949 if key not in self.data:
950 950 raise KeyError(log.error('Missing key: {}'.format(key)))
951 951 if 'spc' in key or not self.buffering:
952 952 ret = self.data[key][self.tm]
953 953 elif 'scope' in key:
954 954 ret = numpy.array(self.data[key][float(self.tm)])
955 955 else:
956 956 ret = numpy.array([self.data[key][x] for x in self.times])
957 957 if ret.ndim > 1:
958 958 ret = numpy.swapaxes(ret, 0, 1)
959 959 return ret
960 960
961 961 def __contains__(self, key):
962 962 return key in self.data
963 963
964 964 def setup(self):
965 965 '''
966 966 Configure object
967 967 '''
968 968 self.type = ''
969 969 self.ready = False
970 970 del self.data
971 971 self.data = {}
972 972 self.__heights = []
973 973 self.__all_heights = set()
974 974 for plot in self.plottypes:
975 975 if 'snr' in plot:
976 976 plot = 'snr'
977 977 elif 'spc_moments' == plot:
978 978 plot = 'moments'
979 979 self.data[plot] = {}
980 980
981 981 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
982 982 self.data['noise'] = {}
983 983 self.data['rti'] = {}
984 984 if 'noise' not in self.plottypes:
985 985 self.plottypes.append('noise')
986 986 if 'rti' not in self.plottypes:
987 987 self.plottypes.append('rti')
988 988
989 989 def shape(self, key):
990 990 '''
991 991 Get the shape of the one-element data for the given key
992 992 '''
993 993
994 994 if len(self.data[key]):
995 995 if 'spc' in key or not self.buffering:
996 996 return self.data[key].shape
997 997 return self.data[key][self.times[0]].shape
998 998 return (0,)
999 999
1000 1000 def update(self, dataOut, tm):
1001 1001 '''
1002 1002 Update data object with new dataOut
1003 1003 '''
1004 1004
1005 1005 self.profileIndex = dataOut.profileIndex
1006 1006 self.tm = tm
1007 1007 self.type = dataOut.type
1008 1008 self.parameters = getattr(dataOut, 'parameters', [])
1009 1009
1010 1010 if hasattr(dataOut, 'meta'):
1011 1011 self.meta.update(dataOut.meta)
1012 1012
1013 1013 if hasattr(dataOut, 'pairsList'):
1014 1014 self.pairs = dataOut.pairsList
1015 1015
1016 1016 self.interval = dataOut.timeInterval
1017 1017 if True in ['spc' in ptype for ptype in self.plottypes]:
1018 1018 self.xrange = (dataOut.getFreqRange(1)/1000.,
1019 1019 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1020 1020 self.__heights.append(dataOut.heightList)
1021 1021 self.__all_heights.update(dataOut.heightList)
1022 1022
1023 1023 for plot in self.plottypes:
1024 1024 if plot in ('spc', 'spc_moments', 'spc_cut'):
1025 1025 z = dataOut.data_spc/dataOut.normFactor
1026 1026 buffer = 10*numpy.log10(z)
1027 1027 if plot == 'cspc':
1028 1028 buffer = (dataOut.data_spc, dataOut.data_cspc)
1029 1029 if plot == 'noise':
1030 1030 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1031 1031 if plot in ('rti', 'spcprofile'):
1032 1032 buffer = dataOut.getPower()
1033 1033 if plot == 'snr_db':
1034 1034 buffer = dataOut.data_SNR
1035 1035 if plot == 'snr':
1036 1036 buffer = 10*numpy.log10(dataOut.data_SNR)
1037 1037 if plot == 'dop':
1038 1038 buffer = dataOut.data_DOP
1039 1039 if plot == 'pow':
1040 1040 buffer = 10*numpy.log10(dataOut.data_POW)
1041 1041 if plot == 'width':
1042 1042 buffer = dataOut.data_WIDTH
1043 1043 if plot == 'coh':
1044 1044 buffer = dataOut.getCoherence()
1045 1045 if plot == 'phase':
1046 1046 buffer = dataOut.getCoherence(phase=True)
1047 1047 if plot == 'output':
1048 1048 buffer = dataOut.data_output
1049 1049 if plot == 'param':
1050 1050 buffer = dataOut.data_param
1051 1051 if plot == 'scope':
1052 1052 buffer = dataOut.data
1053 1053 self.flagDataAsBlock = dataOut.flagDataAsBlock
1054 1054 self.nProfiles = dataOut.nProfiles
1055 1055 if plot == 'pp_power':
1056 1056 buffer = dataOut.dataPP_POWER
1057 1057 self.flagDataAsBlock = dataOut.flagDataAsBlock
1058 1058 self.nProfiles = dataOut.nProfiles
1059 1059 if plot == 'pp_signal':
1060 1060 buffer = dataOut.dataPP_POW
1061 1061 self.flagDataAsBlock = dataOut.flagDataAsBlock
1062 1062 self.nProfiles = dataOut.nProfiles
1063 1063 if plot == 'pp_velocity':
1064 1064 buffer = dataOut.dataPP_DOP
1065 1065 self.flagDataAsBlock = dataOut.flagDataAsBlock
1066 1066 self.nProfiles = dataOut.nProfiles
1067 1067 if plot == 'pp_specwidth':
1068 1068 buffer = dataOut.dataPP_WIDTH
1069 1069 self.flagDataAsBlock = dataOut.flagDataAsBlock
1070 1070 self.nProfiles = dataOut.nProfiles
1071 1071
1072 1072 if plot == 'spc':
1073 1073 self.data['spc'][tm] = buffer
1074 1074 elif plot == 'cspc':
1075 1075 self.data['cspc'][tm] = buffer
1076 1076 elif plot == 'spc_moments':
1077 1077 self.data['spc'][tm] = buffer
1078 1078 self.data['moments'][tm] = dataOut.moments
1079 1079 else:
1080 1080 if self.buffering:
1081 1081 self.data[plot][tm] = buffer
1082 1082 else:
1083 1083 self.data[plot][tm] = buffer
1084 1084
1085 1085 if dataOut.channelList is None:
1086 1086 self.channels = range(buffer.shape[0])
1087 1087 else:
1088 1088 self.channels = dataOut.channelList
1089 1089
1090 1090 if buffer is None:
1091 1091 self.flagNoData = True
1092 1092 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1093 1093
1094 1094 def normalize_heights(self):
1095 1095 '''
1096 1096 Ensure same-dimension of the data for different heighList
1097 1097 '''
1098 1098
1099 1099 H = numpy.array(list(self.__all_heights))
1100 1100 H.sort()
1101 1101 for key in self.data:
1102 1102 shape = self.shape(key)[:-1] + H.shape
1103 1103 for tm, obj in list(self.data[key].items()):
1104 1104 h = self.__heights[self.times.tolist().index(tm)]
1105 1105 if H.size == h.size:
1106 1106 continue
1107 1107 index = numpy.where(numpy.in1d(H, h))[0]
1108 1108 dummy = numpy.zeros(shape) + numpy.nan
1109 1109 if len(shape) == 2:
1110 1110 dummy[:, index] = obj
1111 1111 else:
1112 1112 dummy[index] = obj
1113 1113 self.data[key][tm] = dummy
1114 1114
1115 1115 self.__heights = [H for tm in self.times]
1116 1116
1117 1117 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1118 1118 '''
1119 1119 Convert data to json
1120 1120 '''
1121 1121
1122 1122 dy = int(self.heights.size/self.MAXNUMY) + 1
1123 1123 if self.key in ('spc', 'cspc'):
1124 1124 dx = int(self.data[self.key][tm].shape[1]/self.MAXNUMX) + 1
1125 1125 data = self.roundFloats(
1126 1126 self.data[self.key][tm][::, ::dx, ::dy].tolist())
1127 1127 else:
1128 1128 if self.key is 'noise':
1129 1129 data = [[x] for x in self.roundFloats(self.data[self.key][tm].tolist())]
1130 1130 else:
1131 1131 data = self.roundFloats(self.data[self.key][tm][::, ::dy].tolist())
1132 1132
1133 1133 meta = {}
1134 1134 ret = {
1135 1135 'plot': plot_name,
1136 1136 'code': self.exp_code,
1137 1137 'time': float(tm),
1138 1138 'data': data,
1139 1139 }
1140 1140 meta['type'] = plot_type
1141 1141 meta['interval'] = float(self.interval)
1142 1142 meta['localtime'] = self.localtime
1143 1143 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1144 1144 if 'spc' in self.data or 'cspc' in self.data:
1145 1145 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1146 1146 else:
1147 1147 meta['xrange'] = []
1148 1148
1149 1149 meta.update(self.meta)
1150 1150 ret['metadata'] = meta
1151 1151 return json.dumps(ret)
1152 1152
1153 1153 @property
1154 1154 def times(self):
1155 1155 '''
1156 1156 Return the list of times of the current data
1157 1157 '''
1158 1158
1159 1159 ret = numpy.array([t for t in self.data[self.key]])
1160 1160 if self:
1161 1161 ret.sort()
1162 1162 return ret
1163 1163
1164 1164 @property
1165 1165 def min_time(self):
1166 1166 '''
1167 1167 Return the minimun time value
1168 1168 '''
1169 1169
1170 1170 return self.times[0]
1171 1171
1172 1172 @property
1173 1173 def max_time(self):
1174 1174 '''
1175 1175 Return the maximun time value
1176 1176 '''
1177 1177
1178 1178 return self.times[-1]
1179 1179
1180 1180 @property
1181 1181 def heights(self):
1182 1182 '''
1183 1183 Return the list of heights of the current data
1184 1184 '''
1185 1185
1186 1186 return numpy.array(self.__heights[-1])
1187 1187
1188 1188 @staticmethod
1189 1189 def roundFloats(obj):
1190 1190 if isinstance(obj, list):
1191 1191 return list(map(PlotterData.roundFloats, obj))
1192 1192 elif isinstance(obj, float):
1193 1193 return round(obj, 2)
General Comments 0
You need to be logged in to leave comments. Login now