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