##// END OF EJS Templates
schain
RSS
Add h0 attribute to dataOut
Juan C. Espinoza -
r1481:f29ad632dda8
parent child
Show More
Side by Side Unified
Context file:
r1481:f29ad632dda8 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1068 +1,1069
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 h0 = 0
193 194 metadata_list = ['heightList', 'timeZone', 'type']
194 195
195 196 def __str__(self):
196 197
197 198 return '{} - {}'.format(self.type, self.datatime())
198 199
199 200 def getNoise(self):
200 201
201 202 raise NotImplementedError
202 203
203 204 @property
204 205 def nChannels(self):
205 206
206 207 return len(self.channelList)
207 208
208 209 @property
209 210 def channelIndexList(self):
210 211
211 212 return list(range(self.nChannels))
212 213
213 214 @property
214 215 def nHeights(self):
215 216
216 217 return len(self.heightList)
217 218
218 219 def getDeltaH(self):
219 220
220 221 return self.heightList[1] - self.heightList[0]
221 222
222 223 @property
223 224 def ltctime(self):
224 225
225 226 if self.useLocalTime:
226 227 return self.utctime - self.timeZone * 60
227 228
228 229 return self.utctime
229 230
230 231 @property
231 232 def datatime(self):
232 233
233 234 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
234 235 return datatimeValue
235 236
236 237 def getTimeRange(self):
237 238
238 239 datatime = []
239 240
240 241 datatime.append(self.ltctime)
241 242 datatime.append(self.ltctime + self.timeInterval + 1)
242 243
243 244 datatime = numpy.array(datatime)
244 245
245 246 return datatime
246 247
247 248 def getFmaxTimeResponse(self):
248 249
249 250 period = (10**-6) * self.getDeltaH() / (0.15)
250 251
251 252 PRF = 1. / (period * self.nCohInt)
252 253
253 254 fmax = PRF
254 255
255 256 return fmax
256 257
257 258 def getFmax(self):
258 259 PRF = 1. / (self.ippSeconds * self.nCohInt)
259 260
260 261 fmax = PRF
261 262 return fmax
262 263
263 264 def getVmax(self):
264 265
265 266 _lambda = self.C / self.frequency
266 267
267 268 vmax = self.getFmax() * _lambda / 2
268 269
269 270 return vmax
270 271
271 272 @property
272 273 def ippSeconds(self):
273 274 '''
274 275 '''
275 276 return self.radarControllerHeaderObj.ippSeconds
276 277
277 278 @ippSeconds.setter
278 279 def ippSeconds(self, ippSeconds):
279 280 '''
280 281 '''
281 282 self.radarControllerHeaderObj.ippSeconds = ippSeconds
282 283
283 284 @property
284 285 def code(self):
285 286 '''
286 287 '''
287 288 return self.radarControllerHeaderObj.code
288 289
289 290 @code.setter
290 291 def code(self, code):
291 292 '''
292 293 '''
293 294 self.radarControllerHeaderObj.code = code
294 295
295 296 @property
296 297 def nCode(self):
297 298 '''
298 299 '''
299 300 return self.radarControllerHeaderObj.nCode
300 301
301 302 @nCode.setter
302 303 def nCode(self, ncode):
303 304 '''
304 305 '''
305 306 self.radarControllerHeaderObj.nCode = ncode
306 307
307 308 @property
308 309 def nBaud(self):
309 310 '''
310 311 '''
311 312 return self.radarControllerHeaderObj.nBaud
312 313
313 314 @nBaud.setter
314 315 def nBaud(self, nbaud):
315 316 '''
316 317 '''
317 318 self.radarControllerHeaderObj.nBaud = nbaud
318 319
319 320 @property
320 321 def ipp(self):
321 322 '''
322 323 '''
323 324 return self.radarControllerHeaderObj.ipp
324 325
325 326 @ipp.setter
326 327 def ipp(self, ipp):
327 328 '''
328 329 '''
329 330 self.radarControllerHeaderObj.ipp = ipp
330 331
331 332 @property
332 333 def metadata(self):
333 334 '''
334 335 '''
335 336
336 337 return {attr: getattr(self, attr) for attr in self.metadata_list}
337 338
338 339
339 340 class Voltage(JROData):
340 341
341 342 dataPP_POW = None
342 343 dataPP_DOP = None
343 344 dataPP_WIDTH = None
344 345 dataPP_SNR = None
345 346
346 347 def __init__(self):
347 348 '''
348 349 Constructor
349 350 '''
350 351
351 352 self.useLocalTime = True
352 353 self.radarControllerHeaderObj = RadarControllerHeader()
353 354 self.systemHeaderObj = SystemHeader()
354 355 self.type = "Voltage"
355 356 self.data = None
356 357 self.nProfiles = None
357 358 self.heightList = None
358 359 self.channelList = None
359 360 self.flagNoData = True
360 361 self.flagDiscontinuousBlock = False
361 362 self.utctime = None
362 363 self.timeZone = 0
363 364 self.dstFlag = None
364 365 self.errorCount = None
365 366 self.nCohInt = None
366 367 self.blocksize = None
367 368 self.flagCohInt = False
368 369 self.flagDecodeData = False # asumo q la data no esta decodificada
369 370 self.flagDeflipData = False # asumo q la data no esta sin flip
370 371 self.flagShiftFFT = False
371 372 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
372 373 self.profileIndex = 0
373 374 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
374 375 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
375 376
376 377 def getNoisebyHildebrand(self, channel=None):
377 378 """
378 379 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
379 380
380 381 Return:
381 382 noiselevel
382 383 """
383 384
384 385 if channel != None:
385 386 data = self.data[channel]
386 387 nChannels = 1
387 388 else:
388 389 data = self.data
389 390 nChannels = self.nChannels
390 391
391 392 noise = numpy.zeros(nChannels)
392 393 power = data * numpy.conjugate(data)
393 394
394 395 for thisChannel in range(nChannels):
395 396 if nChannels == 1:
396 397 daux = power[:].real
397 398 else:
398 399 daux = power[thisChannel, :].real
399 400 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
400 401
401 402 return noise
402 403
403 404 def getNoise(self, type=1, channel=None):
404 405
405 406 if type == 1:
406 407 noise = self.getNoisebyHildebrand(channel)
407 408
408 409 return noise
409 410
410 411 def getPower(self, channel=None):
411 412
412 413 if channel != None:
413 414 data = self.data[channel]
414 415 else:
415 416 data = self.data
416 417
417 418 power = data * numpy.conjugate(data)
418 419 powerdB = 10 * numpy.log10(power.real)
419 420 powerdB = numpy.squeeze(powerdB)
420 421
421 422 return powerdB
422 423
423 424 @property
424 425 def timeInterval(self):
425 426
426 427 return self.ippSeconds * self.nCohInt
427 428
428 429 noise = property(getNoise, "I'm the 'nHeights' property.")
429 430
430 431
431 432 class Spectra(JROData):
432 433
433 434 def __init__(self):
434 435 '''
435 436 Constructor
436 437 '''
437 438
438 439 self.data_dc = None
439 440 self.data_spc = None
440 441 self.data_cspc = None
441 442 self.useLocalTime = True
442 443 self.radarControllerHeaderObj = RadarControllerHeader()
443 444 self.systemHeaderObj = SystemHeader()
444 445 self.type = "Spectra"
445 446 self.timeZone = 0
446 447 self.nProfiles = None
447 448 self.heightList = None
448 449 self.channelList = None
449 450 self.pairsList = None
450 451 self.flagNoData = True
451 452 self.flagDiscontinuousBlock = False
452 453 self.utctime = None
453 454 self.nCohInt = None
454 455 self.nIncohInt = None
455 456 self.blocksize = None
456 457 self.nFFTPoints = None
457 458 self.wavelength = None
458 459 self.flagDecodeData = False # asumo q la data no esta decodificada
459 460 self.flagDeflipData = False # asumo q la data no esta sin flip
460 461 self.flagShiftFFT = False
461 462 self.ippFactor = 1
462 463 self.beacon_heiIndexList = []
463 464 self.noise_estimation = None
464 465 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
465 466 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
466 467
467 468 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
468 469 """
469 470 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
470 471
471 472 Return:
472 473 noiselevel
473 474 """
474 475
475 476 noise = numpy.zeros(self.nChannels)
476 477
477 478 for channel in range(self.nChannels):
478 479 daux = self.data_spc[channel,
479 480 xmin_index:xmax_index, ymin_index:ymax_index]
480 481 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
481 482
482 483 return noise
483 484
484 485 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
485 486
486 487 if self.noise_estimation is not None:
487 488 # this was estimated by getNoise Operation defined in jroproc_spectra.py
488 489 return self.noise_estimation
489 490 else:
490 491 noise = self.getNoisebyHildebrand(
491 492 xmin_index, xmax_index, ymin_index, ymax_index)
492 493 return noise
493 494
494 495 def getFreqRangeTimeResponse(self, extrapoints=0):
495 496
496 497 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
497 498 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
498 499
499 500 return freqrange
500 501
501 502 def getAcfRange(self, extrapoints=0):
502 503
503 504 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
504 505 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
505 506
506 507 return freqrange
507 508
508 509 def getFreqRange(self, extrapoints=0):
509 510
510 511 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
511 512 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
512 513
513 514 return freqrange
514 515
515 516 def getVelRange(self, extrapoints=0):
516 517
517 518 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
518 519 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
519 520
520 521 if self.nmodes:
521 522 return velrange/self.nmodes
522 523 else:
523 524 return velrange
524 525
525 526 @property
526 527 def nPairs(self):
527 528
528 529 return len(self.pairsList)
529 530
530 531 @property
531 532 def pairsIndexList(self):
532 533
533 534 return list(range(self.nPairs))
534 535
535 536 @property
536 537 def normFactor(self):
537 538
538 539 pwcode = 1
539 540
540 541 if self.flagDecodeData:
541 542 pwcode = numpy.sum(self.code[0]**2)
542 543 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
543 544 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
544 545
545 546 return normFactor
546 547
547 548 @property
548 549 def flag_cspc(self):
549 550
550 551 if self.data_cspc is None:
551 552 return True
552 553
553 554 return False
554 555
555 556 @property
556 557 def flag_dc(self):
557 558
558 559 if self.data_dc is None:
559 560 return True
560 561
561 562 return False
562 563
563 564 @property
564 565 def timeInterval(self):
565 566
566 567 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
567 568 if self.nmodes:
568 569 return self.nmodes*timeInterval
569 570 else:
570 571 return timeInterval
571 572
572 573 def getPower(self):
573 574
574 575 factor = self.normFactor
575 576 z = self.data_spc / factor
576 577 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
577 578 avg = numpy.average(z, axis=1)
578 579 return 10 * numpy.log10(avg)
579 580
580 581 def getCoherence(self, pairsList=None, phase=False):
581 582
582 583 z = []
583 584 if pairsList is None:
584 585 pairsIndexList = self.pairsIndexList
585 586 else:
586 587 pairsIndexList = []
587 588 for pair in pairsList:
588 589 if pair not in self.pairsList:
589 590 raise ValueError("Pair %s is not in dataOut.pairsList" % (
590 591 pair))
591 592 pairsIndexList.append(self.pairsList.index(pair))
592 593 for i in range(len(pairsIndexList)):
593 594 pair = self.pairsList[pairsIndexList[i]]
594 595 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
595 596 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
596 597 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
597 598 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
598 599 if phase:
599 600 data = numpy.arctan2(avgcoherenceComplex.imag,
600 601 avgcoherenceComplex.real) * 180 / numpy.pi
601 602 else:
602 603 data = numpy.abs(avgcoherenceComplex)
603 604
604 605 z.append(data)
605 606
606 607 return numpy.array(z)
607 608
608 609 def setValue(self, value):
609 610
610 611 print("This property should not be initialized")
611 612
612 613 return
613 614
614 615 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
615 616
616 617
617 618 class SpectraHeis(Spectra):
618 619
619 620 def __init__(self):
620 621
621 622 self.radarControllerHeaderObj = RadarControllerHeader()
622 623 self.systemHeaderObj = SystemHeader()
623 624 self.type = "SpectraHeis"
624 625 self.nProfiles = None
625 626 self.heightList = None
626 627 self.channelList = None
627 628 self.flagNoData = True
628 629 self.flagDiscontinuousBlock = False
629 630 self.utctime = None
630 631 self.blocksize = None
631 632 self.profileIndex = 0
632 633 self.nCohInt = 1
633 634 self.nIncohInt = 1
634 635
635 636 @property
636 637 def normFactor(self):
637 638 pwcode = 1
638 639 if self.flagDecodeData:
639 640 pwcode = numpy.sum(self.code[0]**2)
640 641
641 642 normFactor = self.nIncohInt * self.nCohInt * pwcode
642 643
643 644 return normFactor
644 645
645 646 @property
646 647 def timeInterval(self):
647 648
648 649 return self.ippSeconds * self.nCohInt * self.nIncohInt
649 650
650 651
651 652 class Fits(JROData):
652 653
653 654 def __init__(self):
654 655
655 656 self.type = "Fits"
656 657 self.nProfiles = None
657 658 self.heightList = None
658 659 self.channelList = None
659 660 self.flagNoData = True
660 661 self.utctime = None
661 662 self.nCohInt = 1
662 663 self.nIncohInt = 1
663 664 self.useLocalTime = True
664 665 self.profileIndex = 0
665 666 self.timeZone = 0
666 667
667 668 def getTimeRange(self):
668 669
669 670 datatime = []
670 671
671 672 datatime.append(self.ltctime)
672 673 datatime.append(self.ltctime + self.timeInterval)
673 674
674 675 datatime = numpy.array(datatime)
675 676
676 677 return datatime
677 678
678 679 def getChannelIndexList(self):
679 680
680 681 return list(range(self.nChannels))
681 682
682 683 def getNoise(self, type=1):
683 684
684 685
685 686 if type == 1:
686 687 noise = self.getNoisebyHildebrand()
687 688
688 689 if type == 2:
689 690 noise = self.getNoisebySort()
690 691
691 692 if type == 3:
692 693 noise = self.getNoisebyWindow()
693 694
694 695 return noise
695 696
696 697 @property
697 698 def timeInterval(self):
698 699
699 700 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
700 701
701 702 return timeInterval
702 703
703 704 @property
704 705 def ippSeconds(self):
705 706 '''
706 707 '''
707 708 return self.ipp_sec
708 709
709 710 noise = property(getNoise, "I'm the 'nHeights' property.")
710 711
711 712
712 713 class Correlation(JROData):
713 714
714 715 def __init__(self):
715 716 '''
716 717 Constructor
717 718 '''
718 719 self.radarControllerHeaderObj = RadarControllerHeader()
719 720 self.systemHeaderObj = SystemHeader()
720 721 self.type = "Correlation"
721 722 self.data = None
722 723 self.dtype = None
723 724 self.nProfiles = None
724 725 self.heightList = None
725 726 self.channelList = None
726 727 self.flagNoData = True
727 728 self.flagDiscontinuousBlock = False
728 729 self.utctime = None
729 730 self.timeZone = 0
730 731 self.dstFlag = None
731 732 self.errorCount = None
732 733 self.blocksize = None
733 734 self.flagDecodeData = False # asumo q la data no esta decodificada
734 735 self.flagDeflipData = False # asumo q la data no esta sin flip
735 736 self.pairsList = None
736 737 self.nPoints = None
737 738
738 739 def getPairsList(self):
739 740
740 741 return self.pairsList
741 742
742 743 def getNoise(self, mode=2):
743 744
744 745 indR = numpy.where(self.lagR == 0)[0][0]
745 746 indT = numpy.where(self.lagT == 0)[0][0]
746 747
747 748 jspectra0 = self.data_corr[:, :, indR, :]
748 749 jspectra = copy.copy(jspectra0)
749 750
750 751 num_chan = jspectra.shape[0]
751 752 num_hei = jspectra.shape[2]
752 753
753 754 freq_dc = jspectra.shape[1] / 2
754 755 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
755 756
756 757 if ind_vel[0] < 0:
757 758 ind_vel[list(range(0, 1))] = ind_vel[list(
758 759 range(0, 1))] + self.num_prof
759 760
760 761 if mode == 1:
761 762 jspectra[:, freq_dc, :] = (
762 763 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
763 764
764 765 if mode == 2:
765 766
766 767 vel = numpy.array([-2, -1, 1, 2])
767 768 xx = numpy.zeros([4, 4])
768 769
769 770 for fil in range(4):
770 771 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
771 772
772 773 xx_inv = numpy.linalg.inv(xx)
773 774 xx_aux = xx_inv[0, :]
774 775
775 776 for ich in range(num_chan):
776 777 yy = jspectra[ich, ind_vel, :]
777 778 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
778 779
779 780 junkid = jspectra[ich, freq_dc, :] <= 0
780 781 cjunkid = sum(junkid)
781 782
782 783 if cjunkid.any():
783 784 jspectra[ich, freq_dc, junkid.nonzero()] = (
784 785 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
785 786
786 787 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
787 788
788 789 return noise
789 790
790 791 @property
791 792 def timeInterval(self):
792 793
793 794 return self.ippSeconds * self.nCohInt * self.nProfiles
794 795
795 796 def splitFunctions(self):
796 797
797 798 pairsList = self.pairsList
798 799 ccf_pairs = []
799 800 acf_pairs = []
800 801 ccf_ind = []
801 802 acf_ind = []
802 803 for l in range(len(pairsList)):
803 804 chan0 = pairsList[l][0]
804 805 chan1 = pairsList[l][1]
805 806
806 807 # Obteniendo pares de Autocorrelacion
807 808 if chan0 == chan1:
808 809 acf_pairs.append(chan0)
809 810 acf_ind.append(l)
810 811 else:
811 812 ccf_pairs.append(pairsList[l])
812 813 ccf_ind.append(l)
813 814
814 815 data_acf = self.data_cf[acf_ind]
815 816 data_ccf = self.data_cf[ccf_ind]
816 817
817 818 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
818 819
819 820 @property
820 821 def normFactor(self):
821 822 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
822 823 acf_pairs = numpy.array(acf_pairs)
823 824 normFactor = numpy.zeros((self.nPairs, self.nHeights))
824 825
825 826 for p in range(self.nPairs):
826 827 pair = self.pairsList[p]
827 828
828 829 ch0 = pair[0]
829 830 ch1 = pair[1]
830 831
831 832 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
832 833 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
833 834 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
834 835
835 836 return normFactor
836 837
837 838
838 839 class Parameters(Spectra):
839 840
840 841 groupList = None # List of Pairs, Groups, etc
841 842 data_param = None # Parameters obtained
842 843 data_pre = None # Data Pre Parametrization
843 844 data_SNR = None # Signal to Noise Ratio
844 845 abscissaList = None # Abscissa, can be velocities, lags or time
845 846 utctimeInit = None # Initial UTC time
846 847 paramInterval = None # Time interval to calculate Parameters in seconds
847 848 useLocalTime = True
848 849 # Fitting
849 850 data_error = None # Error of the estimation
850 851 constants = None
851 852 library = None
852 853 # Output signal
853 854 outputInterval = None # Time interval to calculate output signal in seconds
854 855 data_output = None # Out signal
855 856 nAvg = None
856 857 noise_estimation = None
857 858 GauSPC = None # Fit gaussian SPC
858 859
859 860 def __init__(self):
860 861 '''
861 862 Constructor
862 863 '''
863 864 self.radarControllerHeaderObj = RadarControllerHeader()
864 865 self.systemHeaderObj = SystemHeader()
865 866 self.type = "Parameters"
866 867 self.timeZone = 0
867 868
868 869 def getTimeRange1(self, interval):
869 870
870 871 datatime = []
871 872
872 873 if self.useLocalTime:
873 874 time1 = self.utctimeInit - self.timeZone * 60
874 875 else:
875 876 time1 = self.utctimeInit
876 877
877 878 datatime.append(time1)
878 879 datatime.append(time1 + interval)
879 880 datatime = numpy.array(datatime)
880 881
881 882 return datatime
882 883
883 884 @property
884 885 def timeInterval(self):
885 886
886 887 if hasattr(self, 'timeInterval1'):
887 888 return self.timeInterval1
888 889 else:
889 890 return self.paramInterval
890 891
891 892 def setValue(self, value):
892 893
893 894 print("This property should not be initialized")
894 895
895 896 return
896 897
897 898 def getNoise(self):
898 899
899 900 return self.spc_noise
900 901
901 902 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
902 903
903 904
904 905 class PlotterData(object):
905 906 '''
906 907 Object to hold data to be plotted
907 908 '''
908 909
909 910 MAXNUMX = 200
910 911 MAXNUMY = 200
911 912
912 913 def __init__(self, code, exp_code, localtime=True):
913 914
914 915 self.key = code
915 916 self.exp_code = exp_code
916 917 self.ready = False
917 918 self.flagNoData = False
918 919 self.localtime = localtime
919 920 self.data = {}
920 921 self.meta = {}
921 922 self.__heights = []
922 923
923 924 def __str__(self):
924 925 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
925 926 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
926 927
927 928 def __len__(self):
928 929 return len(self.data)
929 930
930 931 def __getitem__(self, key):
931 932 if isinstance(key, int):
932 933 return self.data[self.times[key]]
933 934 elif isinstance(key, str):
934 935 ret = numpy.array([self.data[x][key] for x in self.times])
935 936 if ret.ndim > 1:
936 937 ret = numpy.swapaxes(ret, 0, 1)
937 938 return ret
938 939
939 940 def __contains__(self, key):
940 941 return key in self.data[self.min_time]
941 942
942 943 def setup(self):
943 944 '''
944 945 Configure object
945 946 '''
946 947 self.type = ''
947 948 self.ready = False
948 949 del self.data
949 950 self.data = {}
950 951 self.__heights = []
951 952 self.__all_heights = set()
952 953
953 954 def shape(self, key):
954 955 '''
955 956 Get the shape of the one-element data for the given key
956 957 '''
957 958
958 959 if len(self.data[self.min_time][key]):
959 960 return self.data[self.min_time][key].shape
960 961 return (0,)
961 962
962 963 def update(self, data, tm, meta={}):
963 964 '''
964 965 Update data object with new dataOut
965 966 '''
966 967
967 968 self.data[tm] = data
968 969
969 970 for key, value in meta.items():
970 971 setattr(self, key, value)
971 972
972 973 def normalize_heights(self):
973 974 '''
974 975 Ensure same-dimension of the data for different heighList
975 976 '''
976 977
977 978 H = numpy.array(list(self.__all_heights))
978 979 H.sort()
979 980 for key in self.data:
980 981 shape = self.shape(key)[:-1] + H.shape
981 982 for tm, obj in list(self.data[key].items()):
982 983 h = self.__heights[self.times.tolist().index(tm)]
983 984 if H.size == h.size:
984 985 continue
985 986 index = numpy.where(numpy.in1d(H, h))[0]
986 987 dummy = numpy.zeros(shape) + numpy.nan
987 988 if len(shape) == 2:
988 989 dummy[:, index] = obj
989 990 else:
990 991 dummy[index] = obj
991 992 self.data[key][tm] = dummy
992 993
993 994 self.__heights = [H for tm in self.times]
994 995
995 996 def jsonify(self, tm, plot_name, plot_type, decimate=False):
996 997 '''
997 998 Convert data to json
998 999 '''
999 1000
1000 1001 meta = {}
1001 1002 meta['xrange'] = []
1002 1003 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1003 1004 tmp = self.data[tm][self.key]
1004 1005 shape = tmp.shape
1005 1006 if len(shape) == 2:
1006 1007 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1007 1008 elif len(shape) == 3:
1008 1009 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1009 1010 data = self.roundFloats(
1010 1011 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1011 1012 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1012 1013 else:
1013 1014 data = self.roundFloats(self.data[tm][self.key].tolist())
1014 1015
1015 1016 ret = {
1016 1017 'plot': plot_name,
1017 1018 'code': self.exp_code,
1018 1019 'time': float(tm),
1019 1020 'data': data,
1020 1021 }
1021 1022 meta['type'] = plot_type
1022 1023 meta['interval'] = float(self.interval)
1023 1024 meta['localtime'] = self.localtime
1024 1025 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1025 1026 meta.update(self.meta)
1026 1027 ret['metadata'] = meta
1027 1028 return json.dumps(ret)
1028 1029
1029 1030 @property
1030 1031 def times(self):
1031 1032 '''
1032 1033 Return the list of times of the current data
1033 1034 '''
1034 1035
1035 1036 ret = [t for t in self.data]
1036 1037 ret.sort()
1037 1038 return numpy.array(ret)
1038 1039
1039 1040 @property
1040 1041 def min_time(self):
1041 1042 '''
1042 1043 Return the minimun time value
1043 1044 '''
1044 1045
1045 1046 return self.times[0]
1046 1047
1047 1048 @property
1048 1049 def max_time(self):
1049 1050 '''
1050 1051 Return the maximun time value
1051 1052 '''
1052 1053
1053 1054 return self.times[-1]
1054 1055
1055 1056 # @property
1056 1057 # def heights(self):
1057 1058 # '''
1058 1059 # Return the list of heights of the current data
1059 1060 # '''
1060 1061
1061 1062 # return numpy.array(self.__heights[-1])
1062 1063
1063 1064 @staticmethod
1064 1065 def roundFloats(obj):
1065 1066 if isinstance(obj, list):
1066 1067 return list(map(PlotterData.roundFloats, obj))
1067 1068 elif isinstance(obj, float):
1068 1069 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,1885 +1,1886
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62 62 for channel in channelList:
63 63 if channel not in self.dataOut.channelList:
64 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
65 65
66 66 index = self.dataOut.channelList.index(channel)
67 67 channelIndexList.append(index)
68 68 self.selectChannelsByIndex(channelIndexList)
69 69 return self.dataOut
70 70
71 71 def selectChannelsByIndex(self, channelIndexList):
72 72 """
73 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
74 74
75 75 Input:
76 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
77 77
78 78 Affected:
79 79 self.dataOut.data
80 80 self.dataOut.channelIndexList
81 81 self.dataOut.nChannels
82 82 self.dataOut.m_ProcessingHeader.totalSpectra
83 83 self.dataOut.systemHeaderObj.numChannels
84 84 self.dataOut.m_ProcessingHeader.blockSize
85 85
86 86 Return:
87 87 None
88 88 """
89 89
90 90 for channelIndex in channelIndexList:
91 91 if channelIndex not in self.dataOut.channelIndexList:
92 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
93 93
94 94 if self.dataOut.type == 'Voltage':
95 95 if self.dataOut.flagDataAsBlock:
96 96 """
97 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
98 98 """
99 99 data = self.dataOut.data[channelIndexList,:,:]
100 100 else:
101 101 data = self.dataOut.data[channelIndexList,:]
102 102
103 103 self.dataOut.data = data
104 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
105 105 self.dataOut.channelList = range(len(channelIndexList))
106 106
107 107 elif self.dataOut.type == 'Spectra':
108 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
109 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
110 110
111 111 self.dataOut.data_spc = data_spc
112 112 self.dataOut.data_dc = data_dc
113 113
114 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 self.dataOut.channelList = range(len(channelIndexList))
116 116 self.__selectPairsByChannel(channelIndexList)
117 117
118 118 return 1
119 119
120 120 def __selectPairsByChannel(self, channelList=None):
121 121
122 122 if channelList == None:
123 123 return
124 124
125 125 pairsIndexListSelected = []
126 126 for pairIndex in self.dataOut.pairsIndexList:
127 127 # First pair
128 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
129 129 continue
130 130 # Second pair
131 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
132 132 continue
133 133
134 134 pairsIndexListSelected.append(pairIndex)
135 135
136 136 if not pairsIndexListSelected:
137 137 self.dataOut.data_cspc = None
138 138 self.dataOut.pairsList = []
139 139 return
140 140
141 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
142 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
143 143 for i in pairsIndexListSelected]
144 144
145 145 return
146 146
147 147 class selectHeights(Operation):
148 148
149 149 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
150 150 """
151 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
152 152 minHei <= height <= maxHei
153 153
154 154 Input:
155 155 minHei : valor minimo de altura a considerar
156 156 maxHei : valor maximo de altura a considerar
157 157
158 158 Affected:
159 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
160 160
161 161 Return:
162 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
163 163 """
164 164
165 165 self.dataOut = dataOut
166 166
167 167 if minHei and maxHei:
168 168
169 169 if (minHei < self.dataOut.heightList[0]):
170 170 minHei = self.dataOut.heightList[0]
171 171
172 172 if (maxHei > self.dataOut.heightList[-1]):
173 173 maxHei = self.dataOut.heightList[-1]
174 174
175 175 minIndex = 0
176 176 maxIndex = 0
177 177 heights = self.dataOut.heightList
178 178
179 179 inda = numpy.where(heights >= minHei)
180 180 indb = numpy.where(heights <= maxHei)
181 181
182 182 try:
183 183 minIndex = inda[0][0]
184 184 except:
185 185 minIndex = 0
186 186
187 187 try:
188 188 maxIndex = indb[0][-1]
189 189 except:
190 190 maxIndex = len(heights)
191 191
192 192 self.selectHeightsByIndex(minIndex, maxIndex)
193 193
194 194 return self.dataOut
195 195
196 196 def selectHeightsByIndex(self, minIndex, maxIndex):
197 197 """
198 198 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
199 199 minIndex <= index <= maxIndex
200 200
201 201 Input:
202 202 minIndex : valor de indice minimo de altura a considerar
203 203 maxIndex : valor de indice maximo de altura a considerar
204 204
205 205 Affected:
206 206 self.dataOut.data
207 207 self.dataOut.heightList
208 208
209 209 Return:
210 210 1 si el metodo se ejecuto con exito caso contrario devuelve 0
211 211 """
212 212
213 213 if self.dataOut.type == 'Voltage':
214 214 if (minIndex < 0) or (minIndex > maxIndex):
215 215 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
216 216
217 217 if (maxIndex >= self.dataOut.nHeights):
218 218 maxIndex = self.dataOut.nHeights
219 219 #print("shapeeee",self.dataOut.data.shape)
220 220 #voltage
221 221 if self.dataOut.flagDataAsBlock:
222 222 """
223 223 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
224 224 """
225 225 data = self.dataOut.data[:,:, minIndex:maxIndex]
226 226 else:
227 227 data = self.dataOut.data[:, minIndex:maxIndex]
228 228
229 229 # firstHeight = self.dataOut.heightList[minIndex]
230 230
231 231 self.dataOut.data = data
232 232 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
233 233
234 234 if self.dataOut.nHeights <= 1:
235 235 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
236 236 elif self.dataOut.type == 'Spectra':
237 237 if (minIndex < 0) or (minIndex > maxIndex):
238 238 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
239 239 minIndex, maxIndex))
240 240
241 241 if (maxIndex >= self.dataOut.nHeights):
242 242 maxIndex = self.dataOut.nHeights - 1
243 243
244 244 # Spectra
245 245 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
246 246
247 247 data_cspc = None
248 248 if self.dataOut.data_cspc is not None:
249 249 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
250 250
251 251 data_dc = None
252 252 if self.dataOut.data_dc is not None:
253 253 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
254 254
255 255 self.dataOut.data_spc = data_spc
256 256 self.dataOut.data_cspc = data_cspc
257 257 self.dataOut.data_dc = data_dc
258 258
259 259 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
260 260
261 261 return 1
262 262
263 263
264 264 class filterByHeights(Operation):
265 265
266 266 def run(self, dataOut, window):
267 267
268 268 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
269 269
270 270 if window == None:
271 271 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
272 272
273 273 newdelta = deltaHeight * window
274 274 r = dataOut.nHeights % window
275 275 newheights = (dataOut.nHeights-r)/window
276 276
277 277 if newheights <= 1:
278 278 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
279 279
280 280 if dataOut.flagDataAsBlock:
281 281 """
282 282 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
283 283 """
284 284 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
285 285 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
286 286 buffer = numpy.sum(buffer,3)
287 287
288 288 else:
289 289 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
290 290 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
291 291 buffer = numpy.sum(buffer,2)
292 292
293 293 dataOut.data = buffer
294 294 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
295 295 dataOut.windowOfFilter = window
296 296
297 297 return dataOut
298 298
299 299
300 300 class setH0(Operation):
301 301
302 302 def run(self, dataOut, h0, deltaHeight = None):
303 303
304 304 if not deltaHeight:
305 305 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
306 306
307 307 nHeights = dataOut.nHeights
308 308
309 309 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
310 310
311 311 dataOut.heightList = newHeiRange
312 dataOut.h0 = h0
312 313
313 314 return dataOut
314 315
315 316
316 317 class deFlip(Operation):
317 318
318 319 def run(self, dataOut, channelList = []):
319 320
320 321 data = dataOut.data.copy()
321 322
322 323 if dataOut.flagDataAsBlock:
323 324 flip = self.flip
324 325 profileList = list(range(dataOut.nProfiles))
325 326
326 327 if not channelList:
327 328 for thisProfile in profileList:
328 329 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
329 330 flip *= -1.0
330 331 else:
331 332 for thisChannel in channelList:
332 333 if thisChannel not in dataOut.channelList:
333 334 continue
334 335
335 336 for thisProfile in profileList:
336 337 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
337 338 flip *= -1.0
338 339
339 340 self.flip = flip
340 341
341 342 else:
342 343 if not channelList:
343 344 data[:,:] = data[:,:]*self.flip
344 345 else:
345 346 for thisChannel in channelList:
346 347 if thisChannel not in dataOut.channelList:
347 348 continue
348 349
349 350 data[thisChannel,:] = data[thisChannel,:]*self.flip
350 351
351 352 self.flip *= -1.
352 353
353 354 dataOut.data = data
354 355
355 356 return dataOut
356 357
357 358
358 359 class setAttribute(Operation):
359 360 '''
360 361 Set an arbitrary attribute(s) to dataOut
361 362 '''
362 363
363 364 def __init__(self):
364 365
365 366 Operation.__init__(self)
366 367 self._ready = False
367 368
368 369 def run(self, dataOut, **kwargs):
369 370
370 371 for key, value in kwargs.items():
371 372 setattr(dataOut, key, value)
372 373
373 374 return dataOut
374 375
375 376
376 377 @MPDecorator
377 378 class printAttribute(Operation):
378 379 '''
379 380 Print an arbitrary attribute of dataOut
380 381 '''
381 382
382 383 def __init__(self):
383 384
384 385 Operation.__init__(self)
385 386
386 387 def run(self, dataOut, attributes):
387 388
388 389 if isinstance(attributes, str):
389 390 attributes = [attributes]
390 391 for attr in attributes:
391 392 if hasattr(dataOut, attr):
392 393 log.log(getattr(dataOut, attr), attr)
393 394
394 395
395 396 class interpolateHeights(Operation):
396 397
397 398 def run(self, dataOut, topLim, botLim):
398 399 #69 al 72 para julia
399 400 #82-84 para meteoros
400 401 if len(numpy.shape(dataOut.data))==2:
401 402 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
402 403 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
403 404 #dataOut.data[:,botLim:limSup+1] = sampInterp
404 405 dataOut.data[:,botLim:topLim+1] = sampInterp
405 406 else:
406 407 nHeights = dataOut.data.shape[2]
407 408 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
408 409 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
409 410 f = interpolate.interp1d(x, y, axis = 2)
410 411 xnew = numpy.arange(botLim,topLim+1)
411 412 ynew = f(xnew)
412 413 dataOut.data[:,:,botLim:topLim+1] = ynew
413 414
414 415 return dataOut
415 416
416 417
417 418 class CohInt(Operation):
418 419
419 420 isConfig = False
420 421 __profIndex = 0
421 422 __byTime = False
422 423 __initime = None
423 424 __lastdatatime = None
424 425 __integrationtime = None
425 426 __buffer = None
426 427 __bufferStride = []
427 428 __dataReady = False
428 429 __profIndexStride = 0
429 430 __dataToPutStride = False
430 431 n = None
431 432
432 433 def __init__(self, **kwargs):
433 434
434 435 Operation.__init__(self, **kwargs)
435 436
436 437 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
437 438 """
438 439 Set the parameters of the integration class.
439 440
440 441 Inputs:
441 442
442 443 n : Number of coherent integrations
443 444 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
444 445 overlapping :
445 446 """
446 447
447 448 self.__initime = None
448 449 self.__lastdatatime = 0
449 450 self.__buffer = None
450 451 self.__dataReady = False
451 452 self.byblock = byblock
452 453 self.stride = stride
453 454
454 455 if n == None and timeInterval == None:
455 456 raise ValueError("n or timeInterval should be specified ...")
456 457
457 458 if n != None:
458 459 self.n = n
459 460 self.__byTime = False
460 461 else:
461 462 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
462 463 self.n = 9999
463 464 self.__byTime = True
464 465
465 466 if overlapping:
466 467 self.__withOverlapping = True
467 468 self.__buffer = None
468 469 else:
469 470 self.__withOverlapping = False
470 471 self.__buffer = 0
471 472
472 473 self.__profIndex = 0
473 474
474 475 def putData(self, data):
475 476
476 477 """
477 478 Add a profile to the __buffer and increase in one the __profileIndex
478 479
479 480 """
480 481
481 482 if not self.__withOverlapping:
482 483 self.__buffer += data.copy()
483 484 self.__profIndex += 1
484 485 return
485 486
486 487 #Overlapping data
487 488 nChannels, nHeis = data.shape
488 489 data = numpy.reshape(data, (1, nChannels, nHeis))
489 490
490 491 #If the buffer is empty then it takes the data value
491 492 if self.__buffer is None:
492 493 self.__buffer = data
493 494 self.__profIndex += 1
494 495 return
495 496
496 497 #If the buffer length is lower than n then stakcing the data value
497 498 if self.__profIndex < self.n:
498 499 self.__buffer = numpy.vstack((self.__buffer, data))
499 500 self.__profIndex += 1
500 501 return
501 502
502 503 #If the buffer length is equal to n then replacing the last buffer value with the data value
503 504 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
504 505 self.__buffer[self.n-1] = data
505 506 self.__profIndex = self.n
506 507 return
507 508
508 509
509 510 def pushData(self):
510 511 """
511 512 Return the sum of the last profiles and the profiles used in the sum.
512 513
513 514 Affected:
514 515
515 516 self.__profileIndex
516 517
517 518 """
518 519
519 520 if not self.__withOverlapping:
520 521 data = self.__buffer
521 522 n = self.__profIndex
522 523
523 524 self.__buffer = 0
524 525 self.__profIndex = 0
525 526
526 527 return data, n
527 528
528 529 #Integration with Overlapping
529 530 data = numpy.sum(self.__buffer, axis=0)
530 531 # print data
531 532 # raise
532 533 n = self.__profIndex
533 534
534 535 return data, n
535 536
536 537 def byProfiles(self, data):
537 538
538 539 self.__dataReady = False
539 540 avgdata = None
540 541 # n = None
541 542 # print data
542 543 # raise
543 544 self.putData(data)
544 545
545 546 if self.__profIndex == self.n:
546 547 avgdata, n = self.pushData()
547 548 self.__dataReady = True
548 549
549 550 return avgdata
550 551
551 552 def byTime(self, data, datatime):
552 553
553 554 self.__dataReady = False
554 555 avgdata = None
555 556 n = None
556 557
557 558 self.putData(data)
558 559
559 560 if (datatime - self.__initime) >= self.__integrationtime:
560 561 avgdata, n = self.pushData()
561 562 self.n = n
562 563 self.__dataReady = True
563 564
564 565 return avgdata
565 566
566 567 def integrateByStride(self, data, datatime):
567 568 # print data
568 569 if self.__profIndex == 0:
569 570 self.__buffer = [[data.copy(), datatime]]
570 571 else:
571 572 self.__buffer.append([data.copy(),datatime])
572 573 self.__profIndex += 1
573 574 self.__dataReady = False
574 575
575 576 if self.__profIndex == self.n * self.stride :
576 577 self.__dataToPutStride = True
577 578 self.__profIndexStride = 0
578 579 self.__profIndex = 0
579 580 self.__bufferStride = []
580 581 for i in range(self.stride):
581 582 current = self.__buffer[i::self.stride]
582 583 data = numpy.sum([t[0] for t in current], axis=0)
583 584 avgdatatime = numpy.average([t[1] for t in current])
584 585 # print data
585 586 self.__bufferStride.append((data, avgdatatime))
586 587
587 588 if self.__dataToPutStride:
588 589 self.__dataReady = True
589 590 self.__profIndexStride += 1
590 591 if self.__profIndexStride == self.stride:
591 592 self.__dataToPutStride = False
592 593 # print self.__bufferStride[self.__profIndexStride - 1]
593 594 # raise
594 595 return self.__bufferStride[self.__profIndexStride - 1]
595 596
596 597
597 598 return None, None
598 599
599 600 def integrate(self, data, datatime=None):
600 601
601 602 if self.__initime == None:
602 603 self.__initime = datatime
603 604
604 605 if self.__byTime:
605 606 avgdata = self.byTime(data, datatime)
606 607 else:
607 608 avgdata = self.byProfiles(data)
608 609
609 610
610 611 self.__lastdatatime = datatime
611 612
612 613 if avgdata is None:
613 614 return None, None
614 615
615 616 avgdatatime = self.__initime
616 617
617 618 deltatime = datatime - self.__lastdatatime
618 619
619 620 if not self.__withOverlapping:
620 621 self.__initime = datatime
621 622 else:
622 623 self.__initime += deltatime
623 624
624 625 return avgdata, avgdatatime
625 626
626 627 def integrateByBlock(self, dataOut):
627 628
628 629 times = int(dataOut.data.shape[1]/self.n)
629 630 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
630 631
631 632 id_min = 0
632 633 id_max = self.n
633 634
634 635 for i in range(times):
635 636 junk = dataOut.data[:,id_min:id_max,:]
636 637 avgdata[:,i,:] = junk.sum(axis=1)
637 638 id_min += self.n
638 639 id_max += self.n
639 640
640 641 timeInterval = dataOut.ippSeconds*self.n
641 642 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
642 643 self.__dataReady = True
643 644 return avgdata, avgdatatime
644 645
645 646 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
646 647
647 648 if not self.isConfig:
648 649 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
649 650 self.isConfig = True
650 651
651 652 if dataOut.flagDataAsBlock:
652 653 """
653 654 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
654 655 """
655 656 avgdata, avgdatatime = self.integrateByBlock(dataOut)
656 657 dataOut.nProfiles /= self.n
657 658 else:
658 659 if stride is None:
659 660 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
660 661 else:
661 662 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
662 663
663 664
664 665 # dataOut.timeInterval *= n
665 666 dataOut.flagNoData = True
666 667
667 668 if self.__dataReady:
668 669 dataOut.data = avgdata
669 670 if not dataOut.flagCohInt:
670 671 dataOut.nCohInt *= self.n
671 672 dataOut.flagCohInt = True
672 673 ####################################dataOut.utctime = avgdatatime
673 674 # print avgdata, avgdatatime
674 675 # raise
675 676 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
676 677 dataOut.flagNoData = False
677 678 return dataOut
678 679
679 680 class Decoder(Operation):
680 681
681 682 isConfig = False
682 683 __profIndex = 0
683 684
684 685 code = None
685 686
686 687 nCode = None
687 688 nBaud = None
688 689
689 690 def __init__(self, **kwargs):
690 691
691 692 Operation.__init__(self, **kwargs)
692 693
693 694 self.times = None
694 695 self.osamp = None
695 696 # self.__setValues = False
696 697 self.isConfig = False
697 698 self.setupReq = False
698 699 def setup(self, code, osamp, dataOut):
699 700
700 701 self.__profIndex = 0
701 702
702 703 self.code = code
703 704
704 705 self.nCode = len(code)
705 706 self.nBaud = len(code[0])
706 707
707 708 if (osamp != None) and (osamp >1):
708 709 self.osamp = osamp
709 710 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
710 711 self.nBaud = self.nBaud*self.osamp
711 712
712 713 self.__nChannels = dataOut.nChannels
713 714 self.__nProfiles = dataOut.nProfiles
714 715 self.__nHeis = dataOut.nHeights
715 716
716 717 if self.__nHeis < self.nBaud:
717 718 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
718 719
719 720 #Frequency
720 721 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
721 722
722 723 __codeBuffer[:,0:self.nBaud] = self.code
723 724
724 725 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
725 726
726 727 if dataOut.flagDataAsBlock:
727 728
728 729 self.ndatadec = self.__nHeis #- self.nBaud + 1
729 730
730 731 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
731 732
732 733 else:
733 734
734 735 #Time
735 736 self.ndatadec = self.__nHeis #- self.nBaud + 1
736 737
737 738 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
738 739
739 740 def __convolutionInFreq(self, data):
740 741
741 742 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
742 743
743 744 fft_data = numpy.fft.fft(data, axis=1)
744 745
745 746 conv = fft_data*fft_code
746 747
747 748 data = numpy.fft.ifft(conv,axis=1)
748 749
749 750 return data
750 751
751 752 def __convolutionInFreqOpt(self, data):
752 753
753 754 raise NotImplementedError
754 755
755 756 def __convolutionInTime(self, data):
756 757
757 758 code = self.code[self.__profIndex]
758 759 for i in range(self.__nChannels):
759 760 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
760 761
761 762 return self.datadecTime
762 763
763 764 def __convolutionByBlockInTime(self, data):
764 765
765 766 repetitions = int(self.__nProfiles / self.nCode)
766 767 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
767 768 junk = junk.flatten()
768 769 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
769 770 profilesList = range(self.__nProfiles)
770 771
771 772 for i in range(self.__nChannels):
772 773 for j in profilesList:
773 774 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
774 775 return self.datadecTime
775 776
776 777 def __convolutionByBlockInFreq(self, data):
777 778
778 779 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
779 780
780 781
781 782 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
782 783
783 784 fft_data = numpy.fft.fft(data, axis=2)
784 785
785 786 conv = fft_data*fft_code
786 787
787 788 data = numpy.fft.ifft(conv,axis=2)
788 789
789 790 return data
790 791
791 792
792 793 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
793 794
794 795 if dataOut.flagDecodeData:
795 796 print("This data is already decoded, recoding again ...")
796 797
797 798 if not self.isConfig:
798 799
799 800 if code is None:
800 801 if dataOut.code is None:
801 802 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
802 803
803 804 code = dataOut.code
804 805 else:
805 806 code = numpy.array(code).reshape(nCode,nBaud)
806 807 self.setup(code, osamp, dataOut)
807 808
808 809 self.isConfig = True
809 810
810 811 if mode == 3:
811 812 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
812 813
813 814 if times != None:
814 815 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
815 816
816 817 if self.code is None:
817 818 print("Fail decoding: Code is not defined.")
818 819 return
819 820
820 821 self.__nProfiles = dataOut.nProfiles
821 822 datadec = None
822 823
823 824 if mode == 3:
824 825 mode = 0
825 826
826 827 if dataOut.flagDataAsBlock:
827 828 """
828 829 Decoding when data have been read as block,
829 830 """
830 831
831 832 if mode == 0:
832 833 datadec = self.__convolutionByBlockInTime(dataOut.data)
833 834 if mode == 1:
834 835 datadec = self.__convolutionByBlockInFreq(dataOut.data)
835 836 else:
836 837 """
837 838 Decoding when data have been read profile by profile
838 839 """
839 840 if mode == 0:
840 841 datadec = self.__convolutionInTime(dataOut.data)
841 842
842 843 if mode == 1:
843 844 datadec = self.__convolutionInFreq(dataOut.data)
844 845
845 846 if mode == 2:
846 847 datadec = self.__convolutionInFreqOpt(dataOut.data)
847 848
848 849 if datadec is None:
849 850 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
850 851
851 852 dataOut.code = self.code
852 853 dataOut.nCode = self.nCode
853 854 dataOut.nBaud = self.nBaud
854 855
855 856 dataOut.data = datadec
856 857
857 858 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
858 859
859 860 dataOut.flagDecodeData = True #asumo q la data esta decodificada
860 861
861 862 if self.__profIndex == self.nCode-1:
862 863 self.__profIndex = 0
863 864 return dataOut
864 865
865 866 self.__profIndex += 1
866 867
867 868 return dataOut
868 869 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
869 870
870 871
871 872 class ProfileConcat(Operation):
872 873
873 874 isConfig = False
874 875 buffer = None
875 876
876 877 def __init__(self, **kwargs):
877 878
878 879 Operation.__init__(self, **kwargs)
879 880 self.profileIndex = 0
880 881
881 882 def reset(self):
882 883 self.buffer = numpy.zeros_like(self.buffer)
883 884 self.start_index = 0
884 885 self.times = 1
885 886
886 887 def setup(self, data, m, n=1):
887 888 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
888 889 self.nHeights = data.shape[1]#.nHeights
889 890 self.start_index = 0
890 891 self.times = 1
891 892
892 893 def concat(self, data):
893 894
894 895 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
895 896 self.start_index = self.start_index + self.nHeights
896 897
897 898 def run(self, dataOut, m):
898 899 dataOut.flagNoData = True
899 900
900 901 if not self.isConfig:
901 902 self.setup(dataOut.data, m, 1)
902 903 self.isConfig = True
903 904
904 905 if dataOut.flagDataAsBlock:
905 906 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
906 907
907 908 else:
908 909 self.concat(dataOut.data)
909 910 self.times += 1
910 911 if self.times > m:
911 912 dataOut.data = self.buffer
912 913 self.reset()
913 914 dataOut.flagNoData = False
914 915 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
915 916 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
916 917 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
917 918 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
918 919 dataOut.ippSeconds *= m
919 920 return dataOut
920 921
921 922 class ProfileSelector(Operation):
922 923
923 924 profileIndex = None
924 925 # Tamanho total de los perfiles
925 926 nProfiles = None
926 927
927 928 def __init__(self, **kwargs):
928 929
929 930 Operation.__init__(self, **kwargs)
930 931 self.profileIndex = 0
931 932
932 933 def incProfileIndex(self):
933 934
934 935 self.profileIndex += 1
935 936
936 937 if self.profileIndex >= self.nProfiles:
937 938 self.profileIndex = 0
938 939
939 940 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
940 941
941 942 if profileIndex < minIndex:
942 943 return False
943 944
944 945 if profileIndex > maxIndex:
945 946 return False
946 947
947 948 return True
948 949
949 950 def isThisProfileInList(self, profileIndex, profileList):
950 951
951 952 if profileIndex not in profileList:
952 953 return False
953 954
954 955 return True
955 956
956 957 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
957 958 #print("before",dataOut.data.shape)
958 959 """
959 960 ProfileSelector:
960 961
961 962 Inputs:
962 963 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
963 964
964 965 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
965 966
966 967 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
967 968
968 969 """
969 970
970 971 if rangeList is not None:
971 972 if type(rangeList[0]) not in (tuple, list):
972 973 rangeList = [rangeList]
973 974
974 975 dataOut.flagNoData = True
975 976
976 977 if dataOut.flagDataAsBlock:
977 978 """
978 979 data dimension = [nChannels, nProfiles, nHeis]
979 980 """
980 981 if profileList != None:
981 982 dataOut.data = dataOut.data[:,profileList,:]
982 983
983 984 if profileRangeList != None:
984 985 minIndex = profileRangeList[0]
985 986 maxIndex = profileRangeList[1]
986 987 profileList = list(range(minIndex, maxIndex+1))
987 988
988 989 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
989 990
990 991 if rangeList != None:
991 992
992 993 profileList = []
993 994
994 995 for thisRange in rangeList:
995 996 minIndex = thisRange[0]
996 997 maxIndex = thisRange[1]
997 998
998 999 profileList.extend(list(range(minIndex, maxIndex+1)))
999 1000
1000 1001 dataOut.data = dataOut.data[:,profileList,:]
1001 1002
1002 1003 dataOut.nProfiles = len(profileList)
1003 1004 dataOut.profileIndex = dataOut.nProfiles - 1
1004 1005 dataOut.flagNoData = False
1005 1006 #print(dataOut.data.shape)
1006 1007 return dataOut
1007 1008
1008 1009 """
1009 1010 data dimension = [nChannels, nHeis]
1010 1011 """
1011 1012
1012 1013 if profileList != None:
1013 1014
1014 1015 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1015 1016
1016 1017 self.nProfiles = len(profileList)
1017 1018 dataOut.nProfiles = self.nProfiles
1018 1019 dataOut.profileIndex = self.profileIndex
1019 1020 dataOut.flagNoData = False
1020 1021
1021 1022 self.incProfileIndex()
1022 1023 return dataOut
1023 1024
1024 1025 if profileRangeList != None:
1025 1026
1026 1027 minIndex = profileRangeList[0]
1027 1028 maxIndex = profileRangeList[1]
1028 1029
1029 1030 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1030 1031
1031 1032 self.nProfiles = maxIndex - minIndex + 1
1032 1033 dataOut.nProfiles = self.nProfiles
1033 1034 dataOut.profileIndex = self.profileIndex
1034 1035 dataOut.flagNoData = False
1035 1036
1036 1037 self.incProfileIndex()
1037 1038 return dataOut
1038 1039
1039 1040 if rangeList != None:
1040 1041
1041 1042 nProfiles = 0
1042 1043
1043 1044 for thisRange in rangeList:
1044 1045 minIndex = thisRange[0]
1045 1046 maxIndex = thisRange[1]
1046 1047
1047 1048 nProfiles += maxIndex - minIndex + 1
1048 1049
1049 1050 for thisRange in rangeList:
1050 1051
1051 1052 minIndex = thisRange[0]
1052 1053 maxIndex = thisRange[1]
1053 1054
1054 1055 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1055 1056
1056 1057 self.nProfiles = nProfiles
1057 1058 dataOut.nProfiles = self.nProfiles
1058 1059 dataOut.profileIndex = self.profileIndex
1059 1060 dataOut.flagNoData = False
1060 1061
1061 1062 self.incProfileIndex()
1062 1063
1063 1064 break
1064 1065
1065 1066 return dataOut
1066 1067
1067 1068
1068 1069 if beam != None: #beam is only for AMISR data
1069 1070 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1070 1071 dataOut.flagNoData = False
1071 1072 dataOut.profileIndex = self.profileIndex
1072 1073
1073 1074 self.incProfileIndex()
1074 1075
1075 1076 return dataOut
1076 1077
1077 1078 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1078 1079
1079 1080
1080 1081 class Reshaper(Operation):
1081 1082
1082 1083 def __init__(self, **kwargs):
1083 1084
1084 1085 Operation.__init__(self, **kwargs)
1085 1086
1086 1087 self.__buffer = None
1087 1088 self.__nitems = 0
1088 1089
1089 1090 def __appendProfile(self, dataOut, nTxs):
1090 1091
1091 1092 if self.__buffer is None:
1092 1093 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1093 1094 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1094 1095
1095 1096 ini = dataOut.nHeights * self.__nitems
1096 1097 end = ini + dataOut.nHeights
1097 1098
1098 1099 self.__buffer[:, ini:end] = dataOut.data
1099 1100
1100 1101 self.__nitems += 1
1101 1102
1102 1103 return int(self.__nitems*nTxs)
1103 1104
1104 1105 def __getBuffer(self):
1105 1106
1106 1107 if self.__nitems == int(1./self.__nTxs):
1107 1108
1108 1109 self.__nitems = 0
1109 1110
1110 1111 return self.__buffer.copy()
1111 1112
1112 1113 return None
1113 1114
1114 1115 def __checkInputs(self, dataOut, shape, nTxs):
1115 1116
1116 1117 if shape is None and nTxs is None:
1117 1118 raise ValueError("Reshaper: shape of factor should be defined")
1118 1119
1119 1120 if nTxs:
1120 1121 if nTxs < 0:
1121 1122 raise ValueError("nTxs should be greater than 0")
1122 1123
1123 1124 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1124 1125 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1125 1126
1126 1127 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1127 1128
1128 1129 return shape, nTxs
1129 1130
1130 1131 if len(shape) != 2 and len(shape) != 3:
1131 1132 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1132 1133
1133 1134 if len(shape) == 2:
1134 1135 shape_tuple = [dataOut.nChannels]
1135 1136 shape_tuple.extend(shape)
1136 1137 else:
1137 1138 shape_tuple = list(shape)
1138 1139
1139 1140 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1140 1141
1141 1142 return shape_tuple, nTxs
1142 1143
1143 1144 def run(self, dataOut, shape=None, nTxs=None):
1144 1145
1145 1146 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1146 1147
1147 1148 dataOut.flagNoData = True
1148 1149 profileIndex = None
1149 1150
1150 1151 if dataOut.flagDataAsBlock:
1151 1152
1152 1153 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1153 1154 dataOut.flagNoData = False
1154 1155
1155 1156 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1156 1157
1157 1158 else:
1158 1159
1159 1160 if self.__nTxs < 1:
1160 1161
1161 1162 self.__appendProfile(dataOut, self.__nTxs)
1162 1163 new_data = self.__getBuffer()
1163 1164
1164 1165 if new_data is not None:
1165 1166 dataOut.data = new_data
1166 1167 dataOut.flagNoData = False
1167 1168
1168 1169 profileIndex = dataOut.profileIndex*nTxs
1169 1170
1170 1171 else:
1171 1172 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1172 1173
1173 1174 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1174 1175
1175 1176 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1176 1177
1177 1178 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1178 1179
1179 1180 dataOut.profileIndex = profileIndex
1180 1181
1181 1182 dataOut.ippSeconds /= self.__nTxs
1182 1183
1183 1184 return dataOut
1184 1185
1185 1186 class SplitProfiles(Operation):
1186 1187
1187 1188 def __init__(self, **kwargs):
1188 1189
1189 1190 Operation.__init__(self, **kwargs)
1190 1191
1191 1192 def run(self, dataOut, n):
1192 1193
1193 1194 dataOut.flagNoData = True
1194 1195 profileIndex = None
1195 1196
1196 1197 if dataOut.flagDataAsBlock:
1197 1198
1198 1199 #nchannels, nprofiles, nsamples
1199 1200 shape = dataOut.data.shape
1200 1201
1201 1202 if shape[2] % n != 0:
1202 1203 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1203 1204
1204 1205 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1205 1206
1206 1207 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1207 1208 dataOut.flagNoData = False
1208 1209
1209 1210 profileIndex = int(dataOut.nProfiles/n) - 1
1210 1211
1211 1212 else:
1212 1213
1213 1214 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1214 1215
1215 1216 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1216 1217
1217 1218 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1218 1219
1219 1220 dataOut.nProfiles = int(dataOut.nProfiles*n)
1220 1221
1221 1222 dataOut.profileIndex = profileIndex
1222 1223
1223 1224 dataOut.ippSeconds /= n
1224 1225
1225 1226 return dataOut
1226 1227
1227 1228 class CombineProfiles(Operation):
1228 1229 def __init__(self, **kwargs):
1229 1230
1230 1231 Operation.__init__(self, **kwargs)
1231 1232
1232 1233 self.__remData = None
1233 1234 self.__profileIndex = 0
1234 1235
1235 1236 def run(self, dataOut, n):
1236 1237
1237 1238 dataOut.flagNoData = True
1238 1239 profileIndex = None
1239 1240
1240 1241 if dataOut.flagDataAsBlock:
1241 1242
1242 1243 #nchannels, nprofiles, nsamples
1243 1244 shape = dataOut.data.shape
1244 1245 new_shape = shape[0], shape[1]/n, shape[2]*n
1245 1246
1246 1247 if shape[1] % n != 0:
1247 1248 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1248 1249
1249 1250 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1250 1251 dataOut.flagNoData = False
1251 1252
1252 1253 profileIndex = int(dataOut.nProfiles*n) - 1
1253 1254
1254 1255 else:
1255 1256
1256 1257 #nchannels, nsamples
1257 1258 if self.__remData is None:
1258 1259 newData = dataOut.data
1259 1260 else:
1260 1261 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1261 1262
1262 1263 self.__profileIndex += 1
1263 1264
1264 1265 if self.__profileIndex < n:
1265 1266 self.__remData = newData
1266 1267 #continue
1267 1268 return
1268 1269
1269 1270 self.__profileIndex = 0
1270 1271 self.__remData = None
1271 1272
1272 1273 dataOut.data = newData
1273 1274 dataOut.flagNoData = False
1274 1275
1275 1276 profileIndex = dataOut.profileIndex/n
1276 1277
1277 1278
1278 1279 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1279 1280
1280 1281 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1281 1282
1282 1283 dataOut.nProfiles = int(dataOut.nProfiles/n)
1283 1284
1284 1285 dataOut.profileIndex = profileIndex
1285 1286
1286 1287 dataOut.ippSeconds *= n
1287 1288
1288 1289 return dataOut
1289 1290
1290 1291 class PulsePair(Operation):
1291 1292 '''
1292 1293 Function PulsePair(Signal Power, Velocity)
1293 1294 The real component of Lag[0] provides Intensity Information
1294 1295 The imag component of Lag[1] Phase provides Velocity Information
1295 1296
1296 1297 Configuration Parameters:
1297 1298 nPRF = Number of Several PRF
1298 1299 theta = Degree Azimuth angel Boundaries
1299 1300
1300 1301 Input:
1301 1302 self.dataOut
1302 1303 lag[N]
1303 1304 Affected:
1304 1305 self.dataOut.spc
1305 1306 '''
1306 1307 isConfig = False
1307 1308 __profIndex = 0
1308 1309 __initime = None
1309 1310 __lastdatatime = None
1310 1311 __buffer = None
1311 1312 noise = None
1312 1313 __dataReady = False
1313 1314 n = None
1314 1315 __nch = 0
1315 1316 __nHeis = 0
1316 1317 removeDC = False
1317 1318 ipp = None
1318 1319 lambda_ = 0
1319 1320
1320 1321 def __init__(self,**kwargs):
1321 1322 Operation.__init__(self,**kwargs)
1322 1323
1323 1324 def setup(self, dataOut, n = None, removeDC=False):
1324 1325 '''
1325 1326 n= Numero de PRF's de entrada
1326 1327 '''
1327 1328 self.__initime = None
1328 1329 ####print("[INICIO]-setup del METODO PULSE PAIR")
1329 1330 self.__lastdatatime = 0
1330 1331 self.__dataReady = False
1331 1332 self.__buffer = 0
1332 1333 self.__profIndex = 0
1333 1334 self.noise = None
1334 1335 self.__nch = dataOut.nChannels
1335 1336 self.__nHeis = dataOut.nHeights
1336 1337 self.removeDC = removeDC
1337 1338 self.lambda_ = 3.0e8/(9345.0e6)
1338 1339 self.ippSec = dataOut.ippSeconds
1339 1340 self.nCohInt = dataOut.nCohInt
1340 1341 ####print("IPPseconds",dataOut.ippSeconds)
1341 1342 ####print("ELVALOR DE n es:", n)
1342 1343 if n == None:
1343 1344 raise ValueError("n should be specified.")
1344 1345
1345 1346 if n != None:
1346 1347 if n<2:
1347 1348 raise ValueError("n should be greater than 2")
1348 1349
1349 1350 self.n = n
1350 1351 self.__nProf = n
1351 1352
1352 1353 self.__buffer = numpy.zeros((dataOut.nChannels,
1353 1354 n,
1354 1355 dataOut.nHeights),
1355 1356 dtype='complex')
1356 1357
1357 1358 def putData(self,data):
1358 1359 '''
1359 1360 Add a profile to he __buffer and increase in one the __profiel Index
1360 1361 '''
1361 1362 self.__buffer[:,self.__profIndex,:]= data
1362 1363 self.__profIndex += 1
1363 1364 return
1364 1365
1365 1366 def pushData(self,dataOut):
1366 1367 '''
1367 1368 Return the PULSEPAIR and the profiles used in the operation
1368 1369 Affected : self.__profileIndex
1369 1370 '''
1370 1371 #----------------- Remove DC-----------------------------------
1371 1372 if self.removeDC==True:
1372 1373 mean = numpy.mean(self.__buffer,1)
1373 1374 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1374 1375 dc= numpy.tile(tmp,[1,self.__nProf,1])
1375 1376 self.__buffer = self.__buffer - dc
1376 1377 #------------------Calculo de Potencia ------------------------
1377 1378 pair0 = self.__buffer*numpy.conj(self.__buffer)
1378 1379 pair0 = pair0.real
1379 1380 lag_0 = numpy.sum(pair0,1)
1380 1381 #-----------------Calculo de Cscp------------------------------ New
1381 1382 cspc_pair01 = self.__buffer[0]*self.__buffer[1]
1382 1383 #------------------Calculo de Ruido x canal--------------------
1383 1384 self.noise = numpy.zeros(self.__nch)
1384 1385 for i in range(self.__nch):
1385 1386 daux = numpy.sort(pair0[i,:,:],axis= None)
1386 1387 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1387 1388
1388 1389 self.noise = self.noise.reshape(self.__nch,1)
1389 1390 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1390 1391 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1391 1392 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1392 1393 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1393 1394 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1394 1395 #-------------------- Power --------------------------------------------------
1395 1396 data_power = lag_0/(self.n*self.nCohInt)
1396 1397 #--------------------CCF------------------------------------------------------
1397 1398 data_ccf =numpy.sum(cspc_pair01,axis=0)/(self.n*self.nCohInt)
1398 1399 #------------------ Senal --------------------------------------------------
1399 1400 data_intensity = pair0 - noise_buffer
1400 1401 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1401 1402 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1402 1403 for i in range(self.__nch):
1403 1404 for j in range(self.__nHeis):
1404 1405 if data_intensity[i][j] < 0:
1405 1406 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1406 1407
1407 1408 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1408 1409 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1409 1410 lag_1 = numpy.sum(pair1,1)
1410 1411 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1411 1412 data_velocity = (self.lambda_/2.0)*data_freq
1412 1413
1413 1414 #---------------- Potencia promedio estimada de la Senal-----------
1414 1415 lag_0 = lag_0/self.n
1415 1416 S = lag_0-self.noise
1416 1417
1417 1418 #---------------- Frecuencia Doppler promedio ---------------------
1418 1419 lag_1 = lag_1/(self.n-1)
1419 1420 R1 = numpy.abs(lag_1)
1420 1421
1421 1422 #---------------- Calculo del SNR----------------------------------
1422 1423 data_snrPP = S/self.noise
1423 1424 for i in range(self.__nch):
1424 1425 for j in range(self.__nHeis):
1425 1426 if data_snrPP[i][j] < 1.e-20:
1426 1427 data_snrPP[i][j] = 1.e-20
1427 1428
1428 1429 #----------------- Calculo del ancho espectral ----------------------
1429 1430 L = S/R1
1430 1431 L = numpy.where(L<0,1,L)
1431 1432 L = numpy.log(L)
1432 1433 tmp = numpy.sqrt(numpy.absolute(L))
1433 1434 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1434 1435 n = self.__profIndex
1435 1436
1436 1437 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1437 1438 self.__profIndex = 0
1438 1439 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,data_ccf,n
1439 1440
1440 1441
1441 1442 def pulsePairbyProfiles(self,dataOut):
1442 1443
1443 1444 self.__dataReady = False
1444 1445 data_power = None
1445 1446 data_intensity = None
1446 1447 data_velocity = None
1447 1448 data_specwidth = None
1448 1449 data_snrPP = None
1449 1450 data_ccf = None
1450 1451 self.putData(data=dataOut.data)
1451 1452 if self.__profIndex == self.n:
1452 1453 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, n = self.pushData(dataOut=dataOut)
1453 1454 self.__dataReady = True
1454 1455
1455 1456 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf
1456 1457
1457 1458
1458 1459 def pulsePairOp(self, dataOut, datatime= None):
1459 1460
1460 1461 if self.__initime == None:
1461 1462 self.__initime = datatime
1462 1463 data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf = self.pulsePairbyProfiles(dataOut)
1463 1464 self.__lastdatatime = datatime
1464 1465
1465 1466 if data_power is None:
1466 1467 return None, None, None,None,None,None,None
1467 1468
1468 1469 avgdatatime = self.__initime
1469 1470 deltatime = datatime - self.__lastdatatime
1470 1471 self.__initime = datatime
1471 1472
1472 1473 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf, avgdatatime
1473 1474
1474 1475 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1475 1476 #print("hey")
1476 1477 #print(dataOut.data.shape)
1477 1478 #exit(1)
1478 1479 #print(self.__profIndex)
1479 1480 if not self.isConfig:
1480 1481 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1481 1482 self.isConfig = True
1482 1483 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1483 1484 dataOut.flagNoData = True
1484 1485
1485 1486 if self.__dataReady:
1486 1487 ###print("READY ----------------------------------")
1487 1488 dataOut.nCohInt *= self.n
1488 1489 dataOut.dataPP_POW = data_intensity # S
1489 1490 dataOut.dataPP_POWER = data_power # P valor que corresponde a POTENCIA MOMENTO
1490 1491 dataOut.dataPP_DOP = data_velocity
1491 1492 dataOut.dataPP_SNR = data_snrPP
1492 1493 dataOut.dataPP_WIDTH = data_specwidth
1493 1494 dataOut.dataPP_CCF = data_ccf
1494 1495 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1495 1496 dataOut.nProfiles = int(dataOut.nProfiles/n)
1496 1497 dataOut.utctime = avgdatatime
1497 1498 dataOut.flagNoData = False
1498 1499 return dataOut
1499 1500
1500 1501 class PulsePair_vRF(Operation):
1501 1502 '''
1502 1503 Function PulsePair(Signal Power, Velocity)
1503 1504 The real component of Lag[0] provides Intensity Information
1504 1505 The imag component of Lag[1] Phase provides Velocity Information
1505 1506
1506 1507 Configuration Parameters:
1507 1508 nPRF = Number of Several PRF
1508 1509 theta = Degree Azimuth angel Boundaries
1509 1510
1510 1511 Input:
1511 1512 self.dataOut
1512 1513 lag[N]
1513 1514 Affected:
1514 1515 self.dataOut.spc
1515 1516 '''
1516 1517 isConfig = False
1517 1518 __profIndex = 0
1518 1519 __initime = None
1519 1520 __lastdatatime = None
1520 1521 __buffer = None
1521 1522 noise = None
1522 1523 __dataReady = False
1523 1524 n = None
1524 1525 __nch = 0
1525 1526 __nHeis = 0
1526 1527 removeDC = False
1527 1528 ipp = None
1528 1529 lambda_ = 0
1529 1530
1530 1531 def __init__(self,**kwargs):
1531 1532 Operation.__init__(self,**kwargs)
1532 1533
1533 1534 def setup(self, dataOut, n = None, removeDC=False):
1534 1535 '''
1535 1536 n= Numero de PRF's de entrada
1536 1537 '''
1537 1538 self.__initime = None
1538 1539 ####print("[INICIO]-setup del METODO PULSE PAIR")
1539 1540 self.__lastdatatime = 0
1540 1541 self.__dataReady = False
1541 1542 self.__buffer = 0
1542 1543 self.__profIndex = 0
1543 1544 self.noise = None
1544 1545 self.__nch = dataOut.nChannels
1545 1546 self.__nHeis = dataOut.nHeights
1546 1547 self.removeDC = removeDC
1547 1548 self.lambda_ = 3.0e8/(9345.0e6)
1548 1549 self.ippSec = dataOut.ippSeconds
1549 1550 self.nCohInt = dataOut.nCohInt
1550 1551 ####print("IPPseconds",dataOut.ippSeconds)
1551 1552 ####print("ELVALOR DE n es:", n)
1552 1553 if n == None:
1553 1554 raise ValueError("n should be specified.")
1554 1555
1555 1556 if n != None:
1556 1557 if n<2:
1557 1558 raise ValueError("n should be greater than 2")
1558 1559
1559 1560 self.n = n
1560 1561 self.__nProf = n
1561 1562
1562 1563 self.__buffer = numpy.zeros((dataOut.nChannels,
1563 1564 n,
1564 1565 dataOut.nHeights),
1565 1566 dtype='complex')
1566 1567
1567 1568 def putData(self,data):
1568 1569 '''
1569 1570 Add a profile to he __buffer and increase in one the __profiel Index
1570 1571 '''
1571 1572 self.__buffer[:,self.__profIndex,:]= data
1572 1573 self.__profIndex += 1
1573 1574 return
1574 1575
1575 1576 def putDataByBlock(self,data,n):
1576 1577 '''
1577 1578 Add a profile to he __buffer and increase in one the __profiel Index
1578 1579 '''
1579 1580 self.__buffer[:]= data
1580 1581 self.__profIndex = n
1581 1582 return
1582 1583
1583 1584 def pushData(self,dataOut):
1584 1585 '''
1585 1586 Return the PULSEPAIR and the profiles used in the operation
1586 1587 Affected : self.__profileIndex
1587 1588 NOTA:
1588 1589 1.) Calculo de Potencia
1589 1590 PdBm = 10 *log10(10*(I**2 + Q**2)) Unidades dBm
1590 1591 self.__buffer = I + Qj
1591 1592
1592 1593 2.) Data decodificada
1593 1594 Se toma como referencia el factor estimado en jrodata.py y se adiciona
1594 1595 en PulsePair solo pwcode.
1595 1596 if self.flagDecodeData:
1596 1597 pwcode = numpy.sum(self.code[0]**2)
1597 1598 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
1598 1599 3.) hildebrand_sekhon
1599 1600 Se pasa el arreglo de Potencia pair0 que contiene canales perfiles y altura dividiendole entre el
1600 1601 factor pwcode.
1601 1602 4.) data_power
1602 1603 Este parametro esta dividido por los factores: nro. perfiles, nro intCoh y pwcode
1603 1604 5.) lag_0
1604 1605 Este parametro esta dividido por los factores: nro. perfiles, nro intCoh y pwcode
1605 1606 Igual a data_power
1606 1607
1607 1608 '''
1608 1609 #----------------- Remove DC-----------------------------------
1609 1610 if self.removeDC==True:
1610 1611 mean = numpy.mean(self.__buffer,1)
1611 1612 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1612 1613 dc= numpy.tile(tmp,[1,self.__nProf,1])
1613 1614 self.__buffer = self.__buffer - dc
1614 1615 #------------------Calculo de Potencia ------------------------
1615 1616 pair0 = self.__buffer*numpy.conj(self.__buffer) * 10.0
1616 1617 pair0 = pair0.real
1617 1618 lag_0 = numpy.sum(pair0,1)
1618 1619 #-----------------Calculo de Cscp------------------------------ New
1619 1620 cspc_pair01 = self.__buffer[0]*self.__buffer[1]
1620 1621 #------------------ Data Decodificada------------------------
1621 1622 pwcode = 1
1622 1623 if dataOut.flagDecodeData == True:
1623 1624 pwcode = numpy.sum(dataOut.code[0]**2)
1624 1625 #------------------Calculo de Ruido x canal--------------------
1625 1626 self.noise = numpy.zeros(self.__nch)
1626 1627 for i in range(self.__nch):
1627 1628 daux = numpy.sort(pair0[i,:,:],axis= None)
1628 1629 self.noise[i]=hildebrand_sekhon( daux/pwcode ,self.nCohInt)
1629 1630
1630 1631 self.noise = self.noise.reshape(self.__nch,1)
1631 1632 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1632 1633 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1633 1634 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1634 1635 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1635 1636 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1636 1637 #-------------------- Power --------------------------------------------------
1637 1638 data_power = lag_0/(self.n*self.nCohInt*pwcode)
1638 1639 #--------------------CCF------------------------------------------------------
1639 1640 data_ccf =numpy.sum(cspc_pair01,axis=0)/(self.n*self.nCohInt)
1640 1641 #------------------ Senal --------------------------------------------------
1641 1642 data_intensity = pair0 - noise_buffer
1642 1643 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1643 1644 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1644 1645 for i in range(self.__nch):
1645 1646 for j in range(self.__nHeis):
1646 1647 if data_intensity[i][j] < 0:
1647 1648 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1648 1649
1649 1650 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1650 1651 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1651 1652 lag_1 = numpy.sum(pair1,1)
1652 1653 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1653 1654 data_velocity = (self.lambda_/2.0)*data_freq
1654 1655
1655 1656 #---------------- Potencia promedio estimada de la Senal-----------
1656 1657 lag_0 = data_power
1657 1658 S = lag_0-self.noise
1658 1659
1659 1660 #---------------- Frecuencia Doppler promedio ---------------------
1660 1661 lag_1 = lag_1/(self.n-1)
1661 1662 R1 = numpy.abs(lag_1)
1662 1663
1663 1664 #---------------- Calculo del SNR----------------------------------
1664 1665 data_snrPP = S/self.noise
1665 1666 for i in range(self.__nch):
1666 1667 for j in range(self.__nHeis):
1667 1668 if data_snrPP[i][j] < 1.e-20:
1668 1669 data_snrPP[i][j] = 1.e-20
1669 1670
1670 1671 #----------------- Calculo del ancho espectral ----------------------
1671 1672 L = S/R1
1672 1673 L = numpy.where(L<0,1,L)
1673 1674 L = numpy.log(L)
1674 1675 tmp = numpy.sqrt(numpy.absolute(L))
1675 1676 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1676 1677 n = self.__profIndex
1677 1678
1678 1679 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1679 1680 self.__profIndex = 0
1680 1681 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,data_ccf,n
1681 1682
1682 1683
1683 1684 def pulsePairbyProfiles(self,dataOut,n):
1684 1685
1685 1686 self.__dataReady = False
1686 1687 data_power = None
1687 1688 data_intensity = None
1688 1689 data_velocity = None
1689 1690 data_specwidth = None
1690 1691 data_snrPP = None
1691 1692 data_ccf = None
1692 1693
1693 1694 if dataOut.flagDataAsBlock:
1694 1695 self.putDataByBlock(data=dataOut.data,n=n)
1695 1696 else:
1696 1697 self.putData(data=dataOut.data)
1697 1698 if self.__profIndex == self.n:
1698 1699 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, n = self.pushData(dataOut=dataOut)
1699 1700 self.__dataReady = True
1700 1701
1701 1702 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf
1702 1703
1703 1704
1704 1705 def pulsePairOp(self, dataOut, n, datatime= None):
1705 1706
1706 1707 if self.__initime == None:
1707 1708 self.__initime = datatime
1708 1709 data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf = self.pulsePairbyProfiles(dataOut,n)
1709 1710 self.__lastdatatime = datatime
1710 1711
1711 1712 if data_power is None:
1712 1713 return None, None, None,None,None,None,None
1713 1714
1714 1715 avgdatatime = self.__initime
1715 1716 deltatime = datatime - self.__lastdatatime
1716 1717 self.__initime = datatime
1717 1718
1718 1719 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf, avgdatatime
1719 1720
1720 1721 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1721 1722
1722 1723 if dataOut.flagDataAsBlock:
1723 1724 n = int(dataOut.nProfiles)
1724 1725 #print("n",n)
1725 1726
1726 1727 if not self.isConfig:
1727 1728 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1728 1729 self.isConfig = True
1729 1730
1730 1731
1731 1732 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, avgdatatime = self.pulsePairOp(dataOut, n, dataOut.utctime)
1732 1733
1733 1734
1734 1735 dataOut.flagNoData = True
1735 1736
1736 1737 if self.__dataReady:
1737 1738 ###print("READY ----------------------------------")
1738 1739 dataOut.nCohInt *= self.n
1739 1740 dataOut.dataPP_POW = data_intensity # S
1740 1741 dataOut.dataPP_POWER = data_power # P valor que corresponde a POTENCIA MOMENTO
1741 1742 dataOut.dataPP_DOP = data_velocity
1742 1743 dataOut.dataPP_SNR = data_snrPP
1743 1744 dataOut.dataPP_WIDTH = data_specwidth
1744 1745 dataOut.dataPP_CCF = data_ccf
1745 1746 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1746 1747 dataOut.nProfiles = int(dataOut.nProfiles/n)
1747 1748 dataOut.utctime = avgdatatime
1748 1749 dataOut.flagNoData = False
1749 1750 return dataOut
1750 1751
1751 1752 # import collections
1752 1753 # from scipy.stats import mode
1753 1754 #
1754 1755 # class Synchronize(Operation):
1755 1756 #
1756 1757 # isConfig = False
1757 1758 # __profIndex = 0
1758 1759 #
1759 1760 # def __init__(self, **kwargs):
1760 1761 #
1761 1762 # Operation.__init__(self, **kwargs)
1762 1763 # # self.isConfig = False
1763 1764 # self.__powBuffer = None
1764 1765 # self.__startIndex = 0
1765 1766 # self.__pulseFound = False
1766 1767 #
1767 1768 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1768 1769 #
1769 1770 # #Read data
1770 1771 #
1771 1772 # powerdB = dataOut.getPower(channel = channel)
1772 1773 # noisedB = dataOut.getNoise(channel = channel)[0]
1773 1774 #
1774 1775 # self.__powBuffer.extend(powerdB.flatten())
1775 1776 #
1776 1777 # dataArray = numpy.array(self.__powBuffer)
1777 1778 #
1778 1779 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1779 1780 #
1780 1781 # maxValue = numpy.nanmax(filteredPower)
1781 1782 #
1782 1783 # if maxValue < noisedB + 10:
1783 1784 # #No se encuentra ningun pulso de transmision
1784 1785 # return None
1785 1786 #
1786 1787 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1787 1788 #
1788 1789 # if len(maxValuesIndex) < 2:
1789 1790 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1790 1791 # return None
1791 1792 #
1792 1793 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1793 1794 #
1794 1795 # #Seleccionar solo valores con un espaciamiento de nSamples
1795 1796 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1796 1797 #
1797 1798 # if len(pulseIndex) < 2:
1798 1799 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1799 1800 # return None
1800 1801 #
1801 1802 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1802 1803 #
1803 1804 # #remover senales que se distancien menos de 10 unidades o muestras
1804 1805 # #(No deberian existir IPP menor a 10 unidades)
1805 1806 #
1806 1807 # realIndex = numpy.where(spacing > 10 )[0]
1807 1808 #
1808 1809 # if len(realIndex) < 2:
1809 1810 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1810 1811 # return None
1811 1812 #
1812 1813 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1813 1814 # realPulseIndex = pulseIndex[realIndex]
1814 1815 #
1815 1816 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1816 1817 #
1817 1818 # print "IPP = %d samples" %period
1818 1819 #
1819 1820 # self.__newNSamples = dataOut.nHeights #int(period)
1820 1821 # self.__startIndex = int(realPulseIndex[0])
1821 1822 #
1822 1823 # return 1
1823 1824 #
1824 1825 #
1825 1826 # def setup(self, nSamples, nChannels, buffer_size = 4):
1826 1827 #
1827 1828 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1828 1829 # maxlen = buffer_size*nSamples)
1829 1830 #
1830 1831 # bufferList = []
1831 1832 #
1832 1833 # for i in range(nChannels):
1833 1834 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1834 1835 # maxlen = buffer_size*nSamples)
1835 1836 #
1836 1837 # bufferList.append(bufferByChannel)
1837 1838 #
1838 1839 # self.__nSamples = nSamples
1839 1840 # self.__nChannels = nChannels
1840 1841 # self.__bufferList = bufferList
1841 1842 #
1842 1843 # def run(self, dataOut, channel = 0):
1843 1844 #
1844 1845 # if not self.isConfig:
1845 1846 # nSamples = dataOut.nHeights
1846 1847 # nChannels = dataOut.nChannels
1847 1848 # self.setup(nSamples, nChannels)
1848 1849 # self.isConfig = True
1849 1850 #
1850 1851 # #Append new data to internal buffer
1851 1852 # for thisChannel in range(self.__nChannels):
1852 1853 # bufferByChannel = self.__bufferList[thisChannel]
1853 1854 # bufferByChannel.extend(dataOut.data[thisChannel])
1854 1855 #
1855 1856 # if self.__pulseFound:
1856 1857 # self.__startIndex -= self.__nSamples
1857 1858 #
1858 1859 # #Finding Tx Pulse
1859 1860 # if not self.__pulseFound:
1860 1861 # indexFound = self.__findTxPulse(dataOut, channel)
1861 1862 #
1862 1863 # if indexFound == None:
1863 1864 # dataOut.flagNoData = True
1864 1865 # return
1865 1866 #
1866 1867 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1867 1868 # self.__pulseFound = True
1868 1869 # self.__startIndex = indexFound
1869 1870 #
1870 1871 # #If pulse was found ...
1871 1872 # for thisChannel in range(self.__nChannels):
1872 1873 # bufferByChannel = self.__bufferList[thisChannel]
1873 1874 # #print self.__startIndex
1874 1875 # x = numpy.array(bufferByChannel)
1875 1876 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1876 1877 #
1877 1878 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1878 1879 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1879 1880 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1880 1881 #
1881 1882 # dataOut.data = self.__arrayBuffer
1882 1883 #
1883 1884 # self.__startIndex += self.__newNSamples
1884 1885 #
1885 1886 # return
General Comments 0
You need to be logged in to leave comments. Login now