##// END OF EJS Templates
schain-jc
RSS
Remove C extension
Juan C. Espinoza -
r1176:d4dec4c4d8ae
parent child
Show More
Side by Side Unified
Context file:
r1176:d4dec4c4d8ae Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1251 +1,1234
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
11 11 from .jroheaderIO import SystemHeader, RadarControllerHeader
12 # from schainpy import cSchain
13 12
14 13
15 14 def getNumpyDtype(dataTypeCode):
16 15
17 16 if dataTypeCode == 0:
18 17 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
19 18 elif dataTypeCode == 1:
20 19 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
21 20 elif dataTypeCode == 2:
22 21 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
23 22 elif dataTypeCode == 3:
24 23 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
25 24 elif dataTypeCode == 4:
26 25 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
27 26 elif dataTypeCode == 5:
28 27 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
29 28 else:
30 29 raise ValueError('dataTypeCode was not defined')
31 30
32 31 return numpyDtype
33 32
34 33
35 34 def getDataTypeCode(numpyDtype):
36 35
37 36 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
38 37 datatype = 0
39 38 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
40 39 datatype = 1
41 40 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
42 41 datatype = 2
43 42 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
44 43 datatype = 3
45 44 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
46 45 datatype = 4
47 46 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
48 47 datatype = 5
49 48 else:
50 49 datatype = None
51 50
52 51 return datatype
53 52
54 53
55 54 def hildebrand_sekhon(data, navg):
56 55 """
57 56 This method is for the objective determination of the noise level in Doppler spectra. This
58 57 implementation technique is based on the fact that the standard deviation of the spectral
59 58 densities is equal to the mean spectral density for white Gaussian noise
60 59
61 60 Inputs:
62 61 Data : heights
63 62 navg : numbers of averages
64 63
65 64 Return:
66 -1 : any error
67 anoise : noise's level
65 mean : noise's level
68 66 """
69 67
70 sortdata = numpy.sort(data, axis=None)
71 lenOfData = len(sortdata)
72 nums_min = lenOfData*0.2
73
74 if nums_min <= 5:
75 nums_min = 5
76
77 sump = 0.
78
79 sumq = 0.
80
81 j = 0
82
83 cont = 1
84
85 while((cont==1)and(j<lenOfData)):
86
87 sump += sortdata[j]
88
89 sumq += sortdata[j]**2
90
91 if j > nums_min:
92 rtest = float(j)/(j-1) + 1.0/navg
93 if ((sumq*j) > (rtest*sump**2)):
94 j = j - 1
95 sump = sump - sortdata[j]
96 sumq = sumq - sortdata[j]**2
97 cont = 0
98
99 j += 1
100
101 lnoise = sump /j
102
103 return lnoise
68 sorted_spectrum = numpy.sort(data, axis=None)
69 nnoise = len(sorted_spectrum) # default to all points in the spectrum as noise
70 for npts in range(1, len(sorted_spectrum)+1):
71 partial = sorted_spectrum[:npts]
72 mean = partial.mean()
73 var = partial.var()
74 if var * navg < mean**2.:
75 nnoise = npts
76 else:
77 # partial spectrum no longer has characteristics of white noise
78 break
104 79
105 # return cSchain.hildebrand_sekhon(sortdata, navg)
80 noise_spectrum = sorted_spectrum[:nnoise]
81 mean = noise_spectrum.mean()
82 return mean
106 83
107 84
108 85 class Beam:
109 86
110 87 def __init__(self):
111 88 self.codeList = []
112 89 self.azimuthList = []
113 90 self.zenithList = []
114 91
115 92
116 93 class GenericData(object):
117 94
118 95 flagNoData = True
119 96
120 97 def copy(self, inputObj=None):
121 98
122 99 if inputObj == None:
123 100 return copy.deepcopy(self)
124 101
125 102 for key in list(inputObj.__dict__.keys()):
126 103
127 104 attribute = inputObj.__dict__[key]
128 105
129 106 # If this attribute is a tuple or list
130 107 if type(inputObj.__dict__[key]) in (tuple, list):
131 108 self.__dict__[key] = attribute[:]
132 109 continue
133 110
134 111 # If this attribute is another object or instance
135 112 if hasattr(attribute, '__dict__'):
136 113 self.__dict__[key] = attribute.copy()
137 114 continue
138 115
139 116 self.__dict__[key] = inputObj.__dict__[key]
140 117
141 118 def deepcopy(self):
142 119
143 120 return copy.deepcopy(self)
144 121
145 122 def isEmpty(self):
146 123
147 124 return self.flagNoData
148 125
149 126
150 127 class JROData(GenericData):
151 128
152 129 # m_BasicHeader = BasicHeader()
153 130 # m_ProcessingHeader = ProcessingHeader()
154 131
155 132 systemHeaderObj = SystemHeader()
156 133
157 134 radarControllerHeaderObj = RadarControllerHeader()
158 135
159 136 # data = None
160 137
161 138 type = None
162 139
163 140 datatype = None # dtype but in string
164 141
165 142 # dtype = None
166 143
167 144 # nChannels = None
168 145
169 146 # nHeights = None
170 147
171 148 nProfiles = None
172 149
173 150 heightList = None
174 151
175 152 channelList = None
176 153
177 154 flagDiscontinuousBlock = False
178 155
179 156 useLocalTime = False
180 157
181 158 utctime = None
182 159
183 160 timeZone = None
184 161
185 162 dstFlag = None
186 163
187 164 errorCount = None
188 165
189 166 blocksize = None
190 167
191 168 # nCode = None
192 169 #
193 170 # nBaud = None
194 171 #
195 172 # code = None
196 173
197 174 flagDecodeData = False # asumo q la data no esta decodificada
198 175
199 176 flagDeflipData = False # asumo q la data no esta sin flip
200 177
201 178 flagShiftFFT = False
202 179
203 180 # ippSeconds = None
204 181
205 182 # timeInterval = None
206 183
207 184 nCohInt = None
208 185
209 186 # noise = None
210 187
211 188 windowOfFilter = 1
212 189
213 190 # Speed of ligth
214 191 C = 3e8
215 192
216 193 frequency = 49.92e6
217 194
218 195 realtime = False
219 196
220 197 beacon_heiIndexList = None
221 198
222 199 last_block = None
223 200
224 201 blocknow = None
225 202
226 203 azimuth = None
227 204
228 205 zenith = None
229 206
230 207 beam = Beam()
231 208
232 209 profileIndex = None
233 210
211 error = (0, '')
212
213 def __str__(self):
214
215 return '{} - {}'.format(self.type, self.getDatatime())
216
234 217 def getNoise(self):
235 218
236 219 raise NotImplementedError
237 220
238 221 def getNChannels(self):
239 222
240 223 return len(self.channelList)
241 224
242 225 def getChannelIndexList(self):
243 226
244 227 return list(range(self.nChannels))
245 228
246 229 def getNHeights(self):
247 230
248 231 return len(self.heightList)
249 232
250 233 def getHeiRange(self, extrapoints=0):
251 234
252 235 heis = self.heightList
253 236 # deltah = self.heightList[1] - self.heightList[0]
254 237 #
255 238 # heis.append(self.heightList[-1])
256 239
257 240 return heis
258 241
259 242 def getDeltaH(self):
260 243
261 244 delta = self.heightList[1] - self.heightList[0]
262 245
263 246 return delta
264 247
265 248 def getltctime(self):
266 249
267 250 if self.useLocalTime:
268 251 return self.utctime - self.timeZone * 60
269 252
270 253 return self.utctime
271 254
272 255 def getDatatime(self):
273 256
274 257 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
275 258 return datatimeValue
276 259
277 260 def getTimeRange(self):
278 261
279 262 datatime = []
280 263
281 264 datatime.append(self.ltctime)
282 265 datatime.append(self.ltctime + self.timeInterval + 1)
283 266
284 267 datatime = numpy.array(datatime)
285 268
286 269 return datatime
287 270
288 271 def getFmaxTimeResponse(self):
289 272
290 273 period = (10**-6) * self.getDeltaH() / (0.15)
291 274
292 275 PRF = 1. / (period * self.nCohInt)
293 276
294 277 fmax = PRF
295 278
296 279 return fmax
297 280
298 281 def getFmax(self):
299 282 PRF = 1. / (self.ippSeconds * self.nCohInt)
300 283
301 284 fmax = PRF
302 285 return fmax
303 286
304 287 def getVmax(self):
305 288
306 289 _lambda = self.C / self.frequency
307 290
308 291 vmax = self.getFmax() * _lambda / 2
309 292
310 293 return vmax
311 294
312 295 def get_ippSeconds(self):
313 296 '''
314 297 '''
315 298 return self.radarControllerHeaderObj.ippSeconds
316 299
317 300 def set_ippSeconds(self, ippSeconds):
318 301 '''
319 302 '''
320 303
321 304 self.radarControllerHeaderObj.ippSeconds = ippSeconds
322 305
323 306 return
324 307
325 308 def get_dtype(self):
326 309 '''
327 310 '''
328 311 return getNumpyDtype(self.datatype)
329 312
330 313 def set_dtype(self, numpyDtype):
331 314 '''
332 315 '''
333 316
334 317 self.datatype = getDataTypeCode(numpyDtype)
335 318
336 319 def get_code(self):
337 320 '''
338 321 '''
339 322 return self.radarControllerHeaderObj.code
340 323
341 324 def set_code(self, code):
342 325 '''
343 326 '''
344 327 self.radarControllerHeaderObj.code = code
345 328
346 329 return
347 330
348 331 def get_ncode(self):
349 332 '''
350 333 '''
351 334 return self.radarControllerHeaderObj.nCode
352 335
353 336 def set_ncode(self, nCode):
354 337 '''
355 338 '''
356 339 self.radarControllerHeaderObj.nCode = nCode
357 340
358 341 return
359 342
360 343 def get_nbaud(self):
361 344 '''
362 345 '''
363 346 return self.radarControllerHeaderObj.nBaud
364 347
365 348 def set_nbaud(self, nBaud):
366 349 '''
367 350 '''
368 351 self.radarControllerHeaderObj.nBaud = nBaud
369 352
370 353 return
371 354
372 355 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
373 356 channelIndexList = property(
374 357 getChannelIndexList, "I'm the 'channelIndexList' property.")
375 358 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
376 359 #noise = property(getNoise, "I'm the 'nHeights' property.")
377 360 datatime = property(getDatatime, "I'm the 'datatime' property")
378 361 ltctime = property(getltctime, "I'm the 'ltctime' property")
379 362 ippSeconds = property(get_ippSeconds, set_ippSeconds)
380 363 dtype = property(get_dtype, set_dtype)
381 364 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
382 365 code = property(get_code, set_code)
383 366 nCode = property(get_ncode, set_ncode)
384 367 nBaud = property(get_nbaud, set_nbaud)
385 368
386 369
387 370 class Voltage(JROData):
388 371
389 372 # data es un numpy array de 2 dmensiones (canales, alturas)
390 373 data = None
391 374
392 375 def __init__(self):
393 376 '''
394 377 Constructor
395 378 '''
396 379
397 380 self.useLocalTime = True
398 381
399 382 self.radarControllerHeaderObj = RadarControllerHeader()
400 383
401 384 self.systemHeaderObj = SystemHeader()
402 385
403 386 self.type = "Voltage"
404 387
405 388 self.data = None
406 389
407 390 # self.dtype = None
408 391
409 392 # self.nChannels = 0
410 393
411 394 # self.nHeights = 0
412 395
413 396 self.nProfiles = None
414 397
415 398 self.heightList = None
416 399
417 400 self.channelList = None
418 401
419 402 # self.channelIndexList = None
420 403
421 404 self.flagNoData = True
422 405
423 406 self.flagDiscontinuousBlock = False
424 407
425 408 self.utctime = None
426 409
427 410 self.timeZone = None
428 411
429 412 self.dstFlag = None
430 413
431 414 self.errorCount = None
432 415
433 416 self.nCohInt = None
434 417
435 418 self.blocksize = None
436 419
437 420 self.flagDecodeData = False # asumo q la data no esta decodificada
438 421
439 422 self.flagDeflipData = False # asumo q la data no esta sin flip
440 423
441 424 self.flagShiftFFT = False
442 425
443 426 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
444 427
445 428 self.profileIndex = 0
446 429
447 430 def getNoisebyHildebrand(self, channel=None):
448 431 """
449 432 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
450 433
451 434 Return:
452 435 noiselevel
453 436 """
454 437
455 438 if channel != None:
456 439 data = self.data[channel]
457 440 nChannels = 1
458 441 else:
459 442 data = self.data
460 443 nChannels = self.nChannels
461 444
462 445 noise = numpy.zeros(nChannels)
463 446 power = data * numpy.conjugate(data)
464 447
465 448 for thisChannel in range(nChannels):
466 449 if nChannels == 1:
467 450 daux = power[:].real
468 451 else:
469 452 daux = power[thisChannel, :].real
470 453 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
471 454
472 455 return noise
473 456
474 457 def getNoise(self, type=1, channel=None):
475 458
476 459 if type == 1:
477 460 noise = self.getNoisebyHildebrand(channel)
478 461
479 462 return noise
480 463
481 464 def getPower(self, channel=None):
482 465
483 466 if channel != None:
484 467 data = self.data[channel]
485 468 else:
486 469 data = self.data
487 470
488 471 power = data * numpy.conjugate(data)
489 472 powerdB = 10 * numpy.log10(power.real)
490 473 powerdB = numpy.squeeze(powerdB)
491 474
492 475 return powerdB
493 476
494 477 def getTimeInterval(self):
495 478
496 479 timeInterval = self.ippSeconds * self.nCohInt
497 480
498 481 return timeInterval
499 482
500 483 noise = property(getNoise, "I'm the 'nHeights' property.")
501 484 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
502 485
503 486
504 487 class Spectra(JROData):
505 488
506 489 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
507 490 data_spc = None
508 491
509 492 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
510 493 data_cspc = None
511 494
512 495 # data dc es un numpy array de 2 dmensiones (canales, alturas)
513 496 data_dc = None
514 497
515 498 # data power
516 499 data_pwr = None
517 500
518 501 nFFTPoints = None
519 502
520 503 # nPairs = None
521 504
522 505 pairsList = None
523 506
524 507 nIncohInt = None
525 508
526 509 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
527 510
528 511 nCohInt = None # se requiere para determinar el valor de timeInterval
529 512
530 513 ippFactor = None
531 514
532 515 profileIndex = 0
533 516
534 517 plotting = "spectra"
535 518
536 519 def __init__(self):
537 520 '''
538 521 Constructor
539 522 '''
540 523
541 524 self.useLocalTime = True
542 525
543 526 self.radarControllerHeaderObj = RadarControllerHeader()
544 527
545 528 self.systemHeaderObj = SystemHeader()
546 529
547 530 self.type = "Spectra"
548 531
549 532 # self.data = None
550 533
551 534 # self.dtype = None
552 535
553 536 # self.nChannels = 0
554 537
555 538 # self.nHeights = 0
556 539
557 540 self.nProfiles = None
558 541
559 542 self.heightList = None
560 543
561 544 self.channelList = None
562 545
563 546 # self.channelIndexList = None
564 547
565 548 self.pairsList = None
566 549
567 550 self.flagNoData = True
568 551
569 552 self.flagDiscontinuousBlock = False
570 553
571 554 self.utctime = None
572 555
573 556 self.nCohInt = None
574 557
575 558 self.nIncohInt = None
576 559
577 560 self.blocksize = None
578 561
579 562 self.nFFTPoints = None
580 563
581 564 self.wavelength = None
582 565
583 566 self.flagDecodeData = False # asumo q la data no esta decodificada
584 567
585 568 self.flagDeflipData = False # asumo q la data no esta sin flip
586 569
587 570 self.flagShiftFFT = False
588 571
589 572 self.ippFactor = 1
590 573
591 574 #self.noise = None
592 575
593 576 self.beacon_heiIndexList = []
594 577
595 578 self.noise_estimation = None
596 579
597 580 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
598 581 """
599 582 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
600 583
601 584 Return:
602 585 noiselevel
603 586 """
604 587
605 588 noise = numpy.zeros(self.nChannels)
606 589
607 590 for channel in range(self.nChannels):
608 591 daux = self.data_spc[channel,
609 592 xmin_index:xmax_index, ymin_index:ymax_index]
610 593 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
611 594
612 595 return noise
613 596
614 597 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
615 598
616 599 if self.noise_estimation is not None:
617 600 # this was estimated by getNoise Operation defined in jroproc_spectra.py
618 601 return self.noise_estimation
619 602 else:
620 603 noise = self.getNoisebyHildebrand(
621 604 xmin_index, xmax_index, ymin_index, ymax_index)
622 605 return noise
623 606
624 607 def getFreqRangeTimeResponse(self, extrapoints=0):
625 608
626 609 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
627 610 freqrange = deltafreq * \
628 611 (numpy.arange(self.nFFTPoints + extrapoints) -
629 612 self.nFFTPoints / 2.) - deltafreq / 2
630 613
631 614 return freqrange
632 615
633 616 def getAcfRange(self, extrapoints=0):
634 617
635 618 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
636 619 freqrange = deltafreq * \
637 620 (numpy.arange(self.nFFTPoints + extrapoints) -
638 621 self.nFFTPoints / 2.) - deltafreq / 2
639 622
640 623 return freqrange
641 624
642 625 def getFreqRange(self, extrapoints=0):
643 626
644 627 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
645 628 freqrange = deltafreq * \
646 629 (numpy.arange(self.nFFTPoints + extrapoints) -
647 630 self.nFFTPoints / 2.) - deltafreq / 2
648 631
649 632 return freqrange
650 633
651 634 def getVelRange(self, extrapoints=0):
652 635
653 636 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
654 637 velrange = deltav * (numpy.arange(self.nFFTPoints +
655 638 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
656 639
657 640 return velrange
658 641
659 642 def getNPairs(self):
660 643
661 644 return len(self.pairsList)
662 645
663 646 def getPairsIndexList(self):
664 647
665 648 return list(range(self.nPairs))
666 649
667 650 def getNormFactor(self):
668 651
669 652 pwcode = 1
670 653
671 654 if self.flagDecodeData:
672 655 pwcode = numpy.sum(self.code[0]**2)
673 656 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
674 657 normFactor = self.nProfiles * self.nIncohInt * \
675 658 self.nCohInt * pwcode * self.windowOfFilter
676 659
677 660 return normFactor
678 661
679 662 def getFlagCspc(self):
680 663
681 664 if self.data_cspc is None:
682 665 return True
683 666
684 667 return False
685 668
686 669 def getFlagDc(self):
687 670
688 671 if self.data_dc is None:
689 672 return True
690 673
691 674 return False
692 675
693 676 def getTimeInterval(self):
694 677
695 678 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
696 679
697 680 return timeInterval
698 681
699 682 def getPower(self):
700 683
701 684 factor = self.normFactor
702 685 z = self.data_spc / factor
703 686 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
704 687 avg = numpy.average(z, axis=1)
705 688
706 689 return 10 * numpy.log10(avg)
707 690
708 691 def getCoherence(self, pairsList=None, phase=False):
709 692
710 693 z = []
711 694 if pairsList is None:
712 695 pairsIndexList = self.pairsIndexList
713 696 else:
714 697 pairsIndexList = []
715 698 for pair in pairsList:
716 699 if pair not in self.pairsList:
717 700 raise ValueError("Pair %s is not in dataOut.pairsList" % (
718 701 pair))
719 702 pairsIndexList.append(self.pairsList.index(pair))
720 703 for i in range(len(pairsIndexList)):
721 704 pair = self.pairsList[pairsIndexList[i]]
722 705 ccf = numpy.average(
723 706 self.data_cspc[pairsIndexList[i], :, :], axis=0)
724 707 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
725 708 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
726 709 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
727 710 if phase:
728 711 data = numpy.arctan2(avgcoherenceComplex.imag,
729 712 avgcoherenceComplex.real) * 180 / numpy.pi
730 713 else:
731 714 data = numpy.abs(avgcoherenceComplex)
732 715
733 716 z.append(data)
734 717
735 718 return numpy.array(z)
736 719
737 720 def setValue(self, value):
738 721
739 722 print("This property should not be initialized")
740 723
741 724 return
742 725
743 726 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
744 727 pairsIndexList = property(
745 728 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
746 729 normFactor = property(getNormFactor, setValue,
747 730 "I'm the 'getNormFactor' property.")
748 731 flag_cspc = property(getFlagCspc, setValue)
749 732 flag_dc = property(getFlagDc, setValue)
750 733 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
751 734 timeInterval = property(getTimeInterval, setValue,
752 735 "I'm the 'timeInterval' property")
753 736
754 737
755 738 class SpectraHeis(Spectra):
756 739
757 740 data_spc = None
758 741
759 742 data_cspc = None
760 743
761 744 data_dc = None
762 745
763 746 nFFTPoints = None
764 747
765 748 # nPairs = None
766 749
767 750 pairsList = None
768 751
769 752 nCohInt = None
770 753
771 754 nIncohInt = None
772 755
773 756 def __init__(self):
774 757
775 758 self.radarControllerHeaderObj = RadarControllerHeader()
776 759
777 760 self.systemHeaderObj = SystemHeader()
778 761
779 762 self.type = "SpectraHeis"
780 763
781 764 # self.dtype = None
782 765
783 766 # self.nChannels = 0
784 767
785 768 # self.nHeights = 0
786 769
787 770 self.nProfiles = None
788 771
789 772 self.heightList = None
790 773
791 774 self.channelList = None
792 775
793 776 # self.channelIndexList = None
794 777
795 778 self.flagNoData = True
796 779
797 780 self.flagDiscontinuousBlock = False
798 781
799 782 # self.nPairs = 0
800 783
801 784 self.utctime = None
802 785
803 786 self.blocksize = None
804 787
805 788 self.profileIndex = 0
806 789
807 790 self.nCohInt = 1
808 791
809 792 self.nIncohInt = 1
810 793
811 794 def getNormFactor(self):
812 795 pwcode = 1
813 796 if self.flagDecodeData:
814 797 pwcode = numpy.sum(self.code[0]**2)
815 798
816 799 normFactor = self.nIncohInt * self.nCohInt * pwcode
817 800
818 801 return normFactor
819 802
820 803 def getTimeInterval(self):
821 804
822 805 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
823 806
824 807 return timeInterval
825 808
826 809 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
827 810 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
828 811
829 812
830 813 class Fits(JROData):
831 814
832 815 heightList = None
833 816
834 817 channelList = None
835 818
836 819 flagNoData = True
837 820
838 821 flagDiscontinuousBlock = False
839 822
840 823 useLocalTime = False
841 824
842 825 utctime = None
843 826
844 827 timeZone = None
845 828
846 829 # ippSeconds = None
847 830
848 831 # timeInterval = None
849 832
850 833 nCohInt = None
851 834
852 835 nIncohInt = None
853 836
854 837 noise = None
855 838
856 839 windowOfFilter = 1
857 840
858 841 # Speed of ligth
859 842 C = 3e8
860 843
861 844 frequency = 49.92e6
862 845
863 846 realtime = False
864 847
865 848 def __init__(self):
866 849
867 850 self.type = "Fits"
868 851
869 852 self.nProfiles = None
870 853
871 854 self.heightList = None
872 855
873 856 self.channelList = None
874 857
875 858 # self.channelIndexList = None
876 859
877 860 self.flagNoData = True
878 861
879 862 self.utctime = None
880 863
881 864 self.nCohInt = 1
882 865
883 866 self.nIncohInt = 1
884 867
885 868 self.useLocalTime = True
886 869
887 870 self.profileIndex = 0
888 871
889 872 # self.utctime = None
890 873 # self.timeZone = None
891 874 # self.ltctime = None
892 875 # self.timeInterval = None
893 876 # self.header = None
894 877 # self.data_header = None
895 878 # self.data = None
896 879 # self.datatime = None
897 880 # self.flagNoData = False
898 881 # self.expName = ''
899 882 # self.nChannels = None
900 883 # self.nSamples = None
901 884 # self.dataBlocksPerFile = None
902 885 # self.comments = ''
903 886 #
904 887
905 888 def getltctime(self):
906 889
907 890 if self.useLocalTime:
908 891 return self.utctime - self.timeZone * 60
909 892
910 893 return self.utctime
911 894
912 895 def getDatatime(self):
913 896
914 897 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
915 898 return datatime
916 899
917 900 def getTimeRange(self):
918 901
919 902 datatime = []
920 903
921 904 datatime.append(self.ltctime)
922 905 datatime.append(self.ltctime + self.timeInterval)
923 906
924 907 datatime = numpy.array(datatime)
925 908
926 909 return datatime
927 910
928 911 def getHeiRange(self):
929 912
930 913 heis = self.heightList
931 914
932 915 return heis
933 916
934 917 def getNHeights(self):
935 918
936 919 return len(self.heightList)
937 920
938 921 def getNChannels(self):
939 922
940 923 return len(self.channelList)
941 924
942 925 def getChannelIndexList(self):
943 926
944 927 return list(range(self.nChannels))
945 928
946 929 def getNoise(self, type=1):
947 930
948 931 #noise = numpy.zeros(self.nChannels)
949 932
950 933 if type == 1:
951 934 noise = self.getNoisebyHildebrand()
952 935
953 936 if type == 2:
954 937 noise = self.getNoisebySort()
955 938
956 939 if type == 3:
957 940 noise = self.getNoisebyWindow()
958 941
959 942 return noise
960 943
961 944 def getTimeInterval(self):
962 945
963 946 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
964 947
965 948 return timeInterval
966 949
967 950 datatime = property(getDatatime, "I'm the 'datatime' property")
968 951 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
969 952 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
970 953 channelIndexList = property(
971 954 getChannelIndexList, "I'm the 'channelIndexList' property.")
972 955 noise = property(getNoise, "I'm the 'nHeights' property.")
973 956
974 957 ltctime = property(getltctime, "I'm the 'ltctime' property")
975 958 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
976 959
977 960
978 961 class Correlation(JROData):
979 962
980 963 noise = None
981 964
982 965 SNR = None
983 966
984 967 #--------------------------------------------------
985 968
986 969 mode = None
987 970
988 971 split = False
989 972
990 973 data_cf = None
991 974
992 975 lags = None
993 976
994 977 lagRange = None
995 978
996 979 pairsList = None
997 980
998 981 normFactor = None
999 982
1000 983 #--------------------------------------------------
1001 984
1002 985 # calculateVelocity = None
1003 986
1004 987 nLags = None
1005 988
1006 989 nPairs = None
1007 990
1008 991 nAvg = None
1009 992
1010 993 def __init__(self):
1011 994 '''
1012 995 Constructor
1013 996 '''
1014 997 self.radarControllerHeaderObj = RadarControllerHeader()
1015 998
1016 999 self.systemHeaderObj = SystemHeader()
1017 1000
1018 1001 self.type = "Correlation"
1019 1002
1020 1003 self.data = None
1021 1004
1022 1005 self.dtype = None
1023 1006
1024 1007 self.nProfiles = None
1025 1008
1026 1009 self.heightList = None
1027 1010
1028 1011 self.channelList = None
1029 1012
1030 1013 self.flagNoData = True
1031 1014
1032 1015 self.flagDiscontinuousBlock = False
1033 1016
1034 1017 self.utctime = None
1035 1018
1036 1019 self.timeZone = None
1037 1020
1038 1021 self.dstFlag = None
1039 1022
1040 1023 self.errorCount = None
1041 1024
1042 1025 self.blocksize = None
1043 1026
1044 1027 self.flagDecodeData = False # asumo q la data no esta decodificada
1045 1028
1046 1029 self.flagDeflipData = False # asumo q la data no esta sin flip
1047 1030
1048 1031 self.pairsList = None
1049 1032
1050 1033 self.nPoints = None
1051 1034
1052 1035 def getPairsList(self):
1053 1036
1054 1037 return self.pairsList
1055 1038
1056 1039 def getNoise(self, mode=2):
1057 1040
1058 1041 indR = numpy.where(self.lagR == 0)[0][0]
1059 1042 indT = numpy.where(self.lagT == 0)[0][0]
1060 1043
1061 1044 jspectra0 = self.data_corr[:, :, indR, :]
1062 1045 jspectra = copy.copy(jspectra0)
1063 1046
1064 1047 num_chan = jspectra.shape[0]
1065 1048 num_hei = jspectra.shape[2]
1066 1049
1067 1050 freq_dc = jspectra.shape[1] / 2
1068 1051 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
1069 1052
1070 1053 if ind_vel[0] < 0:
1071 1054 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
1072 1055
1073 1056 if mode == 1:
1074 1057 jspectra[:, freq_dc, :] = (
1075 1058 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
1076 1059
1077 1060 if mode == 2:
1078 1061
1079 1062 vel = numpy.array([-2, -1, 1, 2])
1080 1063 xx = numpy.zeros([4, 4])
1081 1064
1082 1065 for fil in range(4):
1083 1066 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
1084 1067
1085 1068 xx_inv = numpy.linalg.inv(xx)
1086 1069 xx_aux = xx_inv[0, :]
1087 1070
1088 1071 for ich in range(num_chan):
1089 1072 yy = jspectra[ich, ind_vel, :]
1090 1073 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
1091 1074
1092 1075 junkid = jspectra[ich, freq_dc, :] <= 0
1093 1076 cjunkid = sum(junkid)
1094 1077
1095 1078 if cjunkid.any():
1096 1079 jspectra[ich, freq_dc, junkid.nonzero()] = (
1097 1080 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
1098 1081
1099 1082 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
1100 1083
1101 1084 return noise
1102 1085
1103 1086 def getTimeInterval(self):
1104 1087
1105 1088 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1106 1089
1107 1090 return timeInterval
1108 1091
1109 1092 def splitFunctions(self):
1110 1093
1111 1094 pairsList = self.pairsList
1112 1095 ccf_pairs = []
1113 1096 acf_pairs = []
1114 1097 ccf_ind = []
1115 1098 acf_ind = []
1116 1099 for l in range(len(pairsList)):
1117 1100 chan0 = pairsList[l][0]
1118 1101 chan1 = pairsList[l][1]
1119 1102
1120 1103 # Obteniendo pares de Autocorrelacion
1121 1104 if chan0 == chan1:
1122 1105 acf_pairs.append(chan0)
1123 1106 acf_ind.append(l)
1124 1107 else:
1125 1108 ccf_pairs.append(pairsList[l])
1126 1109 ccf_ind.append(l)
1127 1110
1128 1111 data_acf = self.data_cf[acf_ind]
1129 1112 data_ccf = self.data_cf[ccf_ind]
1130 1113
1131 1114 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1132 1115
1133 1116 def getNormFactor(self):
1134 1117 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1135 1118 acf_pairs = numpy.array(acf_pairs)
1136 1119 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1137 1120
1138 1121 for p in range(self.nPairs):
1139 1122 pair = self.pairsList[p]
1140 1123
1141 1124 ch0 = pair[0]
1142 1125 ch1 = pair[1]
1143 1126
1144 1127 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1145 1128 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1146 1129 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1147 1130
1148 1131 return normFactor
1149 1132
1150 1133 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1151 1134 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1152 1135
1153 1136
1154 1137 class Parameters(Spectra):
1155 1138
1156 1139 experimentInfo = None # Information about the experiment
1157 1140
1158 1141 # Information from previous data
1159 1142
1160 1143 inputUnit = None # Type of data to be processed
1161 1144
1162 1145 operation = None # Type of operation to parametrize
1163 1146
1164 1147 # normFactor = None #Normalization Factor
1165 1148
1166 1149 groupList = None # List of Pairs, Groups, etc
1167 1150
1168 1151 # Parameters
1169 1152
1170 1153 data_param = None # Parameters obtained
1171 1154
1172 1155 data_pre = None # Data Pre Parametrization
1173 1156
1174 1157 data_SNR = None # Signal to Noise Ratio
1175 1158
1176 1159 # heightRange = None #Heights
1177 1160
1178 1161 abscissaList = None # Abscissa, can be velocities, lags or time
1179 1162
1180 1163 # noise = None #Noise Potency
1181 1164
1182 1165 utctimeInit = None # Initial UTC time
1183 1166
1184 1167 paramInterval = None # Time interval to calculate Parameters in seconds
1185 1168
1186 1169 useLocalTime = True
1187 1170
1188 1171 # Fitting
1189 1172
1190 1173 data_error = None # Error of the estimation
1191 1174
1192 1175 constants = None
1193 1176
1194 1177 library = None
1195 1178
1196 1179 # Output signal
1197 1180
1198 1181 outputInterval = None # Time interval to calculate output signal in seconds
1199 1182
1200 1183 data_output = None # Out signal
1201 1184
1202 1185 nAvg = None
1203 1186
1204 1187 noise_estimation = None
1205 1188
1206 1189 GauSPC = None # Fit gaussian SPC
1207 1190
1208 1191 def __init__(self):
1209 1192 '''
1210 1193 Constructor
1211 1194 '''
1212 1195 self.radarControllerHeaderObj = RadarControllerHeader()
1213 1196
1214 1197 self.systemHeaderObj = SystemHeader()
1215 1198
1216 1199 self.type = "Parameters"
1217 1200
1218 1201 def getTimeRange1(self, interval):
1219 1202
1220 1203 datatime = []
1221 1204
1222 1205 if self.useLocalTime:
1223 1206 time1 = self.utctimeInit - self.timeZone * 60
1224 1207 else:
1225 1208 time1 = self.utctimeInit
1226 1209
1227 1210 datatime.append(time1)
1228 1211 datatime.append(time1 + interval)
1229 1212 datatime = numpy.array(datatime)
1230 1213
1231 1214 return datatime
1232 1215
1233 1216 def getTimeInterval(self):
1234 1217
1235 1218 if hasattr(self, 'timeInterval1'):
1236 1219 return self.timeInterval1
1237 1220 else:
1238 1221 return self.paramInterval
1239 1222
1240 1223 def setValue(self, value):
1241 1224
1242 1225 print("This property should not be initialized")
1243 1226
1244 1227 return
1245 1228
1246 1229 def getNoise(self):
1247 1230
1248 1231 return self.spc_noise
1249 1232
1250 1233 timeInterval = property(getTimeInterval)
1251 1234 noise = property(getNoise, setValue, "I'm the 'Noise' property.") No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1335 +1,1333
1 1 import sys
2 2 import numpy
3 3 from scipy import interpolate
4 #TODO
5 #from schainpy import cSchain
6 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
7 5 from schainpy.model.data.jrodata import Voltage
8 6 from schainpy.utils import log
9 7 from time import time
10 8
11 9
12 10 @MPDecorator
13 11 class VoltageProc(ProcessingUnit):
14 12
15 13 METHODS = {} #yong
16 14
17 15 def __init__(self):#, **kwargs): #yong
18 16
19 17 ProcessingUnit.__init__(self)#, **kwargs)
20 18
21 19 # self.objectDict = {}
22 20 self.dataOut = Voltage()
23 21 self.flip = 1
24 22 self.setupReq = False #yong
25 23
26 24 def run(self):
27 25
28 26 if self.dataIn.type == 'AMISR':
29 27 self.__updateObjFromAmisrInput()
30 28
31 29 if self.dataIn.type == 'Voltage':
32 30 self.dataOut.copy(self.dataIn)
33 31
34 32 # self.dataOut.copy(self.dataIn)
35 33
36 34 def __updateObjFromAmisrInput(self):
37 35
38 36 self.dataOut.timeZone = self.dataIn.timeZone
39 37 self.dataOut.dstFlag = self.dataIn.dstFlag
40 38 self.dataOut.errorCount = self.dataIn.errorCount
41 39 self.dataOut.useLocalTime = self.dataIn.useLocalTime
42 40
43 41 self.dataOut.flagNoData = self.dataIn.flagNoData
44 42 self.dataOut.data = self.dataIn.data
45 43 self.dataOut.utctime = self.dataIn.utctime
46 44 self.dataOut.channelList = self.dataIn.channelList
47 45 # self.dataOut.timeInterval = self.dataIn.timeInterval
48 46 self.dataOut.heightList = self.dataIn.heightList
49 47 self.dataOut.nProfiles = self.dataIn.nProfiles
50 48
51 49 self.dataOut.nCohInt = self.dataIn.nCohInt
52 50 self.dataOut.ippSeconds = self.dataIn.ippSeconds
53 51 self.dataOut.frequency = self.dataIn.frequency
54 52
55 53 self.dataOut.azimuth = self.dataIn.azimuth
56 54 self.dataOut.zenith = self.dataIn.zenith
57 55
58 56 self.dataOut.beam.codeList = self.dataIn.beam.codeList
59 57 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
60 58 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
61 59 #
62 60 # pass#
63 61 #
64 62 # def init(self):
65 63 #
66 64 #
67 65 # if self.dataIn.type == 'AMISR':
68 66 # self.__updateObjFromAmisrInput()
69 67 #
70 68 # if self.dataIn.type == 'Voltage':
71 69 # self.dataOut.copy(self.dataIn)
72 70 # # No necesita copiar en cada init() los atributos de dataIn
73 71 # # la copia deberia hacerse por cada nuevo bloque de datos
74 72
75 73 def selectChannels(self, channelList):
76 74
77 75 channelIndexList = []
78 76
79 77 for channel in channelList:
80 78 if channel not in self.dataOut.channelList:
81 79 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
82 80
83 81 index = self.dataOut.channelList.index(channel)
84 82 channelIndexList.append(index)
85 83
86 84 self.selectChannelsByIndex(channelIndexList)
87 85
88 86 def selectChannelsByIndex(self, channelIndexList):
89 87 """
90 88 Selecciona un bloque de datos en base a canales segun el channelIndexList
91 89
92 90 Input:
93 91 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
94 92
95 93 Affected:
96 94 self.dataOut.data
97 95 self.dataOut.channelIndexList
98 96 self.dataOut.nChannels
99 97 self.dataOut.m_ProcessingHeader.totalSpectra
100 98 self.dataOut.systemHeaderObj.numChannels
101 99 self.dataOut.m_ProcessingHeader.blockSize
102 100
103 101 Return:
104 102 None
105 103 """
106 104
107 105 for channelIndex in channelIndexList:
108 106 if channelIndex not in self.dataOut.channelIndexList:
109 107 print(channelIndexList)
110 108 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
111 109
112 110 if self.dataOut.flagDataAsBlock:
113 111 """
114 112 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
115 113 """
116 114 data = self.dataOut.data[channelIndexList,:,:]
117 115 else:
118 116 data = self.dataOut.data[channelIndexList,:]
119 117
120 118 self.dataOut.data = data
121 119 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
122 120 # self.dataOut.nChannels = nChannels
123 121
124 122 return 1
125 123
126 124 def selectHeights(self, minHei=None, maxHei=None):
127 125 """
128 126 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
129 127 minHei <= height <= maxHei
130 128
131 129 Input:
132 130 minHei : valor minimo de altura a considerar
133 131 maxHei : valor maximo de altura a considerar
134 132
135 133 Affected:
136 134 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
137 135
138 136 Return:
139 137 1 si el metodo se ejecuto con exito caso contrario devuelve 0
140 138 """
141 139
142 140 if minHei == None:
143 141 minHei = self.dataOut.heightList[0]
144 142
145 143 if maxHei == None:
146 144 maxHei = self.dataOut.heightList[-1]
147 145
148 146 if (minHei < self.dataOut.heightList[0]):
149 147 minHei = self.dataOut.heightList[0]
150 148
151 149 if (maxHei > self.dataOut.heightList[-1]):
152 150 maxHei = self.dataOut.heightList[-1]
153 151
154 152 minIndex = 0
155 153 maxIndex = 0
156 154 heights = self.dataOut.heightList
157 155
158 156 inda = numpy.where(heights >= minHei)
159 157 indb = numpy.where(heights <= maxHei)
160 158
161 159 try:
162 160 minIndex = inda[0][0]
163 161 except:
164 162 minIndex = 0
165 163
166 164 try:
167 165 maxIndex = indb[0][-1]
168 166 except:
169 167 maxIndex = len(heights)
170 168
171 169 self.selectHeightsByIndex(minIndex, maxIndex)
172 170
173 171 return 1
174 172
175 173
176 174 def selectHeightsByIndex(self, minIndex, maxIndex):
177 175 """
178 176 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
179 177 minIndex <= index <= maxIndex
180 178
181 179 Input:
182 180 minIndex : valor de indice minimo de altura a considerar
183 181 maxIndex : valor de indice maximo de altura a considerar
184 182
185 183 Affected:
186 184 self.dataOut.data
187 185 self.dataOut.heightList
188 186
189 187 Return:
190 188 1 si el metodo se ejecuto con exito caso contrario devuelve 0
191 189 """
192 190
193 191 if (minIndex < 0) or (minIndex > maxIndex):
194 192 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
195 193
196 194 if (maxIndex >= self.dataOut.nHeights):
197 195 maxIndex = self.dataOut.nHeights
198 196
199 197 #voltage
200 198 if self.dataOut.flagDataAsBlock:
201 199 """
202 200 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
203 201 """
204 202 data = self.dataOut.data[:,:, minIndex:maxIndex]
205 203 else:
206 204 data = self.dataOut.data[:, minIndex:maxIndex]
207 205
208 206 # firstHeight = self.dataOut.heightList[minIndex]
209 207
210 208 self.dataOut.data = data
211 209 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
212 210
213 211 if self.dataOut.nHeights <= 1:
214 212 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
215 213
216 214 return 1
217 215
218 216
219 217 def filterByHeights(self, window):
220 218
221 219 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
222 220
223 221 if window == None:
224 222 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
225 223
226 224 newdelta = deltaHeight * window
227 225 r = self.dataOut.nHeights % window
228 226 newheights = (self.dataOut.nHeights-r)/window
229 227
230 228 if newheights <= 1:
231 229 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window))
232 230
233 231 if self.dataOut.flagDataAsBlock:
234 232 """
235 233 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
236 234 """
237 235 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
238 236 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
239 237 buffer = numpy.sum(buffer,3)
240 238
241 239 else:
242 240 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
243 241 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
244 242 buffer = numpy.sum(buffer,2)
245 243
246 244 self.dataOut.data = buffer
247 245 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
248 246 self.dataOut.windowOfFilter = window
249 247
250 248 def setH0(self, h0, deltaHeight = None):
251 249
252 250 if not deltaHeight:
253 251 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
254 252
255 253 nHeights = self.dataOut.nHeights
256 254
257 255 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
258 256
259 257 self.dataOut.heightList = newHeiRange
260 258
261 259 def deFlip(self, channelList = []):
262 260
263 261 data = self.dataOut.data.copy()
264 262
265 263 if self.dataOut.flagDataAsBlock:
266 264 flip = self.flip
267 265 profileList = list(range(self.dataOut.nProfiles))
268 266
269 267 if not channelList:
270 268 for thisProfile in profileList:
271 269 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
272 270 flip *= -1.0
273 271 else:
274 272 for thisChannel in channelList:
275 273 if thisChannel not in self.dataOut.channelList:
276 274 continue
277 275
278 276 for thisProfile in profileList:
279 277 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
280 278 flip *= -1.0
281 279
282 280 self.flip = flip
283 281
284 282 else:
285 283 if not channelList:
286 284 data[:,:] = data[:,:]*self.flip
287 285 else:
288 286 for thisChannel in channelList:
289 287 if thisChannel not in self.dataOut.channelList:
290 288 continue
291 289
292 290 data[thisChannel,:] = data[thisChannel,:]*self.flip
293 291
294 292 self.flip *= -1.
295 293
296 294 self.dataOut.data = data
297 295
298 296 def setRadarFrequency(self, frequency=None):
299 297
300 298 if frequency != None:
301 299 self.dataOut.frequency = frequency
302 300
303 301 return 1
304 302
305 303 def interpolateHeights(self, topLim, botLim):
306 304 #69 al 72 para julia
307 305 #82-84 para meteoros
308 306 if len(numpy.shape(self.dataOut.data))==2:
309 307 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
310 308 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
311 309 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
312 310 self.dataOut.data[:,botLim:topLim+1] = sampInterp
313 311 else:
314 312 nHeights = self.dataOut.data.shape[2]
315 313 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
316 314 y = self.dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
317 315 f = interpolate.interp1d(x, y, axis = 2)
318 316 xnew = numpy.arange(botLim,topLim+1)
319 317 ynew = f(xnew)
320 318
321 319 self.dataOut.data[:,:,botLim:topLim+1] = ynew
322 320
323 321 # import collections
324 322 @MPDecorator
325 323 class CohInt(Operation):
326 324
327 325 isConfig = False
328 326 __profIndex = 0
329 327 __byTime = False
330 328 __initime = None
331 329 __lastdatatime = None
332 330 __integrationtime = None
333 331 __buffer = None
334 332 __bufferStride = []
335 333 __dataReady = False
336 334 __profIndexStride = 0
337 335 __dataToPutStride = False
338 336 n = None
339 337
340 338 def __init__(self):#, **kwargs):
341 339
342 340 Operation.__init__(self)#, **kwargs)
343 341
344 342 # self.isConfig = False
345 343
346 344 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
347 345 """
348 346 Set the parameters of the integration class.
349 347
350 348 Inputs:
351 349
352 350 n : Number of coherent integrations
353 351 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
354 352 overlapping :
355 353 """
356 354
357 355 self.__initime = None
358 356 self.__lastdatatime = 0
359 357 self.__buffer = None
360 358 self.__dataReady = False
361 359 self.byblock = byblock
362 360 self.stride = stride
363 361
364 362 if n == None and timeInterval == None:
365 363 raise ValueError("n or timeInterval should be specified ...")
366 364
367 365 if n != None:
368 366 self.n = n
369 367 self.__byTime = False
370 368 else:
371 369 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
372 370 self.n = 9999
373 371 self.__byTime = True
374 372
375 373 if overlapping:
376 374 self.__withOverlapping = True
377 375 self.__buffer = None
378 376 else:
379 377 self.__withOverlapping = False
380 378 self.__buffer = 0
381 379
382 380 self.__profIndex = 0
383 381
384 382 def putData(self, data):
385 383
386 384 """
387 385 Add a profile to the __buffer and increase in one the __profileIndex
388 386
389 387 """
390 388
391 389 if not self.__withOverlapping:
392 390 self.__buffer += data.copy()
393 391 self.__profIndex += 1
394 392 return
395 393
396 394 #Overlapping data
397 395 nChannels, nHeis = data.shape
398 396 data = numpy.reshape(data, (1, nChannels, nHeis))
399 397
400 398 #If the buffer is empty then it takes the data value
401 399 if self.__buffer is None:
402 400 self.__buffer = data
403 401 self.__profIndex += 1
404 402 return
405 403
406 404 #If the buffer length is lower than n then stakcing the data value
407 405 if self.__profIndex < self.n:
408 406 self.__buffer = numpy.vstack((self.__buffer, data))
409 407 self.__profIndex += 1
410 408 return
411 409
412 410 #If the buffer length is equal to n then replacing the last buffer value with the data value
413 411 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
414 412 self.__buffer[self.n-1] = data
415 413 self.__profIndex = self.n
416 414 return
417 415
418 416
419 417 def pushData(self):
420 418 """
421 419 Return the sum of the last profiles and the profiles used in the sum.
422 420
423 421 Affected:
424 422
425 423 self.__profileIndex
426 424
427 425 """
428 426
429 427 if not self.__withOverlapping:
430 428 data = self.__buffer
431 429 n = self.__profIndex
432 430
433 431 self.__buffer = 0
434 432 self.__profIndex = 0
435 433
436 434 return data, n
437 435
438 436 #Integration with Overlapping
439 437 data = numpy.sum(self.__buffer, axis=0)
440 438 # print data
441 439 # raise
442 440 n = self.__profIndex
443 441
444 442 return data, n
445 443
446 444 def byProfiles(self, data):
447 445
448 446 self.__dataReady = False
449 447 avgdata = None
450 448 # n = None
451 449 # print data
452 450 # raise
453 451 self.putData(data)
454 452
455 453 if self.__profIndex == self.n:
456 454 avgdata, n = self.pushData()
457 455 self.__dataReady = True
458 456
459 457 return avgdata
460 458
461 459 def byTime(self, data, datatime):
462 460
463 461 self.__dataReady = False
464 462 avgdata = None
465 463 n = None
466 464
467 465 self.putData(data)
468 466
469 467 if (datatime - self.__initime) >= self.__integrationtime:
470 468 avgdata, n = self.pushData()
471 469 self.n = n
472 470 self.__dataReady = True
473 471
474 472 return avgdata
475 473
476 474 def integrateByStride(self, data, datatime):
477 475 # print data
478 476 if self.__profIndex == 0:
479 477 self.__buffer = [[data.copy(), datatime]]
480 478 else:
481 479 self.__buffer.append([data.copy(),datatime])
482 480 self.__profIndex += 1
483 481 self.__dataReady = False
484 482
485 483 if self.__profIndex == self.n * self.stride :
486 484 self.__dataToPutStride = True
487 485 self.__profIndexStride = 0
488 486 self.__profIndex = 0
489 487 self.__bufferStride = []
490 488 for i in range(self.stride):
491 489 current = self.__buffer[i::self.stride]
492 490 data = numpy.sum([t[0] for t in current], axis=0)
493 491 avgdatatime = numpy.average([t[1] for t in current])
494 492 # print data
495 493 self.__bufferStride.append((data, avgdatatime))
496 494
497 495 if self.__dataToPutStride:
498 496 self.__dataReady = True
499 497 self.__profIndexStride += 1
500 498 if self.__profIndexStride == self.stride:
501 499 self.__dataToPutStride = False
502 500 # print self.__bufferStride[self.__profIndexStride - 1]
503 501 # raise
504 502 return self.__bufferStride[self.__profIndexStride - 1]
505 503
506 504
507 505 return None, None
508 506
509 507 def integrate(self, data, datatime=None):
510 508
511 509 if self.__initime == None:
512 510 self.__initime = datatime
513 511
514 512 if self.__byTime:
515 513 avgdata = self.byTime(data, datatime)
516 514 else:
517 515 avgdata = self.byProfiles(data)
518 516
519 517
520 518 self.__lastdatatime = datatime
521 519
522 520 if avgdata is None:
523 521 return None, None
524 522
525 523 avgdatatime = self.__initime
526 524
527 525 deltatime = datatime - self.__lastdatatime
528 526
529 527 if not self.__withOverlapping:
530 528 self.__initime = datatime
531 529 else:
532 530 self.__initime += deltatime
533 531
534 532 return avgdata, avgdatatime
535 533
536 534 def integrateByBlock(self, dataOut):
537 535
538 536 times = int(dataOut.data.shape[1]/self.n)
539 537 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
540 538
541 539 id_min = 0
542 540 id_max = self.n
543 541
544 542 for i in range(times):
545 543 junk = dataOut.data[:,id_min:id_max,:]
546 544 avgdata[:,i,:] = junk.sum(axis=1)
547 545 id_min += self.n
548 546 id_max += self.n
549 547
550 548 timeInterval = dataOut.ippSeconds*self.n
551 549 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
552 550 self.__dataReady = True
553 551 return avgdata, avgdatatime
554 552
555 553 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
556 554
557 555 if not self.isConfig:
558 556 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
559 557 self.isConfig = True
560 558
561 559 if dataOut.flagDataAsBlock:
562 560 """
563 561 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
564 562 """
565 563 avgdata, avgdatatime = self.integrateByBlock(dataOut)
566 564 dataOut.nProfiles /= self.n
567 565 else:
568 566 if stride is None:
569 567 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
570 568 else:
571 569 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
572 570
573 571
574 572 # dataOut.timeInterval *= n
575 573 dataOut.flagNoData = True
576 574
577 575 if self.__dataReady:
578 576 dataOut.data = avgdata
579 577 dataOut.nCohInt *= self.n
580 578 dataOut.utctime = avgdatatime
581 579 # print avgdata, avgdatatime
582 580 # raise
583 581 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
584 582 dataOut.flagNoData = False
585 583 return dataOut
586 584 @MPDecorator
587 585 class Decoder(Operation):
588 586
589 587 isConfig = False
590 588 __profIndex = 0
591 589
592 590 code = None
593 591
594 592 nCode = None
595 593 nBaud = None
596 594
597 595 def __init__(self):#, **kwargs):
598 596
599 597 Operation.__init__(self)#, **kwargs)
600 598
601 599 self.times = None
602 600 self.osamp = None
603 601 # self.__setValues = False
604 602 # self.isConfig = False
605 603 self.setupReq = False
606 604 def setup(self, code, osamp, dataOut):
607 605
608 606 self.__profIndex = 0
609 607
610 608 self.code = code
611 609
612 610 self.nCode = len(code)
613 611 self.nBaud = len(code[0])
614 612
615 613 if (osamp != None) and (osamp >1):
616 614 self.osamp = osamp
617 615 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
618 616 self.nBaud = self.nBaud*self.osamp
619 617
620 618 self.__nChannels = dataOut.nChannels
621 619 self.__nProfiles = dataOut.nProfiles
622 620 self.__nHeis = dataOut.nHeights
623 621
624 622 if self.__nHeis < self.nBaud:
625 623 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
626 624
627 625 #Frequency
628 626 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
629 627
630 628 __codeBuffer[:,0:self.nBaud] = self.code
631 629
632 630 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
633 631
634 632 if dataOut.flagDataAsBlock:
635 633
636 634 self.ndatadec = self.__nHeis #- self.nBaud + 1
637 635
638 636 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
639 637
640 638 else:
641 639
642 640 #Time
643 641 self.ndatadec = self.__nHeis #- self.nBaud + 1
644 642
645 643 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
646 644
647 645 def __convolutionInFreq(self, data):
648 646
649 647 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
650 648
651 649 fft_data = numpy.fft.fft(data, axis=1)
652 650
653 651 conv = fft_data*fft_code
654 652
655 653 data = numpy.fft.ifft(conv,axis=1)
656 654
657 655 return data
658 656
659 657 def __convolutionInFreqOpt(self, data):
660 658
661 659 raise NotImplementedError
662 660
663 661 def __convolutionInTime(self, data):
664 662
665 663 code = self.code[self.__profIndex]
666 664 for i in range(self.__nChannels):
667 665 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
668 666
669 667 return self.datadecTime
670 668
671 669 def __convolutionByBlockInTime(self, data):
672 670
673 671 repetitions = self.__nProfiles / self.nCode
674 672
675 673 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
676 674 junk = junk.flatten()
677 675 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
678 676 profilesList = range(self.__nProfiles)
679 677
680 678 for i in range(self.__nChannels):
681 679 for j in profilesList:
682 680 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
683 681 return self.datadecTime
684 682
685 683 def __convolutionByBlockInFreq(self, data):
686 684
687 685 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
688 686
689 687
690 688 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
691 689
692 690 fft_data = numpy.fft.fft(data, axis=2)
693 691
694 692 conv = fft_data*fft_code
695 693
696 694 data = numpy.fft.ifft(conv,axis=2)
697 695
698 696 return data
699 697
700 698
701 699 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
702 700
703 701 if dataOut.flagDecodeData:
704 702 print("This data is already decoded, recoding again ...")
705 703
706 704 if not self.isConfig:
707 705
708 706 if code is None:
709 707 if dataOut.code is None:
710 708 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
711 709
712 710 code = dataOut.code
713 711 else:
714 712 code = numpy.array(code).reshape(nCode,nBaud)
715 713 self.setup(code, osamp, dataOut)
716 714
717 715 self.isConfig = True
718 716
719 717 if mode == 3:
720 718 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
721 719
722 720 if times != None:
723 721 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
724 722
725 723 if self.code is None:
726 724 print("Fail decoding: Code is not defined.")
727 725 return
728 726
729 727 self.__nProfiles = dataOut.nProfiles
730 728 datadec = None
731 729
732 730 if mode == 3:
733 731 mode = 0
734 732
735 733 if dataOut.flagDataAsBlock:
736 734 """
737 735 Decoding when data have been read as block,
738 736 """
739 737
740 738 if mode == 0:
741 739 datadec = self.__convolutionByBlockInTime(dataOut.data)
742 740 if mode == 1:
743 741 datadec = self.__convolutionByBlockInFreq(dataOut.data)
744 742 else:
745 743 """
746 744 Decoding when data have been read profile by profile
747 745 """
748 746 if mode == 0:
749 747 datadec = self.__convolutionInTime(dataOut.data)
750 748
751 749 if mode == 1:
752 750 datadec = self.__convolutionInFreq(dataOut.data)
753 751
754 752 if mode == 2:
755 753 datadec = self.__convolutionInFreqOpt(dataOut.data)
756 754
757 755 if datadec is None:
758 756 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
759 757
760 758 dataOut.code = self.code
761 759 dataOut.nCode = self.nCode
762 760 dataOut.nBaud = self.nBaud
763 761
764 762 dataOut.data = datadec
765 763
766 764 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
767 765
768 766 dataOut.flagDecodeData = True #asumo q la data esta decodificada
769 767
770 768 if self.__profIndex == self.nCode-1:
771 769 self.__profIndex = 0
772 770 return dataOut
773 771
774 772 self.__profIndex += 1
775 773
776 774 return dataOut
777 775 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
778 776
779 777 @MPDecorator
780 778 class ProfileConcat(Operation):
781 779
782 780 isConfig = False
783 781 buffer = None
784 782
785 783 def __init__(self):#, **kwargs):
786 784
787 785 Operation.__init__(self)#, **kwargs)
788 786 self.profileIndex = 0
789 787
790 788 def reset(self):
791 789 self.buffer = numpy.zeros_like(self.buffer)
792 790 self.start_index = 0
793 791 self.times = 1
794 792
795 793 def setup(self, data, m, n=1):
796 794 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
797 795 self.nHeights = data.shape[1]#.nHeights
798 796 self.start_index = 0
799 797 self.times = 1
800 798
801 799 def concat(self, data):
802 800
803 801 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
804 802 self.start_index = self.start_index + self.nHeights
805 803
806 804 def run(self, dataOut, m):
807 805
808 806 dataOut.flagNoData = True
809 807
810 808 if not self.isConfig:
811 809 self.setup(dataOut.data, m, 1)
812 810 self.isConfig = True
813 811
814 812 if dataOut.flagDataAsBlock:
815 813 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
816 814
817 815 else:
818 816 self.concat(dataOut.data)
819 817 self.times += 1
820 818 if self.times > m:
821 819 dataOut.data = self.buffer
822 820 self.reset()
823 821 dataOut.flagNoData = False
824 822 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
825 823 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
826 824 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
827 825 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
828 826 dataOut.ippSeconds *= m
829 827 return dataOut
830 828 @MPDecorator
831 829 class ProfileSelector(Operation):
832 830
833 831 profileIndex = None
834 832 # Tamanho total de los perfiles
835 833 nProfiles = None
836 834
837 835 def __init__(self):#, **kwargs):
838 836
839 837 Operation.__init__(self)#, **kwargs)
840 838 self.profileIndex = 0
841 839
842 840 def incProfileIndex(self):
843 841
844 842 self.profileIndex += 1
845 843
846 844 if self.profileIndex >= self.nProfiles:
847 845 self.profileIndex = 0
848 846
849 847 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
850 848
851 849 if profileIndex < minIndex:
852 850 return False
853 851
854 852 if profileIndex > maxIndex:
855 853 return False
856 854
857 855 return True
858 856
859 857 def isThisProfileInList(self, profileIndex, profileList):
860 858
861 859 if profileIndex not in profileList:
862 860 return False
863 861
864 862 return True
865 863
866 864 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
867 865
868 866 """
869 867 ProfileSelector:
870 868
871 869 Inputs:
872 870 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
873 871
874 872 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
875 873
876 874 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
877 875
878 876 """
879 877
880 878 if rangeList is not None:
881 879 if type(rangeList[0]) not in (tuple, list):
882 880 rangeList = [rangeList]
883 881
884 882 dataOut.flagNoData = True
885 883
886 884 if dataOut.flagDataAsBlock:
887 885 """
888 886 data dimension = [nChannels, nProfiles, nHeis]
889 887 """
890 888 if profileList != None:
891 889 dataOut.data = dataOut.data[:,profileList,:]
892 890
893 891 if profileRangeList != None:
894 892 minIndex = profileRangeList[0]
895 893 maxIndex = profileRangeList[1]
896 894 profileList = list(range(minIndex, maxIndex+1))
897 895
898 896 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
899 897
900 898 if rangeList != None:
901 899
902 900 profileList = []
903 901
904 902 for thisRange in rangeList:
905 903 minIndex = thisRange[0]
906 904 maxIndex = thisRange[1]
907 905
908 906 profileList.extend(list(range(minIndex, maxIndex+1)))
909 907
910 908 dataOut.data = dataOut.data[:,profileList,:]
911 909
912 910 dataOut.nProfiles = len(profileList)
913 911 dataOut.profileIndex = dataOut.nProfiles - 1
914 912 dataOut.flagNoData = False
915 913
916 914 return True
917 915
918 916 """
919 917 data dimension = [nChannels, nHeis]
920 918 """
921 919
922 920 if profileList != None:
923 921
924 922 if self.isThisProfileInList(dataOut.profileIndex, profileList):
925 923
926 924 self.nProfiles = len(profileList)
927 925 dataOut.nProfiles = self.nProfiles
928 926 dataOut.profileIndex = self.profileIndex
929 927 dataOut.flagNoData = False
930 928
931 929 self.incProfileIndex()
932 930 return True
933 931
934 932 if profileRangeList != None:
935 933
936 934 minIndex = profileRangeList[0]
937 935 maxIndex = profileRangeList[1]
938 936
939 937 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
940 938
941 939 self.nProfiles = maxIndex - minIndex + 1
942 940 dataOut.nProfiles = self.nProfiles
943 941 dataOut.profileIndex = self.profileIndex
944 942 dataOut.flagNoData = False
945 943
946 944 self.incProfileIndex()
947 945 return True
948 946
949 947 if rangeList != None:
950 948
951 949 nProfiles = 0
952 950
953 951 for thisRange in rangeList:
954 952 minIndex = thisRange[0]
955 953 maxIndex = thisRange[1]
956 954
957 955 nProfiles += maxIndex - minIndex + 1
958 956
959 957 for thisRange in rangeList:
960 958
961 959 minIndex = thisRange[0]
962 960 maxIndex = thisRange[1]
963 961
964 962 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
965 963
966 964 self.nProfiles = nProfiles
967 965 dataOut.nProfiles = self.nProfiles
968 966 dataOut.profileIndex = self.profileIndex
969 967 dataOut.flagNoData = False
970 968
971 969 self.incProfileIndex()
972 970
973 971 break
974 972
975 973 return True
976 974
977 975
978 976 if beam != None: #beam is only for AMISR data
979 977 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
980 978 dataOut.flagNoData = False
981 979 dataOut.profileIndex = self.profileIndex
982 980
983 981 self.incProfileIndex()
984 982
985 983 return True
986 984
987 985 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
988 986
989 987 #return False
990 988 return dataOut
991 989 @MPDecorator
992 990 class Reshaper(Operation):
993 991
994 992 def __init__(self):#, **kwargs):
995 993
996 994 Operation.__init__(self)#, **kwargs)
997 995
998 996 self.__buffer = None
999 997 self.__nitems = 0
1000 998
1001 999 def __appendProfile(self, dataOut, nTxs):
1002 1000
1003 1001 if self.__buffer is None:
1004 1002 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1005 1003 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1006 1004
1007 1005 ini = dataOut.nHeights * self.__nitems
1008 1006 end = ini + dataOut.nHeights
1009 1007
1010 1008 self.__buffer[:, ini:end] = dataOut.data
1011 1009
1012 1010 self.__nitems += 1
1013 1011
1014 1012 return int(self.__nitems*nTxs)
1015 1013
1016 1014 def __getBuffer(self):
1017 1015
1018 1016 if self.__nitems == int(1./self.__nTxs):
1019 1017
1020 1018 self.__nitems = 0
1021 1019
1022 1020 return self.__buffer.copy()
1023 1021
1024 1022 return None
1025 1023
1026 1024 def __checkInputs(self, dataOut, shape, nTxs):
1027 1025
1028 1026 if shape is None and nTxs is None:
1029 1027 raise ValueError("Reshaper: shape of factor should be defined")
1030 1028
1031 1029 if nTxs:
1032 1030 if nTxs < 0:
1033 1031 raise ValueError("nTxs should be greater than 0")
1034 1032
1035 1033 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1036 1034 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1037 1035
1038 1036 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1039 1037
1040 1038 return shape, nTxs
1041 1039
1042 1040 if len(shape) != 2 and len(shape) != 3:
1043 1041 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))
1044 1042
1045 1043 if len(shape) == 2:
1046 1044 shape_tuple = [dataOut.nChannels]
1047 1045 shape_tuple.extend(shape)
1048 1046 else:
1049 1047 shape_tuple = list(shape)
1050 1048
1051 1049 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1052 1050
1053 1051 return shape_tuple, nTxs
1054 1052
1055 1053 def run(self, dataOut, shape=None, nTxs=None):
1056 1054
1057 1055 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1058 1056
1059 1057 dataOut.flagNoData = True
1060 1058 profileIndex = None
1061 1059
1062 1060 if dataOut.flagDataAsBlock:
1063 1061
1064 1062 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1065 1063 dataOut.flagNoData = False
1066 1064
1067 1065 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1068 1066
1069 1067 else:
1070 1068
1071 1069 if self.__nTxs < 1:
1072 1070
1073 1071 self.__appendProfile(dataOut, self.__nTxs)
1074 1072 new_data = self.__getBuffer()
1075 1073
1076 1074 if new_data is not None:
1077 1075 dataOut.data = new_data
1078 1076 dataOut.flagNoData = False
1079 1077
1080 1078 profileIndex = dataOut.profileIndex*nTxs
1081 1079
1082 1080 else:
1083 1081 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1084 1082
1085 1083 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1086 1084
1087 1085 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1088 1086
1089 1087 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1090 1088
1091 1089 dataOut.profileIndex = profileIndex
1092 1090
1093 1091 dataOut.ippSeconds /= self.__nTxs
1094 1092
1095 1093 return dataOut
1096 1094 @MPDecorator
1097 1095 class SplitProfiles(Operation):
1098 1096
1099 1097 def __init__(self):#, **kwargs):
1100 1098
1101 1099 Operation.__init__(self)#, **kwargs)
1102 1100
1103 1101 def run(self, dataOut, n):
1104 1102
1105 1103 dataOut.flagNoData = True
1106 1104 profileIndex = None
1107 1105
1108 1106 if dataOut.flagDataAsBlock:
1109 1107
1110 1108 #nchannels, nprofiles, nsamples
1111 1109 shape = dataOut.data.shape
1112 1110
1113 1111 if shape[2] % n != 0:
1114 1112 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1115 1113
1116 1114 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1117 1115
1118 1116 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1119 1117 dataOut.flagNoData = False
1120 1118
1121 1119 profileIndex = int(dataOut.nProfiles/n) - 1
1122 1120
1123 1121 else:
1124 1122
1125 1123 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1126 1124
1127 1125 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1128 1126
1129 1127 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1130 1128
1131 1129 dataOut.nProfiles = int(dataOut.nProfiles*n)
1132 1130
1133 1131 dataOut.profileIndex = profileIndex
1134 1132
1135 1133 dataOut.ippSeconds /= n
1136 1134
1137 1135 return dataOut
1138 1136 @MPDecorator
1139 1137 class CombineProfiles(Operation):
1140 1138 def __init__(self):#, **kwargs):
1141 1139
1142 1140 Operation.__init__(self)#, **kwargs)
1143 1141
1144 1142 self.__remData = None
1145 1143 self.__profileIndex = 0
1146 1144
1147 1145 def run(self, dataOut, n):
1148 1146
1149 1147 dataOut.flagNoData = True
1150 1148 profileIndex = None
1151 1149
1152 1150 if dataOut.flagDataAsBlock:
1153 1151
1154 1152 #nchannels, nprofiles, nsamples
1155 1153 shape = dataOut.data.shape
1156 1154 new_shape = shape[0], shape[1]/n, shape[2]*n
1157 1155
1158 1156 if shape[1] % n != 0:
1159 1157 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1160 1158
1161 1159 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1162 1160 dataOut.flagNoData = False
1163 1161
1164 1162 profileIndex = int(dataOut.nProfiles*n) - 1
1165 1163
1166 1164 else:
1167 1165
1168 1166 #nchannels, nsamples
1169 1167 if self.__remData is None:
1170 1168 newData = dataOut.data
1171 1169 else:
1172 1170 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1173 1171
1174 1172 self.__profileIndex += 1
1175 1173
1176 1174 if self.__profileIndex < n:
1177 1175 self.__remData = newData
1178 1176 #continue
1179 1177 return
1180 1178
1181 1179 self.__profileIndex = 0
1182 1180 self.__remData = None
1183 1181
1184 1182 dataOut.data = newData
1185 1183 dataOut.flagNoData = False
1186 1184
1187 1185 profileIndex = dataOut.profileIndex/n
1188 1186
1189 1187
1190 1188 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1191 1189
1192 1190 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1193 1191
1194 1192 dataOut.nProfiles = int(dataOut.nProfiles/n)
1195 1193
1196 1194 dataOut.profileIndex = profileIndex
1197 1195
1198 1196 dataOut.ippSeconds *= n
1199 1197
1200 1198 return dataOut
1201 1199 # import collections
1202 1200 # from scipy.stats import mode
1203 1201 #
1204 1202 # class Synchronize(Operation):
1205 1203 #
1206 1204 # isConfig = False
1207 1205 # __profIndex = 0
1208 1206 #
1209 1207 # def __init__(self, **kwargs):
1210 1208 #
1211 1209 # Operation.__init__(self, **kwargs)
1212 1210 # # self.isConfig = False
1213 1211 # self.__powBuffer = None
1214 1212 # self.__startIndex = 0
1215 1213 # self.__pulseFound = False
1216 1214 #
1217 1215 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1218 1216 #
1219 1217 # #Read data
1220 1218 #
1221 1219 # powerdB = dataOut.getPower(channel = channel)
1222 1220 # noisedB = dataOut.getNoise(channel = channel)[0]
1223 1221 #
1224 1222 # self.__powBuffer.extend(powerdB.flatten())
1225 1223 #
1226 1224 # dataArray = numpy.array(self.__powBuffer)
1227 1225 #
1228 1226 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1229 1227 #
1230 1228 # maxValue = numpy.nanmax(filteredPower)
1231 1229 #
1232 1230 # if maxValue < noisedB + 10:
1233 1231 # #No se encuentra ningun pulso de transmision
1234 1232 # return None
1235 1233 #
1236 1234 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1237 1235 #
1238 1236 # if len(maxValuesIndex) < 2:
1239 1237 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1240 1238 # return None
1241 1239 #
1242 1240 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1243 1241 #
1244 1242 # #Seleccionar solo valores con un espaciamiento de nSamples
1245 1243 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1246 1244 #
1247 1245 # if len(pulseIndex) < 2:
1248 1246 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1249 1247 # return None
1250 1248 #
1251 1249 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1252 1250 #
1253 1251 # #remover senales que se distancien menos de 10 unidades o muestras
1254 1252 # #(No deberian existir IPP menor a 10 unidades)
1255 1253 #
1256 1254 # realIndex = numpy.where(spacing > 10 )[0]
1257 1255 #
1258 1256 # if len(realIndex) < 2:
1259 1257 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1260 1258 # return None
1261 1259 #
1262 1260 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1263 1261 # realPulseIndex = pulseIndex[realIndex]
1264 1262 #
1265 1263 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1266 1264 #
1267 1265 # print "IPP = %d samples" %period
1268 1266 #
1269 1267 # self.__newNSamples = dataOut.nHeights #int(period)
1270 1268 # self.__startIndex = int(realPulseIndex[0])
1271 1269 #
1272 1270 # return 1
1273 1271 #
1274 1272 #
1275 1273 # def setup(self, nSamples, nChannels, buffer_size = 4):
1276 1274 #
1277 1275 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1278 1276 # maxlen = buffer_size*nSamples)
1279 1277 #
1280 1278 # bufferList = []
1281 1279 #
1282 1280 # for i in range(nChannels):
1283 1281 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1284 1282 # maxlen = buffer_size*nSamples)
1285 1283 #
1286 1284 # bufferList.append(bufferByChannel)
1287 1285 #
1288 1286 # self.__nSamples = nSamples
1289 1287 # self.__nChannels = nChannels
1290 1288 # self.__bufferList = bufferList
1291 1289 #
1292 1290 # def run(self, dataOut, channel = 0):
1293 1291 #
1294 1292 # if not self.isConfig:
1295 1293 # nSamples = dataOut.nHeights
1296 1294 # nChannels = dataOut.nChannels
1297 1295 # self.setup(nSamples, nChannels)
1298 1296 # self.isConfig = True
1299 1297 #
1300 1298 # #Append new data to internal buffer
1301 1299 # for thisChannel in range(self.__nChannels):
1302 1300 # bufferByChannel = self.__bufferList[thisChannel]
1303 1301 # bufferByChannel.extend(dataOut.data[thisChannel])
1304 1302 #
1305 1303 # if self.__pulseFound:
1306 1304 # self.__startIndex -= self.__nSamples
1307 1305 #
1308 1306 # #Finding Tx Pulse
1309 1307 # if not self.__pulseFound:
1310 1308 # indexFound = self.__findTxPulse(dataOut, channel)
1311 1309 #
1312 1310 # if indexFound == None:
1313 1311 # dataOut.flagNoData = True
1314 1312 # return
1315 1313 #
1316 1314 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1317 1315 # self.__pulseFound = True
1318 1316 # self.__startIndex = indexFound
1319 1317 #
1320 1318 # #If pulse was found ...
1321 1319 # for thisChannel in range(self.__nChannels):
1322 1320 # bufferByChannel = self.__bufferList[thisChannel]
1323 1321 # #print self.__startIndex
1324 1322 # x = numpy.array(bufferByChannel)
1325 1323 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1326 1324 #
1327 1325 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1328 1326 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1329 1327 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1330 1328 #
1331 1329 # dataOut.data = self.__arrayBuffer
1332 1330 #
1333 1331 # self.__startIndex += self.__newNSamples
1334 1332 #
1335 1333 # return
Show file before | Show file after | Hide comments
@@ -1,70 +1,64
1 1 '''
2 2 Created on Jul 16, 2014
3 3
4 4 @author: Miguel Urco
5 5 '''
6 6
7 7 import os
8 8 from setuptools import setup, Extension
9 9 from setuptools.command.build_ext import build_ext as _build_ext
10 10 from schainpy import __version__
11 11
12 12 class build_ext(_build_ext):
13 13 def finalize_options(self):
14 14 _build_ext.finalize_options(self)
15 15 # Prevent numpy from thinking it is still in its setup process:
16 16 __builtins__.__NUMPY_SETUP__ = False
17 17 import numpy
18 18 self.include_dirs.append(numpy.get_include())
19 19
20 20 setup(name = "schainpy",
21 21 version = __version__,
22 22 description = "Python tools to read, write and process Jicamarca data",
23 23 author = "Miguel Urco",
24 24 author_email = "miguel.urco@jro.igp.gob.pe",
25 25 url = "http://jro.igp.gob.pe",
26 26 packages = {'schainpy',
27 27 'schainpy.model',
28 28 'schainpy.model.data',
29 29 'schainpy.model.graphics',
30 30 'schainpy.model.io',
31 31 'schainpy.model.proc',
32 32 'schainpy.model.serializer',
33 33 'schainpy.model.utils',
34 34 'schainpy.utils',
35 35 'schainpy.gui',
36 36 'schainpy.gui.figures',
37 37 'schainpy.gui.viewcontroller',
38 38 'schainpy.gui.viewer',
39 39 'schainpy.gui.viewer.windows',
40 40 'schainpy.cli'},
41 41 ext_package = 'schainpy',
42 42 package_data = {'': ['schain.conf.template'],
43 43 'schainpy.gui.figures': ['*.png', '*.jpg'],
44 44 'schainpy.files': ['*.oga']
45 45 },
46 46 include_package_data = False,
47 47 scripts = ['schainpy/gui/schainGUI'],
48 ext_modules = [
49 Extension("cSchain", ["schainpy/model/proc/extensions.c"])
50 ],
51 48 entry_points = {
52 49 'console_scripts': [
53 50 'schain = schainpy.cli.cli:main',
54 51 ],
55 52 },
56 53 cmdclass = {'build_ext': build_ext},
57 54 setup_requires = ["numpy >= 1.11.2"],
58 55 install_requires = [
59 56 "scipy >= 0.14.0",
60 57 "h5py >= 2.2.1",
61 58 "matplotlib >= 2.0.0",
62 "pyfits >= 3.4",
63 "paramiko >= 2.1.2",
64 "paho-mqtt >= 1.2",
65 59 "zmq",
66 60 "fuzzywuzzy",
67 61 "click",
68 62 "python-Levenshtein"
69 63 ],
70 64 )
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
General Comments 0
You need to be logged in to leave comments. Login now