##// END OF EJS Templates
schain
RSS
update repository
avaldez -
r1259:c52310cc7d94
parent child
Show More
Side by Side Unified
Context file:
r1259:c52310cc7d94 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1261 +1,1264
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 12 from schainpy import cSchain
13 13
14 14
15 15 def getNumpyDtype(dataTypeCode):
16 16
17 17 if dataTypeCode == 0:
18 18 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
19 19 elif dataTypeCode == 1:
20 20 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
21 21 elif dataTypeCode == 2:
22 22 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
23 23 elif dataTypeCode == 3:
24 24 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
25 25 elif dataTypeCode == 4:
26 26 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
27 27 elif dataTypeCode == 5:
28 28 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
29 29 else:
30 30 raise ValueError, 'dataTypeCode was not defined'
31 31
32 32 return numpyDtype
33 33
34 34
35 35 def getDataTypeCode(numpyDtype):
36 36
37 37 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
38 38 datatype = 0
39 39 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
40 40 datatype = 1
41 41 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
42 42 datatype = 2
43 43 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
44 44 datatype = 3
45 45 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
46 46 datatype = 4
47 47 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
48 48 datatype = 5
49 49 else:
50 50 datatype = None
51 51
52 52 return datatype
53 53
54 54
55 55 def hildebrand_sekhon(data, navg):
56 56 """
57 57 This method is for the objective determination of the noise level in Doppler spectra. This
58 58 implementation technique is based on the fact that the standard deviation of the spectral
59 59 densities is equal to the mean spectral density for white Gaussian noise
60 60
61 61 Inputs:
62 62 Data : heights
63 63 navg : numbers of averages
64 64
65 65 Return:
66 66 -1 : any error
67 67 anoise : 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 # nums_min = 5
76 76 #
77 77 # sump = 0.
78 78 #
79 79 # sumq = 0.
80 80 #
81 81 # j = 0
82 82 #
83 83 # cont = 1
84 84 #
85 85 # while((cont==1)and(j<lenOfData)):
86 86 #
87 87 # sump += sortdata[j]
88 88 #
89 89 # sumq += sortdata[j]**2
90 90 #
91 91 # if j > nums_min:
92 92 # rtest = float(j)/(j-1) + 1.0/navg
93 93 # if ((sumq*j) > (rtest*sump**2)):
94 94 # j = j - 1
95 95 # sump = sump - sortdata[j]
96 96 # sumq = sumq - sortdata[j]**2
97 97 # cont = 0
98 98 #
99 99 # j += 1
100 100 #
101 101 # lnoise = sump /j
102 102 #
103 103 # return lnoise
104 104
105 105 return cSchain.hildebrand_sekhon(sortdata, navg)
106 106
107 107
108 108 class Beam:
109 109
110 110 def __init__(self):
111 111 self.codeList = []
112 112 self.azimuthList = []
113 113 self.zenithList = []
114 114
115 115
116 116 class GenericData(object):
117 117
118 118 flagNoData = True
119 119
120 120 def copy(self, inputObj=None):
121 121
122 122 if inputObj == None:
123 123 return copy.deepcopy(self)
124 124
125 125 for key in inputObj.__dict__.keys():
126 126
127 127 attribute = inputObj.__dict__[key]
128 128
129 129 # If this attribute is a tuple or list
130 130 if type(inputObj.__dict__[key]) in (tuple, list):
131 131 self.__dict__[key] = attribute[:]
132 132 continue
133 133
134 134 # If this attribute is another object or instance
135 135 if hasattr(attribute, '__dict__'):
136 136 self.__dict__[key] = attribute.copy()
137 137 continue
138 138
139 139 self.__dict__[key] = inputObj.__dict__[key]
140 140
141 141 def deepcopy(self):
142 142
143 143 return copy.deepcopy(self)
144 144
145 145 def isEmpty(self):
146 146
147 147 return self.flagNoData
148 148
149 149
150 150 class JROData(GenericData):
151 151
152 152 # m_BasicHeader = BasicHeader()
153 153 # m_ProcessingHeader = ProcessingHeader()
154 154
155 155 systemHeaderObj = SystemHeader()
156 156
157 157 radarControllerHeaderObj = RadarControllerHeader()
158 158
159 159 # data = None
160 160
161 161 type = None
162 162
163 163 datatype = None # dtype but in string
164 164
165 165 # dtype = None
166 166
167 167 # nChannels = None
168 168
169 169 # nHeights = None
170 170
171 171 nProfiles = None
172 172
173 173 heightList = None
174 174
175 175 channelList = None
176 176
177 177 flagDiscontinuousBlock = False
178 178
179 179 useLocalTime = False
180 180
181 181 utctime = None
182 182
183 183 timeZone = None
184 184
185 185 dstFlag = None
186 186
187 187 errorCount = None
188 188
189 189 blocksize = None
190 190
191 191 # nCode = None
192 192 #
193 193 # nBaud = None
194 194 #
195 195 # code = None
196 196
197 197 flagDecodeData = False # asumo q la data no esta decodificada
198 198
199 199 flagDeflipData = False # asumo q la data no esta sin flip
200 200
201 201 flagShiftFFT = False
202 202
203 203 # ippSeconds = None
204 204
205 205 # timeInterval = None
206 206
207 207 nCohInt = None
208 208
209 209 # noise = None
210 210
211 211 windowOfFilter = 1
212 212
213 213 # Speed of ligth
214 214 C = 3e8
215 215
216 216 frequency = 49.92e6
217 217
218 218 realtime = False
219 219
220 220 beacon_heiIndexList = None
221 221
222 222 last_block = None
223 223
224 224 blocknow = None
225 225
226 226 azimuth = None
227 227
228 228 zenith = None
229 229
230 230 beam = Beam()
231 231
232 232 profileIndex = None
233 233
234 234 def getNoise(self):
235 235
236 236 raise NotImplementedError
237 237
238 238 def getNChannels(self):
239 239
240 240 return len(self.channelList)
241 241
242 242 def getChannelIndexList(self):
243 243
244 244 return range(self.nChannels)
245 245
246 246 def getNHeights(self):
247 247
248 248 return len(self.heightList)
249 249
250 250 def getHeiRange(self, extrapoints=0):
251 251
252 252 heis = self.heightList
253 253 # deltah = self.heightList[1] - self.heightList[0]
254 254 #
255 255 # heis.append(self.heightList[-1])
256 256
257 257 return heis
258 258
259 259 def getDeltaH(self):
260 260
261 261 delta = self.heightList[1] - self.heightList[0]
262 262
263 263 return delta
264 264
265 265 def getltctime(self):
266 266
267 267 if self.useLocalTime:
268 268 return self.utctime - self.timeZone * 60
269 269
270 270 return self.utctime
271 271
272 272 def getDatatime(self):
273 273
274 274 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
275 275 return datatimeValue
276 276
277 277 def getTimeRange(self):
278 278
279 279 datatime = []
280 280
281 281 datatime.append(self.ltctime)
282 282 datatime.append(self.ltctime + self.timeInterval + 1)
283 283
284 284 datatime = numpy.array(datatime)
285 285
286 286 return datatime
287 287
288 288 def getFmaxTimeResponse(self):
289 289
290 290 period = (10**-6) * self.getDeltaH() / (0.15)
291 291
292 292 PRF = 1. / (period * self.nCohInt)
293 293
294 294 fmax = PRF
295 295
296 296 return fmax
297 297
298 298 def getFmax(self):
299 299 PRF = 1. / (self.ippSeconds * self.nCohInt)
300 300 #print "ippSeconds",self.ippSeconds
301 301 #print "nCohInt",self.nIncohInt
302 302 #print PRF
303 303 #import time
304 304 #time.sleep(30)
305 305 fmax = PRF
306 306 return fmax
307 307
308 308 def getVmax(self):
309 309
310 310 _lambda = self.C / self.frequency
311 311
312 312 vmax = self.getFmax() * _lambda / 2
313 313
314 314 return vmax
315 315
316 316 def get_ippSeconds(self):
317 317 '''
318 318 '''
319 319 return self.radarControllerHeaderObj.ippSeconds
320 320
321 321 def set_ippSeconds(self, ippSeconds):
322 322 '''
323 323 '''
324 324
325 325 self.radarControllerHeaderObj.ippSeconds = ippSeconds
326 326
327 327 return
328 328
329 329 def get_dtype(self):
330 330 '''
331 331 '''
332 332 return getNumpyDtype(self.datatype)
333 333
334 334 def set_dtype(self, numpyDtype):
335 335 '''
336 336 '''
337 337
338 338 self.datatype = getDataTypeCode(numpyDtype)
339 339
340 340 def get_code(self):
341 341 '''
342 342 '''
343 343 return self.radarControllerHeaderObj.code
344 344
345 345 def set_code(self, code):
346 346 '''
347 347 '''
348 348 self.radarControllerHeaderObj.code = code
349 349
350 350 return
351 351
352 352 def get_ncode(self):
353 353 '''
354 354 '''
355 355 return self.radarControllerHeaderObj.nCode
356 356
357 357 def set_ncode(self, nCode):
358 358 '''
359 359 '''
360 360 self.radarControllerHeaderObj.nCode = nCode
361 361
362 362 return
363 363
364 364 def get_nbaud(self):
365 365 '''
366 366 '''
367 367 return self.radarControllerHeaderObj.nBaud
368 368
369 369 def set_nbaud(self, nBaud):
370 370 '''
371 371 '''
372 372 self.radarControllerHeaderObj.nBaud = nBaud
373 373
374 374 return
375 375
376 376 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
377 377 channelIndexList = property(
378 378 getChannelIndexList, "I'm the 'channelIndexList' property.")
379 379 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
380 380 #noise = property(getNoise, "I'm the 'nHeights' property.")
381 381 datatime = property(getDatatime, "I'm the 'datatime' property")
382 382 ltctime = property(getltctime, "I'm the 'ltctime' property")
383 383 ippSeconds = property(get_ippSeconds, set_ippSeconds)
384 384 dtype = property(get_dtype, set_dtype)
385 385 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
386 386 code = property(get_code, set_code)
387 387 nCode = property(get_ncode, set_ncode)
388 388 nBaud = property(get_nbaud, set_nbaud)
389 389
390 390
391 391 class Voltage(JROData):
392 392
393 393 # data es un numpy array de 2 dmensiones (canales, alturas)
394 394 data = None
395 395
396 396 def __init__(self):
397 397 '''
398 398 Constructor
399 399 '''
400 400
401 401 self.useLocalTime = True
402 402
403 403 self.radarControllerHeaderObj = RadarControllerHeader()
404 404
405 405 self.systemHeaderObj = SystemHeader()
406 406
407 407 self.type = "Voltage"
408 408
409 409 self.data = None
410 410
411 411 # self.dtype = None
412 412
413 413 # self.nChannels = 0
414 414
415 415 # self.nHeights = 0
416 416
417 417 self.nProfiles = None
418 418
419 419 self.heightList = None
420 420
421 421 self.channelList = None
422 422
423 423 # self.channelIndexList = None
424 424
425 425 self.flagNoData = True
426 426
427 427 self.flagDiscontinuousBlock = False
428 428
429 429 self.utctime = None
430 430
431 431 self.timeZone = None
432 432
433 433 self.dstFlag = None
434 434
435 435 self.errorCount = None
436 436
437 437 self.nCohInt = None
438 438
439 439 self.blocksize = None
440 440
441 441 self.flagDecodeData = False # asumo q la data no esta decodificada
442 442
443 443 self.flagDeflipData = False # asumo q la data no esta sin flip
444 444
445 445 self.flagShiftFFT = False
446 446
447 447 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
448 448
449 449 self.profileIndex = 0
450 450
451 451 def getNoisebyHildebrand(self, channel=None):
452 452 """
453 453 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
454 454
455 455 Return:
456 456 noiselevel
457 457 """
458 458
459 459 if channel != None:
460 460 data = self.data[channel]
461 461 nChannels = 1
462 462 else:
463 463 data = self.data
464 464 nChannels = self.nChannels
465 465
466 466 noise = numpy.zeros(nChannels)
467 467 power = data * numpy.conjugate(data)
468 468
469 469 for thisChannel in range(nChannels):
470 470 if nChannels == 1:
471 471 daux = power[:].real
472 472 else:
473 473 daux = power[thisChannel, :].real
474 474 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
475 475
476 476 return noise
477 477
478 478 def getNoise(self, type=1, channel=None):
479 479
480 480 if type == 1:
481 481 noise = self.getNoisebyHildebrand(channel)
482 482
483 483 return noise
484 484
485 485 def getPower(self, channel=None):
486 486
487 487 if channel != None:
488 488 data = self.data[channel]
489 489 else:
490 490 data = self.data
491 491
492 492 power = data * numpy.conjugate(data)
493 493 powerdB = 10 * numpy.log10(power.real)
494 494 powerdB = numpy.squeeze(powerdB)
495 495
496 496 return powerdB
497 497
498 498 def getTimeInterval(self):
499 499
500 500 timeInterval = self.ippSeconds * self.nCohInt
501 501
502 502 return timeInterval
503 503
504 504 noise = property(getNoise, "I'm the 'nHeights' property.")
505 505 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
506 506
507 507
508 508 class Spectra(JROData):
509 509
510 510 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
511 511 data_spc = None
512 512
513 513 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
514 514 data_cspc = None
515 515
516 516 # data dc es un numpy array de 2 dmensiones (canales, alturas)
517 517 data_dc = None
518 518
519 519 # data power
520 520 data_pwr = None
521 521
522 522 nFFTPoints = None
523 523
524 524 # nPairs = None
525 525
526 526 pairsList = None
527 527
528 528 nIncohInt = None
529 529
530 530 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
531 531
532 532 nCohInt = None # se requiere para determinar el valor de timeInterval
533 533
534 534 ippFactor = None
535 535
536 536 profileIndex = 0
537 537
538 538 plotting = "spectra"
539 539
540 540 def __init__(self):
541 541 '''
542 542 Constructor
543 543 '''
544 544
545 545 self.useLocalTime = True
546 546
547 547 self.radarControllerHeaderObj = RadarControllerHeader()
548 548
549 549 self.systemHeaderObj = SystemHeader()
550 550
551 551 self.type = "Spectra"
552 552
553 553 # self.data = None
554 554
555 555 # self.dtype = None
556 556
557 557 # self.nChannels = 0
558 558
559 559 # self.nHeights = 0
560 560
561 561 self.nProfiles = None
562 562
563 563 self.heightList = None
564 564
565 565 self.channelList = None
566 566
567 567 # self.channelIndexList = None
568 568
569 569 self.pairsList = None
570 570
571 571 self.flagNoData = True
572 572
573 573 self.flagDiscontinuousBlock = False
574 574
575 575 self.utctime = None
576 576
577 577 self.nCohInt = None
578 578
579 579 self.nIncohInt = None
580 580
581 581 self.blocksize = None
582 582
583 583 self.nFFTPoints = None
584 584
585 585 self.wavelength = None
586 586
587 587 self.flagDecodeData = False # asumo q la data no esta decodificada
588 588
589 589 self.flagDeflipData = False # asumo q la data no esta sin flip
590 590
591 591 self.flagShiftFFT = False
592 592
593 593 self.ippFactor = 1
594 594
595 595 #self.noise = None
596 596
597 597 self.beacon_heiIndexList = []
598 598
599 599 self.noise_estimation = None
600 600
601 601 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
602 602 """
603 603 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
604 604
605 605 Return:
606 606 noiselevel
607 607 """
608 608
609 609 noise = numpy.zeros(self.nChannels)
610 610
611 611 for channel in range(self.nChannels):
612 612 daux = self.data_spc[channel,
613 613 xmin_index:xmax_index, ymin_index:ymax_index]
614 614
615 615 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
616 616
617 617 return noise
618 618
619 619 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
620 620
621 621 if self.noise_estimation is not None:
622 622 # this was estimated by getNoise Operation defined in jroproc_spectra.py
623 623 return self.noise_estimation
624 624 else:
625 625 noise = self.getNoisebyHildebrand(
626 626 xmin_index, xmax_index, ymin_index, ymax_index)
627 627 return noise
628 628
629 629 def getFreqRangeTimeResponse(self, extrapoints=0):
630 630
631 631 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
632 632 freqrange = deltafreq * \
633 633 (numpy.arange(self.nFFTPoints + extrapoints) -
634 634 self.nFFTPoints / 2.) - deltafreq / 2
635 635
636 636 return freqrange
637 637
638 638 def getAcfRange(self, extrapoints=0):
639 #print "GET ACF RANGE"
639 640 #print "NFFTPoints",self.nFFTPoints
640 641 #print "IPPFactor", self.ippFactor
641 deltafreq = 10. / ( self.getFmax() / (self.nFFTPoints * self.ippFactor) )
642 #print "deltafreq",deltafreq
642 #deltafreq = 10. / ( self.getFmax() / (self.nFFTPoints * self.ippFactor) )
643 deltatime = 1./(self.getFmax()/ self.ippFactor)
644 #print "getFmax",self.getFmax()
643 645 #import time
644 646 #time.sleep(30)
645 freqrange = deltafreq * \
647 timerange = deltatime * \
646 648 (numpy.arange(self.nFFTPoints + extrapoints) -
647 self.nFFTPoints / 2.) - deltafreq / 2
648
649 return freqrange
649 self.nFFTPoints / 2.)
650 #- deltafreq / 2
651 #print "timerange",timerange
652 return timerange
650 653
651 654 def getFreqRange(self, extrapoints=0):
652 655
653 656 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
654 657 #print "deltafreq", deltafreq
655 658 freqrange = deltafreq * \
656 659 (numpy.arange(self.nFFTPoints + extrapoints) -
657 660 self.nFFTPoints / 2.) - deltafreq / 2
658 661 #print "freqrange",freqrange
659 662 return freqrange
660 663
661 664 def getVelRange(self, extrapoints=0):
662 665
663 666 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
664 667 velrange = deltav * (numpy.arange(self.nFFTPoints +
665 668 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
666 669
667 670 return velrange
668 671
669 672 def getNPairs(self):
670 673
671 674 return len(self.pairsList)
672 675
673 676 def getPairsIndexList(self):
674 677
675 678 return range(self.nPairs)
676 679
677 680 def getNormFactor(self):
678 681
679 682 pwcode = 1
680 683
681 684 if self.flagDecodeData:
682 685 pwcode = numpy.sum(self.code[0]**2)
683 686 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
684 687 normFactor = self.nProfiles * self.nIncohInt * \
685 688 self.nCohInt * pwcode * self.windowOfFilter
686 689
687 690 return normFactor
688 691
689 692 def getFlagCspc(self):
690 693
691 694 if self.data_cspc is None:
692 695 return True
693 696
694 697 return False
695 698
696 699 def getFlagDc(self):
697 700
698 701 if self.data_dc is None:
699 702 return True
700 703
701 704 return False
702 705
703 706 def getTimeInterval(self):
704 707
705 708 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
706 709
707 710 return timeInterval
708 711
709 712 def getPower(self):
710 713
711 714 factor = self.normFactor
712 715 z = self.data_spc / factor
713 716 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
714 717 avg = numpy.average(z, axis=1)
715 718
716 719 return 10 * numpy.log10(avg)
717 720
718 721 def getCoherence(self, pairsList=None, phase=False):
719 722
720 723 z = []
721 724 if pairsList is None:
722 725 pairsIndexList = self.pairsIndexList
723 726 else:
724 727 pairsIndexList = []
725 728 for pair in pairsList:
726 729 if pair not in self.pairsList:
727 730 raise ValueError, "Pair %s is not in dataOut.pairsList" % (
728 731 pair)
729 732 pairsIndexList.append(self.pairsList.index(pair))
730 733 for i in range(len(pairsIndexList)):
731 734 pair = self.pairsList[pairsIndexList[i]]
732 735 ccf = numpy.average(
733 736 self.data_cspc[pairsIndexList[i], :, :], axis=0)
734 737 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
735 738 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
736 739 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
737 740 if phase:
738 741 data = numpy.arctan2(avgcoherenceComplex.imag,
739 742 avgcoherenceComplex.real) * 180 / numpy.pi
740 743 else:
741 744 data = numpy.abs(avgcoherenceComplex)
742 745
743 746 z.append(data)
744 747
745 748 return numpy.array(z)
746 749
747 750 def setValue(self, value):
748 751
749 752 print "This property should not be initialized"
750 753
751 754 return
752 755
753 756 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
754 757 pairsIndexList = property(
755 758 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
756 759 normFactor = property(getNormFactor, setValue,
757 760 "I'm the 'getNormFactor' property.")
758 761 flag_cspc = property(getFlagCspc, setValue)
759 762 flag_dc = property(getFlagDc, setValue)
760 763 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
761 764 timeInterval = property(getTimeInterval, setValue,
762 765 "I'm the 'timeInterval' property")
763 766
764 767
765 768 class SpectraHeis(Spectra):
766 769
767 770 data_spc = None
768 771
769 772 data_cspc = None
770 773
771 774 data_dc = None
772 775
773 776 nFFTPoints = None
774 777
775 778 # nPairs = None
776 779
777 780 pairsList = None
778 781
779 782 nCohInt = None
780 783
781 784 nIncohInt = None
782 785
783 786 def __init__(self):
784 787
785 788 self.radarControllerHeaderObj = RadarControllerHeader()
786 789
787 790 self.systemHeaderObj = SystemHeader()
788 791
789 792 self.type = "SpectraHeis"
790 793
791 794 # self.dtype = None
792 795
793 796 # self.nChannels = 0
794 797
795 798 # self.nHeights = 0
796 799
797 800 self.nProfiles = None
798 801
799 802 self.heightList = None
800 803
801 804 self.channelList = None
802 805
803 806 # self.channelIndexList = None
804 807
805 808 self.flagNoData = True
806 809
807 810 self.flagDiscontinuousBlock = False
808 811
809 812 # self.nPairs = 0
810 813
811 814 self.utctime = None
812 815
813 816 self.blocksize = None
814 817
815 818 self.profileIndex = 0
816 819
817 820 self.nCohInt = 1
818 821
819 822 self.nIncohInt = 1
820 823
821 824 def getNormFactor(self):
822 825 pwcode = 1
823 826 if self.flagDecodeData:
824 827 pwcode = numpy.sum(self.code[0]**2)
825 828
826 829 normFactor = self.nIncohInt * self.nCohInt * pwcode
827 830
828 831 return normFactor
829 832
830 833 def getTimeInterval(self):
831 834
832 835 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
833 836
834 837 return timeInterval
835 838
836 839 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
837 840 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
838 841
839 842
840 843 class Fits(JROData):
841 844
842 845 heightList = None
843 846
844 847 channelList = None
845 848
846 849 flagNoData = True
847 850
848 851 flagDiscontinuousBlock = False
849 852
850 853 useLocalTime = False
851 854
852 855 utctime = None
853 856
854 857 timeZone = None
855 858
856 859 # ippSeconds = None
857 860
858 861 # timeInterval = None
859 862
860 863 nCohInt = None
861 864
862 865 nIncohInt = None
863 866
864 867 noise = None
865 868
866 869 windowOfFilter = 1
867 870
868 871 # Speed of ligth
869 872 C = 3e8
870 873
871 874 frequency = 49.92e6
872 875
873 876 realtime = False
874 877
875 878 def __init__(self):
876 879
877 880 self.type = "Fits"
878 881
879 882 self.nProfiles = None
880 883
881 884 self.heightList = None
882 885
883 886 self.channelList = None
884 887
885 888 # self.channelIndexList = None
886 889
887 890 self.flagNoData = True
888 891
889 892 self.utctime = None
890 893
891 894 self.nCohInt = 1
892 895
893 896 self.nIncohInt = 1
894 897
895 898 self.useLocalTime = True
896 899
897 900 self.profileIndex = 0
898 901
899 902 # self.utctime = None
900 903 # self.timeZone = None
901 904 # self.ltctime = None
902 905 # self.timeInterval = None
903 906 # self.header = None
904 907 # self.data_header = None
905 908 # self.data = None
906 909 # self.datatime = None
907 910 # self.flagNoData = False
908 911 # self.expName = ''
909 912 # self.nChannels = None
910 913 # self.nSamples = None
911 914 # self.dataBlocksPerFile = None
912 915 # self.comments = ''
913 916 #
914 917
915 918 def getltctime(self):
916 919
917 920 if self.useLocalTime:
918 921 return self.utctime - self.timeZone * 60
919 922
920 923 return self.utctime
921 924
922 925 def getDatatime(self):
923 926
924 927 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
925 928 return datatime
926 929
927 930 def getTimeRange(self):
928 931
929 932 datatime = []
930 933
931 934 datatime.append(self.ltctime)
932 935 datatime.append(self.ltctime + self.timeInterval)
933 936
934 937 datatime = numpy.array(datatime)
935 938
936 939 return datatime
937 940
938 941 def getHeiRange(self):
939 942
940 943 heis = self.heightList
941 944
942 945 return heis
943 946
944 947 def getNHeights(self):
945 948
946 949 return len(self.heightList)
947 950
948 951 def getNChannels(self):
949 952
950 953 return len(self.channelList)
951 954
952 955 def getChannelIndexList(self):
953 956
954 957 return range(self.nChannels)
955 958
956 959 def getNoise(self, type=1):
957 960
958 961 #noise = numpy.zeros(self.nChannels)
959 962
960 963 if type == 1:
961 964 noise = self.getNoisebyHildebrand()
962 965
963 966 if type == 2:
964 967 noise = self.getNoisebySort()
965 968
966 969 if type == 3:
967 970 noise = self.getNoisebyWindow()
968 971
969 972 return noise
970 973
971 974 def getTimeInterval(self):
972 975
973 976 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
974 977
975 978 return timeInterval
976 979
977 980 datatime = property(getDatatime, "I'm the 'datatime' property")
978 981 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
979 982 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
980 983 channelIndexList = property(
981 984 getChannelIndexList, "I'm the 'channelIndexList' property.")
982 985 noise = property(getNoise, "I'm the 'nHeights' property.")
983 986
984 987 ltctime = property(getltctime, "I'm the 'ltctime' property")
985 988 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
986 989
987 990
988 991 class Correlation(JROData):
989 992
990 993 noise = None
991 994
992 995 SNR = None
993 996
994 997 #--------------------------------------------------
995 998
996 999 mode = None
997 1000
998 1001 split = False
999 1002
1000 1003 data_cf = None
1001 1004
1002 1005 lags = None
1003 1006
1004 1007 lagRange = None
1005 1008
1006 1009 pairsList = None
1007 1010
1008 1011 normFactor = None
1009 1012
1010 1013 #--------------------------------------------------
1011 1014
1012 1015 # calculateVelocity = None
1013 1016
1014 1017 nLags = None
1015 1018
1016 1019 nPairs = None
1017 1020
1018 1021 nAvg = None
1019 1022
1020 1023 def __init__(self):
1021 1024 '''
1022 1025 Constructor
1023 1026 '''
1024 1027 self.radarControllerHeaderObj = RadarControllerHeader()
1025 1028
1026 1029 self.systemHeaderObj = SystemHeader()
1027 1030
1028 1031 self.type = "Correlation"
1029 1032
1030 1033 self.data = None
1031 1034
1032 1035 self.dtype = None
1033 1036
1034 1037 self.nProfiles = None
1035 1038
1036 1039 self.heightList = None
1037 1040
1038 1041 self.channelList = None
1039 1042
1040 1043 self.flagNoData = True
1041 1044
1042 1045 self.flagDiscontinuousBlock = False
1043 1046
1044 1047 self.utctime = None
1045 1048
1046 1049 self.timeZone = None
1047 1050
1048 1051 self.dstFlag = None
1049 1052
1050 1053 self.errorCount = None
1051 1054
1052 1055 self.blocksize = None
1053 1056
1054 1057 self.flagDecodeData = False # asumo q la data no esta decodificada
1055 1058
1056 1059 self.flagDeflipData = False # asumo q la data no esta sin flip
1057 1060
1058 1061 self.pairsList = None
1059 1062
1060 1063 self.nPoints = None
1061 1064
1062 1065 def getPairsList(self):
1063 1066
1064 1067 return self.pairsList
1065 1068
1066 1069 def getNoise(self, mode=2):
1067 1070
1068 1071 indR = numpy.where(self.lagR == 0)[0][0]
1069 1072 indT = numpy.where(self.lagT == 0)[0][0]
1070 1073
1071 1074 jspectra0 = self.data_corr[:, :, indR, :]
1072 1075 jspectra = copy.copy(jspectra0)
1073 1076
1074 1077 num_chan = jspectra.shape[0]
1075 1078 num_hei = jspectra.shape[2]
1076 1079
1077 1080 freq_dc = jspectra.shape[1] / 2
1078 1081 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
1079 1082
1080 1083 if ind_vel[0] < 0:
1081 1084 ind_vel[range(0, 1)] = ind_vel[range(0, 1)] + self.num_prof
1082 1085
1083 1086 if mode == 1:
1084 1087 jspectra[:, freq_dc, :] = (
1085 1088 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
1086 1089
1087 1090 if mode == 2:
1088 1091
1089 1092 vel = numpy.array([-2, -1, 1, 2])
1090 1093 xx = numpy.zeros([4, 4])
1091 1094
1092 1095 for fil in range(4):
1093 1096 xx[fil, :] = vel[fil]**numpy.asarray(range(4))
1094 1097
1095 1098 xx_inv = numpy.linalg.inv(xx)
1096 1099 xx_aux = xx_inv[0, :]
1097 1100
1098 1101 for ich in range(num_chan):
1099 1102 yy = jspectra[ich, ind_vel, :]
1100 1103 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
1101 1104
1102 1105 junkid = jspectra[ich, freq_dc, :] <= 0
1103 1106 cjunkid = sum(junkid)
1104 1107
1105 1108 if cjunkid.any():
1106 1109 jspectra[ich, freq_dc, junkid.nonzero()] = (
1107 1110 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
1108 1111
1109 1112 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
1110 1113
1111 1114 return noise
1112 1115
1113 1116 def getTimeInterval(self):
1114 1117
1115 1118 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1116 1119
1117 1120 return timeInterval
1118 1121
1119 1122 def splitFunctions(self):
1120 1123
1121 1124 pairsList = self.pairsList
1122 1125 ccf_pairs = []
1123 1126 acf_pairs = []
1124 1127 ccf_ind = []
1125 1128 acf_ind = []
1126 1129 for l in range(len(pairsList)):
1127 1130 chan0 = pairsList[l][0]
1128 1131 chan1 = pairsList[l][1]
1129 1132
1130 1133 # Obteniendo pares de Autocorrelacion
1131 1134 if chan0 == chan1:
1132 1135 acf_pairs.append(chan0)
1133 1136 acf_ind.append(l)
1134 1137 else:
1135 1138 ccf_pairs.append(pairsList[l])
1136 1139 ccf_ind.append(l)
1137 1140
1138 1141 data_acf = self.data_cf[acf_ind]
1139 1142 data_ccf = self.data_cf[ccf_ind]
1140 1143
1141 1144 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1142 1145
1143 1146 def getNormFactor(self):
1144 1147 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1145 1148 acf_pairs = numpy.array(acf_pairs)
1146 1149 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1147 1150
1148 1151 for p in range(self.nPairs):
1149 1152 pair = self.pairsList[p]
1150 1153
1151 1154 ch0 = pair[0]
1152 1155 ch1 = pair[1]
1153 1156
1154 1157 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1155 1158 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1156 1159 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1157 1160
1158 1161 return normFactor
1159 1162
1160 1163 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1161 1164 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1162 1165
1163 1166
1164 1167 class Parameters(Spectra):
1165 1168
1166 1169 experimentInfo = None # Information about the experiment
1167 1170
1168 1171 # Information from previous data
1169 1172
1170 1173 inputUnit = None # Type of data to be processed
1171 1174
1172 1175 operation = None # Type of operation to parametrize
1173 1176
1174 1177 # normFactor = None #Normalization Factor
1175 1178
1176 1179 groupList = None # List of Pairs, Groups, etc
1177 1180
1178 1181 # Parameters
1179 1182
1180 1183 data_param = None # Parameters obtained
1181 1184
1182 1185 data_pre = None # Data Pre Parametrization
1183 1186
1184 1187 data_SNR = None # Signal to Noise Ratio
1185 1188
1186 1189 # heightRange = None #Heights
1187 1190
1188 1191 abscissaList = None # Abscissa, can be velocities, lags or time
1189 1192
1190 1193 # noise = None #Noise Potency
1191 1194
1192 1195 utctimeInit = None # Initial UTC time
1193 1196
1194 1197 paramInterval = None # Time interval to calculate Parameters in seconds
1195 1198
1196 1199 useLocalTime = True
1197 1200
1198 1201 # Fitting
1199 1202
1200 1203 data_error = None # Error of the estimation
1201 1204
1202 1205 constants = None
1203 1206
1204 1207 library = None
1205 1208
1206 1209 # Output signal
1207 1210
1208 1211 outputInterval = None # Time interval to calculate output signal in seconds
1209 1212
1210 1213 data_output = None # Out signal
1211 1214
1212 1215 nAvg = None
1213 1216
1214 1217 noise_estimation = None
1215 1218
1216 1219 GauSPC = None # Fit gaussian SPC
1217 1220
1218 1221 def __init__(self):
1219 1222 '''
1220 1223 Constructor
1221 1224 '''
1222 1225 self.radarControllerHeaderObj = RadarControllerHeader()
1223 1226
1224 1227 self.systemHeaderObj = SystemHeader()
1225 1228
1226 1229 self.type = "Parameters"
1227 1230
1228 1231 def getTimeRange1(self, interval):
1229 1232
1230 1233 datatime = []
1231 1234
1232 1235 if self.useLocalTime:
1233 1236 time1 = self.utctimeInit - self.timeZone * 60
1234 1237 else:
1235 1238 time1 = self.utctimeInit
1236 1239
1237 1240 datatime.append(time1)
1238 1241 datatime.append(time1 + interval)
1239 1242 datatime = numpy.array(datatime)
1240 1243
1241 1244 return datatime
1242 1245
1243 1246 def getTimeInterval(self):
1244 1247
1245 1248 if hasattr(self, 'timeInterval1'):
1246 1249 return self.timeInterval1
1247 1250 else:
1248 1251 return self.paramInterval
1249 1252
1250 1253 def setValue(self, value):
1251 1254
1252 1255 print "This property should not be initialized"
1253 1256
1254 1257 return
1255 1258
1256 1259 def getNoise(self):
1257 1260
1258 1261 return self.spc_noise
1259 1262
1260 1263 timeInterval = property(getTimeInterval)
1261 1264 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
Show file before | Show file after | Hide comments
@@ -1,1737 +1,1738
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
10 10 from figure import Figure, isRealtime, isTimeInHourRange
11 11 from plotting_codes import *
12 12
13 13
14 14 class SpectraPlot(Figure):
15 15
16 16 isConfig = None
17 17 __nsubplots = None
18 18
19 19 WIDTHPROF = None
20 20 HEIGHTPROF = None
21 21 PREFIX = 'spc'
22 22
23 23 def __init__(self, **kwargs):
24 24 Figure.__init__(self, **kwargs)
25 25 self.isConfig = False
26 26 self.__nsubplots = 1
27 27
28 28 self.WIDTH = 250
29 29 self.HEIGHT = 250
30 30 self.WIDTHPROF = 120
31 31 self.HEIGHTPROF = 0
32 32 self.counter_imagwr = 0
33 33
34 34 self.PLOT_CODE = SPEC_CODE
35 35
36 36 self.FTP_WEI = None
37 37 self.EXP_CODE = None
38 38 self.SUB_EXP_CODE = None
39 39 self.PLOT_POS = None
40 40
41 41 self.__xfilter_ena = False
42 42 self.__yfilter_ena = False
43 43
44 44 def getSubplots(self):
45 45
46 46 ncol = int(numpy.sqrt(self.nplots)+0.9)
47 47 nrow = int(self.nplots*1./ncol + 0.9)
48 48
49 49 return nrow, ncol
50 50
51 51 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
52 52
53 53 self.__showprofile = showprofile
54 54 self.nplots = nplots
55 55
56 56 ncolspan = 1
57 57 colspan = 1
58 58 if showprofile:
59 59 ncolspan = 3
60 60 colspan = 2
61 61 self.__nsubplots = 2
62 62
63 63 self.createFigure(id = id,
64 64 wintitle = wintitle,
65 65 widthplot = self.WIDTH + self.WIDTHPROF,
66 66 heightplot = self.HEIGHT + self.HEIGHTPROF,
67 67 show=show)
68 68
69 69 nrow, ncol = self.getSubplots()
70 70
71 71 counter = 0
72 72 for y in range(nrow):
73 73 for x in range(ncol):
74 74
75 75 if counter >= self.nplots:
76 76 break
77 77
78 78 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
79 79
80 80 if showprofile:
81 81 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
82 82
83 83 counter += 1
84 84
85 85 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
86 86 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
87 87 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
88 88 server=None, folder=None, username=None, password=None,
89 89 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
90 90 xaxis="frequency", colormap='jet', normFactor=None):
91 91
92 92 """
93 93
94 94 Input:
95 95 dataOut :
96 96 id :
97 97 wintitle :
98 98 channelList :
99 99 showProfile :
100 100 xmin : None,
101 101 xmax : None,
102 102 ymin : None,
103 103 ymax : None,
104 104 zmin : None,
105 105 zmax : None
106 106 """
107 107 if realtime:
108 108 if not(isRealtime(utcdatatime = dataOut.utctime)):
109 109 print 'Skipping this plot function'
110 110 return
111 111
112 112 if channelList == None:
113 113 channelIndexList = dataOut.channelIndexList
114 114 else:
115 115 channelIndexList = []
116 116 for channel in channelList:
117 117 if channel not in dataOut.channelList:
118 118 raise ValueError, "Channel %d is not in dataOut.channelList" %channel
119 119 channelIndexList.append(dataOut.channelList.index(channel))
120 120
121 121 if normFactor is None:
122 122 factor = dataOut.normFactor
123 123 else:
124 124 factor = normFactor
125 125 if xaxis == "frequency":
126 126 x = dataOut.getFreqRange(1)/1000.
127 127 xlabel = "Frequency (kHz)"
128 128
129 129 elif xaxis == "time":
130 130 x = dataOut.getAcfRange(1)
131 131 xlabel = "Time (ms)"
132 132
133 133 else:
134 134 x = dataOut.getVelRange(1)
135 135 xlabel = "Velocity (m/s)"
136 136
137 137 ylabel = "Range (Km)"
138 138
139 139 y = dataOut.getHeiRange()
140 140
141 141 z = dataOut.data_spc/factor
142 142 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
143 143 zdB = 10*numpy.log10(z)
144 144
145 145 #print "a000",dataOut.data_spc.dtype
146 146 avg = numpy.average(z, axis=1)
147 147 avgdB = 10*numpy.log10(avg)
148 148 #print "before plot"
149 149 noise = dataOut.getNoise()/factor
150 150 noisedB = 10*numpy.log10(noise)
151 151
152 152 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
153 153 title = wintitle + " Spectra"
154 154 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
155 155 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
156 156
157 157 if not self.isConfig:
158 158
159 159 nplots = len(channelIndexList)
160 160
161 161 self.setup(id=id,
162 162 nplots=nplots,
163 163 wintitle=wintitle,
164 164 showprofile=showprofile,
165 165 show=show)
166 166
167 167 if xmin == None: xmin = numpy.nanmin(x)
168 168 if xmax == None: xmax = numpy.nanmax(x)
169 169 if ymin == None: ymin = numpy.nanmin(y)
170 170 if ymax == None: ymax = numpy.nanmax(y)
171 171 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
172 172 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
173 173
174 174 self.FTP_WEI = ftp_wei
175 175 self.EXP_CODE = exp_code
176 176 self.SUB_EXP_CODE = sub_exp_code
177 177 self.PLOT_POS = plot_pos
178 178
179 179 self.isConfig = True
180 180
181 181 self.setWinTitle(title)
182 182
183 183 for i in range(self.nplots):
184 184 index = channelIndexList[i]
185 185 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
186 186 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
187 187 if len(dataOut.beam.codeList) != 0:
188 188 title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime)
189 189
190 190 axes = self.axesList[i*self.__nsubplots]
191 191 axes.pcolor(x, y, zdB[index,:,:],
192 192 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
193 193 xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap,
194 194 ticksize=9, cblabel='')
195 195
196 196 if self.__showprofile:
197 197 axes = self.axesList[i*self.__nsubplots +1]
198 198 axes.pline(avgdB[index,:], y,
199 199 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
200 200 xlabel='dB', ylabel='', title='',
201 201 ytick_visible=False,
202 202 grid='x')
203 203
204 204 noiseline = numpy.repeat(noisedB[index], len(y))
205 205 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
206 206
207 207 self.draw()
208 208
209 209 if figfile == None:
210 210 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
211 211 name = str_datetime
212 212 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
213 213 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
214 214 figfile = self.getFilename(name)
215 215
216 216 self.save(figpath=figpath,
217 217 figfile=figfile,
218 218 save=save,
219 219 ftp=ftp,
220 220 wr_period=wr_period,
221 221 thisDatetime=thisDatetime)
222 222
223 223 class ACFPlot(Figure):
224 224
225 225 isConfig = None
226 226 __nsubplots = None
227 227
228 228 WIDTHPROF = None
229 229 HEIGHTPROF = None
230 230 PREFIX = 'acf'
231 231
232 232 def __init__(self, **kwargs):
233 233 Figure.__init__(self, **kwargs)
234 234 self.isConfig = False
235 235 self.__nsubplots = 1
236 236
237 237 self.PLOT_CODE = ACF_CODE
238 238
239 239 self.WIDTH = 900
240 240 self.HEIGHT = 700
241 241 self.counter_imagwr = 0
242 242
243 243 def getSubplots(self):
244 244 ncol = 1
245 245 nrow = 1
246 246
247 247 return nrow, ncol
248 248
249 249 def setup(self, id, nplots, wintitle, show):
250 250
251 251 self.nplots = nplots
252 252
253 253 ncolspan = 1
254 254 colspan = 1
255 255
256 256 self.createFigure(id = id,
257 257 wintitle = wintitle,
258 258 widthplot = self.WIDTH,
259 259 heightplot = self.HEIGHT,
260 260 show=show)
261 261
262 262 nrow, ncol = self.getSubplots()
263 263
264 264 counter = 0
265 265 for y in range(nrow):
266 266 for x in range(ncol):
267 267 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
268 268
269 269 def run(self, dataOut, id, wintitle="", channelList=None,channel=None,nSamples=None,
270 270 nSampleList= None,resolutionFactor=None,
271 271 xmin=None, xmax=None, ymin=None, ymax=None,
272 272 save=False, figpath='./', figfile=None, show=True,
273 273 ftp=False, wr_period=1, server=None,
274 274 folder=None, username=None, password=None,
275 275 xaxis="frequency"):
276 276
277 277 channel0 = channel
278 278 nSamples = nSamples
279 279 resFactor = resolutionFactor
280 280
281 281 if nSamples == None:
282 282 nSamples = 20
283 283
284 284 if resFactor == None:
285 285 resFactor = 5
286 286 #else:
287 287 # if nSamples not in dataOut.channelList:
288 288 # raise ValueError, "Channel %d is not in %s dataOut.channelList"%(channel0, dataOut.channelList)
289 289
290 290 if channel0 == None:
291 291 channel0 = 0
292 292 else:
293 293 if channel0 not in dataOut.channelList:
294 294 raise ValueError, "Channel %d is not in %s dataOut.channelList"%(channel0, dataOut.channelList)
295 295
296 296 if channelList == None:
297 297 channelIndexList = dataOut.channelIndexList
298 298 channelList = dataOut.channelList
299 299 else:
300 300 channelIndexList = []
301 301 for channel in channelList:
302 302 if channel not in dataOut.channelList:
303 303 raise ValueError, "Channel %d is not in dataOut.channelList"
304 304 channelIndexList.append(dataOut.channelList.index(channel))
305 305
306 306 #z = dataOut.data_spc/factor
307 307 factor = dataOut.normFactor
308 308 y = dataOut.getHeiRange()
309 309 deltaHeight = dataOut.heightList[1]-dataOut.heightList[0]
310 310 z = dataOut.data_acf
311 311
312 312 #import matplotlib.pyplot as plt
313 313 #plt.plot(z[0,:,85])
314 314 #plt.show()
315 315 #import time
316 316 #time.sleep(20)
317 317
318 318
319 319 #z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
320 320 shape = dataOut.data_acf.shape
321 print "shape_plot",shape
321 #print "shape_plot",shape
322 322 hei_index = numpy.arange(shape[2])
323 323 hei_plot = numpy.arange(nSamples)*resFactor
324 324
325 325 if nSampleList is not None:
326 326 for nsample in nSampleList:
327 327 if nsample not in dataOut.heightList/deltaHeight:
328 328 raise ValueError, "nsample %d is not in %s dataOut.heightList"%(nsample,dataOut.heightList)
329 329
330 330 if nSampleList is not None:
331 331 hei_plot = numpy.array(nSampleList)*resFactor
332 332
333 333 if hei_plot[-1] >= hei_index[-1]:
334 334 print ("La cantidad de puntos en altura es %d y la resolucion es %f Km"%(hei_plot.shape[0],deltaHeight*resFactor ))
335 335 raise ValueError, "resFactor %d multiplicado por el valor de %d nSamples es mayor a %d cantidad total de puntos"%(resFactor,nSamples,hei_index[-1])
336 336
337 337 #escalamiento -1 a 1 a resolucion (factor de resolucion en altura)* deltaHeight
338 338 min = numpy.min(z[0,:,0])
339 339 max =numpy.max(z[0,:,0])
340 340 for i in range(shape[0]):
341 341 for j in range(shape[2]):
342 342 z[i,:,j]= (((z[i,:,j]-min)/(max-min))*deltaHeight*resFactor + j*deltaHeight)
343 343
344 344 if xaxis == "frequency":
345 345 x = dataOut.getFreqRange()/1000.
346 346 zdB = 10*numpy.log10(z[channel0,:,hei_plot])
347 347 xlabel = "Frequency (kHz)"
348 348 ylabel = "Power (dB)"
349 349
350 350 elif xaxis == "time":
351 delta= dataOut.getAcfRange()[1]-dataOut.getAcfRange()[0]
352 x = dataOut.getAcfRange()+delta/2.0
351 #delta= dataOut.getAcfRange()[1]-dataOut.getAcfRange()[0]
352 #x = dataOut.getAcfRange()+delta/2.0
353 x = dataOut.getAcfRange()*1000
353 354 zdB = z[channel0,:,hei_plot]
354 355 xlabel = "Time (ms)"
355 356 ylabel = "ACF"
356 357
357 358 else:
358 359 x = dataOut.getVelRange()
359 360 zdB = 10*numpy.log10(z[channel0,:,hei_plot])
360 361 xlabel = "Velocity (m/s)"
361 362 ylabel = "Power (dB)"
362 363
363 364 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
364 365 title = wintitle + " ACF Plot Ch %s %s" %(channel0,thisDatetime.strftime("%d-%b-%Y"))
365 366
366 367 if not self.isConfig:
367 368
368 369 nplots = 1
369 370
370 371 self.setup(id=id,
371 372 nplots=nplots,
372 373 wintitle=wintitle,
373 374 show=show)
374 375
375 if xmin == None: xmin = numpy.nanmin(x)*0.9
376 if xmax == None: xmax = numpy.nanmax(x)*1.1
376 if xmin == None: xmin = numpy.nanmin(x)#*0.9
377 if xmax == None: xmax = numpy.nanmax(x)#*1.1
377 378 if ymin == None: ymin = numpy.nanmin(zdB)
378 379 if ymax == None: ymax = numpy.nanmax(zdB)
379 380
380 381 print ("El parametro resFactor es %d y la resolucion en altura es %f"%(resFactor,deltaHeight ))
381 382 print ("La cantidad de puntos en altura es %d y la nueva resolucion es %f Km"%(hei_plot.shape[0],deltaHeight*resFactor ))
382 383 print ("La altura maxima es %d Km"%(hei_plot[-1]*deltaHeight ))
383 384
384 385 self.isConfig = True
385 386
386 387 self.setWinTitle(title)
387 388
388 389 title = "ACF Plot: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
389 390 axes = self.axesList[0]
390 391
391 392 legendlabels = ["Range = %dKm" %y[i] for i in hei_plot]
392 393
393 394 axes.pmultilineyaxis( x, zdB,
394 395 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
395 396 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
396 397 ytick_visible=True, nxticks=5,
397 398 grid='x')
398 399
399 400 self.draw()
400 401
401 402 if figfile == None:
402 403 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
403 404 name = str_datetime
404 405 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
405 406 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
406 407 figfile = self.getFilename(name)
407 408
408 409 self.save(figpath=figpath,
409 410 figfile=figfile,
410 411 save=save,
411 412 ftp=ftp,
412 413 wr_period=wr_period,
413 414 thisDatetime=thisDatetime)
414 415
415 416
416 417
417 418 class CrossSpectraPlot(Figure):
418 419
419 420 isConfig = None
420 421 __nsubplots = None
421 422
422 423 WIDTH = None
423 424 HEIGHT = None
424 425 WIDTHPROF = None
425 426 HEIGHTPROF = None
426 427 PREFIX = 'cspc'
427 428
428 429 def __init__(self, **kwargs):
429 430 Figure.__init__(self, **kwargs)
430 431 self.isConfig = False
431 432 self.__nsubplots = 4
432 433 self.counter_imagwr = 0
433 434 self.WIDTH = 250
434 435 self.HEIGHT = 250
435 436 self.WIDTHPROF = 0
436 437 self.HEIGHTPROF = 0
437 438
438 439 self.PLOT_CODE = CROSS_CODE
439 440 self.FTP_WEI = None
440 441 self.EXP_CODE = None
441 442 self.SUB_EXP_CODE = None
442 443 self.PLOT_POS = None
443 444
444 445 def getSubplots(self):
445 446
446 447 ncol = 4
447 448 nrow = self.nplots
448 449
449 450 return nrow, ncol
450 451
451 452 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
452 453
453 454 self.__showprofile = showprofile
454 455 self.nplots = nplots
455 456
456 457 ncolspan = 1
457 458 colspan = 1
458 459
459 460 self.createFigure(id = id,
460 461 wintitle = wintitle,
461 462 widthplot = self.WIDTH + self.WIDTHPROF,
462 463 heightplot = self.HEIGHT + self.HEIGHTPROF,
463 464 show=True)
464 465
465 466 nrow, ncol = self.getSubplots()
466 467
467 468 counter = 0
468 469 for y in range(nrow):
469 470 for x in range(ncol):
470 471 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
471 472
472 473 counter += 1
473 474
474 475 def run(self, dataOut, id, wintitle="", pairsList=None,
475 476 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
476 477 coh_min=None, coh_max=None, phase_min=None, phase_max=None,
477 478 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
478 479 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
479 480 server=None, folder=None, username=None, password=None,
480 481 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None,
481 482 xaxis='frequency'):
482 483
483 484 """
484 485
485 486 Input:
486 487 dataOut :
487 488 id :
488 489 wintitle :
489 490 channelList :
490 491 showProfile :
491 492 xmin : None,
492 493 xmax : None,
493 494 ymin : None,
494 495 ymax : None,
495 496 zmin : None,
496 497 zmax : None
497 498 """
498 499
499 500 if pairsList == None:
500 501 pairsIndexList = dataOut.pairsIndexList
501 502 else:
502 503 pairsIndexList = []
503 504 for pair in pairsList:
504 505 if pair not in dataOut.pairsList:
505 506 raise ValueError, "Pair %s is not in dataOut.pairsList" %str(pair)
506 507 pairsIndexList.append(dataOut.pairsList.index(pair))
507 508
508 509 if not pairsIndexList:
509 510 return
510 511
511 512 if len(pairsIndexList) > 4:
512 513 pairsIndexList = pairsIndexList[0:4]
513 514
514 515 if normFactor is None:
515 516 factor = dataOut.normFactor
516 517 else:
517 518 factor = normFactor
518 519 x = dataOut.getVelRange(1)
519 520 y = dataOut.getHeiRange()
520 521 z = dataOut.data_spc[:,:,:]/factor
521 522 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
522 523
523 524 noise = dataOut.noise/factor
524 525
525 526 zdB = 10*numpy.log10(z)
526 527 noisedB = 10*numpy.log10(noise)
527 528
528 529 if coh_min == None:
529 530 coh_min = 0.0
530 531 if coh_max == None:
531 532 coh_max = 1.0
532 533
533 534 if phase_min == None:
534 535 phase_min = -180
535 536 if phase_max == None:
536 537 phase_max = 180
537 538
538 539 #thisDatetime = dataOut.datatime
539 540 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
540 541 title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
541 542 # xlabel = "Velocity (m/s)"
542 543 ylabel = "Range (Km)"
543 544
544 545 if xaxis == "frequency":
545 546 x = dataOut.getFreqRange(1)/1000.
546 547 xlabel = "Frequency (kHz)"
547 548
548 549 elif xaxis == "time":
549 550 x = dataOut.getAcfRange(1)
550 551 xlabel = "Time (ms)"
551 552
552 553 else:
553 554 x = dataOut.getVelRange(1)
554 555 xlabel = "Velocity (m/s)"
555 556
556 557 if not self.isConfig:
557 558
558 559 nplots = len(pairsIndexList)
559 560
560 561 self.setup(id=id,
561 562 nplots=nplots,
562 563 wintitle=wintitle,
563 564 showprofile=False,
564 565 show=show)
565 566
566 567 avg = numpy.abs(numpy.average(z, axis=1))
567 568 avgdB = 10*numpy.log10(avg)
568 569
569 570 if xmin == None: xmin = numpy.nanmin(x)
570 571 if xmax == None: xmax = numpy.nanmax(x)
571 572 if ymin == None: ymin = numpy.nanmin(y)
572 573 if ymax == None: ymax = numpy.nanmax(y)
573 574 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
574 575 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
575 576
576 577 self.FTP_WEI = ftp_wei
577 578 self.EXP_CODE = exp_code
578 579 self.SUB_EXP_CODE = sub_exp_code
579 580 self.PLOT_POS = plot_pos
580 581
581 582 self.isConfig = True
582 583
583 584 self.setWinTitle(title)
584 585
585 586 for i in range(self.nplots):
586 587 pair = dataOut.pairsList[pairsIndexList[i]]
587 588
588 589 chan_index0 = dataOut.channelList.index(pair[0])
589 590 chan_index1 = dataOut.channelList.index(pair[1])
590 591
591 592 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
592 593 title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[chan_index0], str_datetime)
593 594 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index0,:,:]/factor)
594 595 axes0 = self.axesList[i*self.__nsubplots]
595 596 axes0.pcolor(x, y, zdB,
596 597 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
597 598 xlabel=xlabel, ylabel=ylabel, title=title,
598 599 ticksize=9, colormap=power_cmap, cblabel='')
599 600
600 601 title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[chan_index1], str_datetime)
601 602 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index1,:,:]/factor)
602 603 axes0 = self.axesList[i*self.__nsubplots+1]
603 604 axes0.pcolor(x, y, zdB,
604 605 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
605 606 xlabel=xlabel, ylabel=ylabel, title=title,
606 607 ticksize=9, colormap=power_cmap, cblabel='')
607 608
608 609 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:])
609 610 coherence = numpy.abs(coherenceComplex)
610 611 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
611 612 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
612 613
613 614 title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1])
614 615 axes0 = self.axesList[i*self.__nsubplots+2]
615 616 axes0.pcolor(x, y, coherence,
616 617 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=coh_min, zmax=coh_max,
617 618 xlabel=xlabel, ylabel=ylabel, title=title,
618 619 ticksize=9, colormap=coherence_cmap, cblabel='')
619 620
620 621 title = "Phase Ch%d * Ch%d" %(pair[0], pair[1])
621 622 axes0 = self.axesList[i*self.__nsubplots+3]
622 623 axes0.pcolor(x, y, phase,
623 624 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
624 625 xlabel=xlabel, ylabel=ylabel, title=title,
625 626 ticksize=9, colormap=phase_cmap, cblabel='')
626 627
627 628
628 629
629 630 self.draw()
630 631
631 632 self.save(figpath=figpath,
632 633 figfile=figfile,
633 634 save=save,
634 635 ftp=ftp,
635 636 wr_period=wr_period,
636 637 thisDatetime=thisDatetime)
637 638
638 639
639 640 class RTIPlot(Figure):
640 641
641 642 __isConfig = None
642 643 __nsubplots = None
643 644
644 645 WIDTHPROF = None
645 646 HEIGHTPROF = None
646 647 PREFIX = 'rti'
647 648
648 649 def __init__(self, **kwargs):
649 650
650 651 Figure.__init__(self, **kwargs)
651 652 self.timerange = None
652 653 self.isConfig = False
653 654 self.__nsubplots = 1
654 655
655 656 self.WIDTH = 800
656 657 self.HEIGHT = 180
657 658 self.WIDTHPROF = 120
658 659 self.HEIGHTPROF = 0
659 660 self.counter_imagwr = 0
660 661
661 662 self.PLOT_CODE = RTI_CODE
662 663
663 664 self.FTP_WEI = None
664 665 self.EXP_CODE = None
665 666 self.SUB_EXP_CODE = None
666 667 self.PLOT_POS = None
667 668 self.tmin = None
668 669 self.tmax = None
669 670
670 671 self.xmin = None
671 672 self.xmax = None
672 673
673 674 self.figfile = None
674 675
675 676 def getSubplots(self):
676 677
677 678 ncol = 1
678 679 nrow = self.nplots
679 680
680 681 return nrow, ncol
681 682
682 683 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
683 684
684 685 self.__showprofile = showprofile
685 686 self.nplots = nplots
686 687
687 688 ncolspan = 1
688 689 colspan = 1
689 690 if showprofile:
690 691 ncolspan = 7
691 692 colspan = 6
692 693 self.__nsubplots = 2
693 694
694 695 self.createFigure(id = id,
695 696 wintitle = wintitle,
696 697 widthplot = self.WIDTH + self.WIDTHPROF,
697 698 heightplot = self.HEIGHT + self.HEIGHTPROF,
698 699 show=show)
699 700
700 701 nrow, ncol = self.getSubplots()
701 702
702 703 counter = 0
703 704 for y in range(nrow):
704 705 for x in range(ncol):
705 706
706 707 if counter >= self.nplots:
707 708 break
708 709
709 710 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
710 711
711 712 if showprofile:
712 713 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
713 714
714 715 counter += 1
715 716
716 717 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
717 718 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
718 719 timerange=None, colormap='jet',
719 720 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
720 721 server=None, folder=None, username=None, password=None,
721 722 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, HEIGHT=None):
722 723
723 724 """
724 725
725 726 Input:
726 727 dataOut :
727 728 id :
728 729 wintitle :
729 730 channelList :
730 731 showProfile :
731 732 xmin : None,
732 733 xmax : None,
733 734 ymin : None,
734 735 ymax : None,
735 736 zmin : None,
736 737 zmax : None
737 738 """
738 739
739 740 #colormap = kwargs.get('colormap', 'jet')
740 741 if HEIGHT is not None:
741 742 self.HEIGHT = HEIGHT
742 743
743 744 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
744 745 return
745 746
746 747 if channelList == None:
747 748 channelIndexList = dataOut.channelIndexList
748 749 else:
749 750 channelIndexList = []
750 751 for channel in channelList:
751 752 if channel not in dataOut.channelList:
752 753 raise ValueError, "Channel %d is not in dataOut.channelList"
753 754 channelIndexList.append(dataOut.channelList.index(channel))
754 755
755 756 if normFactor is None:
756 757 factor = dataOut.normFactor
757 758 else:
758 759 factor = normFactor
759 760
760 761 # factor = dataOut.normFactor
761 762 x = dataOut.getTimeRange()
762 763 y = dataOut.getHeiRange()
763 764
764 765 z = dataOut.data_spc/factor
765 766 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
766 767 avg = numpy.average(z, axis=1)
767 768 avgdB = 10.*numpy.log10(avg)
768 769 # avgdB = dataOut.getPower()
769 770
770 771
771 772 thisDatetime = dataOut.datatime
772 773 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
773 774 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
774 775 xlabel = ""
775 776 ylabel = "Range (Km)"
776 777
777 778 update_figfile = False
778 779
779 780 if dataOut.ltctime >= self.xmax:
780 781 self.counter_imagwr = wr_period
781 782 self.isConfig = False
782 783 update_figfile = True
783 784
784 785 if not self.isConfig:
785 786
786 787 nplots = len(channelIndexList)
787 788
788 789 self.setup(id=id,
789 790 nplots=nplots,
790 791 wintitle=wintitle,
791 792 showprofile=showprofile,
792 793 show=show)
793 794
794 795 if timerange != None:
795 796 self.timerange = timerange
796 797
797 798 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
798 799
799 800 noise = dataOut.noise/factor
800 801 noisedB = 10*numpy.log10(noise)
801 802
802 803 if ymin == None: ymin = numpy.nanmin(y)
803 804 if ymax == None: ymax = numpy.nanmax(y)
804 805 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
805 806 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
806 807
807 808 self.FTP_WEI = ftp_wei
808 809 self.EXP_CODE = exp_code
809 810 self.SUB_EXP_CODE = sub_exp_code
810 811 self.PLOT_POS = plot_pos
811 812
812 813 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
813 814 self.isConfig = True
814 815 self.figfile = figfile
815 816 update_figfile = True
816 817
817 818 self.setWinTitle(title)
818 819
819 820 for i in range(self.nplots):
820 821 index = channelIndexList[i]
821 822 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
822 823 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
823 824 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
824 825 axes = self.axesList[i*self.__nsubplots]
825 826 zdB = avgdB[index].reshape((1,-1))
826 827 axes.pcolorbuffer(x, y, zdB,
827 828 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
828 829 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
829 830 ticksize=9, cblabel='', cbsize="1%", colormap=colormap)
830 831
831 832 if self.__showprofile:
832 833 axes = self.axesList[i*self.__nsubplots +1]
833 834 axes.pline(avgdB[index], y,
834 835 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
835 836 xlabel='dB', ylabel='', title='',
836 837 ytick_visible=False,
837 838 grid='x')
838 839
839 840 self.draw()
840 841
841 842 self.save(figpath=figpath,
842 843 figfile=figfile,
843 844 save=save,
844 845 ftp=ftp,
845 846 wr_period=wr_period,
846 847 thisDatetime=thisDatetime,
847 848 update_figfile=update_figfile)
848 849
849 850 class CoherenceMap(Figure):
850 851 isConfig = None
851 852 __nsubplots = None
852 853
853 854 WIDTHPROF = None
854 855 HEIGHTPROF = None
855 856 PREFIX = 'cmap'
856 857
857 858 def __init__(self, **kwargs):
858 859 Figure.__init__(self, **kwargs)
859 860 self.timerange = 2*60*60
860 861 self.isConfig = False
861 862 self.__nsubplots = 1
862 863
863 864 self.WIDTH = 800
864 865 self.HEIGHT = 180
865 866 self.WIDTHPROF = 120
866 867 self.HEIGHTPROF = 0
867 868 self.counter_imagwr = 0
868 869
869 870 self.PLOT_CODE = COH_CODE
870 871
871 872 self.FTP_WEI = None
872 873 self.EXP_CODE = None
873 874 self.SUB_EXP_CODE = None
874 875 self.PLOT_POS = None
875 876 self.counter_imagwr = 0
876 877
877 878 self.xmin = None
878 879 self.xmax = None
879 880
880 881 def getSubplots(self):
881 882 ncol = 1
882 883 nrow = self.nplots*2
883 884
884 885 return nrow, ncol
885 886
886 887 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
887 888 self.__showprofile = showprofile
888 889 self.nplots = nplots
889 890
890 891 ncolspan = 1
891 892 colspan = 1
892 893 if showprofile:
893 894 ncolspan = 7
894 895 colspan = 6
895 896 self.__nsubplots = 2
896 897
897 898 self.createFigure(id = id,
898 899 wintitle = wintitle,
899 900 widthplot = self.WIDTH + self.WIDTHPROF,
900 901 heightplot = self.HEIGHT + self.HEIGHTPROF,
901 902 show=True)
902 903
903 904 nrow, ncol = self.getSubplots()
904 905
905 906 for y in range(nrow):
906 907 for x in range(ncol):
907 908
908 909 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
909 910
910 911 if showprofile:
911 912 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
912 913
913 914 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
914 915 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
915 916 timerange=None, phase_min=None, phase_max=None,
916 917 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
917 918 coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
918 919 server=None, folder=None, username=None, password=None,
919 920 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
920 921
921 922 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
922 923 return
923 924
924 925 if pairsList == None:
925 926 pairsIndexList = dataOut.pairsIndexList
926 927 else:
927 928 pairsIndexList = []
928 929 for pair in pairsList:
929 930 if pair not in dataOut.pairsList:
930 931 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
931 932 pairsIndexList.append(dataOut.pairsList.index(pair))
932 933
933 934 if pairsIndexList == []:
934 935 return
935 936
936 937 if len(pairsIndexList) > 4:
937 938 pairsIndexList = pairsIndexList[0:4]
938 939
939 940 if phase_min == None:
940 941 phase_min = -180
941 942 if phase_max == None:
942 943 phase_max = 180
943 944
944 945 x = dataOut.getTimeRange()
945 946 y = dataOut.getHeiRange()
946 947
947 948 thisDatetime = dataOut.datatime
948 949
949 950 title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
950 951 xlabel = ""
951 952 ylabel = "Range (Km)"
952 953 update_figfile = False
953 954
954 955 if not self.isConfig:
955 956 nplots = len(pairsIndexList)
956 957 self.setup(id=id,
957 958 nplots=nplots,
958 959 wintitle=wintitle,
959 960 showprofile=showprofile,
960 961 show=show)
961 962
962 963 if timerange != None:
963 964 self.timerange = timerange
964 965
965 966 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
966 967
967 968 if ymin == None: ymin = numpy.nanmin(y)
968 969 if ymax == None: ymax = numpy.nanmax(y)
969 970 if zmin == None: zmin = 0.
970 971 if zmax == None: zmax = 1.
971 972
972 973 self.FTP_WEI = ftp_wei
973 974 self.EXP_CODE = exp_code
974 975 self.SUB_EXP_CODE = sub_exp_code
975 976 self.PLOT_POS = plot_pos
976 977
977 978 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
978 979
979 980 self.isConfig = True
980 981 update_figfile = True
981 982
982 983 self.setWinTitle(title)
983 984
984 985 for i in range(self.nplots):
985 986
986 987 pair = dataOut.pairsList[pairsIndexList[i]]
987 988
988 989 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
989 990 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
990 991 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
991 992
992 993
993 994 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
994 995 coherence = numpy.abs(avgcoherenceComplex)
995 996
996 997 z = coherence.reshape((1,-1))
997 998
998 999 counter = 0
999 1000
1000 1001 title = "Coherence Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1001 1002 axes = self.axesList[i*self.__nsubplots*2]
1002 1003 axes.pcolorbuffer(x, y, z,
1003 1004 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1004 1005 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1005 1006 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
1006 1007
1007 1008 if self.__showprofile:
1008 1009 counter += 1
1009 1010 axes = self.axesList[i*self.__nsubplots*2 + counter]
1010 1011 axes.pline(coherence, y,
1011 1012 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
1012 1013 xlabel='', ylabel='', title='', ticksize=7,
1013 1014 ytick_visible=False, nxticks=5,
1014 1015 grid='x')
1015 1016
1016 1017 counter += 1
1017 1018
1018 1019 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1019 1020
1020 1021 z = phase.reshape((1,-1))
1021 1022
1022 1023 title = "Phase Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1023 1024 axes = self.axesList[i*self.__nsubplots*2 + counter]
1024 1025 axes.pcolorbuffer(x, y, z,
1025 1026 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
1026 1027 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1027 1028 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
1028 1029
1029 1030 if self.__showprofile:
1030 1031 counter += 1
1031 1032 axes = self.axesList[i*self.__nsubplots*2 + counter]
1032 1033 axes.pline(phase, y,
1033 1034 xmin=phase_min, xmax=phase_max, ymin=ymin, ymax=ymax,
1034 1035 xlabel='', ylabel='', title='', ticksize=7,
1035 1036 ytick_visible=False, nxticks=4,
1036 1037 grid='x')
1037 1038
1038 1039 self.draw()
1039 1040
1040 1041 if dataOut.ltctime >= self.xmax:
1041 1042 self.counter_imagwr = wr_period
1042 1043 self.isConfig = False
1043 1044 update_figfile = True
1044 1045
1045 1046 self.save(figpath=figpath,
1046 1047 figfile=figfile,
1047 1048 save=save,
1048 1049 ftp=ftp,
1049 1050 wr_period=wr_period,
1050 1051 thisDatetime=thisDatetime,
1051 1052 update_figfile=update_figfile)
1052 1053
1053 1054 class PowerProfilePlot(Figure):
1054 1055
1055 1056 isConfig = None
1056 1057 __nsubplots = None
1057 1058
1058 1059 WIDTHPROF = None
1059 1060 HEIGHTPROF = None
1060 1061 PREFIX = 'spcprofile'
1061 1062
1062 1063 def __init__(self, **kwargs):
1063 1064 Figure.__init__(self, **kwargs)
1064 1065 self.isConfig = False
1065 1066 self.__nsubplots = 1
1066 1067
1067 1068 self.PLOT_CODE = POWER_CODE
1068 1069
1069 1070 self.WIDTH = 300
1070 1071 self.HEIGHT = 500
1071 1072 self.counter_imagwr = 0
1072 1073
1073 1074 def getSubplots(self):
1074 1075 ncol = 1
1075 1076 nrow = 1
1076 1077
1077 1078 return nrow, ncol
1078 1079
1079 1080 def setup(self, id, nplots, wintitle, show):
1080 1081
1081 1082 self.nplots = nplots
1082 1083
1083 1084 ncolspan = 1
1084 1085 colspan = 1
1085 1086
1086 1087 self.createFigure(id = id,
1087 1088 wintitle = wintitle,
1088 1089 widthplot = self.WIDTH,
1089 1090 heightplot = self.HEIGHT,
1090 1091 show=show)
1091 1092
1092 1093 nrow, ncol = self.getSubplots()
1093 1094
1094 1095 counter = 0
1095 1096 for y in range(nrow):
1096 1097 for x in range(ncol):
1097 1098 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1098 1099
1099 1100 def run(self, dataOut, id, wintitle="", channelList=None,
1100 1101 xmin=None, xmax=None, ymin=None, ymax=None,
1101 1102 save=False, figpath='./', figfile=None, show=True,
1102 1103 ftp=False, wr_period=1, server=None,
1103 1104 folder=None, username=None, password=None):
1104 1105
1105 1106
1106 1107 if channelList == None:
1107 1108 channelIndexList = dataOut.channelIndexList
1108 1109 channelList = dataOut.channelList
1109 1110 else:
1110 1111 channelIndexList = []
1111 1112 for channel in channelList:
1112 1113 if channel not in dataOut.channelList:
1113 1114 raise ValueError, "Channel %d is not in dataOut.channelList"
1114 1115 channelIndexList.append(dataOut.channelList.index(channel))
1115 1116
1116 1117 factor = dataOut.normFactor
1117 1118
1118 1119 y = dataOut.getHeiRange()
1119 1120
1120 1121 #for voltage
1121 1122 if dataOut.type == 'Voltage':
1122 1123 x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
1123 1124 x = x.real
1124 1125 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
1125 1126
1126 1127 #for spectra
1127 1128 if dataOut.type == 'Spectra':
1128 1129 x = dataOut.data_spc[channelIndexList,:,:]/factor
1129 1130 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
1130 1131 x = numpy.average(x, axis=1)
1131 1132
1132 1133
1133 1134 xdB = 10*numpy.log10(x)
1134 1135
1135 1136 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1136 1137 title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
1137 1138 xlabel = "dB"
1138 1139 ylabel = "Range (Km)"
1139 1140
1140 1141 if not self.isConfig:
1141 1142
1142 1143 nplots = 1
1143 1144
1144 1145 self.setup(id=id,
1145 1146 nplots=nplots,
1146 1147 wintitle=wintitle,
1147 1148 show=show)
1148 1149
1149 1150 if ymin == None: ymin = numpy.nanmin(y)
1150 1151 if ymax == None: ymax = numpy.nanmax(y)
1151 1152 if xmin == None: xmin = numpy.nanmin(xdB)*0.9
1152 1153 if xmax == None: xmax = numpy.nanmax(xdB)*1.1
1153 1154
1154 1155 self.isConfig = True
1155 1156
1156 1157 self.setWinTitle(title)
1157 1158
1158 1159 title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1159 1160 axes = self.axesList[0]
1160 1161
1161 1162 legendlabels = ["channel %d"%x for x in channelList]
1162 1163 axes.pmultiline(xdB, y,
1163 1164 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1164 1165 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1165 1166 ytick_visible=True, nxticks=5,
1166 1167 grid='x')
1167 1168
1168 1169 self.draw()
1169 1170
1170 1171 self.save(figpath=figpath,
1171 1172 figfile=figfile,
1172 1173 save=save,
1173 1174 ftp=ftp,
1174 1175 wr_period=wr_period,
1175 1176 thisDatetime=thisDatetime)
1176 1177
1177 1178 class SpectraCutPlot(Figure):
1178 1179
1179 1180 isConfig = None
1180 1181 __nsubplots = None
1181 1182
1182 1183 WIDTHPROF = None
1183 1184 HEIGHTPROF = None
1184 1185 PREFIX = 'spc_cut'
1185 1186
1186 1187 def __init__(self, **kwargs):
1187 1188 Figure.__init__(self, **kwargs)
1188 1189 self.isConfig = False
1189 1190 self.__nsubplots = 1
1190 1191
1191 1192 self.PLOT_CODE = POWER_CODE
1192 1193
1193 1194 self.WIDTH = 700
1194 1195 self.HEIGHT = 500
1195 1196 self.counter_imagwr = 0
1196 1197
1197 1198 def getSubplots(self):
1198 1199 ncol = 1
1199 1200 nrow = 1
1200 1201
1201 1202 return nrow, ncol
1202 1203
1203 1204 def setup(self, id, nplots, wintitle, show):
1204 1205
1205 1206 self.nplots = nplots
1206 1207
1207 1208 ncolspan = 1
1208 1209 colspan = 1
1209 1210
1210 1211 self.createFigure(id = id,
1211 1212 wintitle = wintitle,
1212 1213 widthplot = self.WIDTH,
1213 1214 heightplot = self.HEIGHT,
1214 1215 show=show)
1215 1216
1216 1217 nrow, ncol = self.getSubplots()
1217 1218
1218 1219 counter = 0
1219 1220 for y in range(nrow):
1220 1221 for x in range(ncol):
1221 1222 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1222 1223
1223 1224 def run(self, dataOut, id, wintitle="", channelList=None,
1224 1225 xmin=None, xmax=None, ymin=None, ymax=None,
1225 1226 save=False, figpath='./', figfile=None, show=True,
1226 1227 ftp=False, wr_period=1, server=None,
1227 1228 folder=None, username=None, password=None,
1228 1229 xaxis="frequency"):
1229 1230
1230 1231
1231 1232 if channelList == None:
1232 1233 channelIndexList = dataOut.channelIndexList
1233 1234 channelList = dataOut.channelList
1234 1235 else:
1235 1236 channelIndexList = []
1236 1237 for channel in channelList:
1237 1238 if channel not in dataOut.channelList:
1238 1239 raise ValueError, "Channel %d is not in dataOut.channelList"
1239 1240 channelIndexList.append(dataOut.channelList.index(channel))
1240 1241
1241 1242 factor = dataOut.normFactor
1242 1243
1243 1244 y = dataOut.getHeiRange()
1244 1245
1245 1246 z = dataOut.data_spc/factor
1246 1247 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1247 1248
1248 1249 hei_index = numpy.arange(25)*3 + 20
1249 1250
1250 1251 if xaxis == "frequency":
1251 1252 x = dataOut.getFreqRange()/1000.
1252 1253 zdB = 10*numpy.log10(z[0,:,hei_index])
1253 1254 xlabel = "Frequency (kHz)"
1254 1255 ylabel = "Power (dB)"
1255 1256
1256 1257 elif xaxis == "time":
1257 1258 x = dataOut.getAcfRange()
1258 1259 zdB = z[0,:,hei_index]
1259 1260 xlabel = "Time (ms)"
1260 1261 ylabel = "ACF"
1261 1262
1262 1263 else:
1263 1264 x = dataOut.getVelRange()
1264 1265 zdB = 10*numpy.log10(z[0,:,hei_index])
1265 1266 xlabel = "Velocity (m/s)"
1266 1267 ylabel = "Power (dB)"
1267 1268
1268 1269 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1269 1270 title = wintitle + " Range Cuts %s" %(thisDatetime.strftime("%d-%b-%Y"))
1270 1271
1271 1272 if not self.isConfig:
1272 1273
1273 1274 nplots = 1
1274 1275
1275 1276 self.setup(id=id,
1276 1277 nplots=nplots,
1277 1278 wintitle=wintitle,
1278 1279 show=show)
1279 1280
1280 1281 if xmin == None: xmin = numpy.nanmin(x)*0.9
1281 1282 if xmax == None: xmax = numpy.nanmax(x)*1.1
1282 1283 if ymin == None: ymin = numpy.nanmin(zdB)
1283 1284 if ymax == None: ymax = numpy.nanmax(zdB)
1284 1285
1285 1286 self.isConfig = True
1286 1287
1287 1288 self.setWinTitle(title)
1288 1289
1289 1290 title = "Spectra Cuts: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1290 1291 axes = self.axesList[0]
1291 1292
1292 1293 legendlabels = ["Range = %dKm" %y[i] for i in hei_index]
1293 1294
1294 1295 axes.pmultilineyaxis( x, zdB,
1295 1296 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1296 1297 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1297 1298 ytick_visible=True, nxticks=5,
1298 1299 grid='x')
1299 1300
1300 1301 self.draw()
1301 1302
1302 1303 self.save(figpath=figpath,
1303 1304 figfile=figfile,
1304 1305 save=save,
1305 1306 ftp=ftp,
1306 1307 wr_period=wr_period,
1307 1308 thisDatetime=thisDatetime)
1308 1309
1309 1310 class Noise(Figure):
1310 1311
1311 1312 isConfig = None
1312 1313 __nsubplots = None
1313 1314
1314 1315 PREFIX = 'noise'
1315 1316
1316 1317
1317 1318 def __init__(self, **kwargs):
1318 1319 Figure.__init__(self, **kwargs)
1319 1320 self.timerange = 24*60*60
1320 1321 self.isConfig = False
1321 1322 self.__nsubplots = 1
1322 1323 self.counter_imagwr = 0
1323 1324 self.WIDTH = 800
1324 1325 self.HEIGHT = 400
1325 1326 self.WIDTHPROF = 120
1326 1327 self.HEIGHTPROF = 0
1327 1328 self.xdata = None
1328 1329 self.ydata = None
1329 1330
1330 1331 self.PLOT_CODE = NOISE_CODE
1331 1332
1332 1333 self.FTP_WEI = None
1333 1334 self.EXP_CODE = None
1334 1335 self.SUB_EXP_CODE = None
1335 1336 self.PLOT_POS = None
1336 1337 self.figfile = None
1337 1338
1338 1339 self.xmin = None
1339 1340 self.xmax = None
1340 1341
1341 1342 def getSubplots(self):
1342 1343
1343 1344 ncol = 1
1344 1345 nrow = 1
1345 1346
1346 1347 return nrow, ncol
1347 1348
1348 1349 def openfile(self, filename):
1349 1350 dirname = os.path.dirname(filename)
1350 1351
1351 1352 if not os.path.exists(dirname):
1352 1353 os.mkdir(dirname)
1353 1354
1354 1355 f = open(filename,'w+')
1355 1356 f.write('\n\n')
1356 1357 f.write('JICAMARCA RADIO OBSERVATORY - Noise \n')
1357 1358 f.write('DD MM YYYY HH MM SS Channel0 Channel1 Channel2 Channel3\n\n' )
1358 1359 f.close()
1359 1360
1360 1361 def save_data(self, filename_phase, data, data_datetime):
1361 1362
1362 1363 f=open(filename_phase,'a')
1363 1364
1364 1365 timetuple_data = data_datetime.timetuple()
1365 1366 day = str(timetuple_data.tm_mday)
1366 1367 month = str(timetuple_data.tm_mon)
1367 1368 year = str(timetuple_data.tm_year)
1368 1369 hour = str(timetuple_data.tm_hour)
1369 1370 minute = str(timetuple_data.tm_min)
1370 1371 second = str(timetuple_data.tm_sec)
1371 1372
1372 1373 data_msg = ''
1373 1374 for i in range(len(data)):
1374 1375 data_msg += str(data[i]) + ' '
1375 1376
1376 1377 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' ' + data_msg + '\n')
1377 1378 f.close()
1378 1379
1379 1380
1380 1381 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1381 1382
1382 1383 self.__showprofile = showprofile
1383 1384 self.nplots = nplots
1384 1385
1385 1386 ncolspan = 7
1386 1387 colspan = 6
1387 1388 self.__nsubplots = 2
1388 1389
1389 1390 self.createFigure(id = id,
1390 1391 wintitle = wintitle,
1391 1392 widthplot = self.WIDTH+self.WIDTHPROF,
1392 1393 heightplot = self.HEIGHT+self.HEIGHTPROF,
1393 1394 show=show)
1394 1395
1395 1396 nrow, ncol = self.getSubplots()
1396 1397
1397 1398 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1398 1399
1399 1400
1400 1401 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1401 1402 xmin=None, xmax=None, ymin=None, ymax=None,
1402 1403 timerange=None,
1403 1404 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1404 1405 server=None, folder=None, username=None, password=None,
1405 1406 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1406 1407
1407 1408 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1408 1409 return
1409 1410
1410 1411 if channelList == None:
1411 1412 channelIndexList = dataOut.channelIndexList
1412 1413 channelList = dataOut.channelList
1413 1414 else:
1414 1415 channelIndexList = []
1415 1416 for channel in channelList:
1416 1417 if channel not in dataOut.channelList:
1417 1418 raise ValueError, "Channel %d is not in dataOut.channelList"
1418 1419 channelIndexList.append(dataOut.channelList.index(channel))
1419 1420
1420 1421 x = dataOut.getTimeRange()
1421 1422 #y = dataOut.getHeiRange()
1422 1423 factor = dataOut.normFactor
1423 1424 noise = dataOut.noise[channelIndexList]/factor
1424 1425 noisedB = 10*numpy.log10(noise)
1425 1426
1426 1427 thisDatetime = dataOut.datatime
1427 1428
1428 1429 title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1429 1430 xlabel = ""
1430 1431 ylabel = "Intensity (dB)"
1431 1432 update_figfile = False
1432 1433
1433 1434 if not self.isConfig:
1434 1435
1435 1436 nplots = 1
1436 1437
1437 1438 self.setup(id=id,
1438 1439 nplots=nplots,
1439 1440 wintitle=wintitle,
1440 1441 showprofile=showprofile,
1441 1442 show=show)
1442 1443
1443 1444 if timerange != None:
1444 1445 self.timerange = timerange
1445 1446
1446 1447 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1447 1448
1448 1449 if ymin == None: ymin = numpy.floor(numpy.nanmin(noisedB)) - 10.0
1449 1450 if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
1450 1451
1451 1452 self.FTP_WEI = ftp_wei
1452 1453 self.EXP_CODE = exp_code
1453 1454 self.SUB_EXP_CODE = sub_exp_code
1454 1455 self.PLOT_POS = plot_pos
1455 1456
1456 1457
1457 1458 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1458 1459 self.isConfig = True
1459 1460 self.figfile = figfile
1460 1461 self.xdata = numpy.array([])
1461 1462 self.ydata = numpy.array([])
1462 1463
1463 1464 update_figfile = True
1464 1465
1465 1466 #open file beacon phase
1466 1467 path = '%s%03d' %(self.PREFIX, self.id)
1467 1468 noise_file = os.path.join(path,'%s.txt'%self.name)
1468 1469 self.filename_noise = os.path.join(figpath,noise_file)
1469 1470
1470 1471 self.setWinTitle(title)
1471 1472
1472 1473 title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1473 1474
1474 1475 legendlabels = ["channel %d"%(idchannel) for idchannel in channelList]
1475 1476 axes = self.axesList[0]
1476 1477
1477 1478 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1478 1479
1479 1480 if len(self.ydata)==0:
1480 1481 self.ydata = noisedB.reshape(-1,1)
1481 1482 else:
1482 1483 self.ydata = numpy.hstack((self.ydata, noisedB.reshape(-1,1)))
1483 1484
1484 1485
1485 1486 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1486 1487 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1487 1488 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1488 1489 XAxisAsTime=True, grid='both'
1489 1490 )
1490 1491
1491 1492 self.draw()
1492 1493
1493 1494 if dataOut.ltctime >= self.xmax:
1494 1495 self.counter_imagwr = wr_period
1495 1496 self.isConfig = False
1496 1497 update_figfile = True
1497 1498
1498 1499 self.save(figpath=figpath,
1499 1500 figfile=figfile,
1500 1501 save=save,
1501 1502 ftp=ftp,
1502 1503 wr_period=wr_period,
1503 1504 thisDatetime=thisDatetime,
1504 1505 update_figfile=update_figfile)
1505 1506
1506 1507 #store data beacon phase
1507 1508 if save:
1508 1509 self.save_data(self.filename_noise, noisedB, thisDatetime)
1509 1510
1510 1511 class BeaconPhase(Figure):
1511 1512
1512 1513 __isConfig = None
1513 1514 __nsubplots = None
1514 1515
1515 1516 PREFIX = 'beacon_phase'
1516 1517
1517 1518 def __init__(self, **kwargs):
1518 1519 Figure.__init__(self, **kwargs)
1519 1520 self.timerange = 24*60*60
1520 1521 self.isConfig = False
1521 1522 self.__nsubplots = 1
1522 1523 self.counter_imagwr = 0
1523 1524 self.WIDTH = 800
1524 1525 self.HEIGHT = 400
1525 1526 self.WIDTHPROF = 120
1526 1527 self.HEIGHTPROF = 0
1527 1528 self.xdata = None
1528 1529 self.ydata = None
1529 1530
1530 1531 self.PLOT_CODE = BEACON_CODE
1531 1532
1532 1533 self.FTP_WEI = None
1533 1534 self.EXP_CODE = None
1534 1535 self.SUB_EXP_CODE = None
1535 1536 self.PLOT_POS = None
1536 1537
1537 1538 self.filename_phase = None
1538 1539
1539 1540 self.figfile = None
1540 1541
1541 1542 self.xmin = None
1542 1543 self.xmax = None
1543 1544
1544 1545 def getSubplots(self):
1545 1546
1546 1547 ncol = 1
1547 1548 nrow = 1
1548 1549
1549 1550 return nrow, ncol
1550 1551
1551 1552 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1552 1553
1553 1554 self.__showprofile = showprofile
1554 1555 self.nplots = nplots
1555 1556
1556 1557 ncolspan = 7
1557 1558 colspan = 6
1558 1559 self.__nsubplots = 2
1559 1560
1560 1561 self.createFigure(id = id,
1561 1562 wintitle = wintitle,
1562 1563 widthplot = self.WIDTH+self.WIDTHPROF,
1563 1564 heightplot = self.HEIGHT+self.HEIGHTPROF,
1564 1565 show=show)
1565 1566
1566 1567 nrow, ncol = self.getSubplots()
1567 1568
1568 1569 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1569 1570
1570 1571 def save_phase(self, filename_phase):
1571 1572 f = open(filename_phase,'w+')
1572 1573 f.write('\n\n')
1573 1574 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1574 1575 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1575 1576 f.close()
1576 1577
1577 1578 def save_data(self, filename_phase, data, data_datetime):
1578 1579 f=open(filename_phase,'a')
1579 1580 timetuple_data = data_datetime.timetuple()
1580 1581 day = str(timetuple_data.tm_mday)
1581 1582 month = str(timetuple_data.tm_mon)
1582 1583 year = str(timetuple_data.tm_year)
1583 1584 hour = str(timetuple_data.tm_hour)
1584 1585 minute = str(timetuple_data.tm_min)
1585 1586 second = str(timetuple_data.tm_sec)
1586 1587 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1587 1588 f.close()
1588 1589
1589 1590
1590 1591 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1591 1592 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1592 1593 timerange=None,
1593 1594 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1594 1595 server=None, folder=None, username=None, password=None,
1595 1596 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1596 1597
1597 1598 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1598 1599 return
1599 1600
1600 1601 if pairsList == None:
1601 1602 pairsIndexList = dataOut.pairsIndexList[:10]
1602 1603 else:
1603 1604 pairsIndexList = []
1604 1605 for pair in pairsList:
1605 1606 if pair not in dataOut.pairsList:
1606 1607 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
1607 1608 pairsIndexList.append(dataOut.pairsList.index(pair))
1608 1609
1609 1610 if pairsIndexList == []:
1610 1611 return
1611 1612
1612 1613 # if len(pairsIndexList) > 4:
1613 1614 # pairsIndexList = pairsIndexList[0:4]
1614 1615
1615 1616 hmin_index = None
1616 1617 hmax_index = None
1617 1618
1618 1619 if hmin != None and hmax != None:
1619 1620 indexes = numpy.arange(dataOut.nHeights)
1620 1621 hmin_list = indexes[dataOut.heightList >= hmin]
1621 1622 hmax_list = indexes[dataOut.heightList <= hmax]
1622 1623
1623 1624 if hmin_list.any():
1624 1625 hmin_index = hmin_list[0]
1625 1626
1626 1627 if hmax_list.any():
1627 1628 hmax_index = hmax_list[-1]+1
1628 1629
1629 1630 x = dataOut.getTimeRange()
1630 1631 #y = dataOut.getHeiRange()
1631 1632
1632 1633
1633 1634 thisDatetime = dataOut.datatime
1634 1635
1635 1636 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1636 1637 xlabel = "Local Time"
1637 1638 ylabel = "Phase (degrees)"
1638 1639
1639 1640 update_figfile = False
1640 1641
1641 1642 nplots = len(pairsIndexList)
1642 1643 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1643 1644 phase_beacon = numpy.zeros(len(pairsIndexList))
1644 1645 for i in range(nplots):
1645 1646 pair = dataOut.pairsList[pairsIndexList[i]]
1646 1647 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1647 1648 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1648 1649 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1649 1650 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1650 1651 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1651 1652
1652 1653 #print "Phase %d%d" %(pair[0], pair[1])
1653 1654 #print phase[dataOut.beacon_heiIndexList]
1654 1655
1655 1656 if dataOut.beacon_heiIndexList:
1656 1657 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1657 1658 else:
1658 1659 phase_beacon[i] = numpy.average(phase)
1659 1660
1660 1661 if not self.isConfig:
1661 1662
1662 1663 nplots = len(pairsIndexList)
1663 1664
1664 1665 self.setup(id=id,
1665 1666 nplots=nplots,
1666 1667 wintitle=wintitle,
1667 1668 showprofile=showprofile,
1668 1669 show=show)
1669 1670
1670 1671 if timerange != None:
1671 1672 self.timerange = timerange
1672 1673
1673 1674 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1674 1675
1675 1676 if ymin == None: ymin = 0
1676 1677 if ymax == None: ymax = 360
1677 1678
1678 1679 self.FTP_WEI = ftp_wei
1679 1680 self.EXP_CODE = exp_code
1680 1681 self.SUB_EXP_CODE = sub_exp_code
1681 1682 self.PLOT_POS = plot_pos
1682 1683
1683 1684 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1684 1685 self.isConfig = True
1685 1686 self.figfile = figfile
1686 1687 self.xdata = numpy.array([])
1687 1688 self.ydata = numpy.array([])
1688 1689
1689 1690 update_figfile = True
1690 1691
1691 1692 #open file beacon phase
1692 1693 path = '%s%03d' %(self.PREFIX, self.id)
1693 1694 beacon_file = os.path.join(path,'%s.txt'%self.name)
1694 1695 self.filename_phase = os.path.join(figpath,beacon_file)
1695 1696 #self.save_phase(self.filename_phase)
1696 1697
1697 1698
1698 1699 #store data beacon phase
1699 1700 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1700 1701
1701 1702 self.setWinTitle(title)
1702 1703
1703 1704
1704 1705 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1705 1706
1706 1707 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1707 1708
1708 1709 axes = self.axesList[0]
1709 1710
1710 1711 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1711 1712
1712 1713 if len(self.ydata)==0:
1713 1714 self.ydata = phase_beacon.reshape(-1,1)
1714 1715 else:
1715 1716 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1716 1717
1717 1718
1718 1719 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1719 1720 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1720 1721 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1721 1722 XAxisAsTime=True, grid='both'
1722 1723 )
1723 1724
1724 1725 self.draw()
1725 1726
1726 1727 if dataOut.ltctime >= self.xmax:
1727 1728 self.counter_imagwr = wr_period
1728 1729 self.isConfig = False
1729 1730 update_figfile = True
1730 1731
1731 1732 self.save(figpath=figpath,
1732 1733 figfile=figfile,
1733 1734 save=save,
1734 1735 ftp=ftp,
1735 1736 wr_period=wr_period,
1736 1737 thisDatetime=thisDatetime,
1737 1738 update_figfile=update_figfile)
Show file before | Show file after | Hide comments
@@ -1,797 +1,760
1 1 import numpy
2 2
3 3 from jroproc_base import ProcessingUnit, Operation
4 4 from schainpy.model.data.jrodata import Spectra
5 5 from schainpy.model.data.jrodata import hildebrand_sekhon
6 6
7 7 class SpectraAFCProc(ProcessingUnit):
8 8
9 9 def __init__(self, **kwargs):
10 10
11 11 ProcessingUnit.__init__(self, **kwargs)
12 12
13 13 self.buffer = None
14 14 self.firstdatatime = None
15 15 self.profIndex = 0
16 16 self.dataOut = Spectra()
17 17 self.id_min = None
18 18 self.id_max = None
19 19
20 20 def __updateSpecFromVoltage(self):
21 21
22 22 self.dataOut.plotting = "spectra_acf"
23 23
24 24 self.dataOut.timeZone = self.dataIn.timeZone
25 25 self.dataOut.dstFlag = self.dataIn.dstFlag
26 26 self.dataOut.errorCount = self.dataIn.errorCount
27 27 self.dataOut.useLocalTime = self.dataIn.useLocalTime
28 28
29 29 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
30 30 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
31 31 self.dataOut.ippSeconds = self.dataIn.getDeltaH()*(10**-6)/0.15
32 32
33 33 self.dataOut.channelList = self.dataIn.channelList
34 34 self.dataOut.heightList = self.dataIn.heightList
35 35 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
36 36
37 37 self.dataOut.nBaud = self.dataIn.nBaud
38 38 self.dataOut.nCode = self.dataIn.nCode
39 39 self.dataOut.code = self.dataIn.code
40 40 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
41 41
42 42 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
43 43 self.dataOut.utctime = self.firstdatatime
44 44 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
45 45 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
46 46 self.dataOut.flagShiftFFT = False
47 47
48 48 self.dataOut.nCohInt = self.dataIn.nCohInt
49 49 self.dataOut.nIncohInt = 1
50 50
51 51 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
52 52
53 53 self.dataOut.frequency = self.dataIn.frequency
54 54 self.dataOut.realtime = self.dataIn.realtime
55 55
56 56 self.dataOut.azimuth = self.dataIn.azimuth
57 57 self.dataOut.zenith = self.dataIn.zenith
58 58
59 59 self.dataOut.beam.codeList = self.dataIn.beam.codeList
60 60 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
61 61 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
62 62
63 63 def __decodeData(self, nProfiles, code):
64 64
65 65 if code is None:
66 66 return
67 67
68 68 for i in range(nProfiles):
69 69 self.buffer[:,i,:] = self.buffer[:,i,:]*code[0][i]
70 70
71 71 def __getFft(self):
72 72 """
73 73 Convierte valores de Voltaje a Spectra
74 74
75 75 Affected:
76 76 self.dataOut.data_spc
77 77 self.dataOut.data_cspc
78 78 self.dataOut.data_dc
79 79 self.dataOut.heightList
80 80 self.profIndex
81 81 self.buffer
82 82 self.dataOut.flagNoData
83 83 """
84 84 nsegments = self.dataOut.nHeights
85 85
86 86 _fft_buffer = numpy.zeros((self.dataOut.nChannels, self.dataOut.nProfiles, nsegments), dtype='complex')
87 87
88 88 for i in range(nsegments):
89 89 try:
90 90 _fft_buffer[:,:,i] = self.buffer[:,i:i+self.dataOut.nProfiles]
91 91
92 92 if self.code is not None:
93 93 _fft_buffer[:,:,i] = _fft_buffer[:,:,i]*self.code[0]
94 94 except:
95 95 pass
96 96
97 97 fft_volt = numpy.fft.fft(_fft_buffer, n=self.dataOut.nFFTPoints, axis=1)
98 98 fft_volt = fft_volt.astype(numpy.dtype('complex'))
99 99 dc = fft_volt[:,0,:]
100 100
101 101 #calculo de self-spectra
102 102 spc = fft_volt * numpy.conjugate(fft_volt)
103 103 data = numpy.fft.ifft(spc, axis=1)
104 104 data = numpy.fft.fftshift(data, axes=(1,))
105 105
106 106 spc = data
107 107
108 108 blocksize = 0
109 109 blocksize += dc.size
110 110 blocksize += spc.size
111 111
112 112 cspc = None
113 113 pairIndex = 0
114 114
115 115 if self.dataOut.pairsList != None:
116 116 #calculo de cross-spectra
117 117 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
118 118 for pair in self.dataOut.pairsList:
119 119 if pair[0] not in self.dataOut.channelList:
120 120 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
121 121 if pair[1] not in self.dataOut.channelList:
122 122 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
123 123
124 124 chan_index0 = self.dataOut.channelList.index(pair[0])
125 125 chan_index1 = self.dataOut.channelList.index(pair[1])
126 126
127 127 cspc[pairIndex,:,:] = fft_volt[chan_index0,:,:] * numpy.conjugate(fft_volt[chan_index1,:,:])
128 128 pairIndex += 1
129 129 blocksize += cspc.size
130 130
131 131 self.dataOut.data_spc = spc
132 132 self.dataOut.data_cspc = cspc
133 133 self.dataOut.data_dc = dc
134 134 self.dataOut.blockSize = blocksize
135 135 self.dataOut.flagShiftFFT = True
136 136
137 137 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], code=None, nCode=1, nBaud=1,real= None, imag=None):
138 138
139 139 self.dataOut.flagNoData = True
140 140
141 141 if self.dataIn.type == "Spectra":
142 142 self.dataOut.copy(self.dataIn)
143 143 spc = self.dataOut.data_spc
144 data = numpy.fft.ifft(spc, axis=1)
145 data = numpy.fft.fftshift(data, axes=(1,))
144 data = numpy.fft.fftshift( spc, axes=(1,))
145 data = numpy.fft.ifft(data, axis=1)
146 data = numpy.fft.fftshift( data, axes=(1,))
146 147 acf = numpy.abs(data) # Autocorrelacion LLAMAR A ESTE VALOR ACF
147 148 if real:
148 149 acf = data.real
149 150 if imag:
150 151 acf = data.imag
151 shape = acf.shape # nchannels, nprofiles, nsamples
152
153 #import matplotlib.pyplot as plt
154 #print "test",acf.shape
155 #acf_tmp=acf[0,:,]
156 #plt.plot(acf_tmp)
157 #plt.show()
158 #import time
159 #time.sleep(10)
160
152 shape = acf.shape # nchannels, nprofiles, nsamples //nchannles, lags , alturas
153
161 154 for j in range(shape[0]):
162 for i in range(shape[1]):
163 tmp = numpy.argmax(acf[j,:,i])
164 if i>30:
165 value = (acf[j,:,i][tmp+3]+acf[j,:,i][tmp+4])/2.0
166 acf[j,:,i][tmp] = value
167 acf[j,:,i][tmp-1] = value
168 acf[j,:,i][tmp+1] = value
169 acf[j,:,i][tmp-2] = value
170 acf[j,:,i][tmp+2] = value
171
172 import scipy as sp
173 from scipy import signal
174 #acf[3,:,i] = sp.signal.medfilt(acf[3,:,i],21)
175
176
177
178 #print numpy.argmax(acf[0,:,85])
179 #import matplotlib.pyplot as plt
180 #plt.plot(acf[0,:,85])
181 #a= acf[0,:,85]
182 #b= acf[0,:,0]
183 #print a[200],a[198],a[199], a[201],a[202],a[203]
184 #plt.show()
185 #import time
186 #time.sleep(10)
187
188
155 for i in range(shape[2]):
156 tmp = shape[1]/2
157 value = (acf[j,:,i][tmp-1]+acf[j,:,i][tmp+1])/2.0
158 acf[j,:,i][tmp] = value
189 159 # Normalizando
190 160 for i in range(shape[0]):
191 161 for j in range(shape[2]):
192 162 acf[i,:,j]= acf[i,:,j] / numpy.max(numpy.abs(acf[i,:,j]))
193 163
194 #import matplotlib.pyplot as plt
195 #plt.plot(acf[0,:,85])
196 #plt.show()
197 #import time
198 #time.sleep(20)
199
200
201 164 self.dataOut.data_acf = acf
202 165 return True
203 166
204 167 if code is not None:
205 168 self.code = numpy.array(code).reshape(nCode,nBaud)
206 169 else:
207 170 self.code = None
208 171
209 172 if self.dataIn.type == "Voltage":
210 173
211 174 if nFFTPoints == None:
212 175 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
213 176
214 177 if nProfiles == None:
215 178 nProfiles = nFFTPoints
216 179
217 180 self.dataOut.ippFactor = 1
218 181
219 182 self.dataOut.nFFTPoints = nFFTPoints
220 183 self.dataOut.nProfiles = nProfiles
221 184 self.dataOut.pairsList = pairsList
222 185
223 186 # if self.buffer is None:
224 187 # self.buffer = numpy.zeros( (self.dataIn.nChannels, nProfiles, self.dataIn.nHeights),
225 188 # dtype='complex')
226 189
227 190 if not self.dataIn.flagDataAsBlock:
228 191 self.buffer = self.dataIn.data.copy()
229 192
230 193 # for i in range(self.dataIn.nHeights):
231 194 # self.buffer[:, self.profIndex, self.profIndex:] = voltage_data[:,:self.dataIn.nHeights - self.profIndex]
232 195 #
233 196 # self.profIndex += 1
234 197
235 198 else:
236 199 raise ValueError, ""
237 200
238 201 self.firstdatatime = self.dataIn.utctime
239 202
240 203 self.profIndex == nProfiles
241 204
242 205 self.__updateSpecFromVoltage()
243 206
244 207 self.__getFft()
245 208
246 209 self.dataOut.flagNoData = False
247 210
248 211 return True
249 212
250 213 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
251 214
252 215 def __selectPairs(self, pairsList):
253 216
254 217 if channelList == None:
255 218 return
256 219
257 220 pairsIndexListSelected = []
258 221
259 222 for thisPair in pairsList:
260 223
261 224 if thisPair not in self.dataOut.pairsList:
262 225 continue
263 226
264 227 pairIndex = self.dataOut.pairsList.index(thisPair)
265 228
266 229 pairsIndexListSelected.append(pairIndex)
267 230
268 231 if not pairsIndexListSelected:
269 232 self.dataOut.data_cspc = None
270 233 self.dataOut.pairsList = []
271 234 return
272 235
273 236 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
274 237 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
275 238
276 239 return
277 240
278 241 def __selectPairsByChannel(self, channelList=None):
279 242
280 243 if channelList == None:
281 244 return
282 245
283 246 pairsIndexListSelected = []
284 247 for pairIndex in self.dataOut.pairsIndexList:
285 248 #First pair
286 249 if self.dataOut.pairsList[pairIndex][0] not in channelList:
287 250 continue
288 251 #Second pair
289 252 if self.dataOut.pairsList[pairIndex][1] not in channelList:
290 253 continue
291 254
292 255 pairsIndexListSelected.append(pairIndex)
293 256
294 257 if not pairsIndexListSelected:
295 258 self.dataOut.data_cspc = None
296 259 self.dataOut.pairsList = []
297 260 return
298 261
299 262 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
300 263 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
301 264
302 265 return
303 266
304 267 def selectChannels(self, channelList):
305 268
306 269 channelIndexList = []
307 270
308 271 for channel in channelList:
309 272 if channel not in self.dataOut.channelList:
310 273 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
311 274
312 275 index = self.dataOut.channelList.index(channel)
313 276 channelIndexList.append(index)
314 277
315 278 self.selectChannelsByIndex(channelIndexList)
316 279
317 280 def selectChannelsByIndex(self, channelIndexList):
318 281 """
319 282 Selecciona un bloque de datos en base a canales segun el channelIndexList
320 283
321 284 Input:
322 285 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
323 286
324 287 Affected:
325 288 self.dataOut.data_spc
326 289 self.dataOut.channelIndexList
327 290 self.dataOut.nChannels
328 291
329 292 Return:
330 293 None
331 294 """
332 295
333 296 for channelIndex in channelIndexList:
334 297 if channelIndex not in self.dataOut.channelIndexList:
335 298 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
336 299
337 300 # nChannels = len(channelIndexList)
338 301
339 302 data_spc = self.dataOut.data_spc[channelIndexList,:]
340 303 data_dc = self.dataOut.data_dc[channelIndexList,:]
341 304
342 305 self.dataOut.data_spc = data_spc
343 306 self.dataOut.data_dc = data_dc
344 307
345 308 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
346 309 # self.dataOut.nChannels = nChannels
347 310
348 311 self.__selectPairsByChannel(self.dataOut.channelList)
349 312
350 313 return 1
351 314
352 315 def selectHeights(self, minHei, maxHei):
353 316 """
354 317 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
355 318 minHei <= height <= maxHei
356 319
357 320 Input:
358 321 minHei : valor minimo de altura a considerar
359 322 maxHei : valor maximo de altura a considerar
360 323
361 324 Affected:
362 325 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
363 326
364 327 Return:
365 328 1 si el metodo se ejecuto con exito caso contrario devuelve 0
366 329 """
367 330
368 331 if (minHei > maxHei):
369 332 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
370 333
371 334 if (minHei < self.dataOut.heightList[0]):
372 335 minHei = self.dataOut.heightList[0]
373 336
374 337 if (maxHei > self.dataOut.heightList[-1]):
375 338 maxHei = self.dataOut.heightList[-1]
376 339
377 340 minIndex = 0
378 341 maxIndex = 0
379 342 heights = self.dataOut.heightList
380 343
381 344 inda = numpy.where(heights >= minHei)
382 345 indb = numpy.where(heights <= maxHei)
383 346
384 347 try:
385 348 minIndex = inda[0][0]
386 349 except:
387 350 minIndex = 0
388 351
389 352 try:
390 353 maxIndex = indb[0][-1]
391 354 except:
392 355 maxIndex = len(heights)
393 356
394 357 self.selectHeightsByIndex(minIndex, maxIndex)
395 358
396 359 return 1
397 360
398 361 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
399 362 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
400 363
401 364 if hei_ref != None:
402 365 newheis = numpy.where(self.dataOut.heightList>hei_ref)
403 366
404 367 minIndex = min(newheis[0])
405 368 maxIndex = max(newheis[0])
406 369 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
407 370 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
408 371
409 372 # determina indices
410 373 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
411 374 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
412 375 beacon_dB = numpy.sort(avg_dB)[-nheis:]
413 376 beacon_heiIndexList = []
414 377 for val in avg_dB.tolist():
415 378 if val >= beacon_dB[0]:
416 379 beacon_heiIndexList.append(avg_dB.tolist().index(val))
417 380
418 381 #data_spc = data_spc[:,:,beacon_heiIndexList]
419 382 data_cspc = None
420 383 if self.dataOut.data_cspc is not None:
421 384 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
422 385 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
423 386
424 387 data_dc = None
425 388 if self.dataOut.data_dc is not None:
426 389 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
427 390 #data_dc = data_dc[:,beacon_heiIndexList]
428 391
429 392 self.dataOut.data_spc = data_spc
430 393 self.dataOut.data_cspc = data_cspc
431 394 self.dataOut.data_dc = data_dc
432 395 self.dataOut.heightList = heightList
433 396 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
434 397
435 398 return 1
436 399
437 400
438 401 def selectHeightsByIndex(self, minIndex, maxIndex):
439 402 """
440 403 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
441 404 minIndex <= index <= maxIndex
442 405
443 406 Input:
444 407 minIndex : valor de indice minimo de altura a considerar
445 408 maxIndex : valor de indice maximo de altura a considerar
446 409
447 410 Affected:
448 411 self.dataOut.data_spc
449 412 self.dataOut.data_cspc
450 413 self.dataOut.data_dc
451 414 self.dataOut.heightList
452 415
453 416 Return:
454 417 1 si el metodo se ejecuto con exito caso contrario devuelve 0
455 418 """
456 419
457 420 if (minIndex < 0) or (minIndex > maxIndex):
458 421 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
459 422
460 423 if (maxIndex >= self.dataOut.nHeights):
461 424 maxIndex = self.dataOut.nHeights-1
462 425
463 426 #Spectra
464 427 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
465 428
466 429 data_cspc = None
467 430 if self.dataOut.data_cspc is not None:
468 431 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
469 432
470 433 data_dc = None
471 434 if self.dataOut.data_dc is not None:
472 435 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
473 436
474 437 self.dataOut.data_spc = data_spc
475 438 self.dataOut.data_cspc = data_cspc
476 439 self.dataOut.data_dc = data_dc
477 440
478 441 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
479 442
480 443 return 1
481 444
482 445 def removeDC(self, mode = 2):
483 446 jspectra = self.dataOut.data_spc
484 447 jcspectra = self.dataOut.data_cspc
485 448
486 449
487 450 num_chan = jspectra.shape[0]
488 451 num_hei = jspectra.shape[2]
489 452
490 453 if jcspectra is not None:
491 454 jcspectraExist = True
492 455 num_pairs = jcspectra.shape[0]
493 456 else: jcspectraExist = False
494 457
495 458 freq_dc = jspectra.shape[1]/2
496 459 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
497 460
498 461 if ind_vel[0]<0:
499 462 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
500 463
501 464 if mode == 1:
502 465 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
503 466
504 467 if jcspectraExist:
505 468 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
506 469
507 470 if mode == 2:
508 471
509 472 vel = numpy.array([-2,-1,1,2])
510 473 xx = numpy.zeros([4,4])
511 474
512 475 for fil in range(4):
513 476 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
514 477
515 478 xx_inv = numpy.linalg.inv(xx)
516 479 xx_aux = xx_inv[0,:]
517 480
518 481 for ich in range(num_chan):
519 482 yy = jspectra[ich,ind_vel,:]
520 483 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
521 484
522 485 junkid = jspectra[ich,freq_dc,:]<=0
523 486 cjunkid = sum(junkid)
524 487
525 488 if cjunkid.any():
526 489 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
527 490
528 491 if jcspectraExist:
529 492 for ip in range(num_pairs):
530 493 yy = jcspectra[ip,ind_vel,:]
531 494 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
532 495
533 496
534 497 self.dataOut.data_spc = jspectra
535 498 self.dataOut.data_cspc = jcspectra
536 499
537 500 return 1
538 501
539 502 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
540 503
541 504 jspectra = self.dataOut.data_spc
542 505 jcspectra = self.dataOut.data_cspc
543 506 jnoise = self.dataOut.getNoise()
544 507 num_incoh = self.dataOut.nIncohInt
545 508
546 509 num_channel = jspectra.shape[0]
547 510 num_prof = jspectra.shape[1]
548 511 num_hei = jspectra.shape[2]
549 512
550 513 #hei_interf
551 514 if hei_interf is None:
552 515 count_hei = num_hei/2 #Como es entero no importa
553 516 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
554 517 hei_interf = numpy.asarray(hei_interf)[0]
555 518 #nhei_interf
556 519 if (nhei_interf == None):
557 520 nhei_interf = 5
558 521 if (nhei_interf < 1):
559 522 nhei_interf = 1
560 523 if (nhei_interf > count_hei):
561 524 nhei_interf = count_hei
562 525 if (offhei_interf == None):
563 526 offhei_interf = 0
564 527
565 528 ind_hei = range(num_hei)
566 529 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
567 530 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
568 531 mask_prof = numpy.asarray(range(num_prof))
569 532 num_mask_prof = mask_prof.size
570 533 comp_mask_prof = [0, num_prof/2]
571 534
572 535
573 536 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
574 537 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
575 538 jnoise = numpy.nan
576 539 noise_exist = jnoise[0] < numpy.Inf
577 540
578 541 #Subrutina de Remocion de la Interferencia
579 542 for ich in range(num_channel):
580 543 #Se ordena los espectros segun su potencia (menor a mayor)
581 544 power = jspectra[ich,mask_prof,:]
582 545 power = power[:,hei_interf]
583 546 power = power.sum(axis = 0)
584 547 psort = power.ravel().argsort()
585 548
586 549 #Se estima la interferencia promedio en los Espectros de Potencia empleando
587 550 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
588 551
589 552 if noise_exist:
590 553 # tmp_noise = jnoise[ich] / num_prof
591 554 tmp_noise = jnoise[ich]
592 555 junkspc_interf = junkspc_interf - tmp_noise
593 556 #junkspc_interf[:,comp_mask_prof] = 0
594 557
595 558 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
596 559 jspc_interf = jspc_interf.transpose()
597 560 #Calculando el espectro de interferencia promedio
598 561 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
599 562 noiseid = noiseid[0]
600 563 cnoiseid = noiseid.size
601 564 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
602 565 interfid = interfid[0]
603 566 cinterfid = interfid.size
604 567
605 568 if (cnoiseid > 0): jspc_interf[noiseid] = 0
606 569
607 570 #Expandiendo los perfiles a limpiar
608 571 if (cinterfid > 0):
609 572 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
610 573 new_interfid = numpy.asarray(new_interfid)
611 574 new_interfid = {x for x in new_interfid}
612 575 new_interfid = numpy.array(list(new_interfid))
613 576 new_cinterfid = new_interfid.size
614 577 else: new_cinterfid = 0
615 578
616 579 for ip in range(new_cinterfid):
617 580 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
618 581 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
619 582
620 583
621 584 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
622 585
623 586 #Removiendo la interferencia del punto de mayor interferencia
624 587 ListAux = jspc_interf[mask_prof].tolist()
625 588 maxid = ListAux.index(max(ListAux))
626 589
627 590
628 591 if cinterfid > 0:
629 592 for ip in range(cinterfid*(interf == 2) - 1):
630 593 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
631 594 cind = len(ind)
632 595
633 596 if (cind > 0):
634 597 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
635 598
636 599 ind = numpy.array([-2,-1,1,2])
637 600 xx = numpy.zeros([4,4])
638 601
639 602 for id1 in range(4):
640 603 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
641 604
642 605 xx_inv = numpy.linalg.inv(xx)
643 606 xx = xx_inv[:,0]
644 607 ind = (ind + maxid + num_mask_prof)%num_mask_prof
645 608 yy = jspectra[ich,mask_prof[ind],:]
646 609 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
647 610
648 611
649 612 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
650 613 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
651 614
652 615 #Remocion de Interferencia en el Cross Spectra
653 616 if jcspectra is None: return jspectra, jcspectra
654 617 num_pairs = jcspectra.size/(num_prof*num_hei)
655 618 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
656 619
657 620 for ip in range(num_pairs):
658 621
659 622 #-------------------------------------------
660 623
661 624 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
662 625 cspower = cspower[:,hei_interf]
663 626 cspower = cspower.sum(axis = 0)
664 627
665 628 cspsort = cspower.ravel().argsort()
666 629 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
667 630 junkcspc_interf = junkcspc_interf.transpose()
668 631 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
669 632
670 633 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
671 634
672 635 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
673 636 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
674 637 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
675 638
676 639 for iprof in range(num_prof):
677 640 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
678 641 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
679 642
680 643 #Removiendo la Interferencia
681 644 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
682 645
683 646 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
684 647 maxid = ListAux.index(max(ListAux))
685 648
686 649 ind = numpy.array([-2,-1,1,2])
687 650 xx = numpy.zeros([4,4])
688 651
689 652 for id1 in range(4):
690 653 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
691 654
692 655 xx_inv = numpy.linalg.inv(xx)
693 656 xx = xx_inv[:,0]
694 657
695 658 ind = (ind + maxid + num_mask_prof)%num_mask_prof
696 659 yy = jcspectra[ip,mask_prof[ind],:]
697 660 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
698 661
699 662 #Guardar Resultados
700 663 self.dataOut.data_spc = jspectra
701 664 self.dataOut.data_cspc = jcspectra
702 665
703 666 return 1
704 667
705 668 def setRadarFrequency(self, frequency=None):
706 669
707 670 if frequency != None:
708 671 self.dataOut.frequency = frequency
709 672
710 673 return 1
711 674
712 675 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
713 676 #validacion de rango
714 677 if minHei == None:
715 678 minHei = self.dataOut.heightList[0]
716 679
717 680 if maxHei == None:
718 681 maxHei = self.dataOut.heightList[-1]
719 682
720 683 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
721 684 print 'minHei: %.2f is out of the heights range'%(minHei)
722 685 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
723 686 minHei = self.dataOut.heightList[0]
724 687
725 688 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
726 689 print 'maxHei: %.2f is out of the heights range'%(maxHei)
727 690 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
728 691 maxHei = self.dataOut.heightList[-1]
729 692
730 693 # validacion de velocidades
731 694 velrange = self.dataOut.getVelRange(1)
732 695
733 696 if minVel == None:
734 697 minVel = velrange[0]
735 698
736 699 if maxVel == None:
737 700 maxVel = velrange[-1]
738 701
739 702 if (minVel < velrange[0]) or (minVel > maxVel):
740 703 print 'minVel: %.2f is out of the velocity range'%(minVel)
741 704 print 'minVel is setting to %.2f'%(velrange[0])
742 705 minVel = velrange[0]
743 706
744 707 if (maxVel > velrange[-1]) or (maxVel < minVel):
745 708 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
746 709 print 'maxVel is setting to %.2f'%(velrange[-1])
747 710 maxVel = velrange[-1]
748 711
749 712 # seleccion de indices para rango
750 713 minIndex = 0
751 714 maxIndex = 0
752 715 heights = self.dataOut.heightList
753 716
754 717 inda = numpy.where(heights >= minHei)
755 718 indb = numpy.where(heights <= maxHei)
756 719
757 720 try:
758 721 minIndex = inda[0][0]
759 722 except:
760 723 minIndex = 0
761 724
762 725 try:
763 726 maxIndex = indb[0][-1]
764 727 except:
765 728 maxIndex = len(heights)
766 729
767 730 if (minIndex < 0) or (minIndex > maxIndex):
768 731 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
769 732
770 733 if (maxIndex >= self.dataOut.nHeights):
771 734 maxIndex = self.dataOut.nHeights-1
772 735
773 736 # seleccion de indices para velocidades
774 737 indminvel = numpy.where(velrange >= minVel)
775 738 indmaxvel = numpy.where(velrange <= maxVel)
776 739 try:
777 740 minIndexVel = indminvel[0][0]
778 741 except:
779 742 minIndexVel = 0
780 743
781 744 try:
782 745 maxIndexVel = indmaxvel[0][-1]
783 746 except:
784 747 maxIndexVel = len(velrange)
785 748
786 749 #seleccion del espectro
787 750 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
788 751 #estimacion de ruido
789 752 noise = numpy.zeros(self.dataOut.nChannels)
790 753
791 754 for channel in range(self.dataOut.nChannels):
792 755 daux = data_spc[channel,:,:]
793 756 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
794 757
795 758 self.dataOut.noise_estimation = noise.copy()
796 759
797 760 return 1
Show file before | Show file after | Hide comments
@@ -1,1540 +1,1540
1 1 import sys
2 2 import numpy
3 3 from scipy import interpolate
4 4 from schainpy import cSchain
5 5 from jroproc_base import ProcessingUnit, Operation
6 6 from schainpy.model.data.jrodata import Voltage
7 7 from time import time
8 8
9 9 import math
10 10
11 11 def rep_seq(x, rep=10):
12 12 L = len(x) * rep
13 13 res = numpy.zeros(L, dtype=x.dtype)
14 14 idx = numpy.arange(len(x)) * rep
15 15 for i in numpy.arange(rep):
16 16 res[idx + i] = x
17 17 return(res)
18 18
19 19
20 20 def create_pseudo_random_code(clen=10000, seed=0):
21 21 """
22 22 seed is a way of reproducing the random code without
23 23 having to store all actual codes. the seed can then
24 24 act as a sort of station_id.
25 25
26 26 """
27 27 numpy.random.seed(seed)
28 28 phases = numpy.array(
29 29 numpy.exp(1.0j * 2.0 * math.pi * numpy.random.random(clen)),
30 30 dtype=numpy.complex64,
31 31 )
32 32 return(phases)
33 33
34 34
35 35 def periodic_convolution_matrix(envelope, rmin=0, rmax=100):
36 36 """
37 37 we imply that the number of measurements is equal to the number of elements
38 38 in code
39 39
40 40 """
41 41 L = len(envelope)
42 42 ridx = numpy.arange(rmin, rmax)
43 43 A = numpy.zeros([L, rmax-rmin], dtype=numpy.complex64)
44 44 for i in numpy.arange(L):
45 45 A[i, :] = envelope[(i-ridx) % L]
46 46 result = {}
47 47 result['A'] = A
48 48 result['ridx'] = ridx
49 49 return(result)
50 50
51 51
52 52 B_cache = 0
53 53 r_cache = 0
54 54 B_cached = False
55 55 def create_estimation_matrix(code, rmin=0, rmax=1000, cache=True):
56 56 global B_cache
57 57 global r_cache
58 58 global B_cached
59 59
60 60 if not cache or not B_cached:
61 61 r_cache = periodic_convolution_matrix(
62 62 envelope=code, rmin=rmin, rmax=rmax,
63 63 )
64 64 A = r_cache['A']
65 65 Ah = numpy.transpose(numpy.conjugate(A))
66 66 B_cache = numpy.dot(numpy.linalg.inv(numpy.dot(Ah, A)), Ah)
67 67 r_cache['B'] = B_cache
68 68 B_cached = True
69 69 return(r_cache)
70 70 else:
71 71 return(r_cache)
72 72
73 73 class VoltageProc(ProcessingUnit):
74 74
75 75
76 76 def __init__(self, **kwargs):
77 77
78 78 ProcessingUnit.__init__(self, **kwargs)
79 79
80 80 # self.objectDict = {}
81 81 self.dataOut = Voltage()
82 82 self.flip = 1
83 83
84 84 def run(self):
85 85 if self.dataIn.type == 'AMISR':
86 86 self.__updateObjFromAmisrInput()
87 87
88 88 if self.dataIn.type == 'Voltage':
89 89 self.dataOut.copy(self.dataIn)
90 90
91 91 # self.dataOut.copy(self.dataIn)
92 92
93 93 def __updateObjFromAmisrInput(self):
94 94
95 95 self.dataOut.timeZone = self.dataIn.timeZone
96 96 self.dataOut.dstFlag = self.dataIn.dstFlag
97 97 self.dataOut.errorCount = self.dataIn.errorCount
98 98 self.dataOut.useLocalTime = self.dataIn.useLocalTime
99 99
100 100 self.dataOut.flagNoData = self.dataIn.flagNoData
101 101 self.dataOut.data = self.dataIn.data
102 102 self.dataOut.utctime = self.dataIn.utctime
103 103 self.dataOut.channelList = self.dataIn.channelList
104 104 # self.dataOut.timeInterval = self.dataIn.timeInterval
105 105 self.dataOut.heightList = self.dataIn.heightList
106 106 self.dataOut.nProfiles = self.dataIn.nProfiles
107 107
108 108 self.dataOut.nCohInt = self.dataIn.nCohInt
109 109 self.dataOut.ippSeconds = self.dataIn.ippSeconds
110 110 self.dataOut.frequency = self.dataIn.frequency
111 111
112 112 self.dataOut.azimuth = self.dataIn.azimuth
113 113 self.dataOut.zenith = self.dataIn.zenith
114 114
115 115 self.dataOut.beam.codeList = self.dataIn.beam.codeList
116 116 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
117 117 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
118 118 #
119 119 # pass#
120 120 #
121 121 # def init(self):
122 122 #
123 123 #
124 124 # if self.dataIn.type == 'AMISR':
125 125 # self.__updateObjFromAmisrInput()
126 126 #
127 127 # if self.dataIn.type == 'Voltage':
128 128 # self.dataOut.copy(self.dataIn)
129 129 # # No necesita copiar en cada init() los atributos de dataIn
130 130 # # la copia deberia hacerse por cada nuevo bloque de datos
131 131
132 132 def selectChannels(self, channelList):
133 133
134 134 channelIndexList = []
135 135
136 136 for channel in channelList:
137 137 if channel not in self.dataOut.channelList:
138 138 raise ValueError, "Channel %d is not in %s" %(channel, str(self.dataOut.channelList))
139 139
140 140 index = self.dataOut.channelList.index(channel)
141 141 channelIndexList.append(index)
142 142
143 143 self.selectChannelsByIndex(channelIndexList)
144 144
145 145 def selectChannelsByIndex(self, channelIndexList):
146 146 """
147 147 Selecciona un bloque de datos en base a canales segun el channelIndexList
148 148
149 149 Input:
150 150 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
151 151
152 152 Affected:
153 153 self.dataOut.data
154 154 self.dataOut.channelIndexList
155 155 self.dataOut.nChannels
156 156 self.dataOut.m_ProcessingHeader.totalSpectra
157 157 self.dataOut.systemHeaderObj.numChannels
158 158 self.dataOut.m_ProcessingHeader.blockSize
159 159
160 160 Return:
161 161 None
162 162 """
163 163
164 164 for channelIndex in channelIndexList:
165 165 if channelIndex not in self.dataOut.channelIndexList:
166 166 print channelIndexList
167 167 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
168 168
169 169 if self.dataOut.flagDataAsBlock:
170 170 """
171 171 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
172 172 """
173 173 data = self.dataOut.data[channelIndexList,:,:]
174 174 else:
175 175 data = self.dataOut.data[channelIndexList,:]
176 176
177 177 self.dataOut.data = data
178 178 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
179 179 # self.dataOut.nChannels = nChannels
180 180
181 181 return 1
182 182
183 183 def selectHeights(self, minHei=None, maxHei=None):
184 184 """
185 185 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
186 186 minHei <= height <= maxHei
187 187
188 188 Input:
189 189 minHei : valor minimo de altura a considerar
190 190 maxHei : valor maximo de altura a considerar
191 191
192 192 Affected:
193 193 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
194 194
195 195 Return:
196 196 1 si el metodo se ejecuto con exito caso contrario devuelve 0
197 197 """
198 198
199 199 if minHei == None:
200 200 minHei = self.dataOut.heightList[0]
201 201
202 202 if maxHei == None:
203 203 maxHei = self.dataOut.heightList[-1]
204 204
205 205 if (minHei < self.dataOut.heightList[0]):
206 206 minHei = self.dataOut.heightList[0]
207 207
208 208 if (maxHei > self.dataOut.heightList[-1]):
209 209 maxHei = self.dataOut.heightList[-1]
210 210
211 211 minIndex = 0
212 212 maxIndex = 0
213 213 heights = self.dataOut.heightList
214 214
215 215 inda = numpy.where(heights >= minHei)
216 216 indb = numpy.where(heights <= maxHei)
217 217
218 218 try:
219 219 minIndex = inda[0][0]
220 220 except:
221 221 minIndex = 0
222 222
223 223 try:
224 224 maxIndex = indb[0][-1]
225 225 except:
226 226 maxIndex = len(heights)
227 227
228 228 self.selectHeightsByIndex(minIndex, maxIndex)
229 229
230 230 return 1
231 231
232 232
233 233 def selectHeightsByIndex(self, minIndex, maxIndex):
234 234 """
235 235 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
236 236 minIndex <= index <= maxIndex
237 237
238 238 Input:
239 239 minIndex : valor de indice minimo de altura a considerar
240 240 maxIndex : valor de indice maximo de altura a considerar
241 241
242 242 Affected:
243 243 self.dataOut.data
244 244 self.dataOut.heightList
245 245
246 246 Return:
247 247 1 si el metodo se ejecuto con exito caso contrario devuelve 0
248 248 """
249 249
250 250 if (minIndex < 0) or (minIndex > maxIndex):
251 251 raise ValueError, "Height index range (%d,%d) is not valid" % (minIndex, maxIndex)
252 252
253 253 if (maxIndex >= self.dataOut.nHeights):
254 254 maxIndex = self.dataOut.nHeights
255 255
256 256 #voltage
257 257 if self.dataOut.flagDataAsBlock:
258 258 """
259 259 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
260 260 """
261 261 data = self.dataOut.data[:,:, minIndex:maxIndex]
262 262 else:
263 263 data = self.dataOut.data[:, minIndex:maxIndex]
264 264
265 265 # firstHeight = self.dataOut.heightList[minIndex]
266 266
267 267 self.dataOut.data = data
268 268 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
269 269
270 270 if self.dataOut.nHeights <= 1:
271 271 raise ValueError, "selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights)
272 272
273 273 return 1
274 274
275 275
276 276 def filterByHeights(self, window):
277 277
278 278 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
279 279
280 280 if window == None:
281 281 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
282 282
283 283 newdelta = deltaHeight * window
284 284 r = self.dataOut.nHeights % window
285 285 newheights = (self.dataOut.nHeights-r)/window
286 286
287 287 if newheights <= 1:
288 288 raise ValueError, "filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window)
289 289
290 290 if self.dataOut.flagDataAsBlock:
291 291 """
292 292 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
293 293 """
294 294 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
295 295 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
296 296 buffer = numpy.sum(buffer,3)
297 297
298 298 else:
299 299 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
300 300 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
301 301 buffer = numpy.sum(buffer,2)
302 302
303 303 self.dataOut.data = buffer
304 304 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
305 305 self.dataOut.windowOfFilter = window
306 306
307 307 def setH0(self, h0, deltaHeight = None):
308 308
309 309 if not deltaHeight:
310 310 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
311 311
312 312 nHeights = self.dataOut.nHeights
313 313
314 314 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
315 315
316 316 self.dataOut.heightList = newHeiRange
317 317
318 318 def deFlip(self, channelList = []):
319 319
320 320 data = self.dataOut.data.copy()
321 321
322 322 if self.dataOut.flagDataAsBlock:
323 323 flip = self.flip
324 324 profileList = range(self.dataOut.nProfiles)
325 325
326 326 if not channelList:
327 327 for thisProfile in profileList:
328 328 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
329 329 flip *= -1.0
330 330 else:
331 331 for thisChannel in channelList:
332 332 if thisChannel not in self.dataOut.channelList:
333 333 continue
334 334
335 335 for thisProfile in profileList:
336 336 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
337 337 flip *= -1.0
338 338
339 339 self.flip = flip
340 340
341 341 else:
342 342 if not channelList:
343 343 data[:,:] = data[:,:]*self.flip
344 344 else:
345 345 for thisChannel in channelList:
346 346 if thisChannel not in self.dataOut.channelList:
347 347 continue
348 348
349 349 data[thisChannel,:] = data[thisChannel,:]*self.flip
350 350
351 351 self.flip *= -1.
352 352
353 353 self.dataOut.data = data
354 354
355 355 def setRadarFrequency(self, frequency=None):
356 356
357 357 if frequency != None:
358 358 self.dataOut.frequency = frequency
359 359
360 360 return 1
361 361
362 362 def interpolateHeights(self, topLim, botLim):
363 363 #69 al 72 para julia
364 364 #82-84 para meteoros
365 365 if len(numpy.shape(self.dataOut.data))==2:
366 366 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
367 367 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
368 368 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
369 369 self.dataOut.data[:,botLim:topLim+1] = sampInterp
370 370 else:
371 371 nHeights = self.dataOut.data.shape[2]
372 372 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
373 373 y = self.dataOut.data[:,:,range(botLim)+range(topLim+1,nHeights)]
374 374 f = interpolate.interp1d(x, y, axis = 2)
375 375 xnew = numpy.arange(botLim,topLim+1)
376 376 ynew = f(xnew)
377 377
378 378 self.dataOut.data[:,:,botLim:topLim+1] = ynew
379 379
380 380 # import collections
381 381
382 382 class CohInt(Operation):
383 383
384 384 isConfig = False
385 385 __profIndex = 0
386 386 __byTime = False
387 387 __initime = None
388 388 __lastdatatime = None
389 389 __integrationtime = None
390 390 __buffer = None
391 391 __bufferStride = []
392 392 __dataReady = False
393 393 __profIndexStride = 0
394 394 __dataToPutStride = False
395 395 n = None
396 396
397 397 def __init__(self, **kwargs):
398 398
399 399 Operation.__init__(self, **kwargs)
400 400
401 401 # self.isConfig = False
402 402
403 403 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
404 404 """
405 405 Set the parameters of the integration class.
406 406
407 407 Inputs:
408 408
409 409 n : Number of coherent integrations
410 410 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
411 411 overlapping :
412 412 """
413 413
414 414 self.__initime = None
415 415 self.__lastdatatime = 0
416 416 self.__buffer = None
417 417 self.__dataReady = False
418 418 self.byblock = byblock
419 419 self.stride = stride
420 420
421 421 if n == None and timeInterval == None:
422 422 raise ValueError, "n or timeInterval should be specified ..."
423 423
424 424 if n != None:
425 425 self.n = n
426 426 self.__byTime = False
427 427 else:
428 428 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
429 429 self.n = 9999
430 430 self.__byTime = True
431 431
432 432 if overlapping:
433 433 self.__withOverlapping = True
434 434 self.__buffer = None
435 435 else:
436 436 self.__withOverlapping = False
437 437 self.__buffer = 0
438 438
439 439 self.__profIndex = 0
440 440
441 441 def putData(self, data):
442 442
443 443 """
444 444 Add a profile to the __buffer and increase in one the __profileIndex
445 445
446 446 """
447 447
448 448 if not self.__withOverlapping:
449 449 self.__buffer += data.copy()
450 450 self.__profIndex += 1
451 451 return
452 452
453 453 #Overlapping data
454 454 nChannels, nHeis = data.shape
455 455 data = numpy.reshape(data, (1, nChannels, nHeis))
456 456
457 457 #If the buffer is empty then it takes the data value
458 458 if self.__buffer is None:
459 459 self.__buffer = data
460 460 self.__profIndex += 1
461 461 return
462 462
463 463 #If the buffer length is lower than n then stakcing the data value
464 464 if self.__profIndex < self.n:
465 465 self.__buffer = numpy.vstack((self.__buffer, data))
466 466 self.__profIndex += 1
467 467 return
468 468
469 469 #If the buffer length is equal to n then replacing the last buffer value with the data value
470 470 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
471 471 self.__buffer[self.n-1] = data
472 472 self.__profIndex = self.n
473 473 return
474 474
475 475
476 476 def pushData(self):
477 477 """
478 478 Return the sum of the last profiles and the profiles used in the sum.
479 479
480 480 Affected:
481 481
482 482 self.__profileIndex
483 483
484 484 """
485 485
486 486 if not self.__withOverlapping:
487 487 data = self.__buffer
488 488 n = self.__profIndex
489 489
490 490 self.__buffer = 0
491 491 self.__profIndex = 0
492 492
493 493 return data, n
494 494
495 495 #Integration with Overlapping
496 496 data = numpy.sum(self.__buffer, axis=0)
497 497 # print data
498 498 # raise
499 499 n = self.__profIndex
500 500
501 501 return data, n
502 502
503 503 def byProfiles(self, data):
504 504
505 505 self.__dataReady = False
506 506 avgdata = None
507 507 # n = None
508 508 # print data
509 509 # raise
510 510 self.putData(data)
511 511
512 512 if self.__profIndex == self.n:
513 513 avgdata, n = self.pushData()
514 514 self.__dataReady = True
515 515
516 516 return avgdata
517 517
518 518 def byTime(self, data, datatime):
519 519
520 520 self.__dataReady = False
521 521 avgdata = None
522 522 n = None
523 523
524 524 self.putData(data)
525 525
526 526 if (datatime - self.__initime) >= self.__integrationtime:
527 527 avgdata, n = self.pushData()
528 528 self.n = n
529 529 self.__dataReady = True
530 530
531 531 return avgdata
532 532
533 533 def integrateByStride(self, data, datatime):
534 534 # print data
535 535 if self.__profIndex == 0:
536 536 self.__buffer = [[data.copy(), datatime]]
537 537 else:
538 538 self.__buffer.append([data.copy(),datatime])
539 539 self.__profIndex += 1
540 540 self.__dataReady = False
541 541
542 542 if self.__profIndex == self.n * self.stride :
543 543 self.__dataToPutStride = True
544 544 self.__profIndexStride = 0
545 545 self.__profIndex = 0
546 546 self.__bufferStride = []
547 547 for i in range(self.stride):
548 548 current = self.__buffer[i::self.stride]
549 549 data = numpy.sum([t[0] for t in current], axis=0)
550 550 avgdatatime = numpy.average([t[1] for t in current])
551 551 # print data
552 552 self.__bufferStride.append((data, avgdatatime))
553 553
554 554 if self.__dataToPutStride:
555 555 self.__dataReady = True
556 556 self.__profIndexStride += 1
557 557 if self.__profIndexStride == self.stride:
558 558 self.__dataToPutStride = False
559 559 # print self.__bufferStride[self.__profIndexStride - 1]
560 560 # raise
561 561 return self.__bufferStride[self.__profIndexStride - 1]
562 562
563 563
564 564 return None, None
565 565
566 566 def integrate(self, data, datatime=None):
567 567
568 568 if self.__initime == None:
569 569 self.__initime = datatime
570 570
571 571 if self.__byTime:
572 572 avgdata = self.byTime(data, datatime)
573 573 else:
574 574 avgdata = self.byProfiles(data)
575 575
576 576
577 577 self.__lastdatatime = datatime
578 578
579 579 if avgdata is None:
580 580 return None, None
581 581
582 582 avgdatatime = self.__initime
583 583
584 584 deltatime = datatime - self.__lastdatatime
585 585
586 586 if not self.__withOverlapping:
587 587 self.__initime = datatime
588 588 else:
589 589 self.__initime += deltatime
590 590
591 591 return avgdata, avgdatatime
592 592
593 593 def integrateByBlock(self, dataOut):
594 594
595 595 times = int(dataOut.data.shape[1]/self.n)
596 596 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
597 597
598 598 id_min = 0
599 599 id_max = self.n
600 600
601 601 for i in range(times):
602 602 junk = dataOut.data[:,id_min:id_max,:]
603 603 avgdata[:,i,:] = junk.sum(axis=1)
604 604 id_min += self.n
605 605 id_max += self.n
606 606
607 607 timeInterval = dataOut.ippSeconds*self.n
608 608 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
609 609 self.__dataReady = True
610 610 return avgdata, avgdatatime
611 611
612 612 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
613 613 if not self.isConfig:
614 614 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
615 615 self.isConfig = True
616 616
617 617 if dataOut.flagDataAsBlock:
618 618 """
619 619 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
620 620 """
621 621 avgdata, avgdatatime = self.integrateByBlock(dataOut)
622 622 dataOut.nProfiles /= self.n
623 623 else:
624 624 if stride is None:
625 625 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
626 626 else:
627 627 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
628 628
629 629
630 630 # dataOut.timeInterval *= n
631 631 dataOut.flagNoData = True
632 632
633 633 if self.__dataReady:
634 634 dataOut.data = avgdata
635 635 dataOut.nCohInt *= self.n
636 636 dataOut.utctime = avgdatatime
637 637 # print avgdata, avgdatatime
638 638 # raise
639 639 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
640 640 dataOut.flagNoData = False
641 641
642 642 class Decoder(Operation):
643 643
644 644 isConfig = False
645 645 __profIndex = 0
646 646
647 647 code = None
648 648
649 649 nCode = None
650 650 nBaud = None
651 651
652 652 def __init__(self, **kwargs):
653 653
654 654 Operation.__init__(self, **kwargs)
655 655
656 656 self.times = None
657 657 self.osamp = None
658 658 # self.__setValues = False
659 659 self.isConfig = False
660 660
661 661 def setup(self, code, osamp, dataOut):
662 662
663 663 self.__profIndex = 0
664 664
665 665 self.code = code
666 666
667 667 self.nCode = len(code)
668 668 self.nBaud = len(code[0])
669 669
670 670 if (osamp != None) and (osamp >1):
671 671 self.osamp = osamp
672 672 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
673 673 self.nBaud = self.nBaud*self.osamp
674 674
675 675 self.__nChannels = dataOut.nChannels
676 676 self.__nProfiles = dataOut.nProfiles
677 677 self.__nHeis = dataOut.nHeights
678 678
679 679 if self.__nHeis < self.nBaud:
680 680 raise ValueError, 'Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)
681 681
682 682 #Frequency
683 683 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
684 684
685 685 __codeBuffer[:,0:self.nBaud] = self.code
686 686
687 687 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
688 688
689 689 if dataOut.flagDataAsBlock:
690 690
691 691 self.ndatadec = self.__nHeis #- self.nBaud + 1
692 692
693 693 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
694 694
695 695 else:
696 696
697 697 #Time
698 698 self.ndatadec = self.__nHeis #- self.nBaud + 1
699 699
700 700 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
701 701
702 702 def __convolutionInFreq(self, data):
703 703
704 704 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
705 705
706 706 fft_data = numpy.fft.fft(data, axis=1)
707 707
708 708 conv = fft_data*fft_code
709 709
710 710 data = numpy.fft.ifft(conv,axis=1)
711 711
712 712 return data
713 713
714 714 def __convolutionInFreqOpt(self, data):
715 715
716 716 raise NotImplementedError
717 717
718 718 def __convolutionInTime(self, data):
719 719
720 720 code = self.code[self.__profIndex]
721 721 for i in range(self.__nChannels):
722 722 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
723 723
724 724 return self.datadecTime
725 725
726 726 def __convolutionByBlockInTime(self, data):
727 727
728 728 repetitions = self.__nProfiles / self.nCode
729 729
730 730 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
731 731 junk = junk.flatten()
732 732 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
733 733 profilesList = xrange(self.__nProfiles)
734 734
735 735 for i in range(self.__nChannels):
736 736 for j in profilesList:
737 737 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
738 738 return self.datadecTime
739 739
740 740 def __convolutionByBlockInFreq(self, data):
741 741
742 742 raise NotImplementedError, "Decoder by frequency fro Blocks not implemented"
743 743
744 744
745 745 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
746 746
747 747 fft_data = numpy.fft.fft(data, axis=2)
748 748
749 749 conv = fft_data*fft_code
750 750
751 751 data = numpy.fft.ifft(conv,axis=2)
752 752
753 753 return data
754 754
755 755
756 756 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
757 757
758 758 if dataOut.flagDecodeData:
759 759 print "This data is already decoded, recoding again ..."
760 760
761 761 if not self.isConfig:
762 762
763 763 if code is None:
764 764 if dataOut.code is None:
765 765 raise ValueError, "Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type
766 766
767 767 code = dataOut.code
768 768 else:
769 769 code = numpy.array(code).reshape(nCode,nBaud)
770 770 self.setup(code, osamp, dataOut)
771 771
772 772 self.isConfig = True
773 773
774 774 if mode == 3:
775 775 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
776 776
777 777 if times != None:
778 778 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
779 779
780 780 if self.code is None:
781 781 print "Fail decoding: Code is not defined."
782 782 return
783 783
784 784 self.__nProfiles = dataOut.nProfiles
785 785 datadec = None
786 786
787 787 if mode == 3:
788 788 mode = 0
789 789
790 790 if dataOut.flagDataAsBlock:
791 791 """
792 792 Decoding when data have been read as block,
793 793 """
794 794
795 795 if mode == 0:
796 796 datadec = self.__convolutionByBlockInTime(dataOut.data)
797 797 if mode == 1:
798 798 datadec = self.__convolutionByBlockInFreq(dataOut.data)
799 799 else:
800 800 """
801 801 Decoding when data have been read profile by profile
802 802 """
803 803 if mode == 0:
804 804 datadec = self.__convolutionInTime(dataOut.data)
805 805
806 806 if mode == 1:
807 807 datadec = self.__convolutionInFreq(dataOut.data)
808 808
809 809 if mode == 2:
810 810 datadec = self.__convolutionInFreqOpt(dataOut.data)
811 811
812 812 if datadec is None:
813 813 raise ValueError, "Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode
814 814
815 815 dataOut.code = self.code
816 816 dataOut.nCode = self.nCode
817 817 dataOut.nBaud = self.nBaud
818 818
819 819 dataOut.data = datadec
820 820
821 821 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
822 822
823 823 dataOut.flagDecodeData = True #asumo q la data esta decodificada
824 824
825 825 if self.__profIndex == self.nCode-1:
826 826 self.__profIndex = 0
827 827 return 1
828 828
829 829 self.__profIndex += 1
830 830
831 831 return 1
832 832 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
833 833
834 834
835 835 class ProfileConcat(Operation):
836 836
837 837 isConfig = False
838 838 buffer = None
839 839 concat_m =None
840 840
841 841 def __init__(self, **kwargs):
842 842
843 843 Operation.__init__(self, **kwargs)
844 844 self.profileIndex = 0
845 845
846 846 def reset(self):
847 847 self.buffer = numpy.zeros_like(self.buffer)
848 848 self.start_index = 0
849 849 self.times = 1
850 850
851 851 def setup(self, data, m, n=1):
852 852 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
853 853 self.nHeights = data.shape[1]#.nHeights
854 854 self.start_index = 0
855 855 self.times = 1
856 856
857 857 def concat(self, data):
858 858
859 859 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
860 860 self.start_index = self.start_index + self.nHeights
861 861
862 862 def run(self, dataOut, m):
863 863
864 864 self.concat_m= m
865 865 dataOut.flagNoData = True
866 866
867 867 if not self.isConfig:
868 868 self.setup(dataOut.data, m, 1)
869 869 self.isConfig = True
870 870
871 871 if dataOut.flagDataAsBlock:
872 872 raise ValueError, "ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False"
873 873
874 874 else:
875 875 self.concat(dataOut.data)
876 876 self.times += 1
877 877 if self.times > m:
878 878 dataOut.data = self.buffer
879 879 self.reset()
880 880 dataOut.flagNoData = False
881 881 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
882 882 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
883 883 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
884 884 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
885 885 dataOut.ippSeconds *= m
886 886 dataOut.concat_m = int(m)
887 887
888 888 class ProfileSelector(Operation):
889 889
890 890 profileIndex = None
891 891 # Tamanho total de los perfiles
892 892 nProfiles = None
893 893
894 894 def __init__(self, **kwargs):
895 895
896 896 Operation.__init__(self, **kwargs)
897 897 self.profileIndex = 0
898 898
899 899 def incProfileIndex(self):
900 900
901 901 self.profileIndex += 1
902 902
903 903 if self.profileIndex >= self.nProfiles:
904 904 self.profileIndex = 0
905 905
906 906 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
907 907
908 908 if profileIndex < minIndex:
909 909 return False
910 910
911 911 if profileIndex > maxIndex:
912 912 return False
913 913
914 914 return True
915 915
916 916 def isThisProfileInList(self, profileIndex, profileList):
917 917
918 918 if profileIndex not in profileList:
919 919 return False
920 920
921 921 return True
922 922
923 923 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
924 924
925 925 """
926 926 ProfileSelector:
927 927
928 928 Inputs:
929 929 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
930 930
931 931 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
932 932
933 933 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
934 934
935 935 """
936 936
937 937 if rangeList is not None:
938 938 if type(rangeList[0]) not in (tuple, list):
939 939 rangeList = [rangeList]
940 940
941 941 dataOut.flagNoData = True
942 942
943 943 if dataOut.flagDataAsBlock:
944 944 """
945 945 data dimension = [nChannels, nProfiles, nHeis]
946 946 """
947 947 if profileList != None:
948 948 dataOut.data = dataOut.data[:,profileList,:]
949 949
950 950 if profileRangeList != None:
951 951 minIndex = profileRangeList[0]
952 952 maxIndex = profileRangeList[1]
953 953 profileList = range(minIndex, maxIndex+1)
954 954
955 955 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
956 956
957 957 if rangeList != None:
958 958
959 959 profileList = []
960 960
961 961 for thisRange in rangeList:
962 962 minIndex = thisRange[0]
963 963 maxIndex = thisRange[1]
964 964
965 965 profileList.extend(range(minIndex, maxIndex+1))
966 966
967 967 dataOut.data = dataOut.data[:,profileList,:]
968 968
969 969 dataOut.nProfiles = len(profileList)
970 970 dataOut.profileIndex = dataOut.nProfiles - 1
971 971 dataOut.flagNoData = False
972 972
973 973 return True
974 974
975 975 """
976 976 data dimension = [nChannels, nHeis]
977 977 """
978 978
979 979 if profileList != None:
980 980
981 981 if self.isThisProfileInList(dataOut.profileIndex, profileList):
982 982
983 983 self.nProfiles = len(profileList)
984 984 dataOut.nProfiles = self.nProfiles
985 985 dataOut.profileIndex = self.profileIndex
986 986 dataOut.flagNoData = False
987 987
988 988 self.incProfileIndex()
989 989 return True
990 990
991 991 if profileRangeList != None:
992 992
993 993 minIndex = profileRangeList[0]
994 994 maxIndex = profileRangeList[1]
995 995
996 996 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
997 997
998 998 self.nProfiles = maxIndex - minIndex + 1
999 999 dataOut.nProfiles = self.nProfiles
1000 1000 dataOut.profileIndex = self.profileIndex
1001 1001 dataOut.flagNoData = False
1002 1002
1003 1003 self.incProfileIndex()
1004 1004 return True
1005 1005
1006 1006 if rangeList != None:
1007 1007
1008 1008 nProfiles = 0
1009 1009
1010 1010 for thisRange in rangeList:
1011 1011 minIndex = thisRange[0]
1012 1012 maxIndex = thisRange[1]
1013 1013
1014 1014 nProfiles += maxIndex - minIndex + 1
1015 1015
1016 1016 for thisRange in rangeList:
1017 1017
1018 1018 minIndex = thisRange[0]
1019 1019 maxIndex = thisRange[1]
1020 1020
1021 1021 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1022 1022
1023 1023 self.nProfiles = nProfiles
1024 1024 dataOut.nProfiles = self.nProfiles
1025 1025 dataOut.profileIndex = self.profileIndex
1026 1026 dataOut.flagNoData = False
1027 1027
1028 1028 self.incProfileIndex()
1029 1029
1030 1030 break
1031 1031
1032 1032 return True
1033 1033
1034 1034
1035 1035 if beam != None: #beam is only for AMISR data
1036 1036 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1037 1037 dataOut.flagNoData = False
1038 1038 dataOut.profileIndex = self.profileIndex
1039 1039
1040 1040 self.incProfileIndex()
1041 1041
1042 1042 return True
1043 1043
1044 1044 raise ValueError, "ProfileSelector needs profileList, profileRangeList or rangeList parameter"
1045 1045
1046 1046 return False
1047 1047
1048 1048 class Reshaper(Operation):
1049 1049
1050 1050 def __init__(self, **kwargs):
1051 1051
1052 1052 Operation.__init__(self, **kwargs)
1053 1053
1054 1054 self.__buffer = None
1055 1055 self.__nitems = 0
1056 1056
1057 1057 def __appendProfile(self, dataOut, nTxs):
1058 1058
1059 1059 if self.__buffer is None:
1060 1060 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1061 1061 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1062 1062
1063 1063 ini = dataOut.nHeights * self.__nitems
1064 1064 end = ini + dataOut.nHeights
1065 1065
1066 1066 self.__buffer[:, ini:end] = dataOut.data
1067 1067
1068 1068 self.__nitems += 1
1069 1069
1070 1070 return int(self.__nitems*nTxs)
1071 1071
1072 1072 def __getBuffer(self):
1073 1073
1074 1074 if self.__nitems == int(1./self.__nTxs):
1075 1075
1076 1076 self.__nitems = 0
1077 1077
1078 1078 return self.__buffer.copy()
1079 1079
1080 1080 return None
1081 1081
1082 1082 def __checkInputs(self, dataOut, shape, nTxs):
1083 1083
1084 1084 if shape is None and nTxs is None:
1085 1085 raise ValueError, "Reshaper: shape of factor should be defined"
1086 1086
1087 1087 if nTxs:
1088 1088 if nTxs < 0:
1089 1089 raise ValueError, "nTxs should be greater than 0"
1090 1090
1091 1091 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1092 1092 raise ValueError, "nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs))
1093 1093
1094 1094 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1095 1095
1096 1096 return shape, nTxs
1097 1097
1098 1098 if len(shape) != 2 and len(shape) != 3:
1099 1099 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)
1100 1100
1101 1101 if len(shape) == 2:
1102 1102 shape_tuple = [dataOut.nChannels]
1103 1103 shape_tuple.extend(shape)
1104 1104 else:
1105 1105 shape_tuple = list(shape)
1106 1106
1107 1107 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1108 1108
1109 1109 return shape_tuple, nTxs
1110 1110
1111 1111 def run(self, dataOut, shape=None, nTxs=None):
1112 1112
1113 1113 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1114 1114
1115 1115 dataOut.flagNoData = True
1116 1116 profileIndex = None
1117 1117
1118 1118 if dataOut.flagDataAsBlock:
1119 1119
1120 1120 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1121 1121 dataOut.flagNoData = False
1122 1122
1123 1123 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1124 1124
1125 1125 else:
1126 1126
1127 1127 if self.__nTxs < 1:
1128 1128
1129 1129 self.__appendProfile(dataOut, self.__nTxs)
1130 1130 new_data = self.__getBuffer()
1131 1131
1132 1132 if new_data is not None:
1133 1133 dataOut.data = new_data
1134 1134 dataOut.flagNoData = False
1135 1135
1136 1136 profileIndex = dataOut.profileIndex*nTxs
1137 1137
1138 1138 else:
1139 1139 raise ValueError, "nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)"
1140 1140
1141 1141 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1142 1142
1143 1143 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1144 1144
1145 1145 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1146 1146
1147 1147 dataOut.profileIndex = profileIndex
1148 1148
1149 1149 dataOut.ippSeconds /= self.__nTxs
1150 1150
1151 1151 class SplitProfiles(Operation):
1152 1152
1153 1153 def __init__(self, **kwargs):
1154 1154
1155 1155 Operation.__init__(self, **kwargs)
1156 1156
1157 1157 def run(self, dataOut, n):
1158 1158
1159 1159 dataOut.flagNoData = True
1160 1160 profileIndex = None
1161 1161
1162 1162 if dataOut.flagDataAsBlock:
1163 1163
1164 1164 #nchannels, nprofiles, nsamples
1165 1165 shape = dataOut.data.shape
1166 1166
1167 1167 if shape[2] % n != 0:
1168 1168 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[2])
1169 1169
1170 1170 new_shape = shape[0], shape[1]*n, shape[2]/n
1171 1171
1172 1172 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1173 1173 dataOut.flagNoData = False
1174 1174
1175 1175 profileIndex = int(dataOut.nProfiles/n) - 1
1176 1176
1177 1177 else:
1178 1178
1179 1179 raise ValueError, "Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)"
1180 1180
1181 1181 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1182 1182
1183 1183 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1184 1184
1185 1185 dataOut.nProfiles = int(dataOut.nProfiles*n)
1186 1186
1187 1187 dataOut.profileIndex = profileIndex
1188 1188
1189 1189 dataOut.ippSeconds /= n
1190 1190
1191 1191 class CombineProfiles(Operation):
1192 1192
1193 1193 def __init__(self, **kwargs):
1194 1194
1195 1195 Operation.__init__(self, **kwargs)
1196 1196
1197 1197 self.__remData = None
1198 1198 self.__profileIndex = 0
1199 1199
1200 1200 def run(self, dataOut, n):
1201 1201
1202 1202 dataOut.flagNoData = True
1203 1203 profileIndex = None
1204 1204
1205 1205 if dataOut.flagDataAsBlock:
1206 1206
1207 1207 #nchannels, nprofiles, nsamples
1208 1208 shape = dataOut.data.shape
1209 1209 new_shape = shape[0], shape[1]/n, shape[2]*n
1210 1210
1211 1211 if shape[1] % n != 0:
1212 1212 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[1])
1213 1213
1214 1214 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1215 1215 dataOut.flagNoData = False
1216 1216
1217 1217 profileIndex = int(dataOut.nProfiles*n) - 1
1218 1218
1219 1219 else:
1220 1220
1221 1221 #nchannels, nsamples
1222 1222 if self.__remData is None:
1223 1223 newData = dataOut.data
1224 1224 else:
1225 1225 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1226 1226
1227 1227 self.__profileIndex += 1
1228 1228
1229 1229 if self.__profileIndex < n:
1230 1230 self.__remData = newData
1231 1231 #continue
1232 1232 return
1233 1233
1234 1234 self.__profileIndex = 0
1235 1235 self.__remData = None
1236 1236
1237 1237 dataOut.data = newData
1238 1238 dataOut.flagNoData = False
1239 1239
1240 1240 profileIndex = dataOut.profileIndex/n
1241 1241
1242 1242
1243 1243 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1244 1244
1245 1245 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1246 1246
1247 1247 dataOut.nProfiles = int(dataOut.nProfiles/n)
1248 1248
1249 1249 dataOut.profileIndex = profileIndex
1250 1250
1251 1251 dataOut.ippSeconds *= n
1252 1252
1253 1253
1254 1254 class SSheightProfiles(Operation):
1255 1255
1256 1256 step = None
1257 1257 nsamples = None
1258 1258 bufferShape = None
1259 1259 profileShape = None
1260 1260 sshProfiles = None
1261 1261 profileIndex = None
1262 1262
1263 1263 def __init__(self, **kwargs):
1264 1264
1265 1265 Operation.__init__(self, **kwargs)
1266 1266 self.isConfig = False
1267 1267
1268 1268 def setup(self,dataOut ,step = None , nsamples = None):
1269 1269
1270 1270 if step == None and nsamples == None:
1271 1271 raise ValueError, "step or nheights should be specified ..."
1272 1272
1273 1273 self.step = step
1274 1274 self.nsamples = nsamples
1275 1275 self.__nChannels = dataOut.nChannels
1276 1276 self.__nProfiles = dataOut.nProfiles
1277 1277 self.__nHeis = dataOut.nHeights
1278 1278 shape = dataOut.data.shape #nchannels, nprofiles, nsamples
1279 1279 print "input nChannels",self.__nChannels
1280 1280 print "input nProfiles",self.__nProfiles
1281 1281 print "input nHeis",self.__nHeis
1282 1282 print "input Shape",shape
1283 1283
1284 1284
1285 1285
1286 1286 residue = (shape[1] - self.nsamples) % self.step
1287 1287 if residue != 0:
1288 1288 print "The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue)
1289 1289
1290 1290 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1291 1291 numberProfile = self.nsamples
1292 1292 numberSamples = (shape[1] - self.nsamples)/self.step
1293 1293
1294 1294 print "new numberProfile",numberProfile
1295 1295 print "new numberSamples",numberSamples
1296 1296
1297 1297 print "New number of profile: %d, number of height: %d, Resolution %f Km"%(numberProfile,numberSamples,deltaHeight*self.step)
1298 1298
1299 1299 self.bufferShape = shape[0], numberSamples, numberProfile # nchannels, nsamples , nprofiles
1300 1300 self.profileShape = shape[0], numberProfile, numberSamples # nchannels, nprofiles, nsamples
1301 1301
1302 1302 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1303 1303 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
1304 1304
1305 1305 def run(self, dataOut, step, nsamples):
1306 1306
1307 1307 dataOut.flagNoData = True
1308 1308 dataOut.flagDataAsBlock = False
1309 1309 profileIndex = None
1310 1310
1311 1311
1312 1312 if not self.isConfig:
1313 1313 self.setup(dataOut, step=step , nsamples=nsamples)
1314 1314 self.isConfig = True
1315 1315
1316 1316 for i in range(self.buffer.shape[1]):
1317 1317 #self.buffer[:,i] = numpy.flip(dataOut.data[:,i*self.step:i*self.step + self.nsamples])
1318 1318 self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]
1319
1320 #self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights])
1319 #self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights])
1321 1320
1322 1321 for j in range(self.buffer.shape[0]):
1323 1322 self.sshProfiles[j] = numpy.transpose(self.buffer[j])
1324 1323
1325 1324 profileIndex = self.nsamples
1326 1325 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1327 1326 ippSeconds = (deltaHeight*1.0e-6)/(0.15)
1327 #print "ippSeconds",ippSeconds
1328 1328 try:
1329 1329 if dataOut.concat_m is not None:
1330 1330 ippSeconds= ippSeconds/float(dataOut.concat_m)
1331 1331 #print "Profile concat %d"%dataOut.concat_m
1332 1332 except:
1333 1333 pass
1334 1334
1335 1335 dataOut.data = self.sshProfiles
1336 1336 dataOut.flagNoData = False
1337 1337 dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0]
1338 1338 dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
1339 1339 dataOut.profileIndex = profileIndex
1340 1340 dataOut.flagDataAsBlock = True
1341 1341 dataOut.ippSeconds = ippSeconds
1342 1342 dataOut.step = self.step
1343 1343
1344 1344
1345 1345 import time
1346 1346 #################################################
1347 1347
1348 1348 class decoPseudorandom(Operation):
1349 1349
1350 1350 nProfiles= 0
1351 1351 buffer= None
1352 1352 isConfig = False
1353 1353
1354 1354 def setup(self, clen= 10000,seed= 0,Nranges= 1000,oversample=1):
1355 1355 #code = create_pseudo_random_code(clen=clen, seed=seed)
1356 1356 code= rep_seq(create_pseudo_random_code(clen=clen, seed=seed),rep=oversample)
1357 1357 #print ("code_rx", code.shape)
1358 1358 #N = int(an_len/clen) # 100
1359 1359 B_cache = 0
1360 1360 r_cache = 0
1361 1361 B_cached = False
1362 1362 r = create_estimation_matrix(code=code, cache=True, rmax=Nranges)
1363 1363 #print ("code shape", code.shape)
1364 1364 #print ("seed",seed)
1365 1365 #print ("Code", code[0:10])
1366 1366 self.B = r['B']
1367 1367
1368 1368
1369 1369 def run (self,dataOut,length_code= 10000,seed= 0,Nranges= 1000,oversample=1):
1370 1370 #print((dataOut.data.shape))
1371 1371 if not self.isConfig:
1372 1372 self.setup(clen= length_code,seed= seed,Nranges= Nranges,oversample=oversample)
1373 1373 self.isConfig = True
1374 1374
1375 1375 dataOut.flagNoData = True
1376 1376 data =dataOut.data
1377 1377 #print "length_CODE",length_code
1378 1378 data_shape = (data.shape[1])
1379 1379 #print "data_shape",data_shape
1380 1380 n = (length_code /data_shape)
1381 1381 #print "we need this number of sample",n
1382 1382
1383 1383 if n>0 and self.buffer is None:
1384 1384 self.buffer = numpy.zeros([1, length_code], dtype=numpy.complex64)
1385 1385 self.buffer[0][0:data_shape] = data[0]
1386 1386 #print "FIRST CREATION",self.buffer.shape
1387 1387
1388 1388 else:
1389 1389 self.buffer[0][self.nProfiles*data_shape:(self.nProfiles+1)*data_shape]=data[0]
1390 1390
1391 1391 #print "buffer_shape",(self.buffer.shape)
1392 1392 self.nProfiles += 1
1393 1393 #print "count",self.nProfiles
1394 1394
1395 1395 if self.nProfiles== n:
1396 1396 temporal = numpy.dot(self.B, numpy.transpose(self.buffer))
1397 1397 #print temporal.shape
1398 1398 #import time
1399 1399 #time.sleep(40)
1400 1400 dataOut.data=numpy.transpose(temporal)
1401 1401
1402 1402 dataOut.flagNoData = False
1403 1403 self.buffer= None
1404 1404 self.nProfiles = 0
1405 1405
1406 1406 # import collections
1407 1407 # from scipy.stats import mode
1408 1408 #
1409 1409 # class Synchronize(Operation):
1410 1410 #
1411 1411 # isConfig = False
1412 1412 # __profIndex = 0
1413 1413 #
1414 1414 # def __init__(self, **kwargs):
1415 1415 #
1416 1416 # Operation.__init__(self, **kwargs)
1417 1417 # # self.isConfig = False
1418 1418 # self.__powBuffer = None
1419 1419 # self.__startIndex = 0
1420 1420 # self.__pulseFound = False
1421 1421 #
1422 1422 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1423 1423 #
1424 1424 # #Read data
1425 1425 #
1426 1426 # powerdB = dataOut.getPower(channel = channel)
1427 1427 # noisedB = dataOut.getNoise(channel = channel)[0]
1428 1428 #
1429 1429 # self.__powBuffer.extend(powerdB.flatten())
1430 1430 #
1431 1431 # dataArray = numpy.array(self.__powBuffer)
1432 1432 #
1433 1433 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1434 1434 #
1435 1435 # maxValue = numpy.nanmax(filteredPower)
1436 1436 #
1437 1437 # if maxValue < noisedB + 10:
1438 1438 # #No se encuentra ningun pulso de transmision
1439 1439 # return None
1440 1440 #
1441 1441 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1442 1442 #
1443 1443 # if len(maxValuesIndex) < 2:
1444 1444 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1445 1445 # return None
1446 1446 #
1447 1447 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1448 1448 #
1449 1449 # #Seleccionar solo valores con un espaciamiento de nSamples
1450 1450 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1451 1451 #
1452 1452 # if len(pulseIndex) < 2:
1453 1453 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1454 1454 # return None
1455 1455 #
1456 1456 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1457 1457 #
1458 1458 # #remover senales que se distancien menos de 10 unidades o muestras
1459 1459 # #(No deberian existir IPP menor a 10 unidades)
1460 1460 #
1461 1461 # realIndex = numpy.where(spacing > 10 )[0]
1462 1462 #
1463 1463 # if len(realIndex) < 2:
1464 1464 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1465 1465 # return None
1466 1466 #
1467 1467 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1468 1468 # realPulseIndex = pulseIndex[realIndex]
1469 1469 #
1470 1470 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1471 1471 #
1472 1472 # print "IPP = %d samples" %period
1473 1473 #
1474 1474 # self.__newNSamples = dataOut.nHeights #int(period)
1475 1475 # self.__startIndex = int(realPulseIndex[0])
1476 1476 #
1477 1477 # return 1
1478 1478 #
1479 1479 #
1480 1480 # def setup(self, nSamples, nChannels, buffer_size = 4):
1481 1481 #
1482 1482 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1483 1483 # maxlen = buffer_size*nSamples)
1484 1484 #
1485 1485 # bufferList = []
1486 1486 #
1487 1487 # for i in range(nChannels):
1488 1488 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1489 1489 # maxlen = buffer_size*nSamples)
1490 1490 #
1491 1491 # bufferList.append(bufferByChannel)
1492 1492 #
1493 1493 # self.__nSamples = nSamples
1494 1494 # self.__nChannels = nChannels
1495 1495 # self.__bufferList = bufferList
1496 1496 #
1497 1497 # def run(self, dataOut, channel = 0):
1498 1498 #
1499 1499 # if not self.isConfig:
1500 1500 # nSamples = dataOut.nHeights
1501 1501 # nChannels = dataOut.nChannels
1502 1502 # self.setup(nSamples, nChannels)
1503 1503 # self.isConfig = True
1504 1504 #
1505 1505 # #Append new data to internal buffer
1506 1506 # for thisChannel in range(self.__nChannels):
1507 1507 # bufferByChannel = self.__bufferList[thisChannel]
1508 1508 # bufferByChannel.extend(dataOut.data[thisChannel])
1509 1509 #
1510 1510 # if self.__pulseFound:
1511 1511 # self.__startIndex -= self.__nSamples
1512 1512 #
1513 1513 # #Finding Tx Pulse
1514 1514 # if not self.__pulseFound:
1515 1515 # indexFound = self.__findTxPulse(dataOut, channel)
1516 1516 #
1517 1517 # if indexFound == None:
1518 1518 # dataOut.flagNoData = True
1519 1519 # return
1520 1520 #
1521 1521 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1522 1522 # self.__pulseFound = True
1523 1523 # self.__startIndex = indexFound
1524 1524 #
1525 1525 # #If pulse was found ...
1526 1526 # for thisChannel in range(self.__nChannels):
1527 1527 # bufferByChannel = self.__bufferList[thisChannel]
1528 1528 # #print self.__startIndex
1529 1529 # x = numpy.array(bufferByChannel)
1530 1530 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1531 1531 #
1532 1532 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1533 1533 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1534 1534 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1535 1535 #
1536 1536 # dataOut.data = self.__arrayBuffer
1537 1537 #
1538 1538 # self.__startIndex += self.__newNSamples
1539 1539 #
1540 1540 # return
Show file before | Show file after | Hide comments
@@ -1,144 +1,145
1 1 #!/usr/bin/env python
2 2 '''
3 3 Created on Jul 7, 2014
4 4
5 5 @author: roj-idl71
6 6 '''
7 7 import os, sys
8 8
9 9 #path = os.path.dirname(os.getcwd())
10 10 #path = os.path.join(path, 'source')
11 11 #sys.path.insert(0, path)
12 12
13 13 from schainpy.controller import Project
14 14
15 15 if __name__ == '__main__':
16 16
17 17 desc = "High altitude experiment SHORT "
18 18 filename = "schain.xml"
19 19 dpath = '/home/soporte/test_avp'
20 20 figpath = "/home/soporte/pics"
21 21 remotefolder = "/home/wmaster/graficos"
22 22 t=['0','24']
23 23 db_range=['15','35']
24 24 db_range=['25','45']
25 25 db_range=['0','35']
26 26 period=60
27 27
28 28 controllerObj = Project()
29 29
30 30 controllerObj.setup(id = '191', name='test01', description=desc)
31 31
32 32 readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage',
33 33 path=dpath,
34 34 startDate='2019/10/14',
35 35 # startDate='2018/06/18',
36 36 endDate='2019/10/14',
37 37 # endDate='2018/06/18',
38 38 startTime='14:00:00',
39 39 endTime='23:59:59',
40 40 online=0,
41 41 walk=0,
42 42 expLabel='',
43 43 delay=20)
44 44
45 45
46 46 opObj00 = readUnitConfObj.addOperation(name='printInfo')
47 47 opObj00 = readUnitConfObj.addOperation(name='printNumberOfBlock')
48 48
49 49 procUnitConfObj1 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId())
50 50
51 51 opObj11 = procUnitConfObj1.addOperation(name='ProfileSelector', optype='other')
52 52 #opObj11.addParameter(name='profileList', value='(0, 1, 2, 3, 20, 21, 22, 23, 40, 41, 42, 43, 60, 61, 62, 63, 80, 81, 82, 83)', format='intlist')
53 53 opObj11.addParameter(name='rangeList', value='(0,127)', format='intlist') # (0,63),(127,149)
54 54
55 55 opObj11 = procUnitConfObj1.addOperation(name='selectHeights')
56 56 opObj11.addParameter(name='minHei', value='0', format='int')
57 57 opObj11.addParameter(name='maxHei', value='3000', format='int')
58 58
59 59
60 60 #opObj11 = procUnitConfObj1.addOperation(name='filterByHeights')
61 61 #opObj11.addParameter(name='window', value='4', format='int')
62 62
63 63 # opObj11 = procUnitConfObj1.addOperation(name='deFlip')
64 64 # opObj11.addParameter(name='channelList', value='1,3,5,7', format='intlist')
65 65
66 66 opObj11 = procUnitConfObj1.addOperation(name='SSheightProfiles', optype='other')
67 67 opObj11.addParameter(name='step', value='1', format='int')
68 68 opObj11.addParameter(name='nsamples', value='64', format='int')
69 69
70 70
71 71 procUnitConfObj1SPC = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj1.getId())
72 72 procUnitConfObj1SPC.addParameter(name='nFFTPoints', value='128', format='int')
73 73 procUnitConfObj1SPC.addParameter(name='nProfiles', value='64', format='int')
74 74 # procUnitConfObj1SPC.addParameter(name='pairsList', value='(1,0),(3,2),(5,4),(7,6)', format='pairsList')
75 75
76 76 opObj11 = procUnitConfObj1SPC.addOperation(name='IncohInt', optype='other')
77 opObj11.addParameter(name='timeInterval', value='5', format='int')
77 opObj11.addParameter(name='timeInterval', value='10', format='int')
78 78
79 79 procUnitConfObj2SPC = controllerObj.addProcUnit(datatype='SpectraAFCProc', inputId=procUnitConfObj1SPC.getId())
80 procUnitConfObj2SPC.addParameter(name='real', value='1', format='int')
80 81
81 82 opObj11 = procUnitConfObj1SPC.addOperation(name='SpectraPlot', optype='other')
82 83 opObj11.addParameter(name='id', value='1', format='int')
83 84 opObj11.addParameter(name='wintitle', value='October 2019', format='str')
84 85 #opObj11.addParameter(name='xaxis', value='time', format='str')
85 86
86 87 #opObj11.addParameter(name='zmin', value=db_range[0], format='int')
87 88 #opObj11.addParameter(name='zmax', value=db_range[1], format='int')
88 89 #opObj11.addParameter(name='ymin', value='0', format='int')
89 90 #opObj11.addParameter(name='ymax', value='1500', format='int')
90 91 #opObj11.addParameter(name='xmin', value='-5', format='int')
91 92 #opObj11.addParameter(name='xmax', value='5', format='int')
92 93 opObj11.addParameter(name='showprofile', value='1', format='int')
93 94 opObj11.addParameter(name='save', value='1', format='int')
94 95 opObj11.addParameter(name='figpath', value=figpath)
95 96
96 97 opObj11 = procUnitConfObj1SPC.addOperation(name='RTIPlot', optype='other')
97 98 opObj11.addParameter(name='id', value='2', format='int')
98 99 opObj11.addParameter(name='wintitle', value='October 2019', format='str')
99 100 opObj11.addParameter(name='xmin', value=t[0], format='float')
100 101 opObj11.addParameter(name='xmax', value=t[1], format='float')
101 102 #opObj11.addParameter(name='ymin', value='0', format='int')
102 103 #opObj11.addParameter(name='ymax', value='1500', format='int')
103 104 #opObj11.addParameter(name='zmin', value=db_range[0], format='int')
104 105 #opObj11.addParameter(name='zmax', value=db_range[1], format='int')
105 106 opObj11.addParameter(name='showprofile', value='1', format='int')
106 107 opObj11.addParameter(name='save', value='1', format='int')
107 108 opObj11.addParameter(name='figpath', value=figpath)
108 109
109 110
110 111
111 112 opObj11 = procUnitConfObj2SPC.addOperation(name='ACFPlot', optype='other')
112 113 opObj11.addParameter(name='id', value='3', format='int')
113 114 opObj11.addParameter(name='wintitle', value = 'October 2019', format='str')
114 115 opObj11.addParameter(name='xaxis', value='time', format='str')
115 116 opObj11.addParameter(name='channel', value='0', format='int')
116 opObj11.addParameter(name='nSampleList', value='(40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70)', format='intList')
117 opObj11.addParameter(name='resolutionFactor', value='2', format='int')
117 #opObj11.addParameter(name='nSampleList', value='(50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70)', format='intList')
118 opObj11.addParameter(name='resolutionFactor', value='15', format='int')
118 119 #opObj11.addParameter(name='zmin', value=0.5, format='int')
119 120 #opObj11.addParameter(name='zmax', value=-0.5, format='int')
120 121 #opObj11.addParameter(name='ymin', value='0', format='int')
121 122 #opObj11.addParameter(name='ymax', value='0.5', format='int')
122 123 #opObj11.addParameter(name='xmin', value='-1.2', format='int')
123 124 #opObj11.addParameter(name='xmax', value='1.2', format='int')
124 125 opObj11.addParameter(name='show', value='1', format='int')
125 126 opObj11.addParameter(name='save', value='1', format='int')
126 127 opObj11.addParameter(name='figpath', value=figpath)
127 128
128 129 '''
129 130 opObj11 = procUnitConfObj1SPC.addOperation(name='Noise', optype='other')
130 131 opObj11.addParameter(name='id', value='3', format='int')
131 132 opObj11.addParameter(name='wintitle', value='Short Valley August 2019', format='str')
132 133 opObj11.addParameter(name='xmin', value=t[0], format='float')
133 134 opObj11.addParameter(name='xmax', value=t[1], format='float')
134 135 opObj11.addParameter(name='ymin', value=db_range[0], format='int')
135 136 opObj11.addParameter(name='ymax', value=db_range[1], format='int')
136 137 opObj11.addParameter(name='showprofile', value='1', format='int')
137 138 opObj11.addParameter(name='save', value='1', format='int')
138 139 opObj11.addParameter(name='figpath', value=figpath)
139 140
140 141
141 142 '''
142 143
143 144
144 145 controllerObj.start()
General Comments 0
You need to be logged in to leave comments. Login now