##// END OF EJS Templates
schain
RSS
Ultimo avance simulador y pruebas de procesaiento Weather Radar
avaldez -
r1294:575339b32691 v3.0-devel-wrc
parent child
Show More
Side by Side Unified
Context file:
r1294:575339b32691 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1372 +1,1372
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10 import json
11 11
12 12 from schainpy.utils import log
13 13 from .jroheaderIO import SystemHeader, RadarControllerHeader
14 14
15 15
16 16 def getNumpyDtype(dataTypeCode):
17 17
18 18 if dataTypeCode == 0:
19 19 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
20 20 elif dataTypeCode == 1:
21 21 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
22 22 elif dataTypeCode == 2:
23 23 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
24 24 elif dataTypeCode == 3:
25 25 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
26 26 elif dataTypeCode == 4:
27 27 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
28 28 elif dataTypeCode == 5:
29 29 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
30 30 else:
31 31 raise ValueError('dataTypeCode was not defined')
32 32
33 33 return numpyDtype
34 34
35 35
36 36 def getDataTypeCode(numpyDtype):
37 37
38 38 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
39 39 datatype = 0
40 40 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
41 41 datatype = 1
42 42 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
43 43 datatype = 2
44 44 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
45 45 datatype = 3
46 46 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
47 47 datatype = 4
48 48 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
49 49 datatype = 5
50 50 else:
51 51 datatype = None
52 52
53 53 return datatype
54 54
55 55
56 56 def hildebrand_sekhon(data, navg):
57 57 """
58 58 This method is for the objective determination of the noise level in Doppler spectra. This
59 59 implementation technique is based on the fact that the standard deviation of the spectral
60 60 densities is equal to the mean spectral density for white Gaussian noise
61 61
62 62 Inputs:
63 63 Data : heights
64 64 navg : numbers of averages
65 65
66 66 Return:
67 67 mean : noise's level
68 68 """
69 69
70 70 sortdata = numpy.sort(data, axis=None)
71 71 lenOfData = len(sortdata)
72 72 nums_min = lenOfData*0.2
73 73
74 74 if nums_min <= 5:
75 75
76 76 nums_min = 5
77 77
78 78 sump = 0.
79 79 sumq = 0.
80 80
81 81 j = 0
82 82 cont = 1
83 83
84 84 while((cont == 1)and(j < lenOfData)):
85 85
86 86 sump += sortdata[j]
87 87 sumq += sortdata[j]**2
88 88
89 89 if j > nums_min:
90 90 rtest = float(j)/(j-1) + 1.0/navg
91 91 if ((sumq*j) > (rtest*sump**2)):
92 92 j = j - 1
93 93 sump = sump - sortdata[j]
94 94 sumq = sumq - sortdata[j]**2
95 95 cont = 0
96 96
97 97 j += 1
98 98
99 99 lnoise = sump / j
100 100
101 101 return lnoise
102 102
103 103
104 104 class Beam:
105 105
106 106 def __init__(self):
107 107 self.codeList = []
108 108 self.azimuthList = []
109 109 self.zenithList = []
110 110
111 111
112 112 class GenericData(object):
113 113
114 114 flagNoData = True
115 115
116 116 def copy(self, inputObj=None):
117 117
118 118 if inputObj == None:
119 119 return copy.deepcopy(self)
120 120
121 121 for key in list(inputObj.__dict__.keys()):
122 122
123 123 attribute = inputObj.__dict__[key]
124 124
125 125 # If this attribute is a tuple or list
126 126 if type(inputObj.__dict__[key]) in (tuple, list):
127 127 self.__dict__[key] = attribute[:]
128 128 continue
129 129
130 130 # If this attribute is another object or instance
131 131 if hasattr(attribute, '__dict__'):
132 132 self.__dict__[key] = attribute.copy()
133 133 continue
134 134
135 135 self.__dict__[key] = inputObj.__dict__[key]
136 136
137 137 def deepcopy(self):
138 138
139 139 return copy.deepcopy(self)
140 140
141 141 def isEmpty(self):
142 142
143 143 return self.flagNoData
144 144
145 145
146 146 class JROData(GenericData):
147 147
148 148 # m_BasicHeader = BasicHeader()
149 149 # m_ProcessingHeader = ProcessingHeader()
150 150
151 151 systemHeaderObj = SystemHeader()
152 152 radarControllerHeaderObj = RadarControllerHeader()
153 153 # data = None
154 154 type = None
155 155 datatype = None # dtype but in string
156 156 # dtype = None
157 157 # nChannels = None
158 158 # nHeights = None
159 159 nProfiles = None
160 160 heightList = None
161 161 channelList = None
162 162 flagDiscontinuousBlock = False
163 163 useLocalTime = False
164 164 utctime = None
165 165 timeZone = None
166 166 dstFlag = None
167 167 errorCount = None
168 168 blocksize = None
169 169 # nCode = None
170 170 # nBaud = None
171 171 # code = None
172 172 flagDecodeData = False # asumo q la data no esta decodificada
173 173 flagDeflipData = False # asumo q la data no esta sin flip
174 174 flagShiftFFT = False
175 175 # ippSeconds = None
176 176 # timeInterval = None
177 177 nCohInt = None
178 178 # noise = None
179 179 windowOfFilter = 1
180 180 # Speed of ligth
181 181 C = 3e8
182 182 frequency = 49.92e6
183 183 realtime = False
184 184 beacon_heiIndexList = None
185 185 last_block = None
186 186 blocknow = None
187 187 azimuth = None
188 188 zenith = None
189 189 beam = Beam()
190 190 profileIndex = None
191 191 error = None
192 192 data = None
193 193 nmodes = None
194 194
195 195 def __str__(self):
196 196
197 197 return '{} - {}'.format(self.type, self.getDatatime())
198 198
199 199 def getNoise(self):
200 200
201 201 raise NotImplementedError
202 202
203 203 def getNChannels(self):
204 204
205 205 return len(self.channelList)
206 206
207 207 def getChannelIndexList(self):
208 208
209 209 return list(range(self.nChannels))
210 210
211 211 def getNHeights(self):
212 212
213 213 return len(self.heightList)
214 214
215 215 def getHeiRange(self, extrapoints=0):
216 216
217 217 heis = self.heightList
218 218 # deltah = self.heightList[1] - self.heightList[0]
219 219 #
220 220 # heis.append(self.heightList[-1])
221 221
222 222 return heis
223 223
224 224 def getDeltaH(self):
225 225
226 226 delta = self.heightList[1] - self.heightList[0]
227 227
228 228 return delta
229 229
230 230 def getltctime(self):
231 231
232 232 if self.useLocalTime:
233 233 return self.utctime - self.timeZone * 60
234 234
235 235 return self.utctime
236 236
237 237 def getDatatime(self):
238 238
239 239 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
240 240 return datatimeValue
241 241
242 242 def getTimeRange(self):
243 243
244 244 datatime = []
245 245
246 246 datatime.append(self.ltctime)
247 247 datatime.append(self.ltctime + self.timeInterval + 1)
248 248
249 249 datatime = numpy.array(datatime)
250 250
251 251 return datatime
252 252
253 253 def getFmaxTimeResponse(self):
254 254
255 255 period = (10**-6) * self.getDeltaH() / (0.15)
256 256
257 257 PRF = 1. / (period * self.nCohInt)
258 258
259 259 fmax = PRF
260 260
261 261 return fmax
262 262
263 263 def getFmax(self):
264 264 PRF = 1. / (self.ippSeconds * self.nCohInt)
265 265
266 266 fmax = PRF
267 267 return fmax
268 268
269 269 def getVmax(self):
270
270 #print("frequency",self.frequency)
271 271 _lambda = self.C / self.frequency
272 272
273 273 vmax = self.getFmax() * _lambda / 2
274 274
275 275 return vmax
276 276
277 277 def get_ippSeconds(self):
278 278 '''
279 279 '''
280 280 return self.radarControllerHeaderObj.ippSeconds
281 281
282 282 def set_ippSeconds(self, ippSeconds):
283 283 '''
284 284 '''
285 285
286 286 self.radarControllerHeaderObj.ippSeconds = ippSeconds
287 287
288 288 return
289 289
290 290 def get_dtype(self):
291 291 '''
292 292 '''
293 293 return getNumpyDtype(self.datatype)
294 294
295 295 def set_dtype(self, numpyDtype):
296 296 '''
297 297 '''
298 298
299 299 self.datatype = getDataTypeCode(numpyDtype)
300 300
301 301 def get_code(self):
302 302 '''
303 303 '''
304 304 return self.radarControllerHeaderObj.code
305 305
306 306 def set_code(self, code):
307 307 '''
308 308 '''
309 309 self.radarControllerHeaderObj.code = code
310 310
311 311 return
312 312
313 313 def get_ncode(self):
314 314 '''
315 315 '''
316 316 return self.radarControllerHeaderObj.nCode
317 317
318 318 def set_ncode(self, nCode):
319 319 '''
320 320 '''
321 321 self.radarControllerHeaderObj.nCode = nCode
322 322
323 323 return
324 324
325 325 def get_nbaud(self):
326 326 '''
327 327 '''
328 328 return self.radarControllerHeaderObj.nBaud
329 329
330 330 def set_nbaud(self, nBaud):
331 331 '''
332 332 '''
333 333 self.radarControllerHeaderObj.nBaud = nBaud
334 334
335 335 return
336 336
337 337 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
338 338 channelIndexList = property(
339 339 getChannelIndexList, "I'm the 'channelIndexList' property.")
340 340 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
341 341 #noise = property(getNoise, "I'm the 'nHeights' property.")
342 342 datatime = property(getDatatime, "I'm the 'datatime' property")
343 343 ltctime = property(getltctime, "I'm the 'ltctime' property")
344 344 ippSeconds = property(get_ippSeconds, set_ippSeconds)
345 345 dtype = property(get_dtype, set_dtype)
346 346 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
347 347 code = property(get_code, set_code)
348 348 nCode = property(get_ncode, set_ncode)
349 349 nBaud = property(get_nbaud, set_nbaud)
350 350
351 351
352 352 class Voltage(JROData):
353 353
354 354 # data es un numpy array de 2 dmensiones (canales, alturas)
355 355 data = None
356 356
357 357 def __init__(self):
358 358 '''
359 359 Constructor
360 360 '''
361 361
362 362 self.useLocalTime = True
363 363 self.radarControllerHeaderObj = RadarControllerHeader()
364 364 self.systemHeaderObj = SystemHeader()
365 365 self.type = "Voltage"
366 366 self.data = None
367 367 # self.dtype = None
368 368 # self.nChannels = 0
369 369 # self.nHeights = 0
370 370 self.nProfiles = None
371 371 self.heightList = None
372 372 self.channelList = None
373 373 # self.channelIndexList = None
374 374 self.flagNoData = True
375 375 self.flagDiscontinuousBlock = False
376 376 self.utctime = None
377 377 self.timeZone = None
378 378 self.dstFlag = None
379 379 self.errorCount = None
380 380 self.nCohInt = None
381 381 self.blocksize = None
382 382 self.flagDecodeData = False # asumo q la data no esta decodificada
383 383 self.flagDeflipData = False # asumo q la data no esta sin flip
384 384 self.flagShiftFFT = False
385 385 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
386 386 self.profileIndex = 0
387 387
388 388 def getNoisebyHildebrand(self, channel=None):
389 389 """
390 390 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
391 391
392 392 Return:
393 393 noiselevel
394 394 """
395 395
396 396 if channel != None:
397 397 data = self.data[channel]
398 398 nChannels = 1
399 399 else:
400 400 data = self.data
401 401 nChannels = self.nChannels
402 402
403 403 noise = numpy.zeros(nChannels)
404 404 power = data * numpy.conjugate(data)
405 405
406 406 for thisChannel in range(nChannels):
407 407 if nChannels == 1:
408 408 daux = power[:].real
409 409 else:
410 410 daux = power[thisChannel, :].real
411 411 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
412 412
413 413 return noise
414 414
415 415 def getNoise(self, type=1, channel=None):
416 416
417 417 if type == 1:
418 418 noise = self.getNoisebyHildebrand(channel)
419 419
420 420 return noise
421 421
422 422 def getPower(self, channel=None):
423 423
424 424 if channel != None:
425 425 data = self.data[channel]
426 426 else:
427 427 data = self.data
428 428
429 429 power = data * numpy.conjugate(data)
430 430 powerdB = 10 * numpy.log10(power.real)
431 431 powerdB = numpy.squeeze(powerdB)
432 432
433 433 return powerdB
434 434
435 435 def getTimeInterval(self):
436 436
437 437 timeInterval = self.ippSeconds * self.nCohInt
438 438
439 439 return timeInterval
440 440
441 441 noise = property(getNoise, "I'm the 'nHeights' property.")
442 442 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
443 443
444 444
445 445 class Spectra(JROData):
446 446
447 447 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
448 448 data_spc = None
449 449 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
450 450 data_cspc = None
451 451 # data dc es un numpy array de 2 dmensiones (canales, alturas)
452 452 data_dc = None
453 453 # data power
454 454 data_pwr = None
455 455 nFFTPoints = None
456 456 # nPairs = None
457 457 pairsList = None
458 458 nIncohInt = None
459 459 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
460 460 nCohInt = None # se requiere para determinar el valor de timeInterval
461 461 ippFactor = None
462 462 profileIndex = 0
463 463 plotting = "spectra"
464 464
465 465 def __init__(self):
466 466 '''
467 467 Constructor
468 468 '''
469 469
470 470 self.useLocalTime = True
471 471 self.radarControllerHeaderObj = RadarControllerHeader()
472 472 self.systemHeaderObj = SystemHeader()
473 473 self.type = "Spectra"
474 474 # self.data = None
475 475 # self.dtype = None
476 476 # self.nChannels = 0
477 477 # self.nHeights = 0
478 478 self.nProfiles = None
479 479 self.heightList = None
480 480 self.channelList = None
481 481 # self.channelIndexList = None
482 482 self.pairsList = None
483 483 self.flagNoData = True
484 484 self.flagDiscontinuousBlock = False
485 485 self.utctime = None
486 486 self.nCohInt = None
487 487 self.nIncohInt = None
488 488 self.blocksize = None
489 489 self.nFFTPoints = None
490 490 self.wavelength = None
491 491 self.flagDecodeData = False # asumo q la data no esta decodificada
492 492 self.flagDeflipData = False # asumo q la data no esta sin flip
493 493 self.flagShiftFFT = False
494 494 self.ippFactor = 1
495 495 #self.noise = None
496 496 self.beacon_heiIndexList = []
497 497 self.noise_estimation = None
498 498
499 499 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
500 500 """
501 501 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
502 502
503 503 Return:
504 504 noiselevel
505 505 """
506 506
507 507 noise = numpy.zeros(self.nChannels)
508 508
509 509 for channel in range(self.nChannels):
510 510 daux = self.data_spc[channel,
511 511 xmin_index:xmax_index, ymin_index:ymax_index]
512 512 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
513 513
514 514 return noise
515 515
516 516 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
517 517
518 518 if self.noise_estimation is not None:
519 519 # this was estimated by getNoise Operation defined in jroproc_spectra.py
520 520 return self.noise_estimation
521 521 else:
522 522 noise = self.getNoisebyHildebrand(
523 523 xmin_index, xmax_index, ymin_index, ymax_index)
524 524 return noise
525 525
526 526 def getFreqRangeTimeResponse(self, extrapoints=0):
527 527
528 528 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
529 529 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
530 530
531 531 return freqrange
532 532
533 533 def getAcfRange(self, extrapoints=0):
534 534
535 535 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
536 536 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
537 537
538 538 return freqrange
539 539
540 540 def getFreqRange(self, extrapoints=0):
541 541
542 542 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
543 543 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
544 544
545 545 return freqrange
546 546
547 547 def getVelRange(self, extrapoints=0):
548 548
549 549 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
550 550 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
551 551
552 552 if self.nmodes:
553 553 return velrange/self.nmodes
554 554 else:
555 555 return velrange
556 556
557 557 def getNPairs(self):
558 558
559 559 return len(self.pairsList)
560 560
561 561 def getPairsIndexList(self):
562 562
563 563 return list(range(self.nPairs))
564 564
565 565 def getNormFactor(self):
566 566
567 567 pwcode = 1
568 568
569 569 if self.flagDecodeData:
570 570 pwcode = numpy.sum(self.code[0]**2)
571 571 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
572 572 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
573 573
574 574 return normFactor
575 575
576 576 def getFlagCspc(self):
577 577
578 578 if self.data_cspc is None:
579 579 return True
580 580
581 581 return False
582 582
583 583 def getFlagDc(self):
584 584
585 585 if self.data_dc is None:
586 586 return True
587 587
588 588 return False
589 589
590 590 def getTimeInterval(self):
591 591
592 592 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
593 593 if self.nmodes:
594 594 return self.nmodes*timeInterval
595 595 else:
596 596 return timeInterval
597 597
598 598 def getPower(self):
599 599
600 600 factor = self.normFactor
601 601 z = self.data_spc / factor
602 602 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
603 603 avg = numpy.average(z, axis=1)
604 604
605 605 return 10 * numpy.log10(avg)
606 606
607 607 def getCoherence(self, pairsList=None, phase=False):
608 608
609 609 z = []
610 610 if pairsList is None:
611 611 pairsIndexList = self.pairsIndexList
612 612 else:
613 613 pairsIndexList = []
614 614 for pair in pairsList:
615 615 if pair not in self.pairsList:
616 616 raise ValueError("Pair %s is not in dataOut.pairsList" % (
617 617 pair))
618 618 pairsIndexList.append(self.pairsList.index(pair))
619 619 for i in range(len(pairsIndexList)):
620 620 pair = self.pairsList[pairsIndexList[i]]
621 621 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
622 622 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
623 623 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
624 624 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
625 625 if phase:
626 626 data = numpy.arctan2(avgcoherenceComplex.imag,
627 627 avgcoherenceComplex.real) * 180 / numpy.pi
628 628 else:
629 629 data = numpy.abs(avgcoherenceComplex)
630 630
631 631 z.append(data)
632 632
633 633 return numpy.array(z)
634 634
635 635 def setValue(self, value):
636 636
637 637 print("This property should not be initialized")
638 638
639 639 return
640 640
641 641 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
642 642 pairsIndexList = property(
643 643 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
644 644 normFactor = property(getNormFactor, setValue,
645 645 "I'm the 'getNormFactor' property.")
646 646 flag_cspc = property(getFlagCspc, setValue)
647 647 flag_dc = property(getFlagDc, setValue)
648 648 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
649 649 timeInterval = property(getTimeInterval, setValue,
650 650 "I'm the 'timeInterval' property")
651 651
652 652
653 653 class SpectraHeis(Spectra):
654 654
655 655 data_spc = None
656 656 data_cspc = None
657 657 data_dc = None
658 658 nFFTPoints = None
659 659 # nPairs = None
660 660 pairsList = None
661 661 nCohInt = None
662 662 nIncohInt = None
663 663
664 664 def __init__(self):
665 665
666 666 self.radarControllerHeaderObj = RadarControllerHeader()
667 667
668 668 self.systemHeaderObj = SystemHeader()
669 669
670 670 self.type = "SpectraHeis"
671 671
672 672 # self.dtype = None
673 673
674 674 # self.nChannels = 0
675 675
676 676 # self.nHeights = 0
677 677
678 678 self.nProfiles = None
679 679
680 680 self.heightList = None
681 681
682 682 self.channelList = None
683 683
684 684 # self.channelIndexList = None
685 685
686 686 self.flagNoData = True
687 687
688 688 self.flagDiscontinuousBlock = False
689 689
690 690 # self.nPairs = 0
691 691
692 692 self.utctime = None
693 693
694 694 self.blocksize = None
695 695
696 696 self.profileIndex = 0
697 697
698 698 self.nCohInt = 1
699 699
700 700 self.nIncohInt = 1
701 701
702 702 def getNormFactor(self):
703 703 pwcode = 1
704 704 if self.flagDecodeData:
705 705 pwcode = numpy.sum(self.code[0]**2)
706 706
707 707 normFactor = self.nIncohInt * self.nCohInt * pwcode
708 708
709 709 return normFactor
710 710
711 711 def getTimeInterval(self):
712 712
713 713 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
714 714
715 715 return timeInterval
716 716
717 717 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
718 718 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
719 719
720 720
721 721 class Fits(JROData):
722 722
723 723 heightList = None
724 724 channelList = None
725 725 flagNoData = True
726 726 flagDiscontinuousBlock = False
727 727 useLocalTime = False
728 728 utctime = None
729 729 timeZone = None
730 730 # ippSeconds = None
731 731 # timeInterval = None
732 732 nCohInt = None
733 733 nIncohInt = None
734 734 noise = None
735 735 windowOfFilter = 1
736 736 # Speed of ligth
737 737 C = 3e8
738 738 frequency = 49.92e6
739 739 realtime = False
740 740
741 741 def __init__(self):
742 742
743 743 self.type = "Fits"
744 744
745 745 self.nProfiles = None
746 746
747 747 self.heightList = None
748 748
749 749 self.channelList = None
750 750
751 751 # self.channelIndexList = None
752 752
753 753 self.flagNoData = True
754 754
755 755 self.utctime = None
756 756
757 757 self.nCohInt = 1
758 758
759 759 self.nIncohInt = 1
760 760
761 761 self.useLocalTime = True
762 762
763 763 self.profileIndex = 0
764 764
765 765 # self.utctime = None
766 766 # self.timeZone = None
767 767 # self.ltctime = None
768 768 # self.timeInterval = None
769 769 # self.header = None
770 770 # self.data_header = None
771 771 # self.data = None
772 772 # self.datatime = None
773 773 # self.flagNoData = False
774 774 # self.expName = ''
775 775 # self.nChannels = None
776 776 # self.nSamples = None
777 777 # self.dataBlocksPerFile = None
778 778 # self.comments = ''
779 779 #
780 780
781 781 def getltctime(self):
782 782
783 783 if self.useLocalTime:
784 784 return self.utctime - self.timeZone * 60
785 785
786 786 return self.utctime
787 787
788 788 def getDatatime(self):
789 789
790 790 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
791 791 return datatime
792 792
793 793 def getTimeRange(self):
794 794
795 795 datatime = []
796 796
797 797 datatime.append(self.ltctime)
798 798 datatime.append(self.ltctime + self.timeInterval)
799 799
800 800 datatime = numpy.array(datatime)
801 801
802 802 return datatime
803 803
804 804 def getHeiRange(self):
805 805
806 806 heis = self.heightList
807 807
808 808 return heis
809 809
810 810 def getNHeights(self):
811 811
812 812 return len(self.heightList)
813 813
814 814 def getNChannels(self):
815 815
816 816 return len(self.channelList)
817 817
818 818 def getChannelIndexList(self):
819 819
820 820 return list(range(self.nChannels))
821 821
822 822 def getNoise(self, type=1):
823 823
824 824 #noise = numpy.zeros(self.nChannels)
825 825
826 826 if type == 1:
827 827 noise = self.getNoisebyHildebrand()
828 828
829 829 if type == 2:
830 830 noise = self.getNoisebySort()
831 831
832 832 if type == 3:
833 833 noise = self.getNoisebyWindow()
834 834
835 835 return noise
836 836
837 837 def getTimeInterval(self):
838 838
839 839 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
840 840
841 841 return timeInterval
842 842
843 843 def get_ippSeconds(self):
844 844 '''
845 845 '''
846 846 return self.ipp_sec
847 847
848 848
849 849 datatime = property(getDatatime, "I'm the 'datatime' property")
850 850 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
851 851 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
852 852 channelIndexList = property(
853 853 getChannelIndexList, "I'm the 'channelIndexList' property.")
854 854 noise = property(getNoise, "I'm the 'nHeights' property.")
855 855
856 856 ltctime = property(getltctime, "I'm the 'ltctime' property")
857 857 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
858 858 ippSeconds = property(get_ippSeconds, '')
859 859
860 860 class Correlation(JROData):
861 861
862 862 noise = None
863 863 SNR = None
864 864 #--------------------------------------------------
865 865 mode = None
866 866 split = False
867 867 data_cf = None
868 868 lags = None
869 869 lagRange = None
870 870 pairsList = None
871 871 normFactor = None
872 872 #--------------------------------------------------
873 873 # calculateVelocity = None
874 874 nLags = None
875 875 nPairs = None
876 876 nAvg = None
877 877
878 878 def __init__(self):
879 879 '''
880 880 Constructor
881 881 '''
882 882 self.radarControllerHeaderObj = RadarControllerHeader()
883 883
884 884 self.systemHeaderObj = SystemHeader()
885 885
886 886 self.type = "Correlation"
887 887
888 888 self.data = None
889 889
890 890 self.dtype = None
891 891
892 892 self.nProfiles = None
893 893
894 894 self.heightList = None
895 895
896 896 self.channelList = None
897 897
898 898 self.flagNoData = True
899 899
900 900 self.flagDiscontinuousBlock = False
901 901
902 902 self.utctime = None
903 903
904 904 self.timeZone = None
905 905
906 906 self.dstFlag = None
907 907
908 908 self.errorCount = None
909 909
910 910 self.blocksize = None
911 911
912 912 self.flagDecodeData = False # asumo q la data no esta decodificada
913 913
914 914 self.flagDeflipData = False # asumo q la data no esta sin flip
915 915
916 916 self.pairsList = None
917 917
918 918 self.nPoints = None
919 919
920 920 def getPairsList(self):
921 921
922 922 return self.pairsList
923 923
924 924 def getNoise(self, mode=2):
925 925
926 926 indR = numpy.where(self.lagR == 0)[0][0]
927 927 indT = numpy.where(self.lagT == 0)[0][0]
928 928
929 929 jspectra0 = self.data_corr[:, :, indR, :]
930 930 jspectra = copy.copy(jspectra0)
931 931
932 932 num_chan = jspectra.shape[0]
933 933 num_hei = jspectra.shape[2]
934 934
935 935 freq_dc = jspectra.shape[1] / 2
936 936 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
937 937
938 938 if ind_vel[0] < 0:
939 939 ind_vel[list(range(0, 1))] = ind_vel[list(
940 940 range(0, 1))] + self.num_prof
941 941
942 942 if mode == 1:
943 943 jspectra[:, freq_dc, :] = (
944 944 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
945 945
946 946 if mode == 2:
947 947
948 948 vel = numpy.array([-2, -1, 1, 2])
949 949 xx = numpy.zeros([4, 4])
950 950
951 951 for fil in range(4):
952 952 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
953 953
954 954 xx_inv = numpy.linalg.inv(xx)
955 955 xx_aux = xx_inv[0, :]
956 956
957 957 for ich in range(num_chan):
958 958 yy = jspectra[ich, ind_vel, :]
959 959 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
960 960
961 961 junkid = jspectra[ich, freq_dc, :] <= 0
962 962 cjunkid = sum(junkid)
963 963
964 964 if cjunkid.any():
965 965 jspectra[ich, freq_dc, junkid.nonzero()] = (
966 966 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
967 967
968 968 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
969 969
970 970 return noise
971 971
972 972 def getTimeInterval(self):
973 973
974 974 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
975 975
976 976 return timeInterval
977 977
978 978 def splitFunctions(self):
979 979
980 980 pairsList = self.pairsList
981 981 ccf_pairs = []
982 982 acf_pairs = []
983 983 ccf_ind = []
984 984 acf_ind = []
985 985 for l in range(len(pairsList)):
986 986 chan0 = pairsList[l][0]
987 987 chan1 = pairsList[l][1]
988 988
989 989 # Obteniendo pares de Autocorrelacion
990 990 if chan0 == chan1:
991 991 acf_pairs.append(chan0)
992 992 acf_ind.append(l)
993 993 else:
994 994 ccf_pairs.append(pairsList[l])
995 995 ccf_ind.append(l)
996 996
997 997 data_acf = self.data_cf[acf_ind]
998 998 data_ccf = self.data_cf[ccf_ind]
999 999
1000 1000 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1001 1001
1002 1002 def getNormFactor(self):
1003 1003 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1004 1004 acf_pairs = numpy.array(acf_pairs)
1005 1005 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1006 1006
1007 1007 for p in range(self.nPairs):
1008 1008 pair = self.pairsList[p]
1009 1009
1010 1010 ch0 = pair[0]
1011 1011 ch1 = pair[1]
1012 1012
1013 1013 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1014 1014 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1015 1015 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1016 1016
1017 1017 return normFactor
1018 1018
1019 1019 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1020 1020 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1021 1021
1022 1022
1023 1023 class Parameters(Spectra):
1024 1024
1025 1025 experimentInfo = None # Information about the experiment
1026 1026 # Information from previous data
1027 1027 inputUnit = None # Type of data to be processed
1028 1028 operation = None # Type of operation to parametrize
1029 1029 # normFactor = None #Normalization Factor
1030 1030 groupList = None # List of Pairs, Groups, etc
1031 1031 # Parameters
1032 1032 data_param = None # Parameters obtained
1033 1033 data_pre = None # Data Pre Parametrization
1034 1034 data_SNR = None # Signal to Noise Ratio
1035 1035 # heightRange = None #Heights
1036 1036 abscissaList = None # Abscissa, can be velocities, lags or time
1037 1037 # noise = None #Noise Potency
1038 1038 utctimeInit = None # Initial UTC time
1039 1039 paramInterval = None # Time interval to calculate Parameters in seconds
1040 1040 useLocalTime = True
1041 1041 # Fitting
1042 1042 data_error = None # Error of the estimation
1043 1043 constants = None
1044 1044 library = None
1045 1045 # Output signal
1046 1046 outputInterval = None # Time interval to calculate output signal in seconds
1047 1047 data_output = None # Out signal
1048 1048 nAvg = None
1049 1049 noise_estimation = None
1050 1050 GauSPC = None # Fit gaussian SPC
1051 1051
1052 1052 def __init__(self):
1053 1053 '''
1054 1054 Constructor
1055 1055 '''
1056 1056 self.radarControllerHeaderObj = RadarControllerHeader()
1057 1057
1058 1058 self.systemHeaderObj = SystemHeader()
1059 1059
1060 1060 self.type = "Parameters"
1061 1061
1062 1062 def getTimeRange1(self, interval):
1063 1063
1064 1064 datatime = []
1065 1065
1066 1066 if self.useLocalTime:
1067 1067 time1 = self.utctimeInit - self.timeZone * 60
1068 1068 else:
1069 1069 time1 = self.utctimeInit
1070 1070
1071 1071 datatime.append(time1)
1072 1072 datatime.append(time1 + interval)
1073 1073 datatime = numpy.array(datatime)
1074 1074
1075 1075 return datatime
1076 1076
1077 1077 def getTimeInterval(self):
1078 1078
1079 1079 if hasattr(self, 'timeInterval1'):
1080 1080 return self.timeInterval1
1081 1081 else:
1082 1082 return self.paramInterval
1083 1083
1084 1084 def setValue(self, value):
1085 1085
1086 1086 print("This property should not be initialized")
1087 1087
1088 1088 return
1089 1089
1090 1090 def getNoise(self):
1091 1091
1092 1092 return self.spc_noise
1093 1093
1094 1094 timeInterval = property(getTimeInterval)
1095 1095 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1096 1096
1097 1097
1098 1098 class PlotterData(object):
1099 1099 '''
1100 1100 Object to hold data to be plotted
1101 1101 '''
1102 1102
1103 1103 MAXNUMX = 100
1104 1104 MAXNUMY = 100
1105 1105
1106 1106 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1107 1107
1108 1108 self.key = code
1109 1109 self.throttle = throttle_value
1110 1110 self.exp_code = exp_code
1111 1111 self.buffering = buffering
1112 1112 self.ready = False
1113 1113 self.localtime = False
1114 1114 self.data = {}
1115 1115 self.meta = {}
1116 1116 self.__times = []
1117 1117 self.__heights = []
1118 1118
1119 1119 if 'snr' in code:
1120 1120 self.plottypes = ['snr']
1121 1121 elif code == 'spc':
1122 1122 self.plottypes = ['spc', 'noise', 'rti']
1123 1123 elif code == 'rti':
1124 1124 self.plottypes = ['noise', 'rti']
1125 1125 else:
1126 1126 self.plottypes = [code]
1127 1127
1128 1128 if 'snr' not in self.plottypes and snr:
1129 1129 self.plottypes.append('snr')
1130 1130
1131 1131 for plot in self.plottypes:
1132 1132 self.data[plot] = {}
1133 1133
1134 1134 def __str__(self):
1135 1135 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1136 1136 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1137 1137
1138 1138 def __len__(self):
1139 1139 return len(self.__times)
1140 1140
1141 1141 def __getitem__(self, key):
1142 1142
1143 1143 if key not in self.data:
1144 1144 raise KeyError(log.error('Missing key: {}'.format(key)))
1145 1145 if 'spc' in key or not self.buffering:
1146 1146 ret = self.data[key]
1147 1147 elif 'scope' in key:
1148 1148 ret = numpy.array(self.data[key][float(self.tm)])
1149 1149 else:
1150 1150 ret = numpy.array([self.data[key][x] for x in self.times])
1151 1151 if ret.ndim > 1:
1152 1152 ret = numpy.swapaxes(ret, 0, 1)
1153 1153 return ret
1154 1154
1155 1155 def __contains__(self, key):
1156 1156 return key in self.data
1157 1157
1158 1158 def setup(self):
1159 1159 '''
1160 1160 Configure object
1161 1161 '''
1162 1162
1163 1163 self.type = ''
1164 1164 self.ready = False
1165 1165 self.data = {}
1166 1166 self.__times = []
1167 1167 self.__heights = []
1168 1168 self.__all_heights = set()
1169 1169 for plot in self.plottypes:
1170 1170 if 'snr' in plot:
1171 1171 plot = 'snr'
1172 1172 elif 'spc_moments' == plot:
1173 1173 plot = 'moments'
1174 1174 self.data[plot] = {}
1175 1175
1176 1176 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1177 1177 self.data['noise'] = {}
1178 1178 self.data['rti'] = {}
1179 1179 if 'noise' not in self.plottypes:
1180 1180 self.plottypes.append('noise')
1181 1181 if 'rti' not in self.plottypes:
1182 1182 self.plottypes.append('rti')
1183 1183
1184 1184 def shape(self, key):
1185 1185 '''
1186 1186 Get the shape of the one-element data for the given key
1187 1187 '''
1188 1188
1189 1189 if len(self.data[key]):
1190 1190 if 'spc' in key or not self.buffering:
1191 1191 return self.data[key].shape
1192 1192 return self.data[key][self.__times[0]].shape
1193 1193 return (0,)
1194 1194
1195 1195 def update(self, dataOut, tm):
1196 1196 '''
1197 1197 Update data object with new dataOut
1198 1198 '''
1199 1199
1200 1200 if tm in self.__times:
1201 1201 return
1202 1202 self.profileIndex = dataOut.profileIndex
1203 1203 self.tm = tm
1204 1204 self.type = dataOut.type
1205 1205 self.parameters = getattr(dataOut, 'parameters', [])
1206 1206
1207 1207 if hasattr(dataOut, 'meta'):
1208 1208 self.meta.update(dataOut.meta)
1209 1209
1210 1210 self.pairs = dataOut.pairsList
1211 1211 self.interval = dataOut.getTimeInterval()
1212 1212 self.localtime = dataOut.useLocalTime
1213 1213 if 'spc' in self.plottypes or 'cspc' in self.plottypes or 'spc_moments' in self.plottypes:
1214 1214 self.xrange = (dataOut.getFreqRange(1)/1000.,
1215 1215 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1216 1216 self.factor = dataOut.normFactor
1217 1217 self.__heights.append(dataOut.heightList)
1218 1218 self.__all_heights.update(dataOut.heightList)
1219 1219 self.__times.append(tm)
1220 1220
1221 1221 for plot in self.plottypes:
1222 1222 if plot in ('spc', 'spc_moments'):
1223 1223 z = dataOut.data_spc/dataOut.normFactor
1224 1224 buffer = 10*numpy.log10(z)
1225 1225 if plot == 'cspc':
1226 1226 z = dataOut.data_spc/dataOut.normFactor
1227 1227 buffer = (dataOut.data_spc, dataOut.data_cspc)
1228 1228 if plot == 'noise':
1229 1229 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1230 1230 if plot == 'rti':
1231 1231 buffer = dataOut.getPower()
1232 1232 if plot == 'snr_db':
1233 1233 buffer = dataOut.data_SNR
1234 1234 if plot == 'snr':
1235 1235 buffer = 10*numpy.log10(dataOut.data_SNR)
1236 1236 if plot == 'dop':
1237 1237 buffer = dataOut.data_DOP
1238 1238 if plot == 'pow':
1239 1239 buffer = 10*numpy.log10(dataOut.data_POW)
1240 1240 if plot == 'width':
1241 1241 buffer = dataOut.data_WIDTH
1242 1242 if plot == 'coh':
1243 1243 buffer = dataOut.getCoherence()
1244 1244 if plot == 'phase':
1245 1245 buffer = dataOut.getCoherence(phase=True)
1246 1246 if plot == 'output':
1247 1247 buffer = dataOut.data_output
1248 1248 if plot == 'param':
1249 1249 buffer = dataOut.data_param
1250 1250 if plot == 'scope':
1251 1251 buffer = dataOut.data
1252 1252 self.flagDataAsBlock = dataOut.flagDataAsBlock
1253 1253 self.nProfiles = dataOut.nProfiles
1254 1254
1255 1255 if plot == 'spc':
1256 1256 self.data['spc'] = buffer
1257 1257 elif plot == 'cspc':
1258 1258 self.data['spc'] = buffer[0]
1259 1259 self.data['cspc'] = buffer[1]
1260 1260 elif plot == 'spc_moments':
1261 1261 self.data['spc'] = buffer
1262 1262 self.data['moments'][tm] = dataOut.moments
1263 1263 else:
1264 1264 if self.buffering:
1265 1265 self.data[plot][tm] = buffer
1266 1266 else:
1267 1267 self.data[plot] = buffer
1268 1268
1269 1269 if dataOut.channelList is None:
1270 1270 self.channels = range(buffer.shape[0])
1271 1271 else:
1272 1272 self.channels = dataOut.channelList
1273 1273
1274 1274 def normalize_heights(self):
1275 1275 '''
1276 1276 Ensure same-dimension of the data for different heighList
1277 1277 '''
1278 1278
1279 1279 H = numpy.array(list(self.__all_heights))
1280 1280 H.sort()
1281 1281 for key in self.data:
1282 1282 shape = self.shape(key)[:-1] + H.shape
1283 1283 for tm, obj in list(self.data[key].items()):
1284 1284 h = self.__heights[self.__times.index(tm)]
1285 1285 if H.size == h.size:
1286 1286 continue
1287 1287 index = numpy.where(numpy.in1d(H, h))[0]
1288 1288 dummy = numpy.zeros(shape) + numpy.nan
1289 1289 if len(shape) == 2:
1290 1290 dummy[:, index] = obj
1291 1291 else:
1292 1292 dummy[index] = obj
1293 1293 self.data[key][tm] = dummy
1294 1294
1295 1295 self.__heights = [H for tm in self.__times]
1296 1296
1297 1297 def jsonify(self, plot_name, plot_type, decimate=False):
1298 1298 '''
1299 1299 Convert data to json
1300 1300 '''
1301 1301
1302 1302 tm = self.times[-1]
1303 1303 dy = int(self.heights.size/self.MAXNUMY) + 1
1304 1304 if self.key in ('spc', 'cspc') or not self.buffering:
1305 1305 dx = int(self.data[self.key].shape[1]/self.MAXNUMX) + 1
1306 1306 data = self.roundFloats(
1307 1307 self.data[self.key][::, ::dx, ::dy].tolist())
1308 1308 else:
1309 1309 data = self.roundFloats(self.data[self.key][tm].tolist())
1310 1310 if self.key is 'noise':
1311 1311 data = [[x] for x in data]
1312 1312
1313 1313 meta = {}
1314 1314 ret = {
1315 1315 'plot': plot_name,
1316 1316 'code': self.exp_code,
1317 1317 'time': float(tm),
1318 1318 'data': data,
1319 1319 }
1320 1320 meta['type'] = plot_type
1321 1321 meta['interval'] = float(self.interval)
1322 1322 meta['localtime'] = self.localtime
1323 1323 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1324 1324 if 'spc' in self.data or 'cspc' in self.data:
1325 1325 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1326 1326 else:
1327 1327 meta['xrange'] = []
1328 1328
1329 1329 meta.update(self.meta)
1330 1330 ret['metadata'] = meta
1331 1331 return json.dumps(ret)
1332 1332
1333 1333 @property
1334 1334 def times(self):
1335 1335 '''
1336 1336 Return the list of times of the current data
1337 1337 '''
1338 1338
1339 1339 ret = numpy.array(self.__times)
1340 1340 ret.sort()
1341 1341 return ret
1342 1342
1343 1343 @property
1344 1344 def min_time(self):
1345 1345 '''
1346 1346 Return the minimun time value
1347 1347 '''
1348 1348
1349 1349 return self.times[0]
1350 1350
1351 1351 @property
1352 1352 def max_time(self):
1353 1353 '''
1354 1354 Return the maximun time value
1355 1355 '''
1356 1356
1357 1357 return self.times[-1]
1358 1358
1359 1359 @property
1360 1360 def heights(self):
1361 1361 '''
1362 1362 Return the list of heights of the current data
1363 1363 '''
1364 1364
1365 1365 return numpy.array(self.__heights[-1])
1366 1366
1367 1367 @staticmethod
1368 1368 def roundFloats(obj):
1369 1369 if isinstance(obj, list):
1370 1370 return list(map(PlotterData.roundFloats, obj))
1371 1371 elif isinstance(obj, float):
1372 1372 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,295 +1,581
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator #YONG
10 10 from schainpy.utils import log
11 11 from .figure import Figure
12 12
13 13
14 14 @MPDecorator
15 15 class Scope_(Figure):
16 16
17 17 isConfig = None
18 18
19 19 def __init__(self):#, **kwargs): #YONG
20 20 Figure.__init__(self)#, **kwargs)
21 21 self.isConfig = False
22 22 self.WIDTH = 300
23 23 self.HEIGHT = 200
24 24 self.counter_imagwr = 0
25 25
26 26 def getSubplots(self):
27 27
28 28 nrow = self.nplots
29 29 ncol = 3
30 30 return nrow, ncol
31 31
32 32 def setup(self, id, nplots, wintitle, show):
33 33
34 34 self.nplots = nplots
35 35
36 36 self.createFigure(id=id,
37 37 wintitle=wintitle,
38 38 show=show)
39 39
40 40 nrow,ncol = self.getSubplots()
41 41 colspan = 3
42 42 rowspan = 1
43 43
44 44 for i in range(nplots):
45 45 self.addAxes(nrow, ncol, i, 0, colspan, rowspan)
46 46
47 47 def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
48 48 yreal = y[channelIndexList,:].real
49 49 yimag = y[channelIndexList,:].imag
50 50
51 51 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
52 52 xlabel = "Range (Km)"
53 53 ylabel = "Intensity - IQ"
54 54
55 55 if not self.isConfig:
56 56 nplots = len(channelIndexList)
57 57
58 58 self.setup(id=id,
59 59 nplots=nplots,
60 60 wintitle='',
61 61 show=show)
62 62
63 63 if xmin == None: xmin = numpy.nanmin(x)
64 64 if xmax == None: xmax = numpy.nanmax(x)
65 65 if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag))
66 66 if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag))
67 67
68 68 self.isConfig = True
69 69
70 70 self.setWinTitle(title)
71 71
72 72 for i in range(len(self.axesList)):
73 73 title = "Channel %d" %(i)
74 74 axes = self.axesList[i]
75 75
76 76 axes.pline(x, yreal[i,:],
77 77 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
78 78 xlabel=xlabel, ylabel=ylabel, title=title)
79 79
80 80 axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2)
81 81
82 82 def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
83 83 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
84 84 yreal = y.real
85 85
86 86 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
87 87 xlabel = "Range (Km)"
88 88 ylabel = "Intensity"
89 89
90 90 if not self.isConfig:
91 91 nplots = len(channelIndexList)
92 92
93 93 self.setup(id=id,
94 94 nplots=nplots,
95 95 wintitle='',
96 96 show=show)
97 97
98 98 if xmin == None: xmin = numpy.nanmin(x)
99 99 if xmax == None: xmax = numpy.nanmax(x)
100 100 if ymin == None: ymin = numpy.nanmin(yreal)
101 101 if ymax == None: ymax = numpy.nanmax(yreal)
102 102
103
103 104 self.isConfig = True
104 105
105 106 self.setWinTitle(title)
106 107
107 108 for i in range(len(self.axesList)):
108 109 title = "Channel %d" %(i)
109 110 axes = self.axesList[i]
110 111 ychannel = yreal[i,:]
111 112 axes.pline(x, ychannel,
112 113 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
113 114 xlabel=xlabel, ylabel=ylabel, title=title)
114 115
115 116 def plot_weatherpower(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
116 y = y[channelIndexList,:]
117 yreal = y
118 117
118 #x = x[channelIndexList,:]
119 y = y[channelIndexList,:].real
120 y = 10*numpy.log10(y)
119 121 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
120 122 xlabel = "Range (Km)"
121 123 ylabel = "Intensity"
122 124
123 125 if not self.isConfig:
124 126 nplots = len(channelIndexList)
125 127
126 128 self.setup(id=id,
127 129 nplots=nplots,
128 130 wintitle='',
129 131 show=show)
130 132
131 133 if xmin == None: xmin = numpy.nanmin(x)
132 134 if xmax == None: xmax = numpy.nanmax(x)
133 if ymin == None: ymin = numpy.nanmin(yreal)
134 if ymax == None: ymax = numpy.nanmax(yreal)
135 if ymin == None: ymin = numpy.nanmin(y)
136 if ymax == None: ymax = numpy.nanmax(y)
137 #print (xmin,xmax)
135 138
136 139 self.isConfig = True
137 140
138 141 self.setWinTitle(title)
139 142
140 143 for i in range(len(self.axesList)):
141 144 title = "Channel %d" %(i)
142 145 axes = self.axesList[i]
143 ychannel = yreal[i,:]
146 #print(numpy.nanmax(x))
147 ychannel = y[i,:]
148 #ychannel = yreal[i,:]
144 149 axes.pline(x, ychannel,
145 150 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
146 151 xlabel=xlabel, ylabel=ylabel, title=title)
147 152
153 def plot_weathervelocity(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
154 #print(channelIndexList)
155 x = x[channelIndexList,:]
156
157 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
158 xlabel = "Velocity (m/s)"
159 ylabel = "Range (Km)"
160
161 if not self.isConfig:
162 nplots = len(channelIndexList)
163
164 self.setup(id=id,
165 nplots=nplots,
166 wintitle='',
167 show=show)
168
169 if xmin == None: xmin = numpy.nanmin(x)
170 if xmax == None: xmax = numpy.nanmax(x)
171 if ymin == None: ymin = numpy.nanmin(y)
172 if ymax == None: ymax = numpy.nanmax(y)
173 print (xmin,xmax)
174
175 self.isConfig = True
176
177 self.setWinTitle(title)
178
179 for i in range(len(self.axesList)):
180 title = "Channel %d" %(i)
181 axes = self.axesList[i]
182 #print(numpy.nanmax(x))
183 xchannel = x[i,:]
184 #ychannel = yreal[i,:]
185 axes.pline(xchannel, y,
186 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
187 xlabel=xlabel, ylabel=ylabel, title=title)
148 188
149 189
150 190 def run(self, dataOut, id, wintitle="", channelList=None,
151 191 xmin=None, xmax=None, ymin=None, ymax=None, save=False,
152 192 figpath='./', figfile=None, show=True, wr_period=1,
153 193 ftp=False, server=None, folder=None, username=None, password=None, type='power', **kwargs):
154 194
155 195 """
156 196
157 197 Input:
158 198 dataOut :
159 199 id :
160 200 wintitle :
161 201 channelList :
162 202 xmin : None,
163 203 xmax : None,
164 204 ymin : None,
165 205 ymax : None,
166 206 """
167 207 if dataOut.flagNoData:
168 208 return dataOut
169 209
170 210 if channelList == None:
171 211 channelIndexList = dataOut.channelIndexList
172 212 else:
173 213 channelIndexList = []
174 214 for channel in channelList:
175 215 if channel not in dataOut.channelList:
176 216 raise ValueError("Channel %d is not in dataOut.channelList")
177 217 channelIndexList.append(dataOut.channelList.index(channel))
178 218
179 219 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
180 220 #print("***************** PLOTEO **************************")
181 221 #print(dataOut.nProfiles)
182 222 #print(dataOut.heightList.shape)
183 223 #print(dataOut.data.shape)
184 224 if dataOut.flagDataAsBlock:
185 225
186 226 for i in range(dataOut.nProfiles):
187 227
188 228 wintitle1 = wintitle + " [Profile = %d] " %i
189 229
190 230 if type == "power":
191 231 self.plot_power(dataOut.heightList,
192 232 dataOut.data[:,i,:],
193 233 id,
194 234 channelIndexList,
195 235 thisDatetime,
196 236 wintitle1,
197 237 show,
198 238 xmin,
199 239 xmax,
200 240 ymin,
201 241 ymax)
202 242
203 243 if type == "weatherpower":
204 244 self.plot_weatherpower(dataOut.heightList,
205 245 dataOut.data[:,i,:],
206 246 id,
207 247 channelIndexList,
208 248 thisDatetime,
209 wintitle1,
249 wintitle,
210 250 show,
211 251 xmin,
212 252 xmax,
213 253 ymin,
214 254 ymax)
215 255
216 256 if type == "weathervelocity":
217 self.plot_weatherpower(dataOut.heightList,
218 dataOut.data_velocity[:,i,:],
257 self.plot_weathervelocity(dataOut.data_velocity[:,i,:],
258 dataOut.heightList,
219 259 id,
220 260 channelIndexList,
221 261 thisDatetime,
222 262 wintitle1,
223 263 show,
224 264 xmin,
225 265 xmax,
226 266 ymin,
227 267 ymax)
228 268
229 269 if type == "iq":
230 270 self.plot_iq(dataOut.heightList,
231 271 dataOut.data[:,i,:],
232 272 id,
233 273 channelIndexList,
234 274 thisDatetime,
235 275 wintitle1,
236 276 show,
237 277 xmin,
238 278 xmax,
239 279 ymin,
240 280 ymax)
241 281
242 282 self.draw()
243 283
244 284 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
245 285 figfile = self.getFilename(name = str_datetime) + "_" + str(i)
246 286
247 287 self.save(figpath=figpath,
248 288 figfile=figfile,
249 289 save=save,
250 290 ftp=ftp,
251 291 wr_period=wr_period,
252 292 thisDatetime=thisDatetime)
253 293
254 294 else:
255 295 wintitle += " [Profile = %d] " %dataOut.profileIndex
256 296
257 297 if type == "power":
258 298 self.plot_power(dataOut.heightList,
259 299 dataOut.data,
260 300 id,
261 301 channelIndexList,
262 302 thisDatetime,
263 303 wintitle,
264 304 show,
265 305 xmin,
266 306 xmax,
267 307 ymin,
268 308 ymax)
269 309
270 310 if type == "iq":
271 311 self.plot_iq(dataOut.heightList,
272 312 dataOut.data,
273 313 id,
274 314 channelIndexList,
275 315 thisDatetime,
276 316 wintitle,
277 317 show,
278 318 xmin,
279 319 xmax,
280 320 ymin,
281 321 ymax)
282 322
323 if type== "weatherpower":
324 self.plot_weatherpower(dataOut.heightList,
325 dataOut.data,
326 id,
327 channelIndexList,
328 thisDatetime,
329 wintitle,
330 show,
331 xmin,
332 xmax,
333 ymin,
334 ymax)
335 if type== "weathervelocity":
336 self.plot_weathervelocity(dataOut.data_velocity,
337 dataOut.heightList,
338 id,
339 channelIndexList,
340 thisDatetime,
341 wintitle,
342 show,
343 xmin,
344 xmax,
345 ymin,
346 ymax)
347
283 348 self.draw()
284 349
285 350 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + "_" + str(dataOut.profileIndex)
286 351 figfile = self.getFilename(name = str_datetime)
287 352
288 353 self.save(figpath=figpath,
289 354 figfile=figfile,
290 355 save=save,
291 356 ftp=ftp,
292 357 wr_period=wr_period,
293 358 thisDatetime=thisDatetime)
294 359
295 360 return dataOut
361
362
363
364 @MPDecorator
365 class TimePlot_(Figure):
366
367 __isConfig = None
368 __nsubplots = None
369
370 WIDTHPROF = None
371 HEIGHTPROF = None
372 PREFIX = 'time'
373
374 def __init__(self):
375
376 Figure.__init__(self)
377 self.timerange = None
378 self.isConfig = False
379 self.__nsubplots = 1
380
381 self.WIDTH = 800
382 self.HEIGHT = 250
383 self.WIDTHPROF = 120
384 self.HEIGHTPROF = 0
385 self.counter_imagwr = 0
386
387 self.PLOT_CODE = RTIVOLT_CODE
388
389 self.FTP_WEI = None
390 self.EXP_CODE = None
391 self.SUB_EXP_CODE = None
392 self.PLOT_POS = None
393 self.tmin = None
394 self.tmax = None
395
396 self.xmin = None
397 self.xmax = None
398
399 self.figfile = None
400
401 def getSubplots(self):
402
403 ncol = 1
404 nrow = self.nplots
405
406 return nrow, ncol
407
408 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
409
410 self.__showprofile = showprofile
411 self.nplots = nplots
412
413 ncolspan = 1
414 colspan = 1
415 if showprofile:
416 ncolspan = 7
417 colspan = 6
418 self.__nsubplots = 2
419
420 self.createFigure(id = id,
421 wintitle = wintitle,
422 widthplot = self.WIDTH + self.WIDTHPROF,
423 heightplot = self.HEIGHT + self.HEIGHTPROF,
424 show=show)
425
426 nrow, ncol = self.getSubplots()
427
428 counter = 0
429 for y in range(nrow):
430 for x in range(ncol):
431
432 if counter >= self.nplots:
433 break
434
435 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
436
437 if showprofile:
438 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
439
440 counter += 1
441
442 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
443 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,type="intensity",
444 timerange=None, colormap='jet',
445 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
446 server=None, folder=None, username=None, password=None,
447 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, HEIGHT=None):
448
449 """
450
451 Input:
452 dataOut :
453 id :
454 wintitle :
455 channelList :
456 showProfile :
457 xmin : None,
458 xmax : None,
459 ymin : None,
460 ymax : None,
461 zmin : None,
462 zmax : None
463 """
464 print("estoy aqui :D")
465 if dataOut.flagNoData:
466 return dataOut
467
468 #colormap = kwargs.get('colormap', 'jet')
469 if HEIGHT is not None:
470 self.HEIGHT = HEIGHT
471
472 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
473 return
474
475 if channelList == None:
476 channelIndexList = dataOut.channelIndexList
477 else:
478 channelIndexList = []
479 for channel in channelList:
480 if channel not in dataOut.channelList:
481 raise ValueError("Channel %d is not in dataOut.channelList")
482 channelIndexList.append(dataOut.channelList.index(channel))
483
484 if normFactor is None:
485 factor = dataOut.normFactor
486 else:
487 factor = normFactor
488
489 #factor = dataOut.normFactor
490 x = dataOut.getTimeRange()
491 y = dataOut.getHeiRange()
492 if type=="intensity":
493 z = dataOut.data_intensity/factor
494 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
495 avgdB = numpy.average(z, axis=1)
496 avgdB = 10.*numpy.log10(avg)
497 else:
498 z= dataOut.data_velocity
499 avgdB = numpy.average(z, axis=1)
500
501 # avgdB = dataOut.getPower()
502
503
504 thisDatetime = dataOut.datatime
505 #thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
506 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
507 xlabel = ""
508 ylabel = "Range (Km)"
509
510 update_figfile = False
511
512 if self.xmax is not None and dataOut.ltctime >= self.xmax: #yong
513 self.counter_imagwr = wr_period
514 self.isConfig = False
515 update_figfile = True
516
517 if not self.isConfig:
518
519 nplots = len(channelIndexList)
520
521 self.setup(id=id,
522 nplots=nplots,
523 wintitle=wintitle,
524 showprofile=showprofile,
525 show=show)
526
527 if timerange != None:
528 self.timerange = timerange
529
530 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
531
532 noise = dataOut.noise/factor
533 noisedB = 10*numpy.log10(noise)
534
535 if ymin == None: ymin = numpy.nanmin(y)
536 if ymax == None: ymax = numpy.nanmax(y)
537 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
538 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
539
540 self.FTP_WEI = ftp_wei
541 self.EXP_CODE = exp_code
542 self.SUB_EXP_CODE = sub_exp_code
543 self.PLOT_POS = plot_pos
544
545 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
546 self.isConfig = True
547 self.figfile = figfile
548 update_figfile = True
549
550 self.setWinTitle(title)
551
552 for i in range(self.nplots):
553 index = channelIndexList[i]
554 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
555 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
556 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
557 axes = self.axesList[i*self.__nsubplots]
558 zdB = avgdB[index].reshape((1,-1))
559 axes.pcolorbuffer(x, y, zdB,
560 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
561 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
562 ticksize=9, cblabel='', cbsize="1%", colormap=colormap)
563
564 if self.__showprofile:
565 axes = self.axesList[i*self.__nsubplots +1]
566 axes.pline(avgdB[index], y,
567 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
568 xlabel='dB', ylabel='', title='',
569 ytick_visible=False,
570 grid='x')
571
572 self.draw()
573
574 self.save(figpath=figpath,
575 figfile=figfile,
576 save=save,
577 ftp=ftp,
578 wr_period=wr_period,
579 thisDatetime=thisDatetime,
580 update_figfile=update_figfile)
581 return dataOut
Show file before | Show file after | Hide comments
@@ -1,30 +1,31
1 1 '''
2 2 @author: roj-idl71
3 3 '''
4 4 #USED IN jroplot_spectra.py
5 5 RTI_CODE = 0 #Range time intensity (RTI).
6 6 SPEC_CODE = 1 #Spectra (and Cross-spectra) information.
7 7 CROSS_CODE = 2 #Cross-Correlation information.
8 8 COH_CODE = 3 #Coherence map.
9 9 BASE_CODE = 4 #Base lines graphic.
10 10 ROW_CODE = 5 #Row Spectra.
11 11 TOTAL_CODE = 6 #Total Power.
12 12 DRIFT_CODE = 7 #Drifts graphics.
13 13 HEIGHT_CODE = 8 #Height profile.
14 14 PHASE_CODE = 9 #Signal Phase.
15 15
16 16 POWER_CODE = 16
17 17 NOISE_CODE = 17
18 18 BEACON_CODE = 18
19 19
20 20 #USED IN jroplot_parameters.py
21 21 WIND_CODE = 22
22 22 MSKYMAP_CODE = 23
23 23 MPHASE_CODE = 24
24 24
25 25 MOMENTS_CODE = 25
26 26 PARMS_CODE = 26
27 27 SPECFIT_CODE = 27
28 28 EWDRIFT_CODE = 28
29 29
30 30 WPO_CODE = 29 #Weather Intensity - Power
31 RTIVOLT_CODE = 30
Show file before | Show file after | Hide comments
@@ -1,480 +1,496
1 1 import numpy,math,random,time
2 2 import zmq
3 3 import tempfile
4 4 from io import StringIO
5 5 ########## 1 Heredamos JRODatareader
6 6 from schainpy.model.io.jroIO_base import *
7 7 ########## 2 Heredamos las propiedades de ProcessingUnit
8 8 from schainpy.model.proc.jroproc_base import ProcessingUnit,Operation,MPDecorator
9 9 ########## 3 Importaremos las clases BascicHeader, SystemHeader, RadarControlHeader, ProcessingHeader
10 10 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader,SystemHeader,RadarControllerHeader, ProcessingHeader
11 11 ########## 4 Importaremos el objeto Voltge
12 12 from schainpy.model.data.jrodata import Voltage
13 13
14 14 @MPDecorator
15 15 class SimulatorReader(JRODataReader, ProcessingUnit):
16 16 incIntFactor = 1
17 17 nFFTPoints = 0
18 18 FixPP_IncInt = 1
19 19 FixRCP_IPP = 1000
20 20 FixPP_CohInt = 1
21 21 Tau_0 = 250
22 22 AcqH0_0 = 70
23 23 H0 = AcqH0_0
24 24 AcqDH_0 = 1.25
25 25 DH0 = AcqDH_0
26 26 Bauds = 32
27 27 BaudWidth = None
28 28 FixRCP_TXA = 40
29 29 FixRCP_TXB = 70
30 30 fAngle = 2.0*math.pi*(1/16)
31 31 DC_level = 500
32 32 stdev = 8
33 33 Num_Codes = 2
34 34 #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1])
35 35 #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0])
36 36 #Dyn_snCode = numpy.array([Num_Codes,Bauds])
37 37 Dyn_snCode = None
38 38 Samples = 200
39 39 channels = 5
40 40 pulses = None
41 41 Reference = None
42 42 pulse_size = None
43 43 prof_gen = None
44 44 Fdoppler = 100
45 45 Hdoppler = 36
46 Adoppler = 300
47 frequency = 9345
46 48 def __init__(self):
47 49 """
48 50 Inicializador de la clases SimulatorReader para
49 51 generar datos de voltage simulados.
50 52 Input:
51 53 dataOut: Objeto de la clase Voltage.
52 54 Este Objeto sera utilizado apra almacenar
53 55 un perfil de datos cada vez qe se haga psiversho
54 56 un requerimiento (getData)
55 57 """
56 58 ProcessingUnit.__init__(self)
57 59 print(" [ START ] init - Metodo Simulator Reader")
58 60
59 61 self.isConfig = False
60 62 self.basicHeaderObj = BasicHeader(LOCALTIME)
61 63 self.systemHeaderObj = SystemHeader()
62 64 self.radarControllerHeaderObj = RadarControllerHeader()
63 65 self.processingHeaderObj = ProcessingHeader()
64 66 self.profileIndex = 2**32-1
65 67 self.dataOut = Voltage()
66 68 #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1])
67 69 code0 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,1,1,1,-1,1,1,-1,1,-1,-1,-1,1,1,1,-1,1])
68 70 #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0])
69 71 code1 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,-1,-1,1,-1,-1,1,-1,1,1,1,-1,-1,-1,1,-1])
70 72 #self.Dyn_snCode = numpy.array([code0,code1])
71 73 self.Dyn_snCode = None
72 74 print(" [ END ] init - Metodo simulator Reader" )
73 75
74 76
75 77 def __hasNotDataInBuffer(self):
76 78
77 79 if self.profileIndex >= self.processingHeaderObj.profilesPerBlock* self.nTxs:
78 80 if self.nReadBlocks>0:
79 81 tmp = self.dataOut.utctime
80 82 tmp_utc = int(self.dataOut.utctime)
81 83 tmp_milisecond = int((tmp-tmp_utc)*1000)
82 84 self.basicHeaderObj.utc = tmp_utc
83 85 self.basicHeaderObj.miliSecond= tmp_milisecond
84 86 return 1
85 87 return 0
86 88
87 89
88 90 def setNextFile(self):
89 91 """Set the next file to be readed open it and parse de file header"""
90 92
91 93 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
92 print('------------------- [Opening file] ------------------------------')
94 self.nReadFiles=self.nReadFiles+1
95 print('------------------- [Opening file] ------------------------------',self.nReadFiles)
93 96 self.nReadBlocks = 0
94 97
95 98 def __setNewBlock(self):
96 99
97 100 self.setNextFile()
98 101 if self.flagIsNewFile:
99 102 return 1
100 103
101 104 def readNextBlock(self):
102 105 while True:
103 106 self.__setNewBlock()
104 107 if not(self.readBlock()):
105 108 return 0
106 109 self.getBasicHeader()
107 110 break
108 111 if self.verbose:
109 112 print("[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
110 113 self.processingHeaderObj.dataBlocksPerFile,
111 114 self.dataOut.datatime.ctime()) )
112 115 return 1
113 116
114 117 def getFirstHeader(self):
115 118 self.getBasicHeader()
116 119 self.dataOut.processingHeaderObj = self.processingHeaderObj.copy()
117 120 self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
118 121 self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
119 122 #ADD NEW
120 123 self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
121 124 self.dataOut.heightList = numpy.arange(self.processingHeaderObj.nHeights) * self.processingHeaderObj.deltaHeight + self.processingHeaderObj.firstHeight
122 125 self.dataOut.channelList = list(range(self.systemHeaderObj.nChannels))
123 126 self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
124 127 # asumo q la data no esta decodificada
125 128 self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode
126 129 # asumo q la data no esta sin flip
127 130 self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip
128 131 self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft
132 #
133 self.dataOut.frequency = self.frequency
129 134
130 135 def getBasicHeader(self):
131 136
132 137 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
133 138 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
134 139
135 140 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
136 141
137 142 self.dataOut.timeZone = self.basicHeaderObj.timeZone
138 143
139 144 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
140 145
141 146 self.dataOut.errorCount = self.basicHeaderObj.errorCount
142 147
143 148 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
144 149
145 150 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
146 151
147 152 def reshapeData(self):
148 153 if self.nTxs==1:
149 154 return
150 155
151 156 def readBlock(self):
152 157
153 158 self.jro_GenerateBlockOfData(Samples= self.samples,DC_level=self.DC_level,
154 159 stdev=self.stdev,Reference= self.Reference,
155 160 pulses = self.pulses,Num_Codes=self.Num_Codes,
156 161 pulse_size=self.pulse_size,prof_gen=self.profiles,
157 162 H0=self.H0,DH0=self.DH0)
158 163
159 164 self.profileIndex = 0
160 165 self.flagIsNewFile = 0
161 166 self.flagIsNewBlock = 1
162 167 self.nTotalBlocks += 1
163 168 self.nReadBlocks += 1
164 169
165 170 return 1
166 171
167 172
168 173 def getData(self): ### metodo propio de VoltageReader
169 174
170 175 if self.flagNoMoreFiles:
171 176 self.dataOut.flagNodata = True
172 177 self.flagDiscontinuousBlock = 0
173 178 self.flagIsNewBlock = 0
174 179 if self.__hasNotDataInBuffer(): # aqui es verdad
175 180 if not(self.readNextBlock()): # return 1 y por eso el if not salta a getBasic Header
176 181 return 0
177 182 self.getFirstHeader() # atributo
178 183 self.reshapeData() # nTxx1 =1 return , n
179 184
180 185 if not self.getByBlock:
181 186 self.dataOut.flagDataAsBlock = False
182 187 self.dataOut.data = self.datablock[:, self.profileIndex, :]
183 188 self.dataOut.profileIndex = self.profileIndex
184 189 self.profileIndex += 1
185 190 else:
186 191 pass
187 192 self.dataOut.flagNoData = False
188 193 self.getBasicHeader()
189 194 self.dataOut.realtime = self.online
190 195 return self.dataOut.data
191 196
192 197 def set_kwargs(self, **kwargs):
193 198 for key, value in kwargs.items():
194 199 setattr(self, key, value)
195 200
196 201 def set_RCH(self, expType=2, nTx=1,ipp=None, txA=0, txB=0,
197 202 nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None,
198 203 numTaus=0, line6Function=0, line5Function=0, fClock=None,
199 204 prePulseBefore=0, prePulseAfter=0,
200 205 codeType=0, nCode=0, nBaud=0, code=None,
201 206 flip1=0, flip2=0):
202 207
203 208 self.radarControllerHeaderObj.expType = expType
204 209 self.radarControllerHeaderObj.nTx = nTx
205 210 self.radarControllerHeaderObj.ipp = float(ipp)
206 211 self.radarControllerHeaderObj.txA = float(txA)
207 212 self.radarControllerHeaderObj.txB = float(txB)
208 213 self.radarControllerHeaderObj.rangeIPP = ipp
209 214 self.radarControllerHeaderObj.rangeTxA = txA
210 215 self.radarControllerHeaderObj.rangeTxB = txB
211 216
212 217 self.radarControllerHeaderObj.nHeights = int(nHeights)
213 218 self.radarControllerHeaderObj.firstHeight = numpy.array([firstHeight])
214 219 self.radarControllerHeaderObj.deltaHeight = numpy.array([deltaHeight])
215 220 self.radarControllerHeaderObj.samplesWin = numpy.array([nHeights])
216 221
217 222
218 223 self.radarControllerHeaderObj.nWindows = nWindows
219 224 self.radarControllerHeaderObj.numTaus = numTaus
220 225 self.radarControllerHeaderObj.codeType = codeType
221 226 self.radarControllerHeaderObj.line6Function = line6Function
222 227 self.radarControllerHeaderObj.line5Function = line5Function
223 228 self.radarControllerHeaderObj.fclock = fClock
224 229 self.radarControllerHeaderObj.prePulseBefore= prePulseBefore
225 230 self.radarControllerHeaderObj.prePulseAfter = prePulseAfter
226 231
227 232 self.radarControllerHeaderObj.nCode = nCode
228 233 self.radarControllerHeaderObj.nBaud = nBaud
229 234 self.radarControllerHeaderObj.code = code
230 235 self.radarControllerHeaderObj.flip1 = flip1
231 236 self.radarControllerHeaderObj.flip2 = flip2
232 237
233 238 self.radarControllerHeaderObj.code_size = int(numpy.ceil(nBaud / 32.)) * nCode * 4
234 239
235 240 if fClock is None and deltaHeight is not None:
236 241 self.fClock = 0.15 / (deltaHeight * 1e-6)
237 242
238 243 def set_PH(self, dtype=0, blockSize=0, profilesPerBlock=0,
239 244 dataBlocksPerFile=0, nWindows=0, processFlags=0, nCohInt=0,
240 245 nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0,
241 246 deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0,
242 247 code=0, nBaud=None, shif_fft=False, flag_dc=False,
243 248 flag_cspc=False, flag_decode=False, flag_deflip=False):
244 249
245 250 self.processingHeaderObj.profilesPerBlock = profilesPerBlock
246 251 self.processingHeaderObj.dataBlocksPerFile = dataBlocksPerFile
247 252 self.processingHeaderObj.nWindows = nWindows
248 253 self.processingHeaderObj.nCohInt = nCohInt
249 254 self.processingHeaderObj.nIncohInt = nIncohInt
250 255 self.processingHeaderObj.totalSpectra = totalSpectra
251 256 self.processingHeaderObj.nHeights = int(nHeights)
252 257 self.processingHeaderObj.firstHeight = firstHeight
253 258 self.processingHeaderObj.deltaHeight = deltaHeight
254 259 self.processingHeaderObj.samplesWin = nHeights
255 260
256 261 def set_BH(self, utc = 0, miliSecond = 0, timeZone = 0):
257 262 self.basicHeaderObj.utc = utc
258 263 self.basicHeaderObj.miliSecond = miliSecond
259 264 self.basicHeaderObj.timeZone = timeZone
260 265
261 266 def set_SH(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWidth=0):
262 267 self.systemHeaderObj.nSamples = nSamples
263 268 self.systemHeaderObj.nProfiles = nProfiles
264 269 self.systemHeaderObj.nChannels = nChannels
265 270 self.systemHeaderObj.adcResolution = adcResolution
266 271 self.systemHeaderObj.pciDioBusWidth = pciDioBusWidth
267 272
268 def setup(self,incIntFactor= 1, nFFTPoints = 0, FixPP_IncInt=1,FixRCP_IPP=1000,
273 def setup(self,frequency=49.92e6,incIntFactor= 1, nFFTPoints = 0, FixPP_IncInt=1,FixRCP_IPP=1000,
269 274 FixPP_CohInt= 1,Tau_0= 250,AcqH0_0 = 70 ,AcqDH_0=1.25, Bauds= 32,
270 FixRCP_TXA = 40, FixRCP_TXB = 50, fAngle = 2.0*math.pi*(1/16),DC_level= 500,
271 stdev= 8,Num_Codes = 1 , Dyn_snCode = None, samples=200,channels=1,Fdoppler=20,Hdoppler=36,
275 FixRCP_TXA = 40, FixRCP_TXB = 50, fAngle = 2.0*math.pi*(1/16),DC_level= 50,
276 stdev= 8,Num_Codes = 1 , Dyn_snCode = None, samples=200,
277 channels=2,Fdoppler=20,Hdoppler=36,Adoppler=500,
272 278 **kwargs):
273 279
274 280 self.set_kwargs(**kwargs)
275 281 self.nReadBlocks = 0
282 self.nReadFiles = 1
283 print('------------------- [Opening file: ] ------------------------------',self.nReadFiles)
284
276 285 tmp = time.time()
277 286 tmp_utc = int(tmp)
278 287 tmp_milisecond = int((tmp-tmp_utc)*1000)
279 288 print(" SETUP -basicHeaderObj.utc",datetime.datetime.utcfromtimestamp(tmp))
280 289 if Dyn_snCode is None:
281 290 Num_Codes=1
282 291 Bauds =1
283 292
284 293
285 294
286 295 self.set_BH(utc= tmp_utc,miliSecond= tmp_milisecond,timeZone=300 )
287 296
288 297 self.set_RCH( expType=0, nTx=150,ipp=FixRCP_IPP, txA=FixRCP_TXA, txB= FixRCP_TXB,
289 298 nWindows=1 , nHeights=samples, firstHeight=AcqH0_0, deltaHeight=AcqDH_0,
290 299 numTaus=1, line6Function=0, line5Function=0, fClock=None,
291 300 prePulseBefore=0, prePulseAfter=0,
292 301 codeType=14, nCode=Num_Codes, nBaud=32, code=Dyn_snCode,
293 302 flip1=0, flip2=0)
294 303
295 304 self.set_PH(dtype=0, blockSize=0, profilesPerBlock=300,
296 305 dataBlocksPerFile=120, nWindows=1, processFlags=0, nCohInt=1,
297 306 nIncohInt=1, totalSpectra=0, nHeights=samples, firstHeight=AcqH0_0,
298 307 deltaHeight=AcqDH_0, samplesWin=samples, spectraComb=0, nCode=0,
299 308 code=0, nBaud=None, shif_fft=False, flag_dc=False,
300 309 flag_cspc=False, flag_decode=False, flag_deflip=False)
301 310
302 311 self.set_SH(nSamples=samples, nProfiles=300, nChannels=channels)
303 312
313
314 self.frequency = frequency
304 315 self.incIntFactor = incIntFactor
305 316 self.nFFTPoints = nFFTPoints
306 317 self.FixPP_IncInt = FixPP_IncInt
307 318 self.FixRCP_IPP = FixRCP_IPP
308 319 self.FixPP_CohInt = FixPP_CohInt
309 320 self.Tau_0 = Tau_0
310 321 self.AcqH0_0 = AcqH0_0
311 322 self.H0 = AcqH0_0
312 323 self.AcqDH_0 = AcqDH_0
313 324 self.DH0 = AcqDH_0
314 325 self.Bauds = Bauds
315 326 self.FixRCP_TXA = FixRCP_TXA
316 327 self.FixRCP_TXB = FixRCP_TXB
317 328 self.fAngle = fAngle
318 329 self.DC_level = DC_level
319 330 self.stdev = stdev
320 331 self.Num_Codes = Num_Codes
321 332 self.Dyn_snCode = Dyn_snCode
322 333 self.samples = samples
323 334 self.channels = channels
324 335 self.profiles = None
325 336 self.m_nReference = None
326 337 self.Baudwidth = None
327 338 self.Fdoppler = Fdoppler
328 339 self.Hdoppler = Hdoppler
340 self.Adoppler = Adoppler
329 341
330 342 print("IPP ", self.FixRCP_IPP)
331 343 print("Tau_0 ",self.Tau_0)
332 344 print("AcqH0_0",self.AcqH0_0)
333 345 print("samples,window ",self.samples)
334 346 print("AcqDH_0",AcqDH_0)
335 347 print("FixRCP_TXA",self.FixRCP_TXA)
336 348 print("FixRCP_TXB",self.FixRCP_TXB)
337 349 print("Dyn_snCode",Dyn_snCode)
338 350 print("Fdoppler", Fdoppler)
339 351 print("Hdoppler",Hdoppler)
352 print("Vdopplermax",Fdoppler*(3.0e8/self.frequency)/2.0)
340 353
341 354 self.init_acquisition()
342 355 self.pulses,self.pulse_size=self.init_pulse(Num_Codes=self.Num_Codes,Bauds=self.Bauds,BaudWidth=self.BaudWidth,Dyn_snCode=Dyn_snCode)
343 356 print(" [ END ] - SETUP metodo")
344 357 return
345 358
346 359 def run(self,**kwargs): # metodo propio
347 360 if not(self.isConfig):
348 361 self.setup(**kwargs)
349 362 self.isConfig = True
350 363 self.getData()
351 364
352 365 ##################################################################
353 366 ###### Aqui ingresamos las clases y metodos propios del simulador
354 367 ##################################################################
355 368
356 369 #############################################
357 370 ############## INIT_ACQUISITION##############
358 371 #############################################
359 372 def init_acquisition(self):
360 373
361 374 if self.nFFTPoints != 0:
362 375 self.incIntFactor = m_nProfilesperBlock/self.nFFTPoints
363 376 if (self.FixPP_IncInt > self.incIntFactor):
364 377 self.incIntFactor = self.FixPP_IncInt/ self.incIntFactor
365 378 elif(self.FixPP_IncInt< self.incIntFactor):
366 379 print("False alert...")
367 380
368 381 ProfilesperBlock = self.processingHeaderObj.profilesPerBlock
369 382
370 383 self.timeperblock =int(((self.FixRCP_IPP
371 384 *ProfilesperBlock
372 385 *self.FixPP_CohInt
373 386 *self.incIntFactor)
374 387 /150.0)
375 388 *0.9
376 389 +0.5)
377 390 # para cada canal
378 391 self.profiles = ProfilesperBlock*self.FixPP_CohInt
379 392 self.profiles = ProfilesperBlock
380 393 self.Reference = int((self.Tau_0-self.AcqH0_0)/(self.AcqDH_0)+0.5)
381 394 self.BaudWidth = int((self.FixRCP_TXA/self.AcqDH_0)/self.Bauds + 0.5 )
382 395
383 396 if (self.BaudWidth==0):
384 397 self.BaudWidth=1
385 398 #################################################################
386 399 ####################### init_pulse ##############################
387 400 ################################################################
388 401
389 402 def init_pulse(self,Num_Codes=Num_Codes,Bauds=Bauds,BaudWidth=BaudWidth,Dyn_snCode=Dyn_snCode):
390 403
391 404 Num_Codes = Num_Codes
392 405 Bauds = Bauds
393 406 BaudWidth = BaudWidth
394 407 Dyn_snCode = Dyn_snCode
395 408
396 409 if Dyn_snCode:
397 410 print("EXISTE")
398 411 else:
399 412 print("No existe")
400 413
401 414 if Dyn_snCode: # if Bauds:
402 415 pulses = list(range(0,Num_Codes))
403 416 num_codes = Num_Codes
404 417 for i in range(num_codes):
405 418 pulse_size = Bauds*BaudWidth
406 419 pulses[i] = numpy.zeros(pulse_size)
407 420 for j in range(Bauds):
408 421 for k in range(BaudWidth):
409 422 pulses[i][j*BaudWidth+k] = int(Dyn_snCode[i][j]*600)
410 423 else:
411 424 print("sin code")
412 425 pulses = list(range(1))
413 426 pulse_size = int(self.FixRCP_TXB/0.15+0.5)
414 427 pulses[0] = numpy.ones(pulse_size)
415 428 pulses = 600*pulses[0]
416 429
417 430 return pulses,pulse_size
418 431
419 432 #################################################################
420 433 ##################### Generate block data
421 434 ################################################################
422 435
423 436 def jro_GenerateBlockOfData(self,Samples=Samples,DC_level= DC_level,stdev=stdev,
424 437 Reference= Reference,pulses= pulses,
425 438 Num_Codes= Num_Codes,pulse_size=pulse_size,
426 prof_gen= prof_gen,H0 = H0,DH0=DH0,Fdoppler= Fdoppler,Hdoppler=Hdoppler):
439 prof_gen= prof_gen,H0 = H0,DH0=DH0,
440 Adoppler=Adoppler,Fdoppler= Fdoppler,Hdoppler=Hdoppler):
427 441 Samples = Samples
428 442 DC_level = DC_level
429 443 stdev = stdev
430 444 m_nR = Reference
431 445 pulses = pulses
432 446 num_codes = Num_Codes
433 447 ps = pulse_size
434 448 prof_gen = prof_gen
435 449 channels = self.channels
436 450 H0 = H0
437 451 DH0 = DH0
438 452 ippSec = self.radarControllerHeaderObj.ippSeconds
439 453 Fdoppler = self.Fdoppler
440 454 Hdoppler = self.Hdoppler
455 Adoppler = self.Adoppler
441 456
442 457 self.datablock = numpy.zeros([channels,prof_gen,Samples],dtype= numpy.complex64)
443 458 for i in range(channels):
444 459 for k in range(prof_gen):
445 460 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·NOISEΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
446 461 Noise_r = numpy.random.normal(DC_level,stdev,Samples)
447 462 Noise_i = numpy.random.normal(DC_level,stdev,Samples)
448 463 Noise = numpy.zeros(Samples,dtype=complex)
449 464 Noise.real = Noise_r
450 465 Noise.imag = Noise_i
451 466 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·PULSOSΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
452 467 Pulso = numpy.zeros(pulse_size,dtype=complex)
453 468 Pulso.real = pulses[k%num_codes]
454 469 Pulso.imag = pulses[k%num_codes]
455 470 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· PULSES+NOISEΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·
456 471 InBuffer = numpy.zeros(Samples,dtype=complex)
457 472 InBuffer[m_nR:m_nR+ps] = Pulso
458 InBuffer = Noise+ InBuffer
473 InBuffer = InBuffer+Noise
459 474 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· ANGLE Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
460 475 InBuffer.real[m_nR:m_nR+ps] = InBuffer.real[m_nR:m_nR+ps]*(math.cos( self.fAngle)*5)
461 476 InBuffer.imag[m_nR:m_nR+ps] = InBuffer.imag[m_nR:m_nR+ps]*(math.sin( self.fAngle)*5)
462 477 InBuffer=InBuffer
463 478 self.datablock[i][k]= InBuffer
464 479 #plot_cts(InBuffer,H0=H0,DH0=DH0)
465 480 #wave_fft(x=InBuffer,plot_show=True)
466 481 #time.sleep(1)
467 482 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·DOPPLER SIGNAL...............................................
468 time_vec = numpy.linspace(0,(prof_gen-1)*ippSec,int(prof_gen))+self.nReadBlocks*ippSec*prof_gen
483 time_vec = numpy.linspace(0,(prof_gen-1)*ippSec,int(prof_gen))+self.nReadBlocks*ippSec*prof_gen+(self.nReadFiles-1)*ippSec*prof_gen
469 484 fd = Fdoppler #+(600.0/120)*self.nReadBlocks
470 d_signal = 650*numpy.array(numpy.exp(1.0j*2.0*math.pi*fd*time_vec),dtype=numpy.complex64)
485 d_signal = Adoppler*numpy.array(numpy.exp(1.0j*2.0*math.pi*fd*time_vec),dtype=numpy.complex64)
471 486 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· DATABLOCK + DOPPLERΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·...........................
472 487 HD=int(Hdoppler/self.AcqDH_0)
473 self.datablock[0,:,HD]=self.datablock[0,:,HD]+ d_signal # RESULT
488 for i in range(12):
489 self.datablock[:,:,HD+i]=self.datablock[:,:,HD+i]+ d_signal # RESULT
474 490 '''
475 491 a= numpy.zeros(10)
476 492 for i in range(10):
477 493 a[i]=i+self.nReadBlocks+20
478 494 for i in a:
479 495 self.datablock[0,:,int(i)]=self.datablock[0,:,int(i)]+ d_signal # RESULT
480 496 '''
Show file before | Show file after | Hide comments
@@ -1,1266 +1,1266
1 1 import itertools
2 2
3 3 import numpy
4 4
5 5 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
6 6 from schainpy.model.data.jrodata import Spectra
7 7 from schainpy.model.data.jrodata import hildebrand_sekhon
8 8 from schainpy.utils import log
9 9
10 10 @MPDecorator
11 11 class SpectraProc(ProcessingUnit):
12 12
13 13
14 14 def __init__(self):
15 15
16 16 ProcessingUnit.__init__(self)
17 17
18 18 self.buffer = None
19 19 self.firstdatatime = None
20 20 self.profIndex = 0
21 21 self.dataOut = Spectra()
22 22 self.id_min = None
23 23 self.id_max = None
24 24 self.setupReq = False #Agregar a todas las unidades de proc
25 25
26 26 def __updateSpecFromVoltage(self):
27 27
28 28 self.dataOut.timeZone = self.dataIn.timeZone
29 29 self.dataOut.dstFlag = self.dataIn.dstFlag
30 30 self.dataOut.errorCount = self.dataIn.errorCount
31 31 self.dataOut.useLocalTime = self.dataIn.useLocalTime
32 32 try:
33 33 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
34 34 except:
35 35 pass
36 36 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
37 37 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
38 38 self.dataOut.channelList = self.dataIn.channelList
39 39 self.dataOut.heightList = self.dataIn.heightList
40 40 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
41 41
42 42 self.dataOut.nBaud = self.dataIn.nBaud
43 43 self.dataOut.nCode = self.dataIn.nCode
44 44 self.dataOut.code = self.dataIn.code
45 45 self.dataOut.nProfiles = self.dataOut.nFFTPoints
46 46
47 47 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
48 48 self.dataOut.utctime = self.firstdatatime
49 49 # asumo q la data esta decodificada
50 50 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
51 51 # asumo q la data esta sin flip
52 52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
53 53 self.dataOut.flagShiftFFT = False
54 54
55 55 self.dataOut.nCohInt = self.dataIn.nCohInt
56 56 self.dataOut.nIncohInt = 1
57 57
58 58 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59
60 60 self.dataOut.frequency = self.dataIn.frequency
61 61 self.dataOut.realtime = self.dataIn.realtime
62 62
63 63 self.dataOut.azimuth = self.dataIn.azimuth
64 64 self.dataOut.zenith = self.dataIn.zenith
65 65
66 66 self.dataOut.beam.codeList = self.dataIn.beam.codeList
67 67 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
68 68 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
69 69
70 70 def __getFft(self):
71 71 """
72 72 Convierte valores de Voltaje a Spectra
73 73
74 74 Affected:
75 75 self.dataOut.data_spc
76 76 self.dataOut.data_cspc
77 77 self.dataOut.data_dc
78 78 self.dataOut.heightList
79 79 self.profIndex
80 80 self.buffer
81 81 self.dataOut.flagNoData
82 82 """
83 83 fft_volt = numpy.fft.fft(
84 84 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
85 85 fft_volt = fft_volt.astype(numpy.dtype('complex'))
86 86 dc = fft_volt[:, 0, :]
87 87
88 88 # calculo de self-spectra
89 89 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
90 90 spc = fft_volt * numpy.conjugate(fft_volt)
91 91 #print("spcch0",spc[0])
92 92 spc = spc.real
93 93
94 94 blocksize = 0
95 95 blocksize += dc.size
96 96 blocksize += spc.size
97 97
98 98 #print("spc :",spc.shape)
99 99 data_wr = None
100 100 if self.dataOut.flagWR:
101 data_wr = fft_volt
101 data_wr = self.buffer
102 102 blocksize = fft_volt.size
103 103
104 104 cspc = None
105 105 pairIndex = 0
106 106 if self.dataOut.pairsList != None:
107 107 # calculo de cross-spectra
108 108 cspc = numpy.zeros(
109 109 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
110 110 for pair in self.dataOut.pairsList:
111 111 if pair[0] not in self.dataOut.channelList:
112 112 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
113 113 str(pair), str(self.dataOut.channelList)))
114 114 if pair[1] not in self.dataOut.channelList:
115 115 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
116 116 str(pair), str(self.dataOut.channelList)))
117 117
118 118 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
119 119 numpy.conjugate(fft_volt[pair[1], :, :])
120 120 pairIndex += 1
121 121 blocksize += cspc.size
122 122
123 123 self.dataOut.data_spc = spc
124 124 self.dataOut.data_cspc = cspc
125 125 self.dataOut.data_wr = data_wr
126 126 self.dataOut.data_dc = dc
127 127 self.dataOut.blockSize = blocksize
128 128 self.dataOut.flagShiftFFT = False
129 129
130 130 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None, shift_fft=False,flagWR= 0):
131 131
132 132 self.dataOut.flagWR = flagWR
133 133
134 134 if self.dataIn.type == "Spectra":
135 135 self.dataOut.copy(self.dataIn)
136 136
137 137 if shift_fft:
138 138 #desplaza a la derecha en el eje 2 determinadas posiciones
139 139 shift = int(self.dataOut.nFFTPoints/2)
140 140 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
141 141
142 142 if self.dataOut.data_cspc is not None:
143 143 #desplaza a la derecha en el eje 2 determinadas posiciones
144 144 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
145 145
146 146 return True
147 147
148 148 if self.dataIn.type == "Voltage":
149 149 #print("VOLTAGE INPUT SPECTRA")
150 150 self.dataOut.flagNoData = True
151 151
152 152 if nFFTPoints == None:
153 153 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
154 154
155 155 if nProfiles == None:
156 156 nProfiles = nFFTPoints
157 157
158 158 if ippFactor == None:
159 159 ippFactor = 1
160 160
161 161 self.dataOut.ippFactor = ippFactor
162 162
163 163 self.dataOut.nFFTPoints = nFFTPoints
164 164 self.dataOut.pairsList = pairsList
165 165
166 166 if self.buffer is None:
167 167 self.buffer = numpy.zeros((self.dataIn.nChannels,
168 168 nProfiles,
169 169 self.dataIn.nHeights),
170 170 dtype='complex')
171 171 #print("buffer :",self.buffer.shape)
172 172
173 173 if self.dataIn.flagDataAsBlock:
174 174 nVoltProfiles = self.dataIn.data.shape[1]
175 175
176 176 if nVoltProfiles == nProfiles:
177 177 self.buffer = self.dataIn.data.copy()
178 178 self.profIndex = nVoltProfiles
179 179
180 180 elif nVoltProfiles < nProfiles:
181 181
182 182 if self.profIndex == 0:
183 183 self.id_min = 0
184 184 self.id_max = nVoltProfiles
185 185
186 186 self.buffer[:, self.id_min:self.id_max,
187 187 :] = self.dataIn.data
188 188 self.profIndex += nVoltProfiles
189 189 self.id_min += nVoltProfiles
190 190 self.id_max += nVoltProfiles
191 191 else:
192 192 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
193 193 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
194 194 self.dataOut.flagNoData = True
195 195 return 0
196 196 else:
197 197 #print("Spectra ",self.profIndex)
198 198 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
199 199 self.profIndex += 1
200 200
201 201 if self.firstdatatime == None:
202 202 self.firstdatatime = self.dataIn.utctime
203 203
204 204 if self.profIndex == nProfiles:
205 205 self.__updateSpecFromVoltage()
206 206 self.__getFft()
207 207 #print(" DATAOUT SHAPE SPEC",self.dataOut.data_spc.shape)
208 208
209 209 self.dataOut.flagNoData = False
210 210 self.firstdatatime = None
211 211 self.profIndex = 0
212 212
213 213 return True
214 214
215 215 raise ValueError("The type of input object '%s' is not valid" % (
216 216 self.dataIn.type))
217 217
218 218 def __selectPairs(self, pairsList):
219 219
220 220 if not pairsList:
221 221 return
222 222
223 223 pairs = []
224 224 pairsIndex = []
225 225
226 226 for pair in pairsList:
227 227 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
228 228 continue
229 229 pairs.append(pair)
230 230 pairsIndex.append(pairs.index(pair))
231 231
232 232 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
233 233 self.dataOut.pairsList = pairs
234 234
235 235 return
236 236
237 237 def __selectPairsByChannel(self, channelList=None):
238 238
239 239 if channelList == None:
240 240 return
241 241
242 242 pairsIndexListSelected = []
243 243 for pairIndex in self.dataOut.pairsIndexList:
244 244 # First pair
245 245 if self.dataOut.pairsList[pairIndex][0] not in channelList:
246 246 continue
247 247 # Second pair
248 248 if self.dataOut.pairsList[pairIndex][1] not in channelList:
249 249 continue
250 250
251 251 pairsIndexListSelected.append(pairIndex)
252 252
253 253 if not pairsIndexListSelected:
254 254 self.dataOut.data_cspc = None
255 255 self.dataOut.pairsList = []
256 256 return
257 257
258 258 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
259 259 self.dataOut.pairsList = [self.dataOut.pairsList[i]
260 260 for i in pairsIndexListSelected]
261 261
262 262 return
263 263
264 264 def selectChannels(self, channelList):
265 265
266 266 channelIndexList = []
267 267
268 268 for channel in channelList:
269 269 if channel not in self.dataOut.channelList:
270 270 raise ValueError("Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" % (
271 271 channel, str(self.dataOut.channelList)))
272 272
273 273 index = self.dataOut.channelList.index(channel)
274 274 channelIndexList.append(index)
275 275
276 276 self.selectChannelsByIndex(channelIndexList)
277 277
278 278 def selectChannelsByIndex(self, channelIndexList):
279 279 """
280 280 Selecciona un bloque de datos en base a canales segun el channelIndexList
281 281
282 282 Input:
283 283 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
284 284
285 285 Affected:
286 286 self.dataOut.data_spc
287 287 self.dataOut.channelIndexList
288 288 self.dataOut.nChannels
289 289
290 290 Return:
291 291 None
292 292 """
293 293
294 294 for channelIndex in channelIndexList:
295 295 if channelIndex not in self.dataOut.channelIndexList:
296 296 raise ValueError("Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " % (
297 297 channelIndex, self.dataOut.channelIndexList))
298 298
299 299 data_spc = self.dataOut.data_spc[channelIndexList, :]
300 300 data_dc = self.dataOut.data_dc[channelIndexList, :]
301 301
302 302 self.dataOut.data_spc = data_spc
303 303 self.dataOut.data_dc = data_dc
304 304
305 305 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
306 306 self.dataOut.channelList = range(len(channelIndexList))
307 307 self.__selectPairsByChannel(channelIndexList)
308 308
309 309 return 1
310 310
311 311
312 312 def selectFFTs(self, minFFT, maxFFT ):
313 313 """
314 314 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
315 315 minFFT<= FFT <= maxFFT
316 316 """
317 317
318 318 if (minFFT > maxFFT):
319 319 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
320 320
321 321 if (minFFT < self.dataOut.getFreqRange()[0]):
322 322 minFFT = self.dataOut.getFreqRange()[0]
323 323
324 324 if (maxFFT > self.dataOut.getFreqRange()[-1]):
325 325 maxFFT = self.dataOut.getFreqRange()[-1]
326 326
327 327 minIndex = 0
328 328 maxIndex = 0
329 329 FFTs = self.dataOut.getFreqRange()
330 330
331 331 inda = numpy.where(FFTs >= minFFT)
332 332 indb = numpy.where(FFTs <= maxFFT)
333 333
334 334 try:
335 335 minIndex = inda[0][0]
336 336 except:
337 337 minIndex = 0
338 338
339 339 try:
340 340 maxIndex = indb[0][-1]
341 341 except:
342 342 maxIndex = len(FFTs)
343 343
344 344 self.selectFFTsByIndex(minIndex, maxIndex)
345 345
346 346 return 1
347 347
348 348
349 349 def setH0(self, h0, deltaHeight = None):
350 350
351 351 if not deltaHeight:
352 352 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
353 353
354 354 nHeights = self.dataOut.nHeights
355 355
356 356 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
357 357
358 358 self.dataOut.heightList = newHeiRange
359 359
360 360
361 361 def selectHeights(self, minHei, maxHei):
362 362 """
363 363 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
364 364 minHei <= height <= maxHei
365 365
366 366 Input:
367 367 minHei : valor minimo de altura a considerar
368 368 maxHei : valor maximo de altura a considerar
369 369
370 370 Affected:
371 371 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
372 372
373 373 Return:
374 374 1 si el metodo se ejecuto con exito caso contrario devuelve 0
375 375 """
376 376
377 377
378 378 if (minHei > maxHei):
379 379 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei))
380 380
381 381 if (minHei < self.dataOut.heightList[0]):
382 382 minHei = self.dataOut.heightList[0]
383 383
384 384 if (maxHei > self.dataOut.heightList[-1]):
385 385 maxHei = self.dataOut.heightList[-1]
386 386
387 387 minIndex = 0
388 388 maxIndex = 0
389 389 heights = self.dataOut.heightList
390 390
391 391 inda = numpy.where(heights >= minHei)
392 392 indb = numpy.where(heights <= maxHei)
393 393
394 394 try:
395 395 minIndex = inda[0][0]
396 396 except:
397 397 minIndex = 0
398 398
399 399 try:
400 400 maxIndex = indb[0][-1]
401 401 except:
402 402 maxIndex = len(heights)
403 403
404 404 self.selectHeightsByIndex(minIndex, maxIndex)
405 405
406 406
407 407 return 1
408 408
409 409 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
410 410 newheis = numpy.where(
411 411 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
412 412
413 413 if hei_ref != None:
414 414 newheis = numpy.where(self.dataOut.heightList > hei_ref)
415 415
416 416 minIndex = min(newheis[0])
417 417 maxIndex = max(newheis[0])
418 418 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
419 419 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
420 420
421 421 # determina indices
422 422 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
423 423 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
424 424 avg_dB = 10 * \
425 425 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
426 426 beacon_dB = numpy.sort(avg_dB)[-nheis:]
427 427 beacon_heiIndexList = []
428 428 for val in avg_dB.tolist():
429 429 if val >= beacon_dB[0]:
430 430 beacon_heiIndexList.append(avg_dB.tolist().index(val))
431 431
432 432 #data_spc = data_spc[:,:,beacon_heiIndexList]
433 433 data_cspc = None
434 434 if self.dataOut.data_cspc is not None:
435 435 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
436 436 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
437 437
438 438 data_dc = None
439 439 if self.dataOut.data_dc is not None:
440 440 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
441 441 #data_dc = data_dc[:,beacon_heiIndexList]
442 442
443 443 self.dataOut.data_spc = data_spc
444 444 self.dataOut.data_cspc = data_cspc
445 445 self.dataOut.data_dc = data_dc
446 446 self.dataOut.heightList = heightList
447 447 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
448 448
449 449 return 1
450 450
451 451 def selectFFTsByIndex(self, minIndex, maxIndex):
452 452 """
453 453
454 454 """
455 455
456 456 if (minIndex < 0) or (minIndex > maxIndex):
457 457 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
458 458
459 459 if (maxIndex >= self.dataOut.nProfiles):
460 460 maxIndex = self.dataOut.nProfiles-1
461 461
462 462 #Spectra
463 463 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
464 464
465 465 data_cspc = None
466 466 if self.dataOut.data_cspc is not None:
467 467 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
468 468
469 469 data_dc = None
470 470 if self.dataOut.data_dc is not None:
471 471 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
472 472
473 473 self.dataOut.data_spc = data_spc
474 474 self.dataOut.data_cspc = data_cspc
475 475 self.dataOut.data_dc = data_dc
476 476
477 477 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
478 478 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
479 479 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
480 480
481 481 return 1
482 482
483 483
484 484
485 485 def selectHeightsByIndex(self, minIndex, maxIndex):
486 486 """
487 487 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
488 488 minIndex <= index <= maxIndex
489 489
490 490 Input:
491 491 minIndex : valor de indice minimo de altura a considerar
492 492 maxIndex : valor de indice maximo de altura a considerar
493 493
494 494 Affected:
495 495 self.dataOut.data_spc
496 496 self.dataOut.data_cspc
497 497 self.dataOut.data_dc
498 498 self.dataOut.heightList
499 499
500 500 Return:
501 501 1 si el metodo se ejecuto con exito caso contrario devuelve 0
502 502 """
503 503
504 504 if (minIndex < 0) or (minIndex > maxIndex):
505 505 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
506 506 minIndex, maxIndex))
507 507
508 508 if (maxIndex >= self.dataOut.nHeights):
509 509 maxIndex = self.dataOut.nHeights - 1
510 510
511 511 # Spectra
512 512 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
513 513
514 514 data_cspc = None
515 515 if self.dataOut.data_cspc is not None:
516 516 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
517 517
518 518 data_dc = None
519 519 if self.dataOut.data_dc is not None:
520 520 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
521 521
522 522 self.dataOut.data_spc = data_spc
523 523 self.dataOut.data_cspc = data_cspc
524 524 self.dataOut.data_dc = data_dc
525 525
526 526 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
527 527
528 528 return 1
529 529
530 530 def removeDC(self, mode=2):
531 531 jspectra = self.dataOut.data_spc
532 532 jcspectra = self.dataOut.data_cspc
533 533
534 534 num_chan = jspectra.shape[0]
535 535 num_hei = jspectra.shape[2]
536 536
537 537 if jcspectra is not None:
538 538 jcspectraExist = True
539 539 num_pairs = jcspectra.shape[0]
540 540 else:
541 541 jcspectraExist = False
542 542
543 543 freq_dc = int(jspectra.shape[1] / 2)
544 544 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
545 545 ind_vel = ind_vel.astype(int)
546 546
547 547 if ind_vel[0] < 0:
548 548 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
549 549
550 550 if mode == 1:
551 551 jspectra[:, freq_dc, :] = (
552 552 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
553 553
554 554 if jcspectraExist:
555 555 jcspectra[:, freq_dc, :] = (
556 556 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
557 557
558 558 if mode == 2:
559 559
560 560 vel = numpy.array([-2, -1, 1, 2])
561 561 xx = numpy.zeros([4, 4])
562 562
563 563 for fil in range(4):
564 564 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
565 565
566 566 xx_inv = numpy.linalg.inv(xx)
567 567 xx_aux = xx_inv[0, :]
568 568
569 569 for ich in range(num_chan):
570 570 yy = jspectra[ich, ind_vel, :]
571 571 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
572 572
573 573 junkid = jspectra[ich, freq_dc, :] <= 0
574 574 cjunkid = sum(junkid)
575 575
576 576 if cjunkid.any():
577 577 jspectra[ich, freq_dc, junkid.nonzero()] = (
578 578 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
579 579
580 580 if jcspectraExist:
581 581 for ip in range(num_pairs):
582 582 yy = jcspectra[ip, ind_vel, :]
583 583 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
584 584
585 585 self.dataOut.data_spc = jspectra
586 586 self.dataOut.data_cspc = jcspectra
587 587
588 588 return 1
589 589
590 590 def removeInterference2(self):
591 591
592 592 cspc = self.dataOut.data_cspc
593 593 spc = self.dataOut.data_spc
594 594 Heights = numpy.arange(cspc.shape[2])
595 595 realCspc = numpy.abs(cspc)
596 596
597 597 for i in range(cspc.shape[0]):
598 598 LinePower= numpy.sum(realCspc[i], axis=0)
599 599 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
600 600 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
601 601 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
602 602 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
603 603 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
604 604
605 605
606 606 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
607 607 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
608 608 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
609 609 cspc[i,InterferenceRange,:] = numpy.NaN
610 610
611 611
612 612
613 613 self.dataOut.data_cspc = cspc
614 614
615 615 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
616 616
617 617 jspectra = self.dataOut.data_spc
618 618 jcspectra = self.dataOut.data_cspc
619 619 jnoise = self.dataOut.getNoise()
620 620 num_incoh = self.dataOut.nIncohInt
621 621
622 622 num_channel = jspectra.shape[0]
623 623 num_prof = jspectra.shape[1]
624 624 num_hei = jspectra.shape[2]
625 625
626 626 # hei_interf
627 627 if hei_interf is None:
628 628 count_hei = int(num_hei / 2)
629 629 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
630 630 hei_interf = numpy.asarray(hei_interf)[0]
631 631 # nhei_interf
632 632 if (nhei_interf == None):
633 633 nhei_interf = 5
634 634 if (nhei_interf < 1):
635 635 nhei_interf = 1
636 636 if (nhei_interf > count_hei):
637 637 nhei_interf = count_hei
638 638 if (offhei_interf == None):
639 639 offhei_interf = 0
640 640
641 641 ind_hei = list(range(num_hei))
642 642 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
643 643 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
644 644 mask_prof = numpy.asarray(list(range(num_prof)))
645 645 num_mask_prof = mask_prof.size
646 646 comp_mask_prof = [0, num_prof / 2]
647 647
648 648 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
649 649 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
650 650 jnoise = numpy.nan
651 651 noise_exist = jnoise[0] < numpy.Inf
652 652
653 653 # Subrutina de Remocion de la Interferencia
654 654 for ich in range(num_channel):
655 655 # Se ordena los espectros segun su potencia (menor a mayor)
656 656 power = jspectra[ich, mask_prof, :]
657 657 power = power[:, hei_interf]
658 658 power = power.sum(axis=0)
659 659 psort = power.ravel().argsort()
660 660
661 661 # Se estima la interferencia promedio en los Espectros de Potencia empleando
662 662 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
663 663 offhei_interf, nhei_interf + offhei_interf))]]]
664 664
665 665 if noise_exist:
666 666 # tmp_noise = jnoise[ich] / num_prof
667 667 tmp_noise = jnoise[ich]
668 668 junkspc_interf = junkspc_interf - tmp_noise
669 669 #junkspc_interf[:,comp_mask_prof] = 0
670 670
671 671 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
672 672 jspc_interf = jspc_interf.transpose()
673 673 # Calculando el espectro de interferencia promedio
674 674 noiseid = numpy.where(
675 675 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
676 676 noiseid = noiseid[0]
677 677 cnoiseid = noiseid.size
678 678 interfid = numpy.where(
679 679 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
680 680 interfid = interfid[0]
681 681 cinterfid = interfid.size
682 682
683 683 if (cnoiseid > 0):
684 684 jspc_interf[noiseid] = 0
685 685
686 686 # Expandiendo los perfiles a limpiar
687 687 if (cinterfid > 0):
688 688 new_interfid = (
689 689 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
690 690 new_interfid = numpy.asarray(new_interfid)
691 691 new_interfid = {x for x in new_interfid}
692 692 new_interfid = numpy.array(list(new_interfid))
693 693 new_cinterfid = new_interfid.size
694 694 else:
695 695 new_cinterfid = 0
696 696
697 697 for ip in range(new_cinterfid):
698 698 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
699 699 jspc_interf[new_interfid[ip]
700 700 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
701 701
702 702 jspectra[ich, :, ind_hei] = jspectra[ich, :,
703 703 ind_hei] - jspc_interf # Corregir indices
704 704
705 705 # Removiendo la interferencia del punto de mayor interferencia
706 706 ListAux = jspc_interf[mask_prof].tolist()
707 707 maxid = ListAux.index(max(ListAux))
708 708
709 709 if cinterfid > 0:
710 710 for ip in range(cinterfid * (interf == 2) - 1):
711 711 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
712 712 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
713 713 cind = len(ind)
714 714
715 715 if (cind > 0):
716 716 jspectra[ich, interfid[ip], ind] = tmp_noise * \
717 717 (1 + (numpy.random.uniform(cind) - 0.5) /
718 718 numpy.sqrt(num_incoh))
719 719
720 720 ind = numpy.array([-2, -1, 1, 2])
721 721 xx = numpy.zeros([4, 4])
722 722
723 723 for id1 in range(4):
724 724 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
725 725
726 726 xx_inv = numpy.linalg.inv(xx)
727 727 xx = xx_inv[:, 0]
728 728 ind = (ind + maxid + num_mask_prof) % num_mask_prof
729 729 yy = jspectra[ich, mask_prof[ind], :]
730 730 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
731 731 yy.transpose(), xx)
732 732
733 733 indAux = (jspectra[ich, :, :] < tmp_noise *
734 734 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
735 735 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
736 736 (1 - 1 / numpy.sqrt(num_incoh))
737 737
738 738 # Remocion de Interferencia en el Cross Spectra
739 739 if jcspectra is None:
740 740 return jspectra, jcspectra
741 741 num_pairs = int(jcspectra.size / (num_prof * num_hei))
742 742 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
743 743
744 744 for ip in range(num_pairs):
745 745
746 746 #-------------------------------------------
747 747
748 748 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
749 749 cspower = cspower[:, hei_interf]
750 750 cspower = cspower.sum(axis=0)
751 751
752 752 cspsort = cspower.ravel().argsort()
753 753 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
754 754 offhei_interf, nhei_interf + offhei_interf))]]]
755 755 junkcspc_interf = junkcspc_interf.transpose()
756 756 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
757 757
758 758 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
759 759
760 760 median_real = int(numpy.median(numpy.real(
761 761 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
762 762 median_imag = int(numpy.median(numpy.imag(
763 763 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
764 764 comp_mask_prof = [int(e) for e in comp_mask_prof]
765 765 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
766 766 median_real, median_imag)
767 767
768 768 for iprof in range(num_prof):
769 769 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
770 770 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
771 771
772 772 # Removiendo la Interferencia
773 773 jcspectra[ip, :, ind_hei] = jcspectra[ip,
774 774 :, ind_hei] - jcspc_interf
775 775
776 776 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
777 777 maxid = ListAux.index(max(ListAux))
778 778
779 779 ind = numpy.array([-2, -1, 1, 2])
780 780 xx = numpy.zeros([4, 4])
781 781
782 782 for id1 in range(4):
783 783 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
784 784
785 785 xx_inv = numpy.linalg.inv(xx)
786 786 xx = xx_inv[:, 0]
787 787
788 788 ind = (ind + maxid + num_mask_prof) % num_mask_prof
789 789 yy = jcspectra[ip, mask_prof[ind], :]
790 790 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
791 791
792 792 # Guardar Resultados
793 793 self.dataOut.data_spc = jspectra
794 794 self.dataOut.data_cspc = jcspectra
795 795
796 796 return 1
797 797
798 798 def setRadarFrequency(self, frequency=None):
799 799
800 800 if frequency != None:
801 801 self.dataOut.frequency = frequency
802 802
803 803 return 1
804 804
805 805 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
806 806 # validacion de rango
807 807 if minHei == None:
808 808 minHei = self.dataOut.heightList[0]
809 809
810 810 if maxHei == None:
811 811 maxHei = self.dataOut.heightList[-1]
812 812
813 813 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
814 814 print('minHei: %.2f is out of the heights range' % (minHei))
815 815 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
816 816 minHei = self.dataOut.heightList[0]
817 817
818 818 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
819 819 print('maxHei: %.2f is out of the heights range' % (maxHei))
820 820 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
821 821 maxHei = self.dataOut.heightList[-1]
822 822
823 823 # validacion de velocidades
824 824 velrange = self.dataOut.getVelRange(1)
825 825
826 826 if minVel == None:
827 827 minVel = velrange[0]
828 828
829 829 if maxVel == None:
830 830 maxVel = velrange[-1]
831 831
832 832 if (minVel < velrange[0]) or (minVel > maxVel):
833 833 print('minVel: %.2f is out of the velocity range' % (minVel))
834 834 print('minVel is setting to %.2f' % (velrange[0]))
835 835 minVel = velrange[0]
836 836
837 837 if (maxVel > velrange[-1]) or (maxVel < minVel):
838 838 print('maxVel: %.2f is out of the velocity range' % (maxVel))
839 839 print('maxVel is setting to %.2f' % (velrange[-1]))
840 840 maxVel = velrange[-1]
841 841
842 842 # seleccion de indices para rango
843 843 minIndex = 0
844 844 maxIndex = 0
845 845 heights = self.dataOut.heightList
846 846
847 847 inda = numpy.where(heights >= minHei)
848 848 indb = numpy.where(heights <= maxHei)
849 849
850 850 try:
851 851 minIndex = inda[0][0]
852 852 except:
853 853 minIndex = 0
854 854
855 855 try:
856 856 maxIndex = indb[0][-1]
857 857 except:
858 858 maxIndex = len(heights)
859 859
860 860 if (minIndex < 0) or (minIndex > maxIndex):
861 861 raise ValueError("some value in (%d,%d) is not valid" % (
862 862 minIndex, maxIndex))
863 863
864 864 if (maxIndex >= self.dataOut.nHeights):
865 865 maxIndex = self.dataOut.nHeights - 1
866 866
867 867 # seleccion de indices para velocidades
868 868 indminvel = numpy.where(velrange >= minVel)
869 869 indmaxvel = numpy.where(velrange <= maxVel)
870 870 try:
871 871 minIndexVel = indminvel[0][0]
872 872 except:
873 873 minIndexVel = 0
874 874
875 875 try:
876 876 maxIndexVel = indmaxvel[0][-1]
877 877 except:
878 878 maxIndexVel = len(velrange)
879 879
880 880 # seleccion del espectro
881 881 data_spc = self.dataOut.data_spc[:,
882 882 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
883 883 # estimacion de ruido
884 884 noise = numpy.zeros(self.dataOut.nChannels)
885 885
886 886 for channel in range(self.dataOut.nChannels):
887 887 daux = data_spc[channel, :, :]
888 888 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
889 889
890 890 self.dataOut.noise_estimation = noise.copy()
891 891
892 892 return 1
893 893
894 894
895 895 class IncohInt(Operation):
896 896
897 897 __profIndex = 0
898 898 __withOverapping = False
899 899
900 900 __byTime = False
901 901 __initime = None
902 902 __lastdatatime = None
903 903 __integrationtime = None
904 904
905 905 __buffer_spc = None
906 906 __buffer_cspc = None
907 907 __buffer_dc = None
908 908
909 909 __dataReady = False
910 910
911 911 __timeInterval = None
912 912
913 913 n = None
914 914
915 915 def __init__(self):
916 916
917 917 Operation.__init__(self)
918 918
919 919 def setup(self, n=None, timeInterval=None, overlapping=False):
920 920 """
921 921 Set the parameters of the integration class.
922 922
923 923 Inputs:
924 924
925 925 n : Number of coherent integrations
926 926 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
927 927 overlapping :
928 928
929 929 """
930 930
931 931 self.__initime = None
932 932 self.__lastdatatime = 0
933 933
934 934 self.__buffer_spc = 0
935 935 self.__buffer_cspc = 0
936 936 self.__buffer_dc = 0
937 937
938 938 self.__profIndex = 0
939 939 self.__dataReady = False
940 940 self.__byTime = False
941 941
942 942 if n is None and timeInterval is None:
943 943 raise ValueError("n or timeInterval should be specified ...")
944 944
945 945 if n is not None:
946 946 self.n = int(n)
947 947 else:
948 948
949 949 self.__integrationtime = int(timeInterval)
950 950 self.n = None
951 951 self.__byTime = True
952 952
953 953 def putData(self, data_spc, data_cspc, data_dc):
954 954 """
955 955 Add a profile to the __buffer_spc and increase in one the __profileIndex
956 956
957 957 """
958 958 print("profIndex: ",self.__profIndex)
959 959 print("data_spc.shape: ",data_spc.shape)
960 960 print("data_spc.shape: ",data_spc[0,0,:])
961 961
962 962 self.__buffer_spc += data_spc
963 963
964 964 if data_cspc is None:
965 965 self.__buffer_cspc = None
966 966 else:
967 967 self.__buffer_cspc += data_cspc
968 968
969 969 if data_dc is None:
970 970 self.__buffer_dc = None
971 971 else:
972 972 self.__buffer_dc += data_dc
973 973
974 974 self.__profIndex += 1
975 975
976 976 return
977 977
978 978 def pushData(self):
979 979 """
980 980 Return the sum of the last profiles and the profiles used in the sum.
981 981
982 982 Affected:
983 983
984 984 self.__profileIndex
985 985
986 986 """
987 987
988 988 data_spc = self.__buffer_spc
989 989 data_cspc = self.__buffer_cspc
990 990 data_dc = self.__buffer_dc
991 991 n = self.__profIndex
992 992
993 993 self.__buffer_spc = 0
994 994 self.__buffer_cspc = 0
995 995 self.__buffer_dc = 0
996 996 self.__profIndex = 0
997 997
998 998 return data_spc, data_cspc, data_dc, n
999 999
1000 1000 def byProfiles(self, *args):
1001 1001
1002 1002 self.__dataReady = False
1003 1003 avgdata_spc = None
1004 1004 avgdata_cspc = None
1005 1005 avgdata_dc = None
1006 1006
1007 1007 self.putData(*args)
1008 1008
1009 1009 if self.__profIndex == self.n:
1010 1010
1011 1011 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1012 1012 self.n = n
1013 1013 self.__dataReady = True
1014 1014
1015 1015 return avgdata_spc, avgdata_cspc, avgdata_dc
1016 1016
1017 1017 def byTime(self, datatime, *args):
1018 1018
1019 1019 self.__dataReady = False
1020 1020 avgdata_spc = None
1021 1021 avgdata_cspc = None
1022 1022 avgdata_dc = None
1023 1023
1024 1024 self.putData(*args)
1025 1025
1026 1026 if (datatime - self.__initime) >= self.__integrationtime:
1027 1027 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1028 1028 self.n = n
1029 1029 self.__dataReady = True
1030 1030
1031 1031 return avgdata_spc, avgdata_cspc, avgdata_dc
1032 1032
1033 1033 def integrate(self, datatime, *args):
1034 1034
1035 1035 if self.__profIndex == 0:
1036 1036 self.__initime = datatime
1037 1037
1038 1038 if self.__byTime:
1039 1039 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1040 1040 datatime, *args)
1041 1041 else:
1042 1042 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1043 1043
1044 1044 if not self.__dataReady:
1045 1045 return None, None, None, None
1046 1046
1047 1047 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1048 1048
1049 1049 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1050 1050 if n == 1:
1051 1051 return
1052 1052
1053 1053 dataOut.flagNoData = True
1054 1054
1055 1055 if not self.isConfig:
1056 1056 self.setup(n, timeInterval, overlapping)
1057 1057 self.isConfig = True
1058 1058
1059 1059 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1060 1060 dataOut.data_spc,
1061 1061 dataOut.data_cspc,
1062 1062 dataOut.data_dc)
1063 1063
1064 1064 if self.__dataReady:
1065 1065
1066 1066 dataOut.data_spc = avgdata_spc
1067 1067 dataOut.data_cspc = avgdata_cspc
1068 1068 dataOut.data_dc = avgdata_dc
1069 1069 dataOut.nIncohInt *= self.n
1070 1070 dataOut.utctime = avgdatatime
1071 1071 dataOut.flagNoData = False
1072 1072
1073 1073 return dataOut
1074 1074
1075 1075
1076 1076 class PulsePair(Operation):
1077 1077 isConfig = False
1078 1078 __profIndex = 0
1079 1079 __profIndex2 = 0
1080 1080 __initime = None
1081 1081 __lastdatatime = None
1082 1082 __buffer = None
1083 1083 __buffer2 = []
1084 1084 __buffer3 = None
1085 1085 __dataReady = False
1086 1086 n = None
1087 1087
1088 1088 __nch =0
1089 1089 __nProf =0
1090 1090 __nHeis =0
1091 1091
1092 1092 def __init__(self,**kwargs):
1093 1093 Operation.__init__(self,**kwargs)
1094 1094
1095 1095 def setup(self,dataOut,n =None, m = None):
1096 1096
1097 1097 self.__initime = None
1098 1098 self.__lastdatatime = 0
1099 1099 self.__buffer = 0
1100 1100 self.__bufferV = 0
1101 1101 #self.__buffer2 = []
1102 1102 self.__buffer3 = 0
1103 1103 self.__dataReady = False
1104 1104 self.__profIndex = 0
1105 1105 self.__profIndex2 = 0
1106 1106 self.count = 0
1107 1107
1108 1108
1109 1109 self.__nch = dataOut.nChannels
1110 1110 self.__nHeis = dataOut.nHeights
1111 1111 self.__nProf = dataOut.nProfiles
1112 1112 self.__nFFT = dataOut.nFFTPoints
1113 1113 #print("Valores de Ch,Samples,Perfiles,nFFT",self.__nch,self.__nHeis,self.__nProf, self.__nFFT)
1114 1114 #print("EL VALOR DE n es:",n)
1115 1115 if n == None:
1116 1116 raise ValueError("n Should be specified.")
1117 1117
1118 1118 if n != None:
1119 1119 if n<2:
1120 1120 raise ValueError("n Should be greather than 2 ")
1121 1121 self.n = n
1122 1122 if m == None:
1123 1123 m = n
1124 1124 if m != None:
1125 1125 if m<2:
1126 1126 raise ValueError("n Should be greather than 2 ")
1127 1127
1128 1128 self.m = m
1129 1129 self.__buffer2 = numpy.zeros((self.__nch,self.m,self.__nHeis))
1130 1130 self.__bufferV2 = numpy.zeros((self.__nch,self.m,self.__nHeis))
1131 1131
1132 1132
1133 1133
1134 1134 def putData(self,data):
1135 1135 #print("###################################################")
1136 1136 '''
1137 1137 data_tmp = numpy.zeros(self.__nch,self.n,self.__nHeis, dtype= complex)
1138 1138 if self.count < self.__nProf:
1139 1139
1140 1140 for i in range(self.n):
1141 1141 data_tmp[:,i,:] = data[:,i+self.count,:]
1142 1142
1143 1143 self.__buffer = data_tmp*numpy.conjugate(data_tmp)
1144 1144
1145 1145
1146 1146 #####self.__buffer = data*numpy.conjugate(data)
1147 1147 #####self.__bufferV = data[:,(self.__nProf-1):,:]*numpy.conjugate(data[:,1:,:])
1148 1148
1149 1149 #self.__buffer2.append(numpy.conjugate(data))
1150 1150
1151 1151 #####self.__profIndex = data.shape[1]
1152 1152 self.count = self.count + self.n -1
1153 1153 self.__profIndex = self.n
1154 1154 '''
1155 1155 self.__buffer = data*numpy.conjugate(data)
1156 1156 self.__bufferV = data[:,(self.__nProf-1):,:]*numpy.conjugate(data[:,1:,:])
1157 1157 self.__profIndex = self.n
1158 1158 #print("spcch0",self.__buffer)
1159 1159 return
1160 1160
1161 1161 def pushData(self):
1162 1162
1163 1163 data_I = numpy.zeros((self.__nch,self.__nHeis))
1164 1164 data_IV = numpy.zeros((self.__nch,self.__nHeis))
1165 1165
1166 1166 for i in range(self.__nch):
1167 1167 data_I[i,:] = numpy.sum(self.__buffer[i],axis=0)/self.n
1168 1168 data_IV[i,:] = numpy.sum(self.__bufferV[i],axis=0)/(self.n-1)
1169 1169 ##print("******")
1170 1170 #print("data_I",data_I[0])
1171 1171 #print(self.__buffer.shape)
1172 1172 #a=numpy.average(self.__buffer,axis=1)
1173 1173 #print("average", a)
1174 1174 n = self.__profIndex
1175 1175 ####data_intensity = numpy.sum(numpy.sum(self.__buffer,axis=0),axis=0)/self.n
1176 1176 #print("data_intensity push data",data_intensity.shape)
1177 1177 #data_velocity = self.__buffer3/(self.n-1)
1178 1178 ####n = self.__profIndex
1179 1179
1180 1180 self.__buffer = 0
1181 1181 self.__buffer3 = 0
1182 1182 self.__profIndex = 0
1183 1183
1184 1184 #return data_intensity,data_velocity,n
1185 1185 return data_I,data_IV,n
1186 1186
1187 1187 def pulsePairbyProfiles(self,data):
1188 1188 self.__dataReady = False
1189 1189 data_intensity = None
1190 1190 data_velocity = None
1191 1191
1192 1192 self.putData(data)
1193 1193
1194 1194 if self.__profIndex == self.n:
1195 1195 #data_intensity,data_velocity,n = self.pushData()
1196 1196 data_intensity,data_velocity,n = self.pushData()
1197 1197 #print(data_intensity.shape)
1198 1198 #print("self.__profIndex2", self.__profIndex2)
1199 1199 if self.__profIndex2 == 0:
1200 1200 #print("PRIMERA VEZ")
1201 1201 #print("self.__buffer2",self.__buffer2)
1202 1202 for i in range(self.__nch):
1203 1203 self.__buffer2[i][self.__profIndex2] = data_intensity[i]
1204 1204 self.__bufferV2[i][self.__profIndex2] = data_velocity[i]
1205 1205 self.__profIndex2 += 1
1206 1206 return None,None
1207 1207
1208 1208 if self.__profIndex2 > 0:
1209 1209 for i in range(self.__nch):
1210 1210 self.__buffer2[i][self.__profIndex2] = data_intensity[i]
1211 1211 self.__bufferV2[i][self.__profIndex2] = data_velocity[i]
1212 1212 #print("Dentro del bucle",self.__buffer2)
1213 1213 self.__profIndex2 += 1
1214 1214 if self.__profIndex2 == self.m :
1215 1215 data_i = self.__buffer2
1216 1216 data_v = self.__bufferV2
1217 1217 #print(data_i.shape)
1218 1218 self.__dataReady = True
1219 1219 self.__profIndex2 = 0
1220 1220 self.__buffer2 = numpy.zeros((self.__nch,self.m,self.__nHeis))
1221 1221 self.__bufferV2 = numpy.zeros((self.__nch,self.m,self.__nHeis))
1222 1222 return data_i,data_v
1223 1223 return None,None
1224 1224
1225 1225 def pulsePairOp(self,data,datatime=None):
1226 1226 if self.__initime == None:
1227 1227 self.__initime = datatime
1228 1228
1229 1229 data_intensity,data_velocity = self.pulsePairbyProfiles(data)
1230 1230 self.__lastdatatime = datatime
1231 1231
1232 1232 if data_intensity is None:
1233 1233 return None,None,None
1234 1234
1235 1235 avgdatatime = self.__initime
1236 1236 self.__initime = datatime
1237 1237
1238 1238 return data_intensity,data_velocity,avgdatatime
1239 1239
1240 1240 def run(self,dataOut,n =None,m=None):
1241 1241
1242 1242 if not self.isConfig:
1243 1243 self.setup(dataOut = dataOut, n = n, m = m)
1244 1244 self.isConfig = True
1245 1245
1246 1246 data_intensity,data_velocity,avgdatatime = self.pulsePairOp(dataOut.data_wr,dataOut.utctime)
1247 1247 dataOut.flagNoData = True
1248 1248
1249 1249 if self.__dataReady:
1250 1250 #print(" DATA " , data_intensity.shape)
1251 1251 #dataOut.data = numpy.array([data_intensity])#aqui amigo revisa
1252 1252 #tmp = numpy.zeros([1,data_intensity.shape[0],data_intensity.shape[1]])
1253 1253 #tmp[0] = data_intensity
1254 1254 dataOut.data = data_intensity
1255 1255 dataOut.data_velocity = data_velocity
1256 1256 #dataOut.data = tmp
1257 1257 #print(" DATA " , dataOut.data.shape)
1258 1258 dataOut.nIncohInt *= self.n
1259 1259 dataOut.nProfiles = self.m
1260 1260 dataOut.nFFTPoints = self.m
1261 1261 #dataOut.data_intensity = data_intensity
1262 1262 dataOut.PRFbyAngle = self.n
1263 1263 dataOut.utctime = avgdatatime
1264 1264 dataOut.flagNoData = False
1265 1265 #####print("TIEMPO: ",dataOut.utctime)
1266 1266 return dataOut
Show file before | Show file after | Hide comments
@@ -1,1623 +1,1431
1 1 import sys
2 import numpy
2 import time
3 import numpy,math
3 4 from scipy import interpolate
4 5 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 6 from schainpy.model.data.jrodata import Voltage
6 7 from schainpy.utils import log
7 8 from time import time
8 9
9 10
10 11 @MPDecorator
11 12 class VoltageProc(ProcessingUnit):
12 13
13 14 def __init__(self):
14 15
15 16 ProcessingUnit.__init__(self)
16 17
17 18 self.dataOut = Voltage()
18 19 self.flip = 1
19 20 self.setupReq = False
20 21
21 22 def run(self):
22 23
23 24 if self.dataIn.type == 'AMISR':
24 25 self.__updateObjFromAmisrInput()
25 26
26 27 if self.dataIn.type == 'Voltage':
27 28 self.dataOut.copy(self.dataIn)
28 29
29 30 # self.dataOut.copy(self.dataIn)
30 31
31 32 def __updateObjFromAmisrInput(self):
32 33
33 34 self.dataOut.timeZone = self.dataIn.timeZone
34 35 self.dataOut.dstFlag = self.dataIn.dstFlag
35 36 self.dataOut.errorCount = self.dataIn.errorCount
36 37 self.dataOut.useLocalTime = self.dataIn.useLocalTime
37 38
38 39 self.dataOut.flagNoData = self.dataIn.flagNoData
39 40 self.dataOut.data = self.dataIn.data
40 41 self.dataOut.utctime = self.dataIn.utctime
41 42 self.dataOut.channelList = self.dataIn.channelList
42 43 #self.dataOut.timeInterval = self.dataIn.timeInterval
43 44 self.dataOut.heightList = self.dataIn.heightList
44 45 self.dataOut.nProfiles = self.dataIn.nProfiles
45 46
46 47 self.dataOut.nCohInt = self.dataIn.nCohInt
47 48 self.dataOut.ippSeconds = self.dataIn.ippSeconds
48 49 self.dataOut.frequency = self.dataIn.frequency
49 50
50 51 self.dataOut.azimuth = self.dataIn.azimuth
51 52 self.dataOut.zenith = self.dataIn.zenith
52 53
53 54 self.dataOut.beam.codeList = self.dataIn.beam.codeList
54 55 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
55 56 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
56 57 #
57 58 # pass#
58 59 #
59 60 # def init(self):
60 61 #
61 62 #
62 63 # if self.dataIn.type == 'AMISR':
63 64 # self.__updateObjFromAmisrInput()
64 65 #
65 66 # if self.dataIn.type == 'Voltage':
66 67 # self.dataOut.copy(self.dataIn)
67 68 # # No necesita copiar en cada init() los atributos de dataIn
68 69 # # la copia deberia hacerse por cada nuevo bloque de datos
69 70
70 71 def selectChannels(self, channelList):
71 72
72 73 channelIndexList = []
73 74
74 75 for channel in channelList:
75 76 if channel not in self.dataOut.channelList:
76 77 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
77 78
78 79 index = self.dataOut.channelList.index(channel)
79 80 channelIndexList.append(index)
80 81
81 82 self.selectChannelsByIndex(channelIndexList)
82 83
83 84 def selectChannelsByIndex(self, channelIndexList):
84 85 """
85 86 Selecciona un bloque de datos en base a canales segun el channelIndexList
86 87
87 88 Input:
88 89 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
89 90
90 91 Affected:
91 92 self.dataOut.data
92 93 self.dataOut.channelIndexList
93 94 self.dataOut.nChannels
94 95 self.dataOut.m_ProcessingHeader.totalSpectra
95 96 self.dataOut.systemHeaderObj.numChannels
96 97 self.dataOut.m_ProcessingHeader.blockSize
97 98
98 99 Return:
99 100 None
100 101 """
101 102
102 103 for channelIndex in channelIndexList:
103 104 if channelIndex not in self.dataOut.channelIndexList:
104 105 print(channelIndexList)
105 106 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
106 107
107 108 if self.dataOut.flagDataAsBlock:
108 109 """
109 110 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
110 111 """
111 112 data = self.dataOut.data[channelIndexList,:,:]
112 113 else:
113 114 data = self.dataOut.data[channelIndexList,:]
114 115
115 116 self.dataOut.data = data
116 117 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
117 118 self.dataOut.channelList = range(len(channelIndexList))
118 119
119 120 return 1
120 121
121 122 def selectHeights(self, minHei=None, maxHei=None):
122 123 """
123 124 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
124 125 minHei <= height <= maxHei
125 126
126 127 Input:
127 128 minHei : valor minimo de altura a considerar
128 129 maxHei : valor maximo de altura a considerar
129 130
130 131 Affected:
131 132 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
132 133
133 134 Return:
134 135 1 si el metodo se ejecuto con exito caso contrario devuelve 0
135 136 """
136 137
137 138 if minHei == None:
138 139 minHei = self.dataOut.heightList[0]
139 140
140 141 if maxHei == None:
141 142 maxHei = self.dataOut.heightList[-1]
142 143
143 144 if (minHei < self.dataOut.heightList[0]):
144 145 minHei = self.dataOut.heightList[0]
145 146
146 147 if (maxHei > self.dataOut.heightList[-1]):
147 148 maxHei = self.dataOut.heightList[-1]
148 149
149 150 minIndex = 0
150 151 maxIndex = 0
151 152 heights = self.dataOut.heightList
152 153
153 154 inda = numpy.where(heights >= minHei)
154 155 indb = numpy.where(heights <= maxHei)
155 156
156 157 try:
157 158 minIndex = inda[0][0]
158 159 except:
159 160 minIndex = 0
160 161
161 162 try:
162 163 maxIndex = indb[0][-1]
163 164 except:
164 165 maxIndex = len(heights)
165 166
166 167 self.selectHeightsByIndex(minIndex, maxIndex)
167 168
168 169 return 1
169 170
170 171
171 172 def selectHeightsByIndex(self, minIndex, maxIndex):
172 173 """
173 174 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
174 175 minIndex <= index <= maxIndex
175 176
176 177 Input:
177 178 minIndex : valor de indice minimo de altura a considerar
178 179 maxIndex : valor de indice maximo de altura a considerar
179 180
180 181 Affected:
181 182 self.dataOut.data
182 183 self.dataOut.heightList
183 184
184 185 Return:
185 186 1 si el metodo se ejecuto con exito caso contrario devuelve 0
186 187 """
187 188
188 189 if (minIndex < 0) or (minIndex > maxIndex):
189 190 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
190 191
191 192 if (maxIndex >= self.dataOut.nHeights):
192 193 maxIndex = self.dataOut.nHeights
193 194
194 195 #voltage
195 196 if self.dataOut.flagDataAsBlock:
196 197 """
197 198 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
198 199 """
199 200 data = self.dataOut.data[:,:, minIndex:maxIndex]
200 201 else:
201 202 data = self.dataOut.data[:, minIndex:maxIndex]
202 203
203 204 # firstHeight = self.dataOut.heightList[minIndex]
204 205
205 206 self.dataOut.data = data
206 207 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
207 208
208 209 if self.dataOut.nHeights <= 1:
209 210 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
210 211
211 212 return 1
212 213
213 214
214 215 def filterByHeights(self, window):
215 216
216 217 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
217 218
218 219 if window == None:
219 220 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
220 221
221 222 newdelta = deltaHeight * window
222 223 r = self.dataOut.nHeights % window
223 224 newheights = (self.dataOut.nHeights-r)/window
224 225
225 226 if newheights <= 1:
226 227 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window))
227 228
228 229 if self.dataOut.flagDataAsBlock:
229 230 """
230 231 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
231 232 """
232 233 buffer = self.dataOut.data[:, :, 0:int(self.dataOut.nHeights-r)]
233 234 buffer = buffer.reshape(self.dataOut.nChannels, self.dataOut.nProfiles, int(self.dataOut.nHeights/window), window)
234 235 buffer = numpy.sum(buffer,3)
235 236
236 237 else:
237 238 buffer = self.dataOut.data[:,0:int(self.dataOut.nHeights-r)]
238 239 buffer = buffer.reshape(self.dataOut.nChannels,int(self.dataOut.nHeights/window),int(window))
239 240 buffer = numpy.sum(buffer,2)
240 241
241 242 self.dataOut.data = buffer
242 243 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
243 244 self.dataOut.windowOfFilter = window
244 245
245 246 def setH0(self, h0, deltaHeight = None):
246 247
247 248 if not deltaHeight:
248 249 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
249 250
250 251 nHeights = self.dataOut.nHeights
251 252
252 253 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
253 254
254 255 self.dataOut.heightList = newHeiRange
255 256
256 257 def deFlip(self, channelList = []):
257 258
258 259 data = self.dataOut.data.copy()
259 260
260 261 if self.dataOut.flagDataAsBlock:
261 262 flip = self.flip
262 263 profileList = list(range(self.dataOut.nProfiles))
263 264
264 265 if not channelList:
265 266 for thisProfile in profileList:
266 267 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
267 268 flip *= -1.0
268 269 else:
269 270 for thisChannel in channelList:
270 271 if thisChannel not in self.dataOut.channelList:
271 272 continue
272 273
273 274 for thisProfile in profileList:
274 275 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
275 276 flip *= -1.0
276 277
277 278 self.flip = flip
278 279
279 280 else:
280 281 if not channelList:
281 282 data[:,:] = data[:,:]*self.flip
282 283 else:
283 284 for thisChannel in channelList:
284 285 if thisChannel not in self.dataOut.channelList:
285 286 continue
286 287
287 288 data[thisChannel,:] = data[thisChannel,:]*self.flip
288 289
289 290 self.flip *= -1.
290 291
291 292 self.dataOut.data = data
292 293
293 294 def setRadarFrequency(self, frequency=None):
294 295
295 296 if frequency != None:
296 297 self.dataOut.frequency = frequency
297 298
298 299 return 1
299 300
300 301 def interpolateHeights(self, topLim, botLim):
301 302 #69 al 72 para julia
302 303 #82-84 para meteoros
303 304 if len(numpy.shape(self.dataOut.data))==2:
304 305 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
305 306 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
306 307 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
307 308 self.dataOut.data[:,botLim:topLim+1] = sampInterp
308 309 else:
309 310 nHeights = self.dataOut.data.shape[2]
310 311 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
311 312 y = self.dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
312 313 f = interpolate.interp1d(x, y, axis = 2)
313 314 xnew = numpy.arange(botLim,topLim+1)
314 315 ynew = f(xnew)
315 316
316 317 self.dataOut.data[:,:,botLim:topLim+1] = ynew
317 318
318 319 # import collections
319 320
320 321 class CohInt(Operation):
321 322
322 323 isConfig = False
323 324 __profIndex = 0
324 325 __byTime = False
325 326 __initime = None
326 327 __lastdatatime = None
327 328 __integrationtime = None
328 329 __buffer = None
329 330 __bufferStride = []
330 331 __dataReady = False
331 332 __profIndexStride = 0
332 333 __dataToPutStride = False
333 334 n = None
334 335
335 336 def __init__(self, **kwargs):
336 337
337 338 Operation.__init__(self, **kwargs)
338 339
339 340 # self.isConfig = False
340 341
341 342 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
342 343 """
343 344 Set the parameters of the integration class.
344 345
345 346 Inputs:
346 347
347 348 n : Number of coherent integrations
348 349 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
349 350 overlapping :
350 351 """
351 352
352 353 self.__initime = None
353 354 self.__lastdatatime = 0
354 355 self.__buffer = None
355 356 self.__dataReady = False
356 357 self.byblock = byblock
357 358 self.stride = stride
358 359
359 360 if n == None and timeInterval == None:
360 361 raise ValueError("n or timeInterval should be specified ...")
361 362
362 363 if n != None:
363 364 self.n = n
364 365 self.__byTime = False
365 366 else:
366 367 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
367 368 self.n = 9999
368 369 self.__byTime = True
369 370
370 371 if overlapping:
371 372 self.__withOverlapping = True
372 373 self.__buffer = None
373 374 else:
374 375 self.__withOverlapping = False
375 376 self.__buffer = 0
376 377
377 378 self.__profIndex = 0
378 379
379 380 def putData(self, data):
380 381
381 382 """
382 383 Add a profile to the __buffer and increase in one the __profileIndex
383 384
384 385 """
385 386
386 387 if not self.__withOverlapping:
387 print("inside over")
388 #print("inside over")
388 389 self.__buffer += data.copy()
389 390 self.__profIndex += 1
390 391 return
391 392
392 393 #Overlapping data
393 394 nChannels, nHeis = data.shape
394 print("show me the light",data.shape)
395 #print("show me the light",data.shape)
395 396 data = numpy.reshape(data, (1, nChannels, nHeis))
396 print(data.shape)
397 #print(data.shape)
397 398 #If the buffer is empty then it takes the data value
398 399 if self.__buffer is None:
399 400 self.__buffer = data
400 401 self.__profIndex += 1
401 402 return
402 403
403 404 #If the buffer length is lower than n then stakcing the data value
404 405 if self.__profIndex < self.n:
405 406 self.__buffer = numpy.vstack((self.__buffer, data))
406 407 self.__profIndex += 1
407 408 return
408 409
409 410 #If the buffer length is equal to n then replacing the last buffer value with the data value
410 411 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
411 412 self.__buffer[self.n-1] = data
412 413 self.__profIndex = self.n
413 414 return
414 415
415 416
416 417 def pushData(self):
417 418 """
418 419 Return the sum of the last profiles and the profiles used in the sum.
419 420
420 421 Affected:
421 422
422 423 self.__profileIndex
423 424
424 425 """
425 426
426 427 if not self.__withOverlapping:
427 #print("ahora que fue")
428 print("ahora que fue")
428 429 data = self.__buffer
429 430 n = self.__profIndex
430 431
431 432 self.__buffer = 0
432 433 self.__profIndex = 0
433 434
434 435 return data, n
435 436
436 437 #print("cual funciona")
437 438 #Integration with Overlapping
438 439 data = numpy.sum(self.__buffer, axis=0)
439 440 # print data
440 441 # raise
441 442 n = self.__profIndex
442 443
443 444 return data, n
444 445
445 446 def byProfiles(self, data):
446 447
447 448 self.__dataReady = False
448 449 avgdata = None
449 450 # n = None
450 451 # print data
451 452 # raise
452 453 #print("beforeputdata")
453 454 self.putData(data)
454 455
455 456 if self.__profIndex == self.n:
456 457 avgdata, n = self.pushData()
457 458 self.__dataReady = True
458 459
459 460 return avgdata
460 461
461 462 def byTime(self, data, datatime):
462 463
463 464 self.__dataReady = False
464 465 avgdata = None
465 466 n = None
466 467
467 468 self.putData(data)
468 469
469 470 if (datatime - self.__initime) >= self.__integrationtime:
470 471 avgdata, n = self.pushData()
471 472 self.n = n
472 473 self.__dataReady = True
473 474
474 475 return avgdata
475 476
476 477 def integrateByStride(self, data, datatime):
477 478 # print data
478 479 if self.__profIndex == 0:
479 480 self.__buffer = [[data.copy(), datatime]]
480 481 else:
481 482 self.__buffer.append([data.copy(),datatime])
482 483 self.__profIndex += 1
483 484 self.__dataReady = False
484 485
485 486 if self.__profIndex == self.n * self.stride :
486 487 self.__dataToPutStride = True
487 488 self.__profIndexStride = 0
488 489 self.__profIndex = 0
489 490 self.__bufferStride = []
490 491 for i in range(self.stride):
491 492 current = self.__buffer[i::self.stride]
492 493 data = numpy.sum([t[0] for t in current], axis=0)
493 494 avgdatatime = numpy.average([t[1] for t in current])
494 495 # print data
495 496 self.__bufferStride.append((data, avgdatatime))
496 497
497 498 if self.__dataToPutStride:
498 499 self.__dataReady = True
499 500 self.__profIndexStride += 1
500 501 if self.__profIndexStride == self.stride:
501 502 self.__dataToPutStride = False
502 503 # print self.__bufferStride[self.__profIndexStride - 1]
503 504 # raise
504 505 return self.__bufferStride[self.__profIndexStride - 1]
505 506
506 507
507 508 return None, None
508 509
509 510 def integrate(self, data, datatime=None):
510 511
511 512 if self.__initime == None:
512 513 self.__initime = datatime
513 514
514 515 if self.__byTime:
515 516 avgdata = self.byTime(data, datatime)
516 517 else:
517 518 avgdata = self.byProfiles(data)
518 519
519 520
520 521 self.__lastdatatime = datatime
521 522
522 523 if avgdata is None:
523 524 return None, None
524 525
525 526 avgdatatime = self.__initime
526 527
527 528 deltatime = datatime - self.__lastdatatime
528 529
529 530 if not self.__withOverlapping:
530 531 self.__initime = datatime
531 532 else:
532 533 self.__initime += deltatime
533 534
534 535 return avgdata, avgdatatime
535 536
536 537 def integrateByBlock(self, dataOut):
537 538
538 539 times = int(dataOut.data.shape[1]/self.n)
539 540 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
540 541
541 542 id_min = 0
542 543 id_max = self.n
543 544
544 545 for i in range(times):
545 546 junk = dataOut.data[:,id_min:id_max,:]
546 547 avgdata[:,i,:] = junk.sum(axis=1)
547 548 id_min += self.n
548 549 id_max += self.n
549 550
550 551 timeInterval = dataOut.ippSeconds*self.n
551 552 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
552 553 self.__dataReady = True
553 554 return avgdata, avgdatatime
554 555
555 556 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
556 557
557 558 if not self.isConfig:
558 559 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
559 560 self.isConfig = True
560 561
561 562 if dataOut.flagDataAsBlock:
562 563 """
563 564 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
564 565 """
565 566 avgdata, avgdatatime = self.integrateByBlock(dataOut)
566 567 dataOut.nProfiles /= self.n
567 568 else:
568 569 if stride is None:
569 570 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
570 571 else:
571 572 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
572 573
573 574
574 575 # dataOut.timeInterval *= n
575 576 dataOut.flagNoData = True
576 577
577 578 if self.__dataReady:
578 579 dataOut.data = avgdata
579 580 dataOut.nCohInt *= self.n
580 581 dataOut.utctime = avgdatatime
581 582 # print avgdata, avgdatatime
582 583 # raise
583 584 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
584 585 dataOut.flagNoData = False
585 586 return dataOut
586 587
587 588 class Decoder(Operation):
588 589
589 590 isConfig = False
590 591 __profIndex = 0
591 592
592 593 code = None
593 594
594 595 nCode = None
595 596 nBaud = None
596 597
597 598 def __init__(self, **kwargs):
598 599
599 600 Operation.__init__(self, **kwargs)
600 601
601 602 self.times = None
602 603 self.osamp = None
603 604 # self.__setValues = False
604 605 self.isConfig = False
605 606 self.setupReq = False
606 607 def setup(self, code, osamp, dataOut):
607 608
608 609 self.__profIndex = 0
609 610
610 611 self.code = code
611 612
612 613 self.nCode = len(code)
613 614 self.nBaud = len(code[0])
614 615
615 616 if (osamp != None) and (osamp >1):
616 617 self.osamp = osamp
617 618 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
618 619 self.nBaud = self.nBaud*self.osamp
619 620
620 621 self.__nChannels = dataOut.nChannels
621 622 self.__nProfiles = dataOut.nProfiles
622 623 self.__nHeis = dataOut.nHeights
623 624
624 625 if self.__nHeis < self.nBaud:
625 626 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
626 627
627 628 #Frequency
628 629 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
629 630
630 631 __codeBuffer[:,0:self.nBaud] = self.code
631 632
632 633 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
633 634
634 635 if dataOut.flagDataAsBlock:
635 636
636 637 self.ndatadec = self.__nHeis #- self.nBaud + 1
637 638
638 639 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
639 640
640 641 else:
641 642
642 643 #Time
643 644 self.ndatadec = self.__nHeis #- self.nBaud + 1
644 645
645 646 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
646 647
647 648 def __convolutionInFreq(self, data):
648 649
649 650 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
650 651
651 652 fft_data = numpy.fft.fft(data, axis=1)
652 653
653 654 conv = fft_data*fft_code
654 655
655 656 data = numpy.fft.ifft(conv,axis=1)
656 657
657 658 return data
658 659
659 660 def __convolutionInFreqOpt(self, data):
660 661
661 662 raise NotImplementedError
662 663
663 664 def __convolutionInTime(self, data):
664 665
665 666 code = self.code[self.__profIndex]
666 667 for i in range(self.__nChannels):
667 668 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
668 669
669 670 return self.datadecTime
670 671
671 672 def __convolutionByBlockInTime(self, data):
672 673
673 674 repetitions = int(self.__nProfiles / self.nCode)
674 675 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
675 676 junk = junk.flatten()
676 677 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
677 678 profilesList = range(self.__nProfiles)
678 679
679 680 for i in range(self.__nChannels):
680 681 for j in profilesList:
681 682 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
682 683 return self.datadecTime
683 684
684 685 def __convolutionByBlockInFreq(self, data):
685 686
686 687 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
687 688
688 689
689 690 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
690 691
691 692 fft_data = numpy.fft.fft(data, axis=2)
692 693
693 694 conv = fft_data*fft_code
694 695
695 696 data = numpy.fft.ifft(conv,axis=2)
696 697
697 698 return data
698 699
699 700
700 701 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
701 702
702 703 if dataOut.flagDecodeData:
703 704 print("This data is already decoded, recoding again ...")
704 705
705 706 if not self.isConfig:
706 707
707 708 if code is None:
708 709 if dataOut.code is None:
709 710 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
710 711
711 712 code = dataOut.code
712 713 else:
713 714 code = numpy.array(code).reshape(nCode,nBaud)
714 715 self.setup(code, osamp, dataOut)
715 716
716 717 self.isConfig = True
717 718
718 719 if mode == 3:
719 720 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
720 721
721 722 if times != None:
722 723 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
723 724
724 725 if self.code is None:
725 726 print("Fail decoding: Code is not defined.")
726 727 return
727 728
728 729 self.__nProfiles = dataOut.nProfiles
729 730 datadec = None
730 731
731 732 if mode == 3:
732 733 mode = 0
733 734
734 735 if dataOut.flagDataAsBlock:
735 736 """
736 737 Decoding when data have been read as block,
737 738 """
738 739
739 740 if mode == 0:
740 741 datadec = self.__convolutionByBlockInTime(dataOut.data)
741 742 if mode == 1:
742 743 datadec = self.__convolutionByBlockInFreq(dataOut.data)
743 744 else:
744 745 """
745 746 Decoding when data have been read profile by profile
746 747 """
747 748 if mode == 0:
748 749 datadec = self.__convolutionInTime(dataOut.data)
749 750
750 751 if mode == 1:
751 752 datadec = self.__convolutionInFreq(dataOut.data)
752 753
753 754 if mode == 2:
754 755 datadec = self.__convolutionInFreqOpt(dataOut.data)
755 756
756 757 if datadec is None:
757 758 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
758 759
759 760 dataOut.code = self.code
760 761 dataOut.nCode = self.nCode
761 762 dataOut.nBaud = self.nBaud
762 763
763 764 dataOut.data = datadec
764 765
765 766 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
766 767
767 768 dataOut.flagDecodeData = True #asumo q la data esta decodificada
768 769
769 770 if self.__profIndex == self.nCode-1:
770 771 self.__profIndex = 0
771 772 return dataOut
772 773
773 774 self.__profIndex += 1
774 775
775 776 return dataOut
776 777 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
777 778
778 779
779 780 class ProfileConcat(Operation):
780 781
781 782 isConfig = False
782 783 buffer = None
783 784
784 785 def __init__(self, **kwargs):
785 786
786 787 Operation.__init__(self, **kwargs)
787 788 self.profileIndex = 0
788 789
789 790 def reset(self):
790 791 self.buffer = numpy.zeros_like(self.buffer)
791 792 self.start_index = 0
792 793 self.times = 1
793 794
794 795 def setup(self, data, m, n=1):
795 796 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
796 797 self.nHeights = data.shape[1]#.nHeights
797 798 self.start_index = 0
798 799 self.times = 1
799 800
800 801 def concat(self, data):
801 802
802 803 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
803 804 self.start_index = self.start_index + self.nHeights
804 805
805 806 def run(self, dataOut, m):
806 807 dataOut.flagNoData = True
807 808
808 809 if not self.isConfig:
809 810 self.setup(dataOut.data, m, 1)
810 811 self.isConfig = True
811 812
812 813 if dataOut.flagDataAsBlock:
813 814 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
814 815
815 816 else:
816 817 self.concat(dataOut.data)
817 818 self.times += 1
818 819 if self.times > m:
819 820 dataOut.data = self.buffer
820 821 self.reset()
821 822 dataOut.flagNoData = False
822 823 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
823 824 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
824 825 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
825 826 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
826 827 dataOut.ippSeconds *= m
827 828 return dataOut
828 829
829 830 class ProfileSelector(Operation):
830 831
831 832 profileIndex = None
832 833 # Tamanho total de los perfiles
833 834 nProfiles = None
834 835
835 836 def __init__(self, **kwargs):
836 837
837 838 Operation.__init__(self, **kwargs)
838 839 self.profileIndex = 0
839 840
840 841 def incProfileIndex(self):
841 842
842 843 self.profileIndex += 1
843 844
844 845 if self.profileIndex >= self.nProfiles:
845 846 self.profileIndex = 0
846 847
847 848 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
848 849
849 850 if profileIndex < minIndex:
850 851 return False
851 852
852 853 if profileIndex > maxIndex:
853 854 return False
854 855
855 856 return True
856 857
857 858 def isThisProfileInList(self, profileIndex, profileList):
858 859
859 860 if profileIndex not in profileList:
860 861 return False
861 862
862 863 return True
863 864
864 865 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
865 866
866 867 """
867 868 ProfileSelector:
868 869
869 870 Inputs:
870 871 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
871 872
872 873 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
873 874
874 875 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
875 876
876 877 """
877 878
878 879 if rangeList is not None:
879 880 if type(rangeList[0]) not in (tuple, list):
880 881 rangeList = [rangeList]
881 882
882 883 dataOut.flagNoData = True
883 884
884 885 if dataOut.flagDataAsBlock:
885 886 """
886 887 data dimension = [nChannels, nProfiles, nHeis]
887 888 """
888 889 if profileList != None:
889 890 dataOut.data = dataOut.data[:,profileList,:]
890 891
891 892 if profileRangeList != None:
892 893 minIndex = profileRangeList[0]
893 894 maxIndex = profileRangeList[1]
894 895 profileList = list(range(minIndex, maxIndex+1))
895 896
896 897 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
897 898
898 899 if rangeList != None:
899 900
900 901 profileList = []
901 902
902 903 for thisRange in rangeList:
903 904 minIndex = thisRange[0]
904 905 maxIndex = thisRange[1]
905 906
906 907 profileList.extend(list(range(minIndex, maxIndex+1)))
907 908
908 909 dataOut.data = dataOut.data[:,profileList,:]
909 910
910 911 dataOut.nProfiles = len(profileList)
911 912 dataOut.profileIndex = dataOut.nProfiles - 1
912 913 dataOut.flagNoData = False
913 914
914 915 return dataOut
915 916
916 917 """
917 918 data dimension = [nChannels, nHeis]
918 919 """
919 920
920 921 if profileList != None:
921 922
922 923 if self.isThisProfileInList(dataOut.profileIndex, profileList):
923 924
924 925 self.nProfiles = len(profileList)
925 926 dataOut.nProfiles = self.nProfiles
926 927 dataOut.profileIndex = self.profileIndex
927 928 dataOut.flagNoData = False
928 929
929 930 self.incProfileIndex()
930 931 return dataOut
931 932
932 933 if profileRangeList != None:
933 934
934 935 minIndex = profileRangeList[0]
935 936 maxIndex = profileRangeList[1]
936 937
937 938 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
938 939
939 940 self.nProfiles = maxIndex - minIndex + 1
940 941 dataOut.nProfiles = self.nProfiles
941 942 dataOut.profileIndex = self.profileIndex
942 943 dataOut.flagNoData = False
943 944
944 945 self.incProfileIndex()
945 946 return dataOut
946 947
947 948 if rangeList != None:
948 949
949 950 nProfiles = 0
950 951
951 952 for thisRange in rangeList:
952 953 minIndex = thisRange[0]
953 954 maxIndex = thisRange[1]
954 955
955 956 nProfiles += maxIndex - minIndex + 1
956 957
957 958 for thisRange in rangeList:
958 959
959 960 minIndex = thisRange[0]
960 961 maxIndex = thisRange[1]
961 962
962 963 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
963 964
964 965 self.nProfiles = nProfiles
965 966 dataOut.nProfiles = self.nProfiles
966 967 dataOut.profileIndex = self.profileIndex
967 968 dataOut.flagNoData = False
968 969
969 970 self.incProfileIndex()
970 971
971 972 break
972 973
973 974 return dataOut
974 975
975 976
976 977 if beam != None: #beam is only for AMISR data
977 978 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
978 979 dataOut.flagNoData = False
979 980 dataOut.profileIndex = self.profileIndex
980 981
981 982 self.incProfileIndex()
982 983
983 984 return dataOut
984 985
985 986 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
986 987
987 988 #return False
988 989 return dataOut
989 990
990 991 class Reshaper(Operation):
991 992
992 993 def __init__(self, **kwargs):
993 994
994 995 Operation.__init__(self, **kwargs)
995 996
996 997 self.__buffer = None
997 998 self.__nitems = 0
998 999
999 1000 def __appendProfile(self, dataOut, nTxs):
1000 1001
1001 1002 if self.__buffer is None:
1002 1003 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1003 1004 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1004 1005
1005 1006 ini = dataOut.nHeights * self.__nitems
1006 1007 end = ini + dataOut.nHeights
1007 1008
1008 1009 self.__buffer[:, ini:end] = dataOut.data
1009 1010
1010 1011 self.__nitems += 1
1011 1012
1012 1013 return int(self.__nitems*nTxs)
1013 1014
1014 1015 def __getBuffer(self):
1015 1016
1016 1017 if self.__nitems == int(1./self.__nTxs):
1017 1018
1018 1019 self.__nitems = 0
1019 1020
1020 1021 return self.__buffer.copy()
1021 1022
1022 1023 return None
1023 1024
1024 1025 def __checkInputs(self, dataOut, shape, nTxs):
1025 1026
1026 1027 if shape is None and nTxs is None:
1027 1028 raise ValueError("Reshaper: shape of factor should be defined")
1028 1029
1029 1030 if nTxs:
1030 1031 if nTxs < 0:
1031 1032 raise ValueError("nTxs should be greater than 0")
1032 1033
1033 1034 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1034 1035 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1035 1036
1036 1037 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1037 1038
1038 1039 return shape, nTxs
1039 1040
1040 1041 if len(shape) != 2 and len(shape) != 3:
1041 1042 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))
1042 1043
1043 1044 if len(shape) == 2:
1044 1045 shape_tuple = [dataOut.nChannels]
1045 1046 shape_tuple.extend(shape)
1046 1047 else:
1047 1048 shape_tuple = list(shape)
1048 1049
1049 1050 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1050 1051
1051 1052 return shape_tuple, nTxs
1052 1053
1053 1054 def run(self, dataOut, shape=None, nTxs=None):
1054 1055
1055 1056 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1056 1057
1057 1058 dataOut.flagNoData = True
1058 1059 profileIndex = None
1059 1060
1060 1061 if dataOut.flagDataAsBlock:
1061 1062
1062 1063 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1063 1064 dataOut.flagNoData = False
1064 1065
1065 1066 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1066 1067
1067 1068 else:
1068 1069
1069 1070 if self.__nTxs < 1:
1070 1071
1071 1072 self.__appendProfile(dataOut, self.__nTxs)
1072 1073 new_data = self.__getBuffer()
1073 1074
1074 1075 if new_data is not None:
1075 1076 dataOut.data = new_data
1076 1077 dataOut.flagNoData = False
1077 1078
1078 1079 profileIndex = dataOut.profileIndex*nTxs
1079 1080
1080 1081 else:
1081 1082 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1082 1083
1083 1084 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1084 1085
1085 1086 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1086 1087
1087 1088 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1088 1089
1089 1090 dataOut.profileIndex = profileIndex
1090 1091
1091 1092 dataOut.ippSeconds /= self.__nTxs
1092 1093
1093 1094 return dataOut
1094 1095
1095 1096 class SplitProfiles(Operation):
1096 1097
1097 1098 def __init__(self, **kwargs):
1098 1099
1099 1100 Operation.__init__(self, **kwargs)
1100 1101
1101 1102 def run(self, dataOut, n):
1102 1103
1103 1104 dataOut.flagNoData = True
1104 1105 profileIndex = None
1105 1106
1106 1107 if dataOut.flagDataAsBlock:
1107 1108
1108 1109 #nchannels, nprofiles, nsamples
1109 1110 shape = dataOut.data.shape
1110 1111
1111 1112 if shape[2] % n != 0:
1112 1113 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1113 1114
1114 1115 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1115 1116
1116 1117 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1117 1118 dataOut.flagNoData = False
1118 1119
1119 1120 profileIndex = int(dataOut.nProfiles/n) - 1
1120 1121
1121 1122 else:
1122 1123
1123 1124 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1124 1125
1125 1126 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1126 1127
1127 1128 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1128 1129
1129 1130 dataOut.nProfiles = int(dataOut.nProfiles*n)
1130 1131
1131 1132 dataOut.profileIndex = profileIndex
1132 1133
1133 1134 dataOut.ippSeconds /= n
1134 1135
1135 1136 return dataOut
1136 1137
1137 1138 class CombineProfiles(Operation):
1138 1139 def __init__(self, **kwargs):
1139 1140
1140 1141 Operation.__init__(self, **kwargs)
1141 1142
1142 1143 self.__remData = None
1143 1144 self.__profileIndex = 0
1144 1145
1145 1146 def run(self, dataOut, n):
1146 1147
1147 1148 dataOut.flagNoData = True
1148 1149 profileIndex = None
1149 1150
1150 1151 if dataOut.flagDataAsBlock:
1151 1152
1152 1153 #nchannels, nprofiles, nsamples
1153 1154 shape = dataOut.data.shape
1154 1155 new_shape = shape[0], shape[1]/n, shape[2]*n
1155 1156
1156 1157 if shape[1] % n != 0:
1157 1158 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1158 1159
1159 1160 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1160 1161 dataOut.flagNoData = False
1161 1162
1162 1163 profileIndex = int(dataOut.nProfiles*n) - 1
1163 1164
1164 1165 else:
1165 1166
1166 1167 #nchannels, nsamples
1167 1168 if self.__remData is None:
1168 1169 newData = dataOut.data
1169 1170 else:
1170 1171 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1171 1172
1172 1173 self.__profileIndex += 1
1173 1174
1174 1175 if self.__profileIndex < n:
1175 1176 self.__remData = newData
1176 1177 #continue
1177 1178 return
1178 1179
1179 1180 self.__profileIndex = 0
1180 1181 self.__remData = None
1181 1182
1182 1183 dataOut.data = newData
1183 1184 dataOut.flagNoData = False
1184 1185
1185 1186 profileIndex = dataOut.profileIndex/n
1186 1187
1187 1188
1188 1189 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1189 1190
1190 1191 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1191 1192
1192 1193 dataOut.nProfiles = int(dataOut.nProfiles/n)
1193 1194
1194 1195 dataOut.profileIndex = profileIndex
1195 1196
1196 1197 dataOut.ippSeconds *= n
1197 1198
1198 1199 return dataOut
1199 1200
1200 1201
1201 1202
1202 1203 class CreateBlockVoltage(Operation):
1203 1204
1204 1205 isConfig = False
1205 1206 __Index = 0
1206 1207 bufferShape = None
1207 1208 buffer = None
1208 1209 firstdatatime = None
1209 1210
1210 1211 def __init__(self,**kwargs):
1211 1212 Operation.__init__(self,**kwargs)
1212 1213 self.isConfig = False
1213 1214 self.__Index = 0
1214 1215 self.firstdatatime = None
1215 1216
1216 1217 def setup(self,dataOut, m = None ):
1217 1218 '''
1218 1219 m= Numero perfiles
1219 1220 '''
1220 1221 #print("CONFIGURANDO CBV")
1221 1222 self.__nChannels = dataOut.nChannels
1222 1223 self.__nHeis = dataOut.nHeights
1223 1224 shape = dataOut.data.shape #nchannels, nprofiles, nsamples
1224 1225 #print("input nChannels",self.__nChannels)
1225 1226 #print("input nHeis",self.__nHeis)
1226 1227 #print("SETUP CREATE BLOCK VOLTAGE")
1227 1228 #print("input Shape",shape)
1228 1229 #print("dataOut.nProfiles",dataOut.nProfiles)
1229 1230 numberSamples = self.__nHeis
1230 1231 numberProfile = int(m)
1231 1232 dataOut.nProfiles = numberProfile
1232 1233 #print("new numberProfile",numberProfile)
1233 1234 #print("new numberSamples",numberSamples)
1234 1235
1235 1236 self.bufferShape = shape[0], numberProfile, numberSamples # nchannels,nprofiles,nsamples
1236 self.buffer = numpy.zeros((self.bufferShape))
1237 self.bufferVel = numpy.zeros((self.bufferShape))
1237 self.buffer = numpy.zeros([shape[0], numberProfile, numberSamples])
1238 self.bufferVel = numpy.zeros([shape[0], numberProfile, numberSamples])
1238 1239
1239 1240 def run(self, dataOut, m=None):
1240 1241 #print("RUN")
1241 1242 dataOut.flagNoData = True
1242 1243 dataOut.flagDataAsBlock = False
1243 1244 #print("BLOCK INDEX ",self.__Index)
1244 1245
1245 1246 if not self.isConfig:
1246 1247 self.setup(dataOut, m= m)
1247 1248 self.isConfig = True
1248 1249 if self.__Index < m:
1249 1250 #print("PROFINDEX BLOCK CBV",self.__Index)
1250 1251 self.buffer[:,self.__Index,:] = dataOut.data
1252 #corregir porque debe tener un perfil menos ojo
1251 1253 self.bufferVel[:,self.__Index,:] = dataOut.data_velocity
1252 1254 self.__Index += 1
1253 1255 dataOut.flagNoData = True
1254 1256
1255 1257 if self.firstdatatime == None:
1256 1258 self.firstdatatime = dataOut.utctime
1257 1259
1258 1260 if self.__Index == m:
1259 1261 #print("**********************************************")
1260 1262 #print("self.buffer.shape ",self.buffer.shape)
1261 1263 #print("##############",self.firstdatatime)
1262 1264 ##print("*********************************************")
1263 1265 ##print("*********************************************")
1264 1266 ##print("******* nProfiles *******", dataOut.nProfiles)
1265 1267 ##print("*********************************************")
1266 1268 ##print("*********************************************")
1267 1269 dataOut.data = self.buffer
1268 1270 dataOut.data_velocity = self.bufferVel
1269 1271 dataOut.utctime = self.firstdatatime
1270 1272 dataOut.nProfiles = m
1271 1273 self.firstdatatime = None
1272 1274 dataOut.flagNoData = False
1273 1275 dataOut.flagDataAsBlock = True
1274 1276 self.__Index = 0
1275 1277 dataOut.identifierWR = True
1276 1278 return dataOut
1277 1279
1278 1280 class PulsePairVoltage(Operation):
1279 1281 '''
1280 1282 Function PulsePair(Signal Power, Velocity)
1281 1283 The real component of Lag[0] provides Intensity Information
1282 1284 The imag component of Lag[1] Phase provides Velocity Information
1283 1285
1284 1286 Configuration Parameters:
1285 1287 nPRF = Number of Several PRF
1286 1288 theta = Degree Azimuth angel Boundaries
1287 1289
1288 1290 Input:
1289 1291 self.dataOut
1290 1292 lag[N]
1291 1293 Affected:
1292 1294 self.dataOut.spc
1293 1295 '''
1294 1296 isConfig = False
1295 1297 __profIndex = 0
1296 1298 __initime = None
1297 1299 __lastdatatime = None
1298 1300 __buffer = None
1299 1301 __buffer2 = []
1300 1302 __buffer3 = None
1301 1303 __dataReady = False
1302 1304 n = None
1303 1305 __nch = 0
1304 1306 __nHeis = 0
1307 removeDC = False
1308 ipp = None
1309 lambda_ = 0
1305 1310
1306 1311 def __init__(self,**kwargs):
1307 1312 Operation.__init__(self,**kwargs)
1308 1313
1309 def setup(self, dataOut, n = None ):
1314 def setup(self, dataOut, n = None, removeDC=False):
1310 1315 '''
1311 1316 n= Numero de PRF's de entrada
1312 1317 '''
1313 1318 self.__initime = None
1314 1319 self.__lastdatatime = 0
1315 1320 self.__dataReady = False
1316 1321 self.__buffer = 0
1317 1322 self.__buffer2 = []
1318 1323 self.__buffer3 = 0
1319 1324 self.__profIndex = 0
1320 1325
1321 1326 self.__nch = dataOut.nChannels
1322 1327 self.__nHeis = dataOut.nHeights
1328 self.removeDC = removeDC
1329 self.lambda_ = 3.0e8/(9345.0e6)
1330 self.ippSec = dataOut.ippSeconds
1331 print("IPPseconds",dataOut.ippSeconds)
1323 1332
1324 1333 print("ELVALOR DE n es:", n)
1325 1334 if n == None:
1326 1335 raise ValueError("n should be specified.")
1327 1336
1328 1337 if n != None:
1329 1338 if n<2:
1330 1339 raise ValueError("n should be greater than 2")
1331 1340
1332 1341 self.n = n
1333 1342 self.__nProf = n
1334 '''
1335 if overlapping:
1336 self.__withOverlapping = True
1337 self.__buffer = None
1338 1343
1339 else:
1340 #print ("estoy sin __withO")
1341 self.__withOverlapping = False
1342 self.__buffer = 0
1343 self.__buffer2 = []
1344 self.__buffer3 = 0
1345 '''
1344 self.__buffer = numpy.zeros((dataOut.nChannels,
1345 n,
1346 dataOut.nHeights),
1347 dtype='complex')
1348
1349
1346 1350
1347 1351 def putData(self,data):
1348 1352 '''
1349 1353 Add a profile to he __buffer and increase in one the __profiel Index
1350 1354 '''
1351 #print("self.__profIndex :",self.__profIndex)
1352 self.__buffer += data*numpy.conjugate(data)
1353 self.__buffer2.append(numpy.conjugate(data))
1354 if self.__profIndex > 0:
1355 self.__buffer3 += self.__buffer2[self.__profIndex-1]*data
1356 self.__profIndex += 1
1357 return
1358 '''
1359 if not self.__withOverlapping:
1360 #print("Putdata inside over")
1361 self.__buffer += data* numpy.conjugate(data)
1362 self.__buffer2.append(numpy.conjugate(data))
1363
1364 if self.__profIndex >0:
1365 self.__buffer3 += self.__buffer2[self.__profIndex-1]*data
1355 self.__buffer[:,self.__profIndex,:]= data
1366 1356 self.__profIndex += 1
1367 1357 return
1368 1358
1369 if self.__buffer is None:
1370 #print("aqui bro")
1371 self.__buffer = data* numpy.conjugate(data)
1372 self.__buffer2.append(numpy.conjugate(data))
1373 self.__profIndex += 1
1374
1375 return
1376
1377 if self.__profIndex < self.n:
1378 self.__buffer = numpy.vstack(self.__buffer,data* numpy.conjugate(data))
1379 self.__buffer2.append(numpy.conjugate(data))
1380
1381 if self.__profIndex == 1:
1382 self.__buffer3 = self.__buffer2[self.__profIndex -1] * data
1383 else:
1384 self.__buffer3 = numpy.vstack(self.__buffer3, self.__buffer2[self.profIndex-1]*data)
1385
1386 self.__profIndex += 1
1387 return
1388 '''
1389
1390 1359 def pushData(self):
1391 1360 '''
1392 1361 Return the PULSEPAIR and the profiles used in the operation
1393 1362 Affected : self.__profileIndex
1394 1363 '''
1395 #print("************************************************")
1396 #print("push data int vel n")
1397 data_intensity = self.__buffer/self.n
1398 data_velocity = self.__buffer3/(self.n-1)
1399 n = self.__profIndex
1400 1364
1401 self.__buffer = 0
1402 self.__buffer2 = []
1403 self.__buffer3 = 0
1404 self.__profIndex = 0
1405
1406 return data_intensity, data_velocity,n
1407 '''
1408 if not self.__withOverlapping:
1409 #print("ahora que fue")
1410 data_intensity = self.__buffer/self.n
1411 data_velocity = self.__buffer3/(self.n-1)
1365 if self.removeDC==True:
1366 mean = numpy.mean(self.__buffer,1)
1367 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1368 dc= numpy.tile(tmp,[1,self.__nProf,1])
1369 self.__buffer = self.__buffer - dc
1370
1371 data_intensity = numpy.sum(self.__buffer*numpy.conj(self.__buffer),1)/self.n
1372 pair1 = self.__buffer[:,1:,:]*numpy.conjugate(self.__buffer[:,:-1,:])
1373 angle=numpy.angle(numpy.sum(pair1,1))*180/(math.pi)
1374 #print(angle.shape)#print("__ANGLE__") #print("angle",angle[:,:10])
1375 data_velocity = (self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(numpy.sum(pair1,1))
1412 1376 n = self.__profIndex
1413 1377
1414 self.__buffer = 0
1415 self.__buffer2 = []
1416 self.__buffer3 = 0
1378 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1417 1379 self.__profIndex = 0
1418 1380 return data_intensity, data_velocity,n
1419 1381
1420 data_intensity = numpy.sum(self.__buffer,axis = 0)
1421 data_velocity = numpy.sum(self.__buffer3,axis = 0)
1422 n = self.__profIndex
1423 #self.__buffer = 0
1424 #self.__buffer2 = []
1425 #self.__buffer3 = 0
1426 #self.__profIndex = 0
1427 return data_intensity, data_velocity,n
1428 '''
1429
1430 1382 def pulsePairbyProfiles(self,data):
1431 1383
1432 1384 self.__dataReady = False
1433 1385 data_intensity = None
1434 1386 data_velocity = None
1435 #print("beforeputada")
1436 1387 self.putData(data)
1437 #print("ProfileIndex:",self.__profIndex)
1438 1388 if self.__profIndex == self.n:
1439 1389 data_intensity, data_velocity, n = self.pushData()
1440 1390 self.__dataReady = True
1441 #print("-----------------------------------------------")
1442 #print("data_intensity",data_intensity.shape,"data_velocity",data_velocity.shape)
1391
1443 1392 return data_intensity, data_velocity
1444 1393
1445 1394 def pulsePairOp(self, data, datatime= None):
1446 1395
1447 1396 if self.__initime == None:
1448 1397 self.__initime = datatime
1449 1398
1450 1399 data_intensity, data_velocity = self.pulsePairbyProfiles(data)
1451 1400 self.__lastdatatime = datatime
1452 1401
1453 1402 if data_intensity is None:
1454 1403 return None, None, None
1455 1404
1456 1405 avgdatatime = self.__initime
1457 1406 deltatime = datatime - self.__lastdatatime
1458 1407 self.__initime = datatime
1459 '''
1460 if not self.__withOverlapping:
1461 self.__initime = datatime
1462 else:
1463 self.__initime += deltatime
1464 '''
1408
1465 1409 return data_intensity, data_velocity, avgdatatime
1466 1410
1467 def run(self, dataOut,n = None, overlapping= False,**kwargs):
1411 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1468 1412
1469 1413 if not self.isConfig:
1470 self.setup(dataOut = dataOut, n = n , **kwargs)
1414 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1471 1415 self.isConfig = True
1472 1416 #print("*******************")
1473 1417 #print("print Shape input data:",dataOut.data.shape)
1474 1418 data_intensity, data_velocity, avgdatatime = self.pulsePairOp(dataOut.data, dataOut.utctime)
1475 1419 dataOut.flagNoData = True
1476 1420
1477 1421 if self.__dataReady:
1478 1422 #print("#------------------------------------------------------")
1479 1423 #print("data_ready",data_intensity.shape)
1480 1424 dataOut.data = data_intensity #valor para plotear RTI
1481 1425 dataOut.nCohInt *= self.n
1482 1426 dataOut.data_intensity = data_intensity #valor para intensidad
1483 1427 dataOut.data_velocity = data_velocity #valor para velocidad
1484 1428 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1485 1429 dataOut.utctime = avgdatatime
1486 1430 dataOut.flagNoData = False
1487 1431 return dataOut
1488
1489 # import collections
1490 # from scipy.stats import mode
1491 #
1492 # class Synchronize(Operation):
1493 #
1494 # isConfig = False
1495 # __profIndex = 0
1496 #
1497 # def __init__(self, **kwargs):
1498 #
1499 # Operation.__init__(self, **kwargs)
1500 # # self.isConfig = False
1501 # self.__powBuffer = None
1502 # self.__startIndex = 0
1503 # self.__pulseFound = False
1504 #
1505 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1506 #
1507 # #Read data
1508 #
1509 # powerdB = dataOut.getPower(channel = channel)
1510 # noisedB = dataOut.getNoise(channel = channel)[0]
1511 #
1512 # self.__powBuffer.extend(powerdB.flatten())
1513 #
1514 # dataArray = numpy.array(self.__powBuffer)
1515 #
1516 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1517 #
1518 # maxValue = numpy.nanmax(filteredPower)
1519 #
1520 # if maxValue < noisedB + 10:
1521 # #No se encuentra ningun pulso de transmision
1522 # return None
1523 #
1524 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1525 #
1526 # if len(maxValuesIndex) < 2:
1527 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1528 # return None
1529 #
1530 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1531 #
1532 # #Seleccionar solo valores con un espaciamiento de nSamples
1533 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1534 #
1535 # if len(pulseIndex) < 2:
1536 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1537 # return None
1538 #
1539 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1540 #
1541 # #remover senales que se distancien menos de 10 unidades o muestras
1542 # #(No deberian existir IPP menor a 10 unidades)
1543 #
1544 # realIndex = numpy.where(spacing > 10 )[0]
1545 #
1546 # if len(realIndex) < 2:
1547 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1548 # return None
1549 #
1550 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1551 # realPulseIndex = pulseIndex[realIndex]
1552 #
1553 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1554 #
1555 # print "IPP = %d samples" %period
1556 #
1557 # self.__newNSamples = dataOut.nHeights #int(period)
1558 # self.__startIndex = int(realPulseIndex[0])
1559 #
1560 # return 1
1561 #
1562 #
1563 # def setup(self, nSamples, nChannels, buffer_size = 4):
1564 #
1565 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1566 # maxlen = buffer_size*nSamples)
1567 #
1568 # bufferList = []
1569 #
1570 # for i in range(nChannels):
1571 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1572 # maxlen = buffer_size*nSamples)
1573 #
1574 # bufferList.append(bufferByChannel)
1575 #
1576 # self.__nSamples = nSamples
1577 # self.__nChannels = nChannels
1578 # self.__bufferList = bufferList
1579 #
1580 # def run(self, dataOut, channel = 0):
1581 #
1582 # if not self.isConfig:
1583 # nSamples = dataOut.nHeights
1584 # nChannels = dataOut.nChannels
1585 # self.setup(nSamples, nChannels)
1586 # self.isConfig = True
1587 #
1588 # #Append new data to internal buffer
1589 # for thisChannel in range(self.__nChannels):
1590 # bufferByChannel = self.__bufferList[thisChannel]
1591 # bufferByChannel.extend(dataOut.data[thisChannel])
1592 #
1593 # if self.__pulseFound:
1594 # self.__startIndex -= self.__nSamples
1595 #
1596 # #Finding Tx Pulse
1597 # if not self.__pulseFound:
1598 # indexFound = self.__findTxPulse(dataOut, channel)
1599 #
1600 # if indexFound == None:
1601 # dataOut.flagNoData = True
1602 # return
1603 #
1604 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1605 # self.__pulseFound = True
1606 # self.__startIndex = indexFound
1607 #
1608 # #If pulse was found ...
1609 # for thisChannel in range(self.__nChannels):
1610 # bufferByChannel = self.__bufferList[thisChannel]
1611 # #print self.__startIndex
1612 # x = numpy.array(bufferByChannel)
1613 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1614 #
1615 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1616 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1617 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1618 #
1619 # dataOut.data = self.__arrayBuffer
1620 #
1621 # self.__startIndex += self.__newNSamples
1622 #
1623 # return
General Comments 0
You need to be logged in to leave comments. Login now