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