##// END OF EJS Templates
Filtering AMISR files for Datetime Range...
Daniel Valdez -
r499:f28411fcbfa0
parent child
Show More
@@ -1,61 +1,65
1 1 import numpy
2 2
3 3 class AMISR:
4 4 def __init__(self):
5 5 self.flagNoData = True
6 6 self.data = None
7 7 self.utctime = None
8 8 self.type = "AMISR"
9 9
10 10 #propiedades para compatibilidad con Voltages
11 11 self.timeZone = 0#timezone like jroheader, difference in minutes between UTC and localtime
12 12 self.dstFlag = 0#self.dataIn.dstFlag
13 13 self.errorCount = 0#self.dataIn.errorCount
14 14 self.useLocalTime = True#self.dataIn.useLocalTime
15 15
16 16 self.radarControllerHeaderObj = None#self.dataIn.radarControllerHeaderObj.copy()
17 17 self.systemHeaderObj = None#self.dataIn.systemHeaderObj.copy()
18 18 self.channelList = [0]#self.dataIn.channelList esto solo aplica para el caso de AMISR
19 19 self.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
20 20
21 21 self.flagTimeBlock = None#self.dataIn.flagTimeBlock
22 22 #self.utctime = #self.firstdatatime
23 23 self.flagDecodeData = None#self.dataIn.flagDecodeData #asumo q la data esta decodificada
24 24 self.flagDeflipData = None#self.dataIn.flagDeflipData #asumo q la data esta sin flip
25 25
26 26 self.nCohInt = 1#self.dataIn.nCohInt
27 27 self.nIncohInt = 1
28 28 self.ippSeconds = None#self.dataIn.ippSeconds, segun el filename/Setup/Tufile
29 29 self.windowOfFilter = None#self.dataIn.windowOfFilter
30 30
31 31 self.timeInterval = None#self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt
32 32 self.frequency = None#self.dataIn.frequency
33 33 self.realtime = 0#self.dataIn.realtime
34 34
35 35 #actualizar en la lectura de datos
36 36 self.heightList = None#self.dataIn.heightList
37 37 self.nProfiles = None#Number of samples or nFFTPoints
38 38 self.nRecords = None
39 39 self.nBeams = None
40 40 self.nBaud = None#self.dataIn.nBaud
41 41 self.nCode = None#self.dataIn.nCode
42 42 self.code = None#self.dataIn.code
43 43
44 44 #consideracion para los Beams
45 45 self.beamCodeDict = None
46 46 self.beamRangeDict = None
47
47 self.beamcode = None
48 self.azimuth = None
49 self.zenith = None
50 self.gain = None
51
48 52 self.npulseByFrame = None
49 53
50 54 def copy(self, inputObj=None):
51 55
52 56 if inputObj == None:
53 57 return copy.deepcopy(self)
54 58
55 59 for key in inputObj.__dict__.keys():
56 60 self.__dict__[key] = inputObj.__dict__[key]
57 61
58 62
59 63 def isEmpty(self):
60 64
61 65 return self.flagNoData No newline at end of file
@@ -1,736 +1,740
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
13 13 def getNumpyDtype(dataTypeCode):
14 14
15 15 if dataTypeCode == 0:
16 16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 17 elif dataTypeCode == 1:
18 18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 19 elif dataTypeCode == 2:
20 20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 21 elif dataTypeCode == 3:
22 22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 23 elif dataTypeCode == 4:
24 24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 25 elif dataTypeCode == 5:
26 26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 27 else:
28 28 raise ValueError, 'dataTypeCode was not defined'
29 29
30 30 return numpyDtype
31 31
32 32 def getDataTypeCode(numpyDtype):
33 33
34 34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 35 datatype = 0
36 36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 37 datatype = 1
38 38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 39 datatype = 2
40 40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 41 datatype = 3
42 42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 43 datatype = 4
44 44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 45 datatype = 5
46 46 else:
47 47 datatype = None
48 48
49 49 return datatype
50 50
51 51 def hildebrand_sekhon(data, navg):
52 52
53 53 data = data.copy()
54 54
55 55 sortdata = numpy.sort(data,axis=None)
56 56 lenOfData = len(sortdata)
57 57 nums_min = lenOfData/10
58 58
59 59 if (lenOfData/10) > 2:
60 60 nums_min = lenOfData/10
61 61 else:
62 62 nums_min = 2
63 63
64 64 sump = 0.
65 65
66 66 sumq = 0.
67 67
68 68 j = 0
69 69
70 70 cont = 1
71 71
72 72 while((cont==1)and(j<lenOfData)):
73 73
74 74 sump += sortdata[j]
75 75
76 76 sumq += sortdata[j]**2
77 77
78 78 j += 1
79 79
80 80 if j > nums_min:
81 81 rtest = float(j)/(j-1) + 1.0/navg
82 82 if ((sumq*j) > (rtest*sump**2)):
83 83 j = j - 1
84 84 sump = sump - sortdata[j]
85 85 sumq = sumq - sortdata[j]**2
86 86 cont = 0
87 87
88 88 lnoise = sump /j
89 89 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
90 90 return lnoise
91 91
92 92 class GenericData(object):
93 93
94 94 flagNoData = True
95 95
96 96 def __init__(self):
97 97
98 98 raise ValueError, "This class has not been implemented"
99 99
100 100 def copy(self, inputObj=None):
101 101
102 102 if inputObj == None:
103 103 return copy.deepcopy(self)
104 104
105 105 for key in inputObj.__dict__.keys():
106 106 self.__dict__[key] = inputObj.__dict__[key]
107 107
108 108 def deepcopy(self):
109 109
110 110 return copy.deepcopy(self)
111 111
112 112 def isEmpty(self):
113 113
114 114 return self.flagNoData
115 115
116 116 class JROData(GenericData):
117 117
118 118 # m_BasicHeader = BasicHeader()
119 119 # m_ProcessingHeader = ProcessingHeader()
120 120
121 121 systemHeaderObj = SystemHeader()
122 122
123 123 radarControllerHeaderObj = RadarControllerHeader()
124 124
125 125 # data = None
126 126
127 127 type = None
128 128
129 129 datatype = None #dtype but in string
130 130
131 131 # dtype = None
132 132
133 133 # nChannels = None
134 134
135 135 # nHeights = None
136 136
137 137 nProfiles = None
138 138
139 139 heightList = None
140 140
141 141 channelList = None
142 142
143 143 flagTimeBlock = False
144 144
145 145 useLocalTime = False
146 146
147 147 utctime = None
148 148
149 149 timeZone = None
150 150
151 151 dstFlag = None
152 152
153 153 errorCount = None
154 154
155 155 blocksize = None
156 156
157 157 nCode = None
158 158
159 159 nBaud = None
160 160
161 161 code = None
162 162
163 163 flagDecodeData = False #asumo q la data no esta decodificada
164 164
165 165 flagDeflipData = False #asumo q la data no esta sin flip
166 166
167 167 flagShiftFFT = False
168 168
169 169 # ippSeconds = None
170 170
171 171 timeInterval = None
172 172
173 173 nCohInt = None
174 174
175 175 noise = None
176 176
177 177 windowOfFilter = 1
178 178
179 179 #Speed of ligth
180 180 C = 3e8
181 181
182 182 frequency = 49.92e6
183 183
184 184 realtime = False
185 185
186 186 beacon_heiIndexList = None
187 187
188 188 last_block = None
189 189
190 190 blocknow = None
191 191
192 azimuth = None
193
194 zenith = None
195
192 196 def __init__(self):
193 197
194 198 raise ValueError, "This class has not been implemented"
195 199
196 200 def getNoise(self):
197 201
198 202 raise ValueError, "Not implemented"
199 203
200 204 def getNChannels(self):
201 205
202 206 return len(self.channelList)
203 207
204 208 def getChannelIndexList(self):
205 209
206 210 return range(self.nChannels)
207 211
208 212 def getNHeights(self):
209 213
210 214 return len(self.heightList)
211 215
212 216 def getHeiRange(self, extrapoints=0):
213 217
214 218 heis = self.heightList
215 219 # deltah = self.heightList[1] - self.heightList[0]
216 220 #
217 221 # heis.append(self.heightList[-1])
218 222
219 223 return heis
220 224
221 225 def getltctime(self):
222 226
223 227 if self.useLocalTime:
224 228 return self.utctime - self.timeZone*60
225 229
226 230 return self.utctime
227 231
228 232 def getDatatime(self):
229 233
230 234 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
231 235 return datatimeValue
232 236
233 237 def getTimeRange(self):
234 238
235 239 datatime = []
236 240
237 241 datatime.append(self.ltctime)
238 242 datatime.append(self.ltctime + self.timeInterval)
239 243
240 244 datatime = numpy.array(datatime)
241 245
242 246 return datatime
243 247
244 248 def getFmax(self):
245 249
246 250 PRF = 1./(self.ippSeconds * self.nCohInt)
247 251
248 252 fmax = PRF/2.
249 253
250 254 return fmax
251 255
252 256 def getVmax(self):
253 257
254 258 _lambda = self.C/self.frequency
255 259
256 260 vmax = self.getFmax() * _lambda
257 261
258 262 return vmax
259 263
260 264 def get_ippSeconds(self):
261 265 '''
262 266 '''
263 267 return self.radarControllerHeaderObj.ippSeconds
264 268
265 269 def set_ippSeconds(self, ippSeconds):
266 270 '''
267 271 '''
268 272
269 273 self.radarControllerHeaderObj.ippSeconds = ippSeconds
270 274
271 275 return
272 276
273 277 def get_dtype(self):
274 278 '''
275 279 '''
276 280 return getNumpyDtype(self.datatype)
277 281
278 282 def set_dtype(self, numpyDtype):
279 283 '''
280 284 '''
281 285
282 286 self.datatype = getDataTypeCode(numpyDtype)
283 287
284 288 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
285 289 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
286 290 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
287 291 #noise = property(getNoise, "I'm the 'nHeights' property.")
288 292 datatime = property(getDatatime, "I'm the 'datatime' property")
289 293 ltctime = property(getltctime, "I'm the 'ltctime' property")
290 294 ippSeconds = property(get_ippSeconds, set_ippSeconds)
291 295 dtype = property(get_dtype, set_dtype)
292 296
293 297 class Voltage(JROData):
294 298
295 299 #data es un numpy array de 2 dmensiones (canales, alturas)
296 300 data = None
297 301
298 302 def __init__(self):
299 303 '''
300 304 Constructor
301 305 '''
302 306
303 307 self.radarControllerHeaderObj = RadarControllerHeader()
304 308
305 309 self.systemHeaderObj = SystemHeader()
306 310
307 311 self.type = "Voltage"
308 312
309 313 self.data = None
310 314
311 315 # self.dtype = None
312 316
313 317 # self.nChannels = 0
314 318
315 319 # self.nHeights = 0
316 320
317 321 self.nProfiles = None
318 322
319 323 self.heightList = None
320 324
321 325 self.channelList = None
322 326
323 327 # self.channelIndexList = None
324 328
325 329 self.flagNoData = True
326 330
327 331 self.flagTimeBlock = False
328 332
329 333 self.utctime = None
330 334
331 335 self.timeZone = None
332 336
333 337 self.dstFlag = None
334 338
335 339 self.errorCount = None
336 340
337 341 self.nCohInt = None
338 342
339 343 self.blocksize = None
340 344
341 345 self.flagDecodeData = False #asumo q la data no esta decodificada
342 346
343 347 self.flagDeflipData = False #asumo q la data no esta sin flip
344 348
345 349 self.flagShiftFFT = False
346 350
347 351
348 352 def getNoisebyHildebrand(self):
349 353 """
350 354 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
351 355
352 356 Return:
353 357 noiselevel
354 358 """
355 359
356 360 for channel in range(self.nChannels):
357 361 daux = self.data_spc[channel,:,:]
358 362 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
359 363
360 364 return self.noise
361 365
362 366 def getNoise(self, type = 1):
363 367
364 368 self.noise = numpy.zeros(self.nChannels)
365 369
366 370 if type == 1:
367 371 noise = self.getNoisebyHildebrand()
368 372
369 373 return 10*numpy.log10(noise)
370 374
371 375 noise = property(getNoise, "I'm the 'nHeights' property.")
372 376
373 377 class Spectra(JROData):
374 378
375 379 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
376 380 data_spc = None
377 381
378 382 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
379 383 data_cspc = None
380 384
381 385 #data es un numpy array de 2 dmensiones (canales, alturas)
382 386 data_dc = None
383 387
384 388 nFFTPoints = None
385 389
386 390 # nPairs = None
387 391
388 392 pairsList = None
389 393
390 394 nIncohInt = None
391 395
392 396 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
393 397
394 398 nCohInt = None #se requiere para determinar el valor de timeInterval
395 399
396 400 ippFactor = None
397 401
398 402 def __init__(self):
399 403 '''
400 404 Constructor
401 405 '''
402 406
403 407 self.radarControllerHeaderObj = RadarControllerHeader()
404 408
405 409 self.systemHeaderObj = SystemHeader()
406 410
407 411 self.type = "Spectra"
408 412
409 413 # self.data = None
410 414
411 415 # self.dtype = None
412 416
413 417 # self.nChannels = 0
414 418
415 419 # self.nHeights = 0
416 420
417 421 self.nProfiles = None
418 422
419 423 self.heightList = None
420 424
421 425 self.channelList = None
422 426
423 427 # self.channelIndexList = None
424 428
425 429 self.pairsList = None
426 430
427 431 self.flagNoData = True
428 432
429 433 self.flagTimeBlock = False
430 434
431 435 self.utctime = None
432 436
433 437 self.nCohInt = None
434 438
435 439 self.nIncohInt = None
436 440
437 441 self.blocksize = None
438 442
439 443 self.nFFTPoints = None
440 444
441 445 self.wavelength = None
442 446
443 447 self.flagDecodeData = False #asumo q la data no esta decodificada
444 448
445 449 self.flagDeflipData = False #asumo q la data no esta sin flip
446 450
447 451 self.flagShiftFFT = False
448 452
449 453 self.ippFactor = 1
450 454
451 455 #self.noise = None
452 456
453 457 self.beacon_heiIndexList = []
454 458
455 459 self.noise_estimation = None
456 460
457 461
458 462 def getNoisebyHildebrand(self):
459 463 """
460 464 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
461 465
462 466 Return:
463 467 noiselevel
464 468 """
465 469
466 470 noise = numpy.zeros(self.nChannels)
467 471 for channel in range(self.nChannels):
468 472 daux = self.data_spc[channel,:,:]
469 473 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
470 474
471 475 return noise
472 476
473 477 def getNoise(self):
474 478 if self.noise_estimation != None:
475 479 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
476 480 else:
477 481 noise = self.getNoisebyHildebrand()
478 482 return noise
479 483
480 484
481 485 def getFreqRange(self, extrapoints=0):
482 486
483 487 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
484 488 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
485 489
486 490 return freqrange
487 491
488 492 def getVelRange(self, extrapoints=0):
489 493
490 494 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
491 495 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
492 496
493 497 return velrange
494 498
495 499 def getNPairs(self):
496 500
497 501 return len(self.pairsList)
498 502
499 503 def getPairsIndexList(self):
500 504
501 505 return range(self.nPairs)
502 506
503 507 def getNormFactor(self):
504 508 pwcode = 1
505 509 if self.flagDecodeData:
506 510 pwcode = numpy.sum(self.code[0]**2)
507 511 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
508 512 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
509 513
510 514 return normFactor
511 515
512 516 def getFlagCspc(self):
513 517
514 518 if self.data_cspc == None:
515 519 return True
516 520
517 521 return False
518 522
519 523 def getFlagDc(self):
520 524
521 525 if self.data_dc == None:
522 526 return True
523 527
524 528 return False
525 529
526 530 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
527 531 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
528 532 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
529 533 flag_cspc = property(getFlagCspc)
530 534 flag_dc = property(getFlagDc)
531 535 noise = property(getNoise, "I'm the 'nHeights' property.")
532 536
533 537 class SpectraHeis(Spectra):
534 538
535 539 data_spc = None
536 540
537 541 data_cspc = None
538 542
539 543 data_dc = None
540 544
541 545 nFFTPoints = None
542 546
543 547 # nPairs = None
544 548
545 549 pairsList = None
546 550
547 551 nIncohInt = None
548 552
549 553 def __init__(self):
550 554
551 555 self.radarControllerHeaderObj = RadarControllerHeader()
552 556
553 557 self.systemHeaderObj = SystemHeader()
554 558
555 559 self.type = "SpectraHeis"
556 560
557 561 # self.dtype = None
558 562
559 563 # self.nChannels = 0
560 564
561 565 # self.nHeights = 0
562 566
563 567 self.nProfiles = None
564 568
565 569 self.heightList = None
566 570
567 571 self.channelList = None
568 572
569 573 # self.channelIndexList = None
570 574
571 575 self.flagNoData = True
572 576
573 577 self.flagTimeBlock = False
574 578
575 579 # self.nPairs = 0
576 580
577 581 self.utctime = None
578 582
579 583 self.blocksize = None
580 584
581 585 def getNormFactor(self):
582 586 pwcode = 1
583 587 if self.flagDecodeData:
584 588 pwcode = numpy.sum(self.code[0]**2)
585 589
586 590 normFactor = self.nIncohInt*self.nCohInt*pwcode
587 591
588 592 return normFactor
589 593
590 594 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
591 595
592 596 class Fits:
593 597
594 598 heightList = None
595 599
596 600 channelList = None
597 601
598 602 flagNoData = True
599 603
600 604 flagTimeBlock = False
601 605
602 606 useLocalTime = False
603 607
604 608 utctime = None
605 609
606 610 timeZone = None
607 611
608 612 # ippSeconds = None
609 613
610 614 timeInterval = None
611 615
612 616 nCohInt = None
613 617
614 618 nIncohInt = None
615 619
616 620 noise = None
617 621
618 622 windowOfFilter = 1
619 623
620 624 #Speed of ligth
621 625 C = 3e8
622 626
623 627 frequency = 49.92e6
624 628
625 629 realtime = False
626 630
627 631
628 632 def __init__(self):
629 633
630 634 self.type = "Fits"
631 635
632 636 self.nProfiles = None
633 637
634 638 self.heightList = None
635 639
636 640 self.channelList = None
637 641
638 642 # self.channelIndexList = None
639 643
640 644 self.flagNoData = True
641 645
642 646 self.utctime = None
643 647
644 648 self.nCohInt = None
645 649
646 650 self.nIncohInt = None
647 651
648 652 self.useLocalTime = True
649 653
650 654 # self.utctime = None
651 655 # self.timeZone = None
652 656 # self.ltctime = None
653 657 # self.timeInterval = None
654 658 # self.header = None
655 659 # self.data_header = None
656 660 # self.data = None
657 661 # self.datatime = None
658 662 # self.flagNoData = False
659 663 # self.expName = ''
660 664 # self.nChannels = None
661 665 # self.nSamples = None
662 666 # self.dataBlocksPerFile = None
663 667 # self.comments = ''
664 668 #
665 669
666 670
667 671 def getltctime(self):
668 672
669 673 if self.useLocalTime:
670 674 return self.utctime - self.timeZone*60
671 675
672 676 return self.utctime
673 677
674 678 def getDatatime(self):
675 679
676 680 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
677 681 return datatime
678 682
679 683 def getTimeRange(self):
680 684
681 685 datatime = []
682 686
683 687 datatime.append(self.ltctime)
684 688 datatime.append(self.ltctime + self.timeInterval)
685 689
686 690 datatime = numpy.array(datatime)
687 691
688 692 return datatime
689 693
690 694 def getHeiRange(self):
691 695
692 696 heis = self.heightList
693 697
694 698 return heis
695 699
696 700 def isEmpty(self):
697 701
698 702 return self.flagNoData
699 703
700 704 def getNHeights(self):
701 705
702 706 return len(self.heightList)
703 707
704 708 def getNChannels(self):
705 709
706 710 return len(self.channelList)
707 711
708 712 def getChannelIndexList(self):
709 713
710 714 return range(self.nChannels)
711 715
712 716 def getNoise(self, type = 1):
713 717
714 718 self.noise = numpy.zeros(self.nChannels)
715 719
716 720 if type == 1:
717 721 noise = self.getNoisebyHildebrand()
718 722
719 723 if type == 2:
720 724 noise = self.getNoisebySort()
721 725
722 726 if type == 3:
723 727 noise = self.getNoisebyWindow()
724 728
725 729 return noise
726 730
727 731 datatime = property(getDatatime, "I'm the 'datatime' property")
728 732 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
729 733 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
730 734 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
731 735 noise = property(getNoise, "I'm the 'nHeights' property.")
732 736 datatime = property(getDatatime, "I'm the 'datatime' property")
733 737 ltctime = property(getltctime, "I'm the 'ltctime' property")
734 738
735 739 ltctime = property(getltctime, "I'm the 'ltctime' property")
736 740
@@ -1,1346 +1,1354
1 1 '''
2 2 @author: Daniel Suarez
3 3 '''
4 4 import os
5 5 import datetime
6 6 import numpy
7 7
8 8 from figure import Figure, isRealtime
9 9
10 10 class SpectraPlot(Figure):
11 11
12 12 isConfig = None
13 13 __nsubplots = None
14 14
15 15 WIDTHPROF = None
16 16 HEIGHTPROF = None
17 17 PREFIX = 'spc'
18 18
19 19 def __init__(self):
20 20
21 21 self.isConfig = False
22 22 self.__nsubplots = 1
23 23
24 24 self.WIDTH = 280
25 25 self.HEIGHT = 250
26 26 self.WIDTHPROF = 120
27 27 self.HEIGHTPROF = 0
28 28 self.counter_imagwr = 0
29 29
30 30 self.PLOT_CODE = 1
31 31 self.FTP_WEI = None
32 32 self.EXP_CODE = None
33 33 self.SUB_EXP_CODE = None
34 34 self.PLOT_POS = None
35 35
36 36 def getSubplots(self):
37 37
38 38 ncol = int(numpy.sqrt(self.nplots)+0.9)
39 39 nrow = int(self.nplots*1./ncol + 0.9)
40 40
41 41 return nrow, ncol
42 42
43 43 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
44 44
45 45 self.__showprofile = showprofile
46 46 self.nplots = nplots
47 47
48 48 ncolspan = 1
49 49 colspan = 1
50 50 if showprofile:
51 51 ncolspan = 3
52 52 colspan = 2
53 53 self.__nsubplots = 2
54 54
55 55 self.createFigure(id = id,
56 56 wintitle = wintitle,
57 57 widthplot = self.WIDTH + self.WIDTHPROF,
58 58 heightplot = self.HEIGHT + self.HEIGHTPROF,
59 59 show=show)
60 60
61 61 nrow, ncol = self.getSubplots()
62 62
63 63 counter = 0
64 64 for y in range(nrow):
65 65 for x in range(ncol):
66 66
67 67 if counter >= self.nplots:
68 68 break
69 69
70 70 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
71 71
72 72 if showprofile:
73 73 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
74 74
75 75 counter += 1
76 76
77 77 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
78 78 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
79 79 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
80 80 server=None, folder=None, username=None, password=None,
81 81 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
82 82
83 83 """
84 84
85 85 Input:
86 86 dataOut :
87 87 id :
88 88 wintitle :
89 89 channelList :
90 90 showProfile :
91 91 xmin : None,
92 92 xmax : None,
93 93 ymin : None,
94 94 ymax : None,
95 95 zmin : None,
96 96 zmax : None
97 97 """
98 98
99 99 if dataOut.flagNoData:
100 100 return None
101 101
102 102 if realtime:
103 103 if not(isRealtime(utcdatatime = dataOut.utctime)):
104 104 print 'Skipping this plot function'
105 105 return
106 106
107 107 if channelList == None:
108 108 channelIndexList = dataOut.channelIndexList
109 109 else:
110 110 channelIndexList = []
111 111 for channel in channelList:
112 112 if channel not in dataOut.channelList:
113 113 raise ValueError, "Channel %d is not in dataOut.channelList"
114 114 channelIndexList.append(dataOut.channelList.index(channel))
115 115
116 116 factor = dataOut.normFactor
117 117
118 118 x = dataOut.getVelRange(1)
119 119 y = dataOut.getHeiRange()
120 120
121 121 z = dataOut.data_spc[channelIndexList,:,:]/factor
122 122 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
123 123 avg = numpy.average(z, axis=1)
124 124 avg = numpy.nanmean(z, axis=1)
125 125 noise = dataOut.noise/factor
126 126
127 127 zdB = 10*numpy.log10(z)
128 128 avgdB = 10*numpy.log10(avg)
129 129 noisedB = 10*numpy.log10(noise)
130 130
131 131 #thisDatetime = dataOut.datatime
132 132 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
133 133 title = wintitle + " Spectra"
134 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
135 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
136
134 137 xlabel = "Velocity (m/s)"
135 138 ylabel = "Range (Km)"
136 139
137 140 if not self.isConfig:
138 141
139 142 nplots = len(channelIndexList)
140 143
141 144 self.setup(id=id,
142 145 nplots=nplots,
143 146 wintitle=wintitle,
144 147 showprofile=showprofile,
145 148 show=show)
146 149
147 150 if xmin == None: xmin = numpy.nanmin(x)
148 151 if xmax == None: xmax = numpy.nanmax(x)
149 152 if ymin == None: ymin = numpy.nanmin(y)
150 153 if ymax == None: ymax = numpy.nanmax(y)
151 154 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
152 155 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
153 156
154 157 self.FTP_WEI = ftp_wei
155 158 self.EXP_CODE = exp_code
156 159 self.SUB_EXP_CODE = sub_exp_code
157 160 self.PLOT_POS = plot_pos
158 161
159 162 self.isConfig = True
160 163
161 164 self.setWinTitle(title)
162 165
163 166 for i in range(self.nplots):
164 167 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
165 168 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
166 169 axes = self.axesList[i*self.__nsubplots]
167 170 axes.pcolor(x, y, zdB[i,:,:],
168 171 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
169 172 xlabel=xlabel, ylabel=ylabel, title=title,
170 173 ticksize=9, cblabel='')
171 174
172 175 if self.__showprofile:
173 176 axes = self.axesList[i*self.__nsubplots +1]
174 177 axes.pline(avgdB[i], y,
175 178 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
176 179 xlabel='dB', ylabel='', title='',
177 180 ytick_visible=False,
178 181 grid='x')
179 182
180 183 noiseline = numpy.repeat(noisedB[i], len(y))
181 184 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
182 185
183 186 self.draw()
184 187
185 188 if figfile == None:
186 189 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
187 190 figfile = self.getFilename(name = str_datetime)
188
191 name = str_datetime
192 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
193 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
194 figfile = self.getFilename(name)
189 195 if figpath != '':
190 196 self.counter_imagwr += 1
191 197 if (self.counter_imagwr>=wr_period):
192 198 # store png plot to local folder
193 199 self.saveFigure(figpath, figfile)
194 200 # store png plot to FTP server according to RT-Web format
195 201 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
196 202 ftp_filename = os.path.join(figpath, name)
197 203 self.saveFigure(figpath, ftp_filename)
198 204 self.counter_imagwr = 0
199 205
200 206
201 207 class CrossSpectraPlot(Figure):
202 208
203 209 isConfig = None
204 210 __nsubplots = None
205 211
206 212 WIDTH = None
207 213 HEIGHT = None
208 214 WIDTHPROF = None
209 215 HEIGHTPROF = None
210 216 PREFIX = 'cspc'
211 217
212 218 def __init__(self):
213 219
214 220 self.isConfig = False
215 221 self.__nsubplots = 4
216 222 self.counter_imagwr = 0
217 223 self.WIDTH = 250
218 224 self.HEIGHT = 250
219 225 self.WIDTHPROF = 0
220 226 self.HEIGHTPROF = 0
221 227
222 228 self.PLOT_CODE = 1
223 229 self.FTP_WEI = None
224 230 self.EXP_CODE = None
225 231 self.SUB_EXP_CODE = None
226 232 self.PLOT_POS = None
227 233
228 234 def getSubplots(self):
229 235
230 236 ncol = 4
231 237 nrow = self.nplots
232 238
233 239 return nrow, ncol
234 240
235 241 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
236 242
237 243 self.__showprofile = showprofile
238 244 self.nplots = nplots
239 245
240 246 ncolspan = 1
241 247 colspan = 1
242 248
243 249 self.createFigure(id = id,
244 250 wintitle = wintitle,
245 251 widthplot = self.WIDTH + self.WIDTHPROF,
246 252 heightplot = self.HEIGHT + self.HEIGHTPROF,
247 253 show=True)
248 254
249 255 nrow, ncol = self.getSubplots()
250 256
251 257 counter = 0
252 258 for y in range(nrow):
253 259 for x in range(ncol):
254 260 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
255 261
256 262 counter += 1
257 263
258 264 def run(self, dataOut, id, wintitle="", pairsList=None,
259 265 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
260 266 save=False, figpath='', figfile=None, ftp=False, wr_period=1,
261 267 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
262 268 server=None, folder=None, username=None, password=None,
263 269 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
264 270
265 271 """
266 272
267 273 Input:
268 274 dataOut :
269 275 id :
270 276 wintitle :
271 277 channelList :
272 278 showProfile :
273 279 xmin : None,
274 280 xmax : None,
275 281 ymin : None,
276 282 ymax : None,
277 283 zmin : None,
278 284 zmax : None
279 285 """
280 286
281 287 if pairsList == None:
282 288 pairsIndexList = dataOut.pairsIndexList
283 289 else:
284 290 pairsIndexList = []
285 291 for pair in pairsList:
286 292 if pair not in dataOut.pairsList:
287 293 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
288 294 pairsIndexList.append(dataOut.pairsList.index(pair))
289 295
290 296 if pairsIndexList == []:
291 297 return
292 298
293 299 if len(pairsIndexList) > 4:
294 300 pairsIndexList = pairsIndexList[0:4]
295 301 factor = dataOut.normFactor
296 302 x = dataOut.getVelRange(1)
297 303 y = dataOut.getHeiRange()
298 304 z = dataOut.data_spc[:,:,:]/factor
299 305 # z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
300 306 avg = numpy.abs(numpy.average(z, axis=1))
301 307 noise = dataOut.noise()/factor
302 308
303 309 zdB = 10*numpy.log10(z)
304 310 avgdB = 10*numpy.log10(avg)
305 311 noisedB = 10*numpy.log10(noise)
306 312
307 313
308 314 #thisDatetime = dataOut.datatime
309 315 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
310 316 title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
311 317 xlabel = "Velocity (m/s)"
312 318 ylabel = "Range (Km)"
313 319
314 320 if not self.isConfig:
315 321
316 322 nplots = len(pairsIndexList)
317 323
318 324 self.setup(id=id,
319 325 nplots=nplots,
320 326 wintitle=wintitle,
321 327 showprofile=False,
322 328 show=show)
323 329
324 330 if xmin == None: xmin = numpy.nanmin(x)
325 331 if xmax == None: xmax = numpy.nanmax(x)
326 332 if ymin == None: ymin = numpy.nanmin(y)
327 333 if ymax == None: ymax = numpy.nanmax(y)
328 334 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
329 335 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
330 336
331 337 self.FTP_WEI = ftp_wei
332 338 self.EXP_CODE = exp_code
333 339 self.SUB_EXP_CODE = sub_exp_code
334 340 self.PLOT_POS = plot_pos
335 341
336 342 self.isConfig = True
337 343
338 344 self.setWinTitle(title)
339 345
340 346 for i in range(self.nplots):
341 347 pair = dataOut.pairsList[pairsIndexList[i]]
342 348 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
343 349 title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[pair[0]], str_datetime)
344 350 zdB = 10.*numpy.log10(dataOut.data_spc[pair[0],:,:]/factor)
345 351 axes0 = self.axesList[i*self.__nsubplots]
346 352 axes0.pcolor(x, y, zdB,
347 353 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
348 354 xlabel=xlabel, ylabel=ylabel, title=title,
349 355 ticksize=9, colormap=power_cmap, cblabel='')
350 356
351 357 title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[pair[1]], str_datetime)
352 358 zdB = 10.*numpy.log10(dataOut.data_spc[pair[1],:,:]/factor)
353 359 axes0 = self.axesList[i*self.__nsubplots+1]
354 360 axes0.pcolor(x, y, zdB,
355 361 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
356 362 xlabel=xlabel, ylabel=ylabel, title=title,
357 363 ticksize=9, colormap=power_cmap, cblabel='')
358 364
359 365 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:])
360 366 coherence = numpy.abs(coherenceComplex)
361 367 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
362 368 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
363 369
364 370 title = "Coherence %d%d" %(pair[0], pair[1])
365 371 axes0 = self.axesList[i*self.__nsubplots+2]
366 372 axes0.pcolor(x, y, coherence,
367 373 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=1,
368 374 xlabel=xlabel, ylabel=ylabel, title=title,
369 375 ticksize=9, colormap=coherence_cmap, cblabel='')
370 376
371 377 title = "Phase %d%d" %(pair[0], pair[1])
372 378 axes0 = self.axesList[i*self.__nsubplots+3]
373 379 axes0.pcolor(x, y, phase,
374 380 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
375 381 xlabel=xlabel, ylabel=ylabel, title=title,
376 382 ticksize=9, colormap=phase_cmap, cblabel='')
377 383
378 384
379 385
380 386 self.draw()
381 387
382 388 if figfile == None:
383 389 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
384 390 figfile = self.getFilename(name = str_datetime)
385 391
386 392 if figpath != '':
387 393 self.counter_imagwr += 1
388 394 if (self.counter_imagwr>=wr_period):
389 395 # store png plot to local folder
390 396 self.saveFigure(figpath, figfile)
391 397 # store png plot to FTP server according to RT-Web format
392 398 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
393 399 ftp_filename = os.path.join(figpath, name)
394 400 self.saveFigure(figpath, ftp_filename)
395 401 self.counter_imagwr = 0
396 402
397 403
398 404 class RTIPlot(Figure):
399 405
400 406 isConfig = None
401 407 __nsubplots = None
402 408
403 409 WIDTHPROF = None
404 410 HEIGHTPROF = None
405 411 PREFIX = 'rti'
406 412
407 413 def __init__(self):
408 414
409 415 self.timerange = 2*60*60
410 416 self.isConfig = False
411 417 self.__nsubplots = 1
412 418
413 419 self.WIDTH = 800
414 420 self.HEIGHT = 150
415 421 self.WIDTHPROF = 120
416 422 self.HEIGHTPROF = 0
417 423 self.counter_imagwr = 0
418 424
419 425 self.PLOT_CODE = 0
420 426 self.FTP_WEI = None
421 427 self.EXP_CODE = None
422 428 self.SUB_EXP_CODE = None
423 429 self.PLOT_POS = None
424 430 self.tmin = None
425 431 self.tmax = None
426 432
427 433 self.xmin = None
428 434 self.xmax = None
429 435
430 436 self.figfile = None
431 437
432 438 def getSubplots(self):
433 439
434 440 ncol = 1
435 441 nrow = self.nplots
436 442
437 443 return nrow, ncol
438 444
439 445 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
440 446
441 447 self.__showprofile = showprofile
442 448 self.nplots = nplots
443 449
444 450 ncolspan = 1
445 451 colspan = 1
446 452 if showprofile:
447 453 ncolspan = 7
448 454 colspan = 6
449 455 self.__nsubplots = 2
450 456
451 457 self.createFigure(id = id,
452 458 wintitle = wintitle,
453 459 widthplot = self.WIDTH + self.WIDTHPROF,
454 460 heightplot = self.HEIGHT + self.HEIGHTPROF,
455 461 show=show)
456 462
457 463 nrow, ncol = self.getSubplots()
458 464
459 465 counter = 0
460 466 for y in range(nrow):
461 467 for x in range(ncol):
462 468
463 469 if counter >= self.nplots:
464 470 break
465 471
466 472 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
467 473
468 474 if showprofile:
469 475 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
470 476
471 477 counter += 1
472 478
473 479 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
474 480 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
475 481 timerange=None,
476 482 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
477 483 server=None, folder=None, username=None, password=None,
478 484 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
479 485
480 486 """
481 487
482 488 Input:
483 489 dataOut :
484 490 id :
485 491 wintitle :
486 492 channelList :
487 493 showProfile :
488 494 xmin : None,
489 495 xmax : None,
490 496 ymin : None,
491 497 ymax : None,
492 498 zmin : None,
493 499 zmax : None
494 500 """
495 501
496 502 if channelList == None:
497 503 channelIndexList = dataOut.channelIndexList
498 504 else:
499 505 channelIndexList = []
500 506 for channel in channelList:
501 507 if channel not in dataOut.channelList:
502 508 raise ValueError, "Channel %d is not in dataOut.channelList"
503 509 channelIndexList.append(dataOut.channelList.index(channel))
504 510
505 511 if timerange != None:
506 512 self.timerange = timerange
507 513
508 514 #tmin = None
509 515 #tmax = None
510 516 factor = dataOut.normFactor
511 517 x = dataOut.getTimeRange()
512 518 y = dataOut.getHeiRange()
513 519
514 520 z = dataOut.data_spc[channelIndexList,:,:]/factor
515 521 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
516 522 avg = numpy.average(z, axis=1)
517 523
518 524 avgdB = 10.*numpy.log10(avg)
519 525
520 526
521 527 # thisDatetime = dataOut.datatime
522 528 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
523 529 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
524 530 xlabel = ""
525 531 ylabel = "Range (Km)"
526 532
527 533 if not self.isConfig:
528 534
529 535 nplots = len(channelIndexList)
530 536
531 537 self.setup(id=id,
532 538 nplots=nplots,
533 539 wintitle=wintitle,
534 540 showprofile=showprofile,
535 541 show=show)
536 542
537 543 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
538 544
539 545 # if timerange != None:
540 546 # self.timerange = timerange
541 547 # self.xmin, self.tmax = self.getTimeLim(x, xmin, xmax, timerange)
542 548
543 549
544 550
545 551 if ymin == None: ymin = numpy.nanmin(y)
546 552 if ymax == None: ymax = numpy.nanmax(y)
547 553 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
548 554 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
549 555
550 556 self.FTP_WEI = ftp_wei
551 557 self.EXP_CODE = exp_code
552 558 self.SUB_EXP_CODE = sub_exp_code
553 559 self.PLOT_POS = plot_pos
554 560
555 561 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
556 562 self.isConfig = True
557 563 self.figfile = figfile
558 564
559 565 self.setWinTitle(title)
560 566
561 567 if ((self.xmax - x[1]) < (x[1]-x[0])):
562 568 x[1] = self.xmax
563 569
564 570 for i in range(self.nplots):
565 571 title = "Channel %d: %s" %(dataOut.channelList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
572 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
573 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
566 574 axes = self.axesList[i*self.__nsubplots]
567 575 zdB = avgdB[i].reshape((1,-1))
568 576 axes.pcolorbuffer(x, y, zdB,
569 577 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
570 578 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
571 579 ticksize=9, cblabel='', cbsize="1%")
572 580
573 581 if self.__showprofile:
574 582 axes = self.axesList[i*self.__nsubplots +1]
575 583 axes.pline(avgdB[i], y,
576 584 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
577 585 xlabel='dB', ylabel='', title='',
578 586 ytick_visible=False,
579 587 grid='x')
580 588
581 589 self.draw()
582 590
583 591 if x[1] >= self.axesList[0].xmax:
584 592 self.counter_imagwr = wr_period
585 593 self.__isConfig = False
586 594
587 595
588 596 if self.figfile == None:
589 597 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
590 598 self.figfile = self.getFilename(name = str_datetime)
591 599
592 600 if figpath != '':
593 601
594 602 self.counter_imagwr += 1
595 603 if (self.counter_imagwr>=wr_period):
596 604 # store png plot to local folder
597 605 self.saveFigure(figpath, self.figfile)
598 606 # store png plot to FTP server according to RT-Web format
599 607 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
600 608 ftp_filename = os.path.join(figpath, name)
601 609 self.saveFigure(figpath, ftp_filename)
602 610
603 611 self.counter_imagwr = 0
604 612
605 613
606 614 class CoherenceMap(Figure):
607 615 isConfig = None
608 616 __nsubplots = None
609 617
610 618 WIDTHPROF = None
611 619 HEIGHTPROF = None
612 620 PREFIX = 'cmap'
613 621
614 622 def __init__(self):
615 623 self.timerange = 2*60*60
616 624 self.isConfig = False
617 625 self.__nsubplots = 1
618 626
619 627 self.WIDTH = 800
620 628 self.HEIGHT = 150
621 629 self.WIDTHPROF = 120
622 630 self.HEIGHTPROF = 0
623 631 self.counter_imagwr = 0
624 632
625 633 self.PLOT_CODE = 3
626 634 self.FTP_WEI = None
627 635 self.EXP_CODE = None
628 636 self.SUB_EXP_CODE = None
629 637 self.PLOT_POS = None
630 638 self.counter_imagwr = 0
631 639
632 640 self.xmin = None
633 641 self.xmax = None
634 642
635 643 def getSubplots(self):
636 644 ncol = 1
637 645 nrow = self.nplots*2
638 646
639 647 return nrow, ncol
640 648
641 649 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
642 650 self.__showprofile = showprofile
643 651 self.nplots = nplots
644 652
645 653 ncolspan = 1
646 654 colspan = 1
647 655 if showprofile:
648 656 ncolspan = 7
649 657 colspan = 6
650 658 self.__nsubplots = 2
651 659
652 660 self.createFigure(id = id,
653 661 wintitle = wintitle,
654 662 widthplot = self.WIDTH + self.WIDTHPROF,
655 663 heightplot = self.HEIGHT + self.HEIGHTPROF,
656 664 show=True)
657 665
658 666 nrow, ncol = self.getSubplots()
659 667
660 668 for y in range(nrow):
661 669 for x in range(ncol):
662 670
663 671 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
664 672
665 673 if showprofile:
666 674 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
667 675
668 676 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
669 677 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
670 678 timerange=None,
671 679 save=False, figpath='', figfile=None, ftp=False, wr_period=1,
672 680 coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
673 681 server=None, folder=None, username=None, password=None,
674 682 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
675 683
676 684 if pairsList == None:
677 685 pairsIndexList = dataOut.pairsIndexList
678 686 else:
679 687 pairsIndexList = []
680 688 for pair in pairsList:
681 689 if pair not in dataOut.pairsList:
682 690 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
683 691 pairsIndexList.append(dataOut.pairsList.index(pair))
684 692
685 693 if timerange != None:
686 694 self.timerange = timerange
687 695
688 696 if pairsIndexList == []:
689 697 return
690 698
691 699 if len(pairsIndexList) > 4:
692 700 pairsIndexList = pairsIndexList[0:4]
693 701
694 702 # tmin = None
695 703 # tmax = None
696 704 x = dataOut.getTimeRange()
697 705 y = dataOut.getHeiRange()
698 706
699 707 #thisDatetime = dataOut.datatime
700 708 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
701 709 title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
702 710 xlabel = ""
703 711 ylabel = "Range (Km)"
704 712
705 713 if not self.isConfig:
706 714 nplots = len(pairsIndexList)
707 715 self.setup(id=id,
708 716 nplots=nplots,
709 717 wintitle=wintitle,
710 718 showprofile=showprofile,
711 719 show=show)
712 720
713 721 #tmin, tmax = self.getTimeLim(x, xmin, xmax)
714 722
715 723 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
716 724
717 725 if ymin == None: ymin = numpy.nanmin(y)
718 726 if ymax == None: ymax = numpy.nanmax(y)
719 727 if zmin == None: zmin = 0.
720 728 if zmax == None: zmax = 1.
721 729
722 730 self.FTP_WEI = ftp_wei
723 731 self.EXP_CODE = exp_code
724 732 self.SUB_EXP_CODE = sub_exp_code
725 733 self.PLOT_POS = plot_pos
726 734
727 735 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
728 736
729 737 self.isConfig = True
730 738
731 739 self.setWinTitle(title)
732 740
733 741 if ((self.xmax - x[1]) < (x[1]-x[0])):
734 742 x[1] = self.xmax
735 743
736 744 for i in range(self.nplots):
737 745
738 746 pair = dataOut.pairsList[pairsIndexList[i]]
739 747
740 748 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
741 749 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
742 750 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
743 751
744 752
745 753 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
746 754 coherence = numpy.abs(avgcoherenceComplex)
747 755
748 756 z = coherence.reshape((1,-1))
749 757
750 758 counter = 0
751 759
752 760 title = "Coherence %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
753 761 axes = self.axesList[i*self.__nsubplots*2]
754 762 axes.pcolorbuffer(x, y, z,
755 763 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
756 764 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
757 765 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
758 766
759 767 if self.__showprofile:
760 768 counter += 1
761 769 axes = self.axesList[i*self.__nsubplots*2 + counter]
762 770 axes.pline(coherence, y,
763 771 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
764 772 xlabel='', ylabel='', title='', ticksize=7,
765 773 ytick_visible=False, nxticks=5,
766 774 grid='x')
767 775
768 776 counter += 1
769 777
770 778 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
771 779
772 780 z = phase.reshape((1,-1))
773 781
774 782 title = "Phase %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
775 783 axes = self.axesList[i*self.__nsubplots*2 + counter]
776 784 axes.pcolorbuffer(x, y, z,
777 785 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
778 786 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
779 787 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
780 788
781 789 if self.__showprofile:
782 790 counter += 1
783 791 axes = self.axesList[i*self.__nsubplots*2 + counter]
784 792 axes.pline(phase, y,
785 793 xmin=-180, xmax=180, ymin=ymin, ymax=ymax,
786 794 xlabel='', ylabel='', title='', ticksize=7,
787 795 ytick_visible=False, nxticks=4,
788 796 grid='x')
789 797
790 798 self.draw()
791 799
792 800 if x[1] >= self.axesList[0].xmax:
793 801 self.counter_imagwr = wr_period
794 802 self.__isConfig = False
795 803
796 804 if figfile == None:
797 805 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
798 806 figfile = self.getFilename(name = str_datetime)
799 807
800 808 if figpath != '':
801 809
802 810 self.counter_imagwr += 1
803 811 if (self.counter_imagwr>=wr_period):
804 812 # store png plot to local folder
805 813 self.saveFigure(figpath, figfile)
806 814 # store png plot to FTP server according to RT-Web format
807 815 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
808 816 ftp_filename = os.path.join(figpath, name)
809 817 self.saveFigure(figpath, ftp_filename)
810 818
811 819 self.counter_imagwr = 0
812 820
813 821 class PowerProfile(Figure):
814 822 isConfig = None
815 823 __nsubplots = None
816 824
817 825 WIDTHPROF = None
818 826 HEIGHTPROF = None
819 827 PREFIX = 'spcprofile'
820 828
821 829 def __init__(self):
822 830 self.isConfig = False
823 831 self.__nsubplots = 1
824 832
825 833 self.WIDTH = 300
826 834 self.HEIGHT = 500
827 835 self.counter_imagwr = 0
828 836
829 837 def getSubplots(self):
830 838 ncol = 1
831 839 nrow = 1
832 840
833 841 return nrow, ncol
834 842
835 843 def setup(self, id, nplots, wintitle, show):
836 844
837 845 self.nplots = nplots
838 846
839 847 ncolspan = 1
840 848 colspan = 1
841 849
842 850 self.createFigure(id = id,
843 851 wintitle = wintitle,
844 852 widthplot = self.WIDTH,
845 853 heightplot = self.HEIGHT,
846 854 show=show)
847 855
848 856 nrow, ncol = self.getSubplots()
849 857
850 858 counter = 0
851 859 for y in range(nrow):
852 860 for x in range(ncol):
853 861 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
854 862
855 863 def run(self, dataOut, id, wintitle="", channelList=None,
856 864 xmin=None, xmax=None, ymin=None, ymax=None,
857 865 save=False, figpath='', figfile=None, show=True, wr_period=1,
858 866 server=None, folder=None, username=None, password=None,):
859 867
860 868 if dataOut.flagNoData:
861 869 return None
862 870
863 871 if channelList == None:
864 872 channelIndexList = dataOut.channelIndexList
865 873 channelList = dataOut.channelList
866 874 else:
867 875 channelIndexList = []
868 876 for channel in channelList:
869 877 if channel not in dataOut.channelList:
870 878 raise ValueError, "Channel %d is not in dataOut.channelList"
871 879 channelIndexList.append(dataOut.channelList.index(channel))
872 880
873 881 try:
874 882 factor = dataOut.normFactor
875 883 except:
876 884 factor = 1
877 885
878 886 y = dataOut.getHeiRange()
879 887
880 888 #for voltage
881 889 if dataOut.type == 'Voltage':
882 890 x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
883 891 x = x.real
884 892 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
885 893
886 894 #for spectra
887 895 if dataOut.type == 'Spectra':
888 896 x = dataOut.data_spc[channelIndexList,:,:]/factor
889 897 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
890 898 x = numpy.average(x, axis=1)
891 899
892 900
893 901 xdB = 10*numpy.log10(x)
894 902
895 903 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
896 904 title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
897 905 xlabel = "dB"
898 906 ylabel = "Range (Km)"
899 907
900 908 if not self.isConfig:
901 909
902 910 nplots = 1
903 911
904 912 self.setup(id=id,
905 913 nplots=nplots,
906 914 wintitle=wintitle,
907 915 show=show)
908 916
909 917 if ymin == None: ymin = numpy.nanmin(y)
910 918 if ymax == None: ymax = numpy.nanmax(y)
911 919 if xmin == None: xmin = numpy.nanmin(xdB)*0.9
912 920 if xmax == None: xmax = numpy.nanmax(xdB)*0.9
913 921
914 922 self.__isConfig = True
915 923
916 924 self.setWinTitle(title)
917 925
918 926 title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
919 927 axes = self.axesList[0]
920 928
921 929 legendlabels = ["channel %d"%x for x in channelList]
922 930 axes.pmultiline(xdB, y,
923 931 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
924 932 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
925 933 ytick_visible=True, nxticks=5,
926 934 grid='x')
927 935
928 936 self.draw()
929 937
930 938 if figfile == None:
931 939 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
932 940 figfile = self.getFilename(name = str_datetime)
933 941
934 942 if figpath != '':
935 943 self.counter_imagwr += 1
936 944 if (self.counter_imagwr>=wr_period):
937 945 # store png plot to local folder
938 946 self.saveFigure(figpath, figfile)
939 947 # store png plot to FTP server according to RT-Web format
940 948 #name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
941 949 #ftp_filename = os.path.join(figpath, name)
942 950 #self.saveFigure(figpath, ftp_filename)
943 951 self.counter_imagwr = 0
944 952
945 953
946 954
947 955 class Noise(Figure):
948 956
949 957 isConfig = None
950 958 __nsubplots = None
951 959
952 960 PREFIX = 'noise'
953 961
954 962 def __init__(self):
955 963
956 964 self.timerange = 24*60*60
957 965 self.isConfig = False
958 966 self.__nsubplots = 1
959 967 self.counter_imagwr = 0
960 968 self.WIDTH = 600
961 969 self.HEIGHT = 300
962 970 self.WIDTHPROF = 120
963 971 self.HEIGHTPROF = 0
964 972 self.xdata = None
965 973 self.ydata = None
966 974
967 975 self.PLOT_CODE = 17
968 976 self.FTP_WEI = None
969 977 self.EXP_CODE = None
970 978 self.SUB_EXP_CODE = None
971 979 self.PLOT_POS = None
972 980 self.figfile = None
973 981
974 982 def getSubplots(self):
975 983
976 984 ncol = 1
977 985 nrow = 1
978 986
979 987 return nrow, ncol
980 988
981 989 def openfile(self, filename):
982 990 f = open(filename,'w+')
983 991 f.write('\n\n')
984 992 f.write('JICAMARCA RADIO OBSERVATORY - Noise \n')
985 993 f.write('DD MM YYYY HH MM SS Channel0 Channel1 Channel2 Channel3\n\n' )
986 994 f.close()
987 995
988 996 def save_data(self, filename_phase, data, data_datetime):
989 997 f=open(filename_phase,'a')
990 998 timetuple_data = data_datetime.timetuple()
991 999 day = str(timetuple_data.tm_mday)
992 1000 month = str(timetuple_data.tm_mon)
993 1001 year = str(timetuple_data.tm_year)
994 1002 hour = str(timetuple_data.tm_hour)
995 1003 minute = str(timetuple_data.tm_min)
996 1004 second = str(timetuple_data.tm_sec)
997 1005 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
998 1006 f.close()
999 1007
1000 1008
1001 1009 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1002 1010
1003 1011 self.__showprofile = showprofile
1004 1012 self.nplots = nplots
1005 1013
1006 1014 ncolspan = 7
1007 1015 colspan = 6
1008 1016 self.__nsubplots = 2
1009 1017
1010 1018 self.createFigure(id = id,
1011 1019 wintitle = wintitle,
1012 1020 widthplot = self.WIDTH+self.WIDTHPROF,
1013 1021 heightplot = self.HEIGHT+self.HEIGHTPROF,
1014 1022 show=show)
1015 1023
1016 1024 nrow, ncol = self.getSubplots()
1017 1025
1018 1026 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1019 1027
1020 1028
1021 1029 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1022 1030 xmin=None, xmax=None, ymin=None, ymax=None,
1023 1031 timerange=None,
1024 1032 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
1025 1033 server=None, folder=None, username=None, password=None,
1026 1034 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1027 1035
1028 1036 if channelList == None:
1029 1037 channelIndexList = dataOut.channelIndexList
1030 1038 channelList = dataOut.channelList
1031 1039 else:
1032 1040 channelIndexList = []
1033 1041 for channel in channelList:
1034 1042 if channel not in dataOut.channelList:
1035 1043 raise ValueError, "Channel %d is not in dataOut.channelList"
1036 1044 channelIndexList.append(dataOut.channelList.index(channel))
1037 1045
1038 1046 if timerange != None:
1039 1047 self.timerange = timerange
1040 1048
1041 1049 tmin = None
1042 1050 tmax = None
1043 1051 x = dataOut.getTimeRange()
1044 1052 y = dataOut.getHeiRange()
1045 1053 factor = dataOut.normFactor
1046 1054 noise = dataOut.noise()/factor
1047 1055 noisedB = 10*numpy.log10(noise)
1048 1056
1049 1057 #thisDatetime = dataOut.datatime
1050 1058 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1051 1059 title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1052 1060 xlabel = ""
1053 1061 ylabel = "Intensity (dB)"
1054 1062
1055 1063 if not self.isConfig:
1056 1064
1057 1065 nplots = 1
1058 1066
1059 1067 self.setup(id=id,
1060 1068 nplots=nplots,
1061 1069 wintitle=wintitle,
1062 1070 showprofile=showprofile,
1063 1071 show=show)
1064 1072
1065 1073 tmin, tmax = self.getTimeLim(x, xmin, xmax)
1066 1074 if ymin == None: ymin = numpy.nanmin(noisedB) - 10.0
1067 1075 if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
1068 1076
1069 1077 self.FTP_WEI = ftp_wei
1070 1078 self.EXP_CODE = exp_code
1071 1079 self.SUB_EXP_CODE = sub_exp_code
1072 1080 self.PLOT_POS = plot_pos
1073 1081
1074 1082
1075 1083 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1076 1084 self.isConfig = True
1077 1085 self.figfile = figfile
1078 1086 self.xdata = numpy.array([])
1079 1087 self.ydata = numpy.array([])
1080 1088
1081 1089 #open file beacon phase
1082 1090 path = '%s%03d' %(self.PREFIX, self.id)
1083 1091 noise_file = os.path.join(path,'%s.txt'%self.name)
1084 1092 self.filename_noise = os.path.join(figpath,noise_file)
1085 1093 self.openfile(self.filename_noise)
1086 1094
1087 1095
1088 1096 #store data beacon phase
1089 1097 self.save_data(self.filename_noise, noisedB, thisDatetime)
1090 1098
1091 1099
1092 1100 self.setWinTitle(title)
1093 1101
1094 1102
1095 1103 title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1096 1104
1097 1105 legendlabels = ["channel %d"%(idchannel+1) for idchannel in channelList]
1098 1106 axes = self.axesList[0]
1099 1107
1100 1108 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1101 1109
1102 1110 if len(self.ydata)==0:
1103 1111 self.ydata = noisedB[channelIndexList].reshape(-1,1)
1104 1112 else:
1105 1113 self.ydata = numpy.hstack((self.ydata, noisedB[channelIndexList].reshape(-1,1)))
1106 1114
1107 1115
1108 1116 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1109 1117 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax,
1110 1118 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1111 1119 XAxisAsTime=True, grid='both'
1112 1120 )
1113 1121
1114 1122 self.draw()
1115 1123
1116 1124 if x[1] >= self.axesList[0].xmax:
1117 1125 self.counter_imagwr = wr_period
1118 1126 del self.xdata
1119 1127 del self.ydata
1120 1128 self.__isConfig = False
1121 1129
1122 1130 if self.figfile == None:
1123 1131 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
1124 1132 self.figfile = self.getFilename(name = str_datetime)
1125 1133
1126 1134 if figpath != '':
1127 1135 self.counter_imagwr += 1
1128 1136 if (self.counter_imagwr>=wr_period):
1129 1137 # store png plot to local folder
1130 1138 self.saveFigure(figpath, self.figfile)
1131 1139 # store png plot to FTP server according to RT-Web format
1132 1140 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
1133 1141 ftp_filename = os.path.join(figpath, name)
1134 1142 self.saveFigure(figpath, ftp_filename)
1135 1143 self.counter_imagwr = 0
1136 1144
1137 1145
1138 1146 class BeaconPhase(Figure):
1139 1147
1140 1148 __isConfig = None
1141 1149 __nsubplots = None
1142 1150
1143 1151 PREFIX = 'beacon_phase'
1144 1152
1145 1153 def __init__(self):
1146 1154
1147 1155 self.timerange = 24*60*60
1148 1156 self.__isConfig = False
1149 1157 self.__nsubplots = 1
1150 1158 self.counter_imagwr = 0
1151 1159 self.WIDTH = 600
1152 1160 self.HEIGHT = 300
1153 1161 self.WIDTHPROF = 120
1154 1162 self.HEIGHTPROF = 0
1155 1163 self.xdata = None
1156 1164 self.ydata = None
1157 1165
1158 1166 self.PLOT_CODE = 18
1159 1167 self.FTP_WEI = None
1160 1168 self.EXP_CODE = None
1161 1169 self.SUB_EXP_CODE = None
1162 1170 self.PLOT_POS = None
1163 1171
1164 1172 self.filename_phase = None
1165 1173
1166 1174 self.figfile = None
1167 1175
1168 1176 def getSubplots(self):
1169 1177
1170 1178 ncol = 1
1171 1179 nrow = 1
1172 1180
1173 1181 return nrow, ncol
1174 1182
1175 1183 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1176 1184
1177 1185 self.__showprofile = showprofile
1178 1186 self.nplots = nplots
1179 1187
1180 1188 ncolspan = 7
1181 1189 colspan = 6
1182 1190 self.__nsubplots = 2
1183 1191
1184 1192 self.createFigure(id = id,
1185 1193 wintitle = wintitle,
1186 1194 widthplot = self.WIDTH+self.WIDTHPROF,
1187 1195 heightplot = self.HEIGHT+self.HEIGHTPROF,
1188 1196 show=show)
1189 1197
1190 1198 nrow, ncol = self.getSubplots()
1191 1199
1192 1200 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1193 1201
1194 1202 def save_phase(self, filename_phase):
1195 1203 f = open(filename_phase,'w+')
1196 1204 f.write('\n\n')
1197 1205 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1198 1206 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1199 1207 f.close()
1200 1208
1201 1209 def save_data(self, filename_phase, data, data_datetime):
1202 1210 f=open(filename_phase,'a')
1203 1211 timetuple_data = data_datetime.timetuple()
1204 1212 day = str(timetuple_data.tm_mday)
1205 1213 month = str(timetuple_data.tm_mon)
1206 1214 year = str(timetuple_data.tm_year)
1207 1215 hour = str(timetuple_data.tm_hour)
1208 1216 minute = str(timetuple_data.tm_min)
1209 1217 second = str(timetuple_data.tm_sec)
1210 1218 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1211 1219 f.close()
1212 1220
1213 1221
1214 1222 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1215 1223 xmin=None, xmax=None, ymin=None, ymax=None,
1216 1224 timerange=None,
1217 1225 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
1218 1226 server=None, folder=None, username=None, password=None,
1219 1227 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1220 1228
1221 1229 if pairsList == None:
1222 1230 pairsIndexList = dataOut.pairsIndexList
1223 1231 else:
1224 1232 pairsIndexList = []
1225 1233 for pair in pairsList:
1226 1234 if pair not in dataOut.pairsList:
1227 1235 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
1228 1236 pairsIndexList.append(dataOut.pairsList.index(pair))
1229 1237
1230 1238 if pairsIndexList == []:
1231 1239 return
1232 1240
1233 1241 # if len(pairsIndexList) > 4:
1234 1242 # pairsIndexList = pairsIndexList[0:4]
1235 1243
1236 1244 if timerange != None:
1237 1245 self.timerange = timerange
1238 1246
1239 1247 tmin = None
1240 1248 tmax = None
1241 1249 x = dataOut.getTimeRange()
1242 1250 y = dataOut.getHeiRange()
1243 1251
1244 1252
1245 1253 #thisDatetime = dataOut.datatime
1246 1254 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1247 1255 title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1248 1256 xlabel = "Local Time"
1249 1257 ylabel = "Phase"
1250 1258
1251 1259 nplots = len(pairsIndexList)
1252 1260 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1253 1261 phase_beacon = numpy.zeros(len(pairsIndexList))
1254 1262 for i in range(nplots):
1255 1263 pair = dataOut.pairsList[pairsIndexList[i]]
1256 1264 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
1257 1265 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
1258 1266 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
1259 1267 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1260 1268 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1261 1269
1262 1270 #print "Phase %d%d" %(pair[0], pair[1])
1263 1271 #print phase[dataOut.beacon_heiIndexList]
1264 1272
1265 1273 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1266 1274
1267 1275 if not self.__isConfig:
1268 1276
1269 1277 nplots = len(pairsIndexList)
1270 1278
1271 1279 self.setup(id=id,
1272 1280 nplots=nplots,
1273 1281 wintitle=wintitle,
1274 1282 showprofile=showprofile,
1275 1283 show=show)
1276 1284
1277 1285 tmin, tmax = self.getTimeLim(x, xmin, xmax)
1278 1286 if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0
1279 1287 if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0
1280 1288
1281 1289 self.FTP_WEI = ftp_wei
1282 1290 self.EXP_CODE = exp_code
1283 1291 self.SUB_EXP_CODE = sub_exp_code
1284 1292 self.PLOT_POS = plot_pos
1285 1293
1286 1294 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1287 1295 self.__isConfig = True
1288 1296 self.figfile = figfile
1289 1297 self.xdata = numpy.array([])
1290 1298 self.ydata = numpy.array([])
1291 1299
1292 1300 #open file beacon phase
1293 1301 path = '%s%03d' %(self.PREFIX, self.id)
1294 1302 beacon_file = os.path.join(path,'%s.txt'%self.name)
1295 1303 self.filename_phase = os.path.join(figpath,beacon_file)
1296 1304 #self.save_phase(self.filename_phase)
1297 1305
1298 1306
1299 1307 #store data beacon phase
1300 1308 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1301 1309
1302 1310 self.setWinTitle(title)
1303 1311
1304 1312
1305 1313 title = "Beacon Signal %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1306 1314
1307 1315 legendlabels = ["pairs %d%d"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1308 1316
1309 1317 axes = self.axesList[0]
1310 1318
1311 1319 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1312 1320
1313 1321 if len(self.ydata)==0:
1314 1322 self.ydata = phase_beacon.reshape(-1,1)
1315 1323 else:
1316 1324 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1317 1325
1318 1326
1319 1327 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1320 1328 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax,
1321 1329 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1322 1330 XAxisAsTime=True, grid='both'
1323 1331 )
1324 1332
1325 1333 self.draw()
1326 1334
1327 1335 if x[1] >= self.axesList[0].xmax:
1328 1336 self.counter_imagwr = wr_period
1329 1337 del self.xdata
1330 1338 del self.ydata
1331 1339 self.__isConfig = False
1332 1340
1333 1341 if self.figfile == None:
1334 1342 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
1335 1343 self.figfile = self.getFilename(name = str_datetime)
1336 1344
1337 1345 if figpath != '':
1338 1346 self.counter_imagwr += 1
1339 1347 if (self.counter_imagwr>=wr_period):
1340 1348 # store png plot to local folder
1341 1349 self.saveFigure(figpath, self.figfile)
1342 1350 # store png plot to FTP server according to RT-Web format
1343 1351 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
1344 1352 ftp_filename = os.path.join(figpath, name)
1345 1353 self.saveFigure(figpath, ftp_filename)
1346 1354 self.counter_imagwr = 0
@@ -1,573 +1,616
1 1 '''
2 2 @author: Daniel Suarez
3 3 '''
4 4
5 5 import os
6 6 import sys
7 7 import glob
8 8 import fnmatch
9 9 import datetime
10 10 import time
11 11 import re
12 12 import h5py
13 13 import numpy
14 14
15 15 from model.proc.jroproc_base import ProcessingUnit, Operation
16 16 from model.data.jroamisr import AMISR
17 17
18 18 class RadacHeader():
19 19 def __init__(self, fp):
20 20 header = 'Raw11/Data/RadacHeader'
21 21 self.beamCodeByPulse = fp.get(header+'/BeamCode')
22 22 self.beamCode = fp.get('Raw11/Data/Beamcodes')
23 23 self.code = fp.get(header+'/Code')
24 24 self.frameCount = fp.get(header+'/FrameCount')
25 25 self.modeGroup = fp.get(header+'/ModeGroup')
26 26 self.nsamplesPulse = fp.get(header+'/NSamplesPulse')
27 27 self.pulseCount = fp.get(header+'/PulseCount')
28 28 self.radacTime = fp.get(header+'/RadacTime')
29 29 self.timeCount = fp.get(header+'/TimeCount')
30 30 self.timeStatus = fp.get(header+'/TimeStatus')
31 31
32 32 self.nrecords = self.pulseCount.shape[0] #nblocks
33 33 self.npulses = self.pulseCount.shape[1] #nprofile
34 34 self.nsamples = self.nsamplesPulse[0,0] #ngates
35 35 self.nbeams = self.beamCode.shape[1]
36 36
37 37
38 38 def getIndexRangeToPulse(self, idrecord=0):
39 39 #indexToZero = numpy.where(self.pulseCount.value[idrecord,:]==0)
40 40 #startPulseCountId = indexToZero[0][0]
41 41 #endPulseCountId = startPulseCountId - 1
42 42 #range1 = numpy.arange(startPulseCountId,self.npulses,1)
43 43 #range2 = numpy.arange(0,startPulseCountId,1)
44 44 #return range1, range2
45 45
46 46 looking_zeros_index = numpy.where(self.pulseCount.value[idrecord,:]==0)[0]
47 47 getLastIndexZero = looking_zeros_index[-1]
48 48 index_data = numpy.arange(0,getLastIndexZero,1)
49 49 index_buffer = numpy.arange(getLastIndexZero,self.npulses,1)
50 50 return index_data, index_buffer
51 51
52 52 class AMISRReader(ProcessingUnit):
53 53
54 54 path = None
55 55 startDate = None
56 56 endDate = None
57 57 startTime = None
58 58 endTime = None
59 59 walk = None
60 60 isConfig = False
61 61
62 62 def __init__(self):
63 63 self.set = None
64 64 self.subset = None
65 65 self.extension_file = '.h5'
66 66 self.dtc_str = 'dtc'
67 67 self.dtc_id = 0
68 68 self.status = True
69 69 self.isConfig = False
70 70 self.dirnameList = []
71 71 self.filenameList = []
72 72 self.fileIndex = None
73 73 self.flagNoMoreFiles = False
74 74 self.flagIsNewFile = 0
75 75 self.filename = ''
76 76 self.amisrFilePointer = None
77 77 self.radacHeaderObj = None
78 78 self.dataOut = self.__createObjByDefault()
79 79 self.datablock = None
80 80 self.rest_datablock = None
81 81 self.range = None
82 82 self.idrecord_count = 0
83 83 self.profileIndex = 0
84 84 self.index_amisr_sample = None
85 85 self.index_amisr_buffer = None
86 86 self.beamCodeByFrame = None
87 87 self.radacTimeByFrame = None
88 88 #atributos originales tal y como esta en el archivo de datos
89 89 self.beamCodesFromFile = None
90 90 self.radacTimeFromFile = None
91 91 self.rangeFromFile = None
92 92 self.dataByFrame = None
93 93 self.dataset = None
94 94
95 95 self.beamCodeDict = {}
96 96 self.beamRangeDict = {}
97 97
98 98 #experiment cgf file
99 99 self.npulsesint_fromfile = None
100 100 self.recordsperfile_fromfile = None
101 101 self.nbeamcodes_fromfile = None
102 102 self.ngates_fromfile = None
103 103 self.ippSeconds_fromfile = None
104 104 self.frequency_h5file = None
105 105
106 106
107 107 self.__firstFile = True
108 108 self.buffer_radactime = None
109 109
110 110 self.index4_schain_datablock = None
111 111 self.index4_buffer = None
112 112 self.schain_datablock = None
113 113 self.buffer = None
114 114 self.linear_pulseCount = None
115 115 self.npulseByFrame = None
116 116 self.profileIndex_offset = None
117 117 self.timezone = 'ut'
118 118
119 119 def __createObjByDefault(self):
120 120
121 121 dataObj = AMISR()
122 122
123 123 return dataObj
124 124
125 125 def __setParameters(self,path,startDate,endDate,startTime,endTime,walk):
126 126 self.path = path
127 127 self.startDate = startDate
128 128 self.endDate = endDate
129 129 self.startTime = startTime
130 130 self.endTime = endTime
131 131 self.walk = walk
132 132
133 133 def __checkPath(self):
134 134 if os.path.exists(self.path):
135 135 self.status = 1
136 136 else:
137 137 self.status = 0
138 138 print 'Path:%s does not exists'%self.path
139 139
140 140 return
141 141
142 142 def __selDates(self, amisr_dirname_format):
143 143 try:
144 144 year = int(amisr_dirname_format[0:4])
145 145 month = int(amisr_dirname_format[4:6])
146 146 dom = int(amisr_dirname_format[6:8])
147 147 thisDate = datetime.date(year,month,dom)
148 148
149 149 if (thisDate>=self.startDate and thisDate <= self.endDate):
150 150 return amisr_dirname_format
151 151 except:
152 152 return None
153 153
154 154 def __findDataForDates(self):
155 155
156 156
157 157
158 158 if not(self.status):
159 159 return None
160 160
161 161 pat = '\d+.\d+'
162 162 dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
163 163 dirnameList = filter(lambda x:x!=None,dirnameList)
164 164 dirnameList = [x.string for x in dirnameList]
165 165 dirnameList = [self.__selDates(x) for x in dirnameList]
166 166 dirnameList = filter(lambda x:x!=None,dirnameList)
167 167 if len(dirnameList)>0:
168 168 self.status = 1
169 169 self.dirnameList = dirnameList
170 170 self.dirnameList.sort()
171 171 else:
172 172 self.status = 0
173 173 return None
174 174
175 175 def __getTimeFromData(self):
176 pass
176 startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
177 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
178
179 print 'Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader)
180 print '........................................'
181 filter_filenameList = []
182 for filename in self.filenameList:
183 fp = h5py.File(filename,'r')
184 time_str = fp.get('Time/RadacTimeString')
185
186 startDateTimeStr_File = time_str[0][0].split('.')[0]
187 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
188 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
189
190 endDateTimeStr_File = time_str[-1][-1].split('.')[0]
191 junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
192 endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
193
194 fp.close()
195
196 if self.timezone == 'lt':
197 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
198 endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
199
200 if (endDateTime_File>=startDateTime_Reader and endDateTime_File<endDateTime_Reader):
201 #self.filenameList.remove(filename)
202 filter_filenameList.append(filename)
203
204 filter_filenameList.sort()
205 self.filenameList = filter_filenameList
206 return 1
177 207
178 208 def __filterByGlob1(self, dirName):
179 209 filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
180 210 filterDict = {}
181 211 filterDict.setdefault(dirName)
182 212 filterDict[dirName] = filter_files
183 213 return filterDict
184 214
185 215 def __getFilenameList(self, fileListInKeys, dirList):
186 216 for value in fileListInKeys:
187 217 dirName = value.keys()[0]
188 218 for file in value[dirName]:
189 219 filename = os.path.join(dirName, file)
190 220 self.filenameList.append(filename)
191 221
192 222
193 223 def __selectDataForTimes(self):
194 224 #aun no esta implementado el filtro for tiempo
195 225 if not(self.status):
196 226 return None
197 227
198 228 dirList = [os.path.join(self.path,x) for x in self.dirnameList]
199 229
200 230 fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
201 231
202 232 self.__getFilenameList(fileListInKeys, dirList)
203
233 #filtro por tiempo
234 if not(self.all):
235 self.__getTimeFromData()
236
237
204 238 if len(self.filenameList)>0:
205 239 self.status = 1
206 240 self.filenameList.sort()
207 241 else:
208 242 self.status = 0
209 243 return None
210 244
211 245
212 246 def __searchFilesOffline(self,
213 247 path,
214 248 startDate,
215 249 endDate,
216 250 startTime=datetime.time(0,0,0),
217 251 endTime=datetime.time(23,59,59),
218 252 walk=True):
219 253
220 254 self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
221 255
222 256 self.__checkPath()
223 257
224 258 self.__findDataForDates()
225 259
226 260 self.__selectDataForTimes()
227 261
228 262 for i in range(len(self.filenameList)):
229 263 print "%s" %(self.filenameList[i])
230 264
231 265 return
232 266
233 267 def __setNextFileOffline(self):
234 268 idFile = self.fileIndex
235 269
236 270 while (True):
237 271 idFile += 1
238 272 if not(idFile < len(self.filenameList)):
239 273 self.flagNoMoreFiles = 1
240 274 print "No more Files"
241 275 return 0
242 276
243 277 filename = self.filenameList[idFile]
244 278
245 279 amisrFilePointer = h5py.File(filename,'r')
246 280
247 281 break
248 282
249 283 self.flagIsNewFile = 1
250 284 self.fileIndex = idFile
251 285 self.filename = filename
252 286
253 287 self.amisrFilePointer = amisrFilePointer
254 288
255 289 print "Setting the file: %s"%self.filename
256 290
257 291 return 1
258 292
259 293 def __readHeader(self):
260 294 self.radacHeaderObj = RadacHeader(self.amisrFilePointer)
261 295
262 296 #update values from experiment cfg file
263 297 if self.radacHeaderObj.nrecords == self.recordsperfile_fromfile:
264 298 self.radacHeaderObj.nrecords = self.recordsperfile_fromfile
265 299 self.radacHeaderObj.nbeams = self.nbeamcodes_fromfile
266 300 self.radacHeaderObj.npulses = self.npulsesint_fromfile
267 301 self.radacHeaderObj.nsamples = self.ngates_fromfile
268 302
269 303 #looking index list for data
270 304 start_index = self.radacHeaderObj.pulseCount[0,:][0]
271 305 end_index = self.radacHeaderObj.npulses
272 306 range4data = range(start_index, end_index)
273 307 self.index4_schain_datablock = numpy.array(range4data)
274 308
275 309 buffer_start_index = 0
276 310 buffer_end_index = self.radacHeaderObj.pulseCount[0,:][0]
277 311 range4buffer = range(buffer_start_index, buffer_end_index)
278 312 self.index4_buffer = numpy.array(range4buffer)
279 313
280 314 self.linear_pulseCount = numpy.array(range4data + range4buffer)
281 315 self.npulseByFrame = max(self.radacHeaderObj.pulseCount[0,:]+1)
282 316
283 317 #get tuning frequency
284 318 frequency_h5file_dataset = self.amisrFilePointer.get('Rx'+'/TuningFrequency')
285 319 self.frequency_h5file = frequency_h5file_dataset[0,0]
286 320
287 321 self.flagIsNewFile = 1
288 322
289 323 def __getBeamCode(self):
290 324 self.beamCodeDict = {}
291 325 self.beamRangeDict = {}
292 326
327 beamCodeMap = self.amisrFilePointer.get('Setup/BeamcodeMap')
328
293 329 for i in range(len(self.radacHeaderObj.beamCode[0,:])):
294 330 self.beamCodeDict.setdefault(i)
295 331 self.beamRangeDict.setdefault(i)
296 self.beamCodeDict[i] = self.radacHeaderObj.beamCode[0,i]
332 beamcodeValue = self.radacHeaderObj.beamCode[0,i]
333 beamcodeIndex = numpy.where(beamCodeMap[:,0] == beamcodeValue)[0][0]
334 x = beamCodeMap[beamcodeIndex][1]
335 y = beamCodeMap[beamcodeIndex][2]
336 z = beamCodeMap[beamcodeIndex][3]
337 self.beamCodeDict[i] = [beamcodeValue, x, y, z]
297 338
298
299 339 just4record0 = self.radacHeaderObj.beamCodeByPulse[0,:]
300 340
301 341 for i in range(len(self.beamCodeDict.values())):
302 xx = numpy.where(just4record0==self.beamCodeDict.values()[i])
342 xx = numpy.where(just4record0==self.beamCodeDict.values()[i][0])
303 343 self.beamRangeDict[i] = xx[0]
304 344
305 345 def __getExpParameters(self):
306 346 if not(self.status):
307 347 return None
308 348
309 349 experimentCfgPath = os.path.join(self.path, self.dirnameList[0], 'Setup')
310 350
311 351 expFinder = glob.glob1(experimentCfgPath,'*.exp')
312 352 if len(expFinder)== 0:
313 353 self.status = 0
314 354 return None
315 355
316 356 experimentFilename = os.path.join(experimentCfgPath,expFinder[0])
317 357
318 358 f = open(experimentFilename)
319 359 lines = f.readlines()
320 360 f.close()
321 361
322 362 parmsList = ['npulsesint*','recordsperfile*','nbeamcodes*','ngates*']
323 363 filterList = [fnmatch.filter(lines, x) for x in parmsList]
324 364
325 365
326 366 values = [re.sub(r'\D',"",x[0]) for x in filterList]
327 367
328 368 self.npulsesint_fromfile = int(values[0])
329 369 self.recordsperfile_fromfile = int(values[1])
330 370 self.nbeamcodes_fromfile = int(values[2])
331 371 self.ngates_fromfile = int(values[3])
332 372
333 373 tufileFinder = fnmatch.filter(lines, 'tufile=*')
334 374 tufile = tufileFinder[0].split('=')[1].split('\n')[0]
375 tufile = tufile.split('\r')[0]
335 376 tufilename = os.path.join(experimentCfgPath,tufile)
336 377
337 378 f = open(tufilename)
338 379 lines = f.readlines()
339 380 f.close()
340 381 self.ippSeconds_fromfile = float(lines[1].split()[2])/1E6
341 382
342 383
343 384 self.status = 1
344 385
345 386 def __setIdsAndArrays(self):
346 387 self.dataByFrame = self.__setDataByFrame()
347 388 self.beamCodeByFrame = self.amisrFilePointer.get('Raw11/Data/RadacHeader/BeamCode').value[0, :]
348 389 self.readRanges()
349 390 self.index_amisr_sample, self.index_amisr_buffer = self.radacHeaderObj.getIndexRangeToPulse(0)
350 391 self.radacTimeByFrame = numpy.zeros(self.radacHeaderObj.npulses)
351 392 self.buffer_radactime = numpy.zeros_like(self.radacTimeByFrame)
352 393
353 394
354 395 def __setNextFile(self):
355 396
356 397 newFile = self.__setNextFileOffline()
357 398
358 399 if not(newFile):
359 400 return 0
360 401
361 402 self.__readHeader()
362 403
363 404 if self.__firstFile:
364 405 self.__setIdsAndArrays()
365 406 self.__firstFile = False
366 407
367 408 self.__getBeamCode()
368 409 self.readDataBlock()
369 410
370 411
371 412 def setup(self,path=None,
372 413 startDate=None,
373 414 endDate=None,
374 415 startTime=datetime.time(0,0,0),
375 416 endTime=datetime.time(23,59,59),
376 417 walk=True,
377 timezone='ut',):
418 timezone='ut',
419 all=0,):
378 420
379 421 self.timezone = timezone
422 self.all = all
380 423 #Busqueda de archivos offline
381 424 self.__searchFilesOffline(path, startDate, endDate, startTime, endTime, walk)
382 425
383 426 if not(self.filenameList):
384 427 print "There is no files into the folder: %s"%(path)
385 428
386 429 sys.exit(-1)
387 430
388 431 self.__getExpParameters()
389 432
390 433 self.fileIndex = -1
391 434
392 435 self.__setNextFile()
393 436
394 437 first_beamcode = self.radacHeaderObj.beamCodeByPulse[0,0]
395 438 index = numpy.where(self.radacHeaderObj.beamCodeByPulse[0,:]!=first_beamcode)[0][0]
396 439 self.profileIndex_offset = self.radacHeaderObj.pulseCount[0,:][index]
397 440 self.profileIndex = self.profileIndex_offset
398 441
399 442 def readRanges(self):
400 443 dataset = self.amisrFilePointer.get('Raw11/Data/Samples/Range')
401 444 #self.rangeFromFile = dataset.value
402 445 self.rangeFromFile = numpy.reshape(dataset.value,(-1))
403 446 return range
404 447
405 448
406 449 def readRadacTime(self,idrecord, range1, range2):
407 450 self.radacTimeFromFile = self.radacHeaderObj.radacTime.value
408 451
409 452 radacTimeByFrame = numpy.zeros((self.radacHeaderObj.npulses))
410 453 #radacTimeByFrame = dataset[idrecord - 1,range1]
411 454 #radacTimeByFrame = dataset[idrecord,range2]
412 455
413 456 return radacTimeByFrame
414 457
415 458 def readBeamCode(self, idrecord, range1, range2):
416 459 dataset = self.amisrFilePointer.get('Raw11/Data/RadacHeader/BeamCode')
417 460 beamcodeByFrame = numpy.zeros((self.radacHeaderObj.npulses))
418 461 self.beamCodesFromFile = dataset.value
419 462
420 463 #beamcodeByFrame[range1] = dataset[idrecord - 1, range1]
421 464 #beamcodeByFrame[range2] = dataset[idrecord, range2]
422 465 beamcodeByFrame[range1] = dataset[idrecord, range1]
423 466 beamcodeByFrame[range2] = dataset[idrecord, range2]
424 467
425 468 return beamcodeByFrame
426 469
427 470
428 471 def __setDataByFrame(self):
429 472 ndata = 2 # porque es complejo
430 473 dataByFrame = numpy.zeros((self.radacHeaderObj.npulses, self.radacHeaderObj.nsamples, ndata))
431 474 return dataByFrame
432 475
433 476 def __readDataSet(self):
434 477 dataset = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
435 478 return dataset
436 479
437 480 def __setDataBlock(self,):
438 481 real = self.dataByFrame[:,:,0] #asumo que 0 es real
439 482 imag = self.dataByFrame[:,:,1] #asumo que 1 es imaginario
440 483 datablock = real + imag*1j #armo el complejo
441 484 return datablock
442 485
443 486 def readSamples_version1(self,idrecord):
444 487 #estas tres primeras lineas solo se deben ejecutar una vez
445 488 if self.flagIsNewFile:
446 489 #reading dataset
447 490 self.dataset = self.__readDataSet()
448 491 self.flagIsNewFile = 0
449 492
450 493 if idrecord == 0:
451 494 self.dataByFrame[self.index4_schain_datablock, : ,:] = self.dataset[0, self.index_amisr_sample,:,:]
452 495 self.radacTimeByFrame[self.index4_schain_datablock] = self.radacHeaderObj.radacTime[0, self.index_amisr_sample]
453 496 datablock = self.__setDataBlock()
454 497
455 498 self.buffer = self.dataset[0, self.index_amisr_buffer,:,:]
456 499 self.buffer_radactime = self.radacHeaderObj.radacTime[0, self.index_amisr_buffer]
457 500
458 501 return datablock
459 502
460 503 self.dataByFrame[self.index4_buffer,:,:] = self.buffer.copy()
461 504 self.radacTimeByFrame[self.index4_buffer] = self.buffer_radactime.copy()
462 505 self.dataByFrame[self.index4_schain_datablock,:,:] = self.dataset[idrecord, self.index_amisr_sample,:,:]
463 506 self.radacTimeByFrame[self.index4_schain_datablock] = self.radacHeaderObj.radacTime[idrecord, self.index_amisr_sample]
464 507 datablock = self.__setDataBlock()
465 508
466 509 self.buffer = self.dataset[idrecord, self.index_amisr_buffer, :, :]
467 510 self.buffer_radactime = self.radacHeaderObj.radacTime[idrecord, self.index_amisr_buffer]
468 511
469 512 return datablock
470 513
471 514
472 515 def readSamples(self,idrecord):
473 516 if self.flagIsNewFile:
474 517 self.dataByFrame = self.__setDataByFrame()
475 518 self.beamCodeByFrame = self.amisrFilePointer.get('Raw11/Data/RadacHeader/BeamCode').value[idrecord, :]
476 519
477 520 #reading ranges
478 521 self.readRanges()
479 522 #reading dataset
480 523 self.dataset = self.__readDataSet()
481 524
482 525 self.flagIsNewFile = 0
483 526 self.radacTimeByFrame = self.radacHeaderObj.radacTime.value[idrecord, :]
484 527 self.dataByFrame = self.dataset[idrecord, :, :, :]
485 528 datablock = self.__setDataBlock()
486 529 return datablock
487 530
488 531
489 532 def readDataBlock(self):
490 533
491 534 self.datablock = self.readSamples_version1(self.idrecord_count)
492 535 #self.datablock = self.readSamples(self.idrecord_count)
493 536 #print 'record:', self.idrecord_count
494 537
495 538 self.idrecord_count += 1
496 539 self.profileIndex = 0
497 540
498 541 if self.idrecord_count >= self.radacHeaderObj.nrecords:
499 542 self.idrecord_count = 0
500 543 self.flagIsNewFile = 1
501 544
502 545 def readNextBlock(self):
503 546
504 547 self.readDataBlock()
505 548
506 549 if self.flagIsNewFile:
507 550 self.__setNextFile()
508 551 pass
509 552
510 553 def __hasNotDataInBuffer(self):
511 554 #self.radacHeaderObj.npulses debe ser otra variable para considerar el numero de pulsos a tomar en el primer y ultimo record
512 555 if self.profileIndex >= self.radacHeaderObj.npulses:
513 556 return 1
514 557 return 0
515 558
516 559 def printUTC(self):
517 560 print self.dataOut.utctime
518 561 print ''
519 562
520 563 def setObjProperties(self):
521 564 self.dataOut.heightList = self.rangeFromFile/1000.0 #km
522 565 self.dataOut.nProfiles = self.radacHeaderObj.npulses
523 566 self.dataOut.nRecords = self.radacHeaderObj.nrecords
524 567 self.dataOut.nBeams = self.radacHeaderObj.nbeams
525 568 self.dataOut.ippSeconds = self.ippSeconds_fromfile
526 569 self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
527 570 self.dataOut.frequency = self.frequency_h5file
528 571 self.dataOut.npulseByFrame = self.npulseByFrame
529 572 self.dataOut.nBaud = None
530 573 self.dataOut.nCode = None
531 574 self.dataOut.code = None
532 575
533 576 self.dataOut.beamCodeDict = self.beamCodeDict
534 577 self.dataOut.beamRangeDict = self.beamRangeDict
535 578
536 579 if self.timezone == 'lt':
537 580 self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
538 581 else:
539 582 self.dataOut.timeZone = 0 #by default time is UTC
540 583
541 584 def getData(self):
542 585
543 586 if self.flagNoMoreFiles:
544 587 self.dataOut.flagNoData = True
545 588 print 'Process finished'
546 589 return 0
547 590
548 591 if self.__hasNotDataInBuffer():
549 592 self.readNextBlock()
550 593
551 594
552 595 if self.datablock == None: # setear esta condicion cuando no hayan datos por leers
553 596 self.dataOut.flagNoData = True
554 597 return 0
555 598
556 599 self.dataOut.data = numpy.reshape(self.datablock[self.profileIndex,:],(1,-1))
557 600
558 601 self.dataOut.utctime = self.radacTimeByFrame[self.profileIndex]
559 602
560 603 self.dataOut.flagNoData = False
561 604
562 605 self.profileIndex += 1
563 606
564 607 return self.dataOut.data
565 608
566 609
567 610 def run(self, **kwargs):
568 611 if not(self.isConfig):
569 612 self.setup(**kwargs)
570 613 self.setObjProperties()
571 614 self.isConfig = True
572 615
573 616 self.getData()
@@ -1,87 +1,91
1 1 '''
2 2 @author: Daniel Suarez
3 3 '''
4 4
5 5 from jroproc_base import ProcessingUnit, Operation
6 6 from model.data.jroamisr import AMISR
7 7
8 8 class AMISRProc(ProcessingUnit):
9 9 def __init__(self):
10 10 ProcessingUnit.__init__(self)
11 11 self.objectDict = {}
12 12 self.dataOut = AMISR()
13 13
14 14 def run(self):
15 15 if self.dataIn.type == 'AMISR':
16 16 self.dataOut.copy(self.dataIn)
17 17
18 18
19 19 class PrintInfo(Operation):
20 20 def __init__(self):
21 21 self.__isPrinted = False
22 22
23 23 def run(self, dataOut):
24 24
25 25 if not self.__isPrinted:
26 26 print 'Number of Records by File: %d'%dataOut.nRecords
27 27 print 'Number of Pulses: %d'%dataOut.nProfiles
28 28 print 'Number of Pulses by Frame: %d'%dataOut.npulseByFrame
29 29 print 'Number of Samples by Pulse: %d'%len(dataOut.heightList)
30 30 print 'Ipp Seconds: %f'%dataOut.ippSeconds
31 31 print 'Number of Beams: %d'%dataOut.nBeams
32 32 print 'BeamCodes:'
33 beamStrList = ['Beam %d -> Code %d'%(k,v) for k,v in dataOut.beamCodeDict.items()]
33 beamStrList = ['Beam %d -> Code=%d, azimuth=%2.2f, zenith=%2.2f, gain=%2.2f'%(k,v[0],v[1],v[2],v[3]) for k,v in dataOut.beamCodeDict.items()]
34 34 for b in beamStrList:
35 35 print b
36 36 self.__isPrinted = True
37 37
38 38 return
39 39
40 40
41 41 class BeamSelector(Operation):
42 42 profileIndex = None
43 43 nProfiles = None
44 44
45 45 def __init__(self):
46 46
47 47 self.profileIndex = 0
48 48
49 49 def incIndex(self):
50 50 self.profileIndex += 1
51 51
52 52 if self.profileIndex >= self.nProfiles:
53 53 self.profileIndex = 0
54 54
55 55 def isProfileInRange(self, minIndex, maxIndex):
56 56
57 57 if self.profileIndex < minIndex:
58 58 return False
59 59
60 60 if self.profileIndex > maxIndex:
61 61 return False
62 62
63 63 return True
64 64
65 65 def isProfileInList(self, profileList):
66 66
67 67 if self.profileIndex not in profileList:
68 68 return False
69 69
70 70 return True
71 71
72 72 def run(self, dataOut, beam=None):
73 73
74 74 dataOut.flagNoData = True
75 75 self.nProfiles = dataOut.nProfiles
76 76
77 77 if beam != None:
78 78 if self.isProfileInList(dataOut.beamRangeDict[beam]):
79 beamInfo = dataOut.beamCodeDict[beam]
80 dataOut.azimuth = beamInfo[1]
81 dataOut.zenith = beamInfo[2]
82 dataOut.gain = beamInfo[3]
79 83 dataOut.flagNoData = False
80 84
81 85 self.incIndex()
82 86 return 1
83 87
84 88 else:
85 89 raise ValueError, "BeamSelector needs beam value"
86 90
87 91 return 0 No newline at end of file
@@ -1,927 +1,930
1 1 import numpy
2 2
3 3 from jroproc_base import ProcessingUnit, Operation
4 4 from model.data.jrodata import Spectra
5 5 from model.data.jrodata import hildebrand_sekhon
6 6
7 7 class SpectraProc(ProcessingUnit):
8 8
9 9 def __init__(self):
10 10
11 11 ProcessingUnit.__init__(self)
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 __updateObjFromInput(self):
21 21
22 22 self.dataOut.timeZone = self.dataIn.timeZone
23 23 self.dataOut.dstFlag = self.dataIn.dstFlag
24 24 self.dataOut.errorCount = self.dataIn.errorCount
25 25 self.dataOut.useLocalTime = self.dataIn.useLocalTime
26 26
27 27 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
28 28 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
29 29 self.dataOut.channelList = self.dataIn.channelList
30 30 self.dataOut.heightList = self.dataIn.heightList
31 31 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
32 32 # self.dataOut.nHeights = self.dataIn.nHeights
33 33 # self.dataOut.nChannels = self.dataIn.nChannels
34 34 self.dataOut.nBaud = self.dataIn.nBaud
35 35 self.dataOut.nCode = self.dataIn.nCode
36 36 self.dataOut.code = self.dataIn.code
37 37 self.dataOut.nProfiles = self.dataOut.nFFTPoints
38 38 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
39 39 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
40 40 self.dataOut.utctime = self.firstdatatime
41 41 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
42 42 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
43 43 # self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.nIncohInt = 1
46 46 # self.dataOut.ippSeconds = self.dataIn.ippSeconds
47 47 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
48 48
49 49 self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt
50 50 self.dataOut.frequency = self.dataIn.frequency
51 51 self.dataOut.realtime = self.dataIn.realtime
52 52
53 self.dataOut.azimuth = self.dataIn.azimuth
54 self.dataOut.zenith = self.dataIn.zenith
55
53 56 def __getFft(self):
54 57 """
55 58 Convierte valores de Voltaje a Spectra
56 59
57 60 Affected:
58 61 self.dataOut.data_spc
59 62 self.dataOut.data_cspc
60 63 self.dataOut.data_dc
61 64 self.dataOut.heightList
62 65 self.profIndex
63 66 self.buffer
64 67 self.dataOut.flagNoData
65 68 """
66 69 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
67 70 fft_volt = fft_volt.astype(numpy.dtype('complex'))
68 71 dc = fft_volt[:,0,:]
69 72
70 73 #calculo de self-spectra
71 74 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
72 75 spc = fft_volt * numpy.conjugate(fft_volt)
73 76 spc = spc.real
74 77
75 78 blocksize = 0
76 79 blocksize += dc.size
77 80 blocksize += spc.size
78 81
79 82 cspc = None
80 83 pairIndex = 0
81 84 if self.dataOut.pairsList != None:
82 85 #calculo de cross-spectra
83 86 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
84 87 for pair in self.dataOut.pairsList:
85 88 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
86 89 pairIndex += 1
87 90 blocksize += cspc.size
88 91
89 92 self.dataOut.data_spc = spc
90 93 self.dataOut.data_cspc = cspc
91 94 self.dataOut.data_dc = dc
92 95 self.dataOut.blockSize = blocksize
93 96 self.dataOut.flagShiftFFT = False
94 97
95 98 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None):
96 99
97 100 self.dataOut.flagNoData = True
98 101
99 102 if self.dataIn.type == "Spectra":
100 103 self.dataOut.copy(self.dataIn)
101 104 return True
102 105
103 106 if self.dataIn.type == "Voltage":
104 107
105 108 if nFFTPoints == None:
106 109 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
107 110
108 111 if nProfiles == None:
109 112 raise ValueError, "This SpectraProc.run() need nProfiles input variable"
110 113
111 114
112 115 if ippFactor == None:
113 116 ippFactor = 1
114 117 self.dataOut.ippFactor = ippFactor
115 118
116 119 self.dataOut.nFFTPoints = nFFTPoints
117 120 self.dataOut.pairsList = pairsList
118 121
119 122 if self.buffer == None:
120 123 self.buffer = numpy.zeros((self.dataIn.nChannels,
121 124 nProfiles,
122 125 self.dataIn.nHeights),
123 126 dtype='complex')
124 127 self.id_min = 0
125 128 self.id_max = self.dataIn.data.shape[1]
126 129
127 130 if len(self.dataIn.data.shape) == 2:
128 131 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
129 132 self.profIndex += 1
130 133 else:
131 134 if self.dataIn.data.shape[1] == nProfiles:
132 135 self.buffer = self.dataIn.data.copy()
133 136 self.profIndex = nProfiles
134 137 elif self.dataIn.data.shape[1] < nProfiles:
135 138 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
136 139 self.profIndex += self.dataIn.data.shape[1]
137 140 self.id_min += self.dataIn.data.shape[1]
138 141 self.id_max += self.dataIn.data.shape[1]
139 142 else:
140 143 raise ValueError, "The type object %s has %d profiles, it should be equal to %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
141 144 self.dataOut.flagNoData = True
142 145 return 0
143 146
144 147
145 148 if self.firstdatatime == None:
146 149 self.firstdatatime = self.dataIn.utctime
147 150
148 151 if self.profIndex == nProfiles:
149 152 self.__updateObjFromInput()
150 153 self.__getFft()
151 154
152 155 self.dataOut.flagNoData = False
153 156
154 157 self.buffer = None
155 158 self.firstdatatime = None
156 159 self.profIndex = 0
157 160
158 161 return True
159 162
160 163 raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
161 164
162 165 def selectChannels(self, channelList):
163 166
164 167 channelIndexList = []
165 168
166 169 for channel in channelList:
167 170 index = self.dataOut.channelList.index(channel)
168 171 channelIndexList.append(index)
169 172
170 173 self.selectChannelsByIndex(channelIndexList)
171 174
172 175 def selectChannelsByIndex(self, channelIndexList):
173 176 """
174 177 Selecciona un bloque de datos en base a canales segun el channelIndexList
175 178
176 179 Input:
177 180 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
178 181
179 182 Affected:
180 183 self.dataOut.data_spc
181 184 self.dataOut.channelIndexList
182 185 self.dataOut.nChannels
183 186
184 187 Return:
185 188 None
186 189 """
187 190
188 191 for channelIndex in channelIndexList:
189 192 if channelIndex not in self.dataOut.channelIndexList:
190 193 print channelIndexList
191 194 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
192 195
193 196 # nChannels = len(channelIndexList)
194 197
195 198 data_spc = self.dataOut.data_spc[channelIndexList,:]
196 199
197 200 self.dataOut.data_spc = data_spc
198 201 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
199 202 # self.dataOut.nChannels = nChannels
200 203
201 204 return 1
202 205
203 206 def selectHeights(self, minHei, maxHei):
204 207 """
205 208 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
206 209 minHei <= height <= maxHei
207 210
208 211 Input:
209 212 minHei : valor minimo de altura a considerar
210 213 maxHei : valor maximo de altura a considerar
211 214
212 215 Affected:
213 216 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
214 217
215 218 Return:
216 219 1 si el metodo se ejecuto con exito caso contrario devuelve 0
217 220 """
218 221 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
219 222 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
220 223
221 224 if (maxHei > self.dataOut.heightList[-1]):
222 225 maxHei = self.dataOut.heightList[-1]
223 226 # raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
224 227
225 228 minIndex = 0
226 229 maxIndex = 0
227 230 heights = self.dataOut.heightList
228 231
229 232 inda = numpy.where(heights >= minHei)
230 233 indb = numpy.where(heights <= maxHei)
231 234
232 235 try:
233 236 minIndex = inda[0][0]
234 237 except:
235 238 minIndex = 0
236 239
237 240 try:
238 241 maxIndex = indb[0][-1]
239 242 except:
240 243 maxIndex = len(heights)
241 244
242 245 self.selectHeightsByIndex(minIndex, maxIndex)
243 246
244 247 return 1
245 248
246 249 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
247 250 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
248 251
249 252 if hei_ref != None:
250 253 newheis = numpy.where(self.dataOut.heightList>hei_ref)
251 254
252 255 minIndex = min(newheis[0])
253 256 maxIndex = max(newheis[0])
254 257 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
255 258 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
256 259
257 260 # determina indices
258 261 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
259 262 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
260 263 beacon_dB = numpy.sort(avg_dB)[-nheis:]
261 264 beacon_heiIndexList = []
262 265 for val in avg_dB.tolist():
263 266 if val >= beacon_dB[0]:
264 267 beacon_heiIndexList.append(avg_dB.tolist().index(val))
265 268
266 269 #data_spc = data_spc[:,:,beacon_heiIndexList]
267 270 data_cspc = None
268 271 if self.dataOut.data_cspc != None:
269 272 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
270 273 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
271 274
272 275 data_dc = None
273 276 if self.dataOut.data_dc != None:
274 277 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
275 278 #data_dc = data_dc[:,beacon_heiIndexList]
276 279
277 280 self.dataOut.data_spc = data_spc
278 281 self.dataOut.data_cspc = data_cspc
279 282 self.dataOut.data_dc = data_dc
280 283 self.dataOut.heightList = heightList
281 284 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
282 285
283 286 return 1
284 287
285 288
286 289 def selectHeightsByIndex(self, minIndex, maxIndex):
287 290 """
288 291 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
289 292 minIndex <= index <= maxIndex
290 293
291 294 Input:
292 295 minIndex : valor de indice minimo de altura a considerar
293 296 maxIndex : valor de indice maximo de altura a considerar
294 297
295 298 Affected:
296 299 self.dataOut.data_spc
297 300 self.dataOut.data_cspc
298 301 self.dataOut.data_dc
299 302 self.dataOut.heightList
300 303
301 304 Return:
302 305 1 si el metodo se ejecuto con exito caso contrario devuelve 0
303 306 """
304 307
305 308 if (minIndex < 0) or (minIndex > maxIndex):
306 309 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
307 310
308 311 if (maxIndex >= self.dataOut.nHeights):
309 312 maxIndex = self.dataOut.nHeights-1
310 313 # raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
311 314
312 315 # nHeights = maxIndex - minIndex + 1
313 316
314 317 #Spectra
315 318 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
316 319
317 320 data_cspc = None
318 321 if self.dataOut.data_cspc != None:
319 322 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
320 323
321 324 data_dc = None
322 325 if self.dataOut.data_dc != None:
323 326 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
324 327
325 328 self.dataOut.data_spc = data_spc
326 329 self.dataOut.data_cspc = data_cspc
327 330 self.dataOut.data_dc = data_dc
328 331
329 332 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
330 333
331 334 return 1
332 335
333 336 def removeDC(self, mode = 2):
334 337 jspectra = self.dataOut.data_spc
335 338 jcspectra = self.dataOut.data_cspc
336 339
337 340
338 341 num_chan = jspectra.shape[0]
339 342 num_hei = jspectra.shape[2]
340 343
341 344 if jcspectra != None:
342 345 jcspectraExist = True
343 346 num_pairs = jcspectra.shape[0]
344 347 else: jcspectraExist = False
345 348
346 349 freq_dc = jspectra.shape[1]/2
347 350 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
348 351
349 352 if ind_vel[0]<0:
350 353 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
351 354
352 355 if mode == 1:
353 356 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
354 357
355 358 if jcspectraExist:
356 359 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
357 360
358 361 if mode == 2:
359 362
360 363 vel = numpy.array([-2,-1,1,2])
361 364 xx = numpy.zeros([4,4])
362 365
363 366 for fil in range(4):
364 367 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
365 368
366 369 xx_inv = numpy.linalg.inv(xx)
367 370 xx_aux = xx_inv[0,:]
368 371
369 372 for ich in range(num_chan):
370 373 yy = jspectra[ich,ind_vel,:]
371 374 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
372 375
373 376 junkid = jspectra[ich,freq_dc,:]<=0
374 377 cjunkid = sum(junkid)
375 378
376 379 if cjunkid.any():
377 380 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
378 381
379 382 if jcspectraExist:
380 383 for ip in range(num_pairs):
381 384 yy = jcspectra[ip,ind_vel,:]
382 385 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
383 386
384 387
385 388 self.dataOut.data_spc = jspectra
386 389 self.dataOut.data_cspc = jcspectra
387 390
388 391 return 1
389 392
390 393 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
391 394
392 395 jspectra = self.dataOut.data_spc
393 396 jcspectra = self.dataOut.data_cspc
394 397 jnoise = self.dataOut.getNoise()
395 398 num_incoh = self.dataOut.nIncohInt
396 399
397 400 num_channel = jspectra.shape[0]
398 401 num_prof = jspectra.shape[1]
399 402 num_hei = jspectra.shape[2]
400 403
401 404 #hei_interf
402 405 if hei_interf == None:
403 406 count_hei = num_hei/2 #Como es entero no importa
404 407 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
405 408 hei_interf = numpy.asarray(hei_interf)[0]
406 409 #nhei_interf
407 410 if (nhei_interf == None):
408 411 nhei_interf = 5
409 412 if (nhei_interf < 1):
410 413 nhei_interf = 1
411 414 if (nhei_interf > count_hei):
412 415 nhei_interf = count_hei
413 416 if (offhei_interf == None):
414 417 offhei_interf = 0
415 418
416 419 ind_hei = range(num_hei)
417 420 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
418 421 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
419 422 mask_prof = numpy.asarray(range(num_prof))
420 423 num_mask_prof = mask_prof.size
421 424 comp_mask_prof = [0, num_prof/2]
422 425
423 426
424 427 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
425 428 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
426 429 jnoise = numpy.nan
427 430 noise_exist = jnoise[0] < numpy.Inf
428 431
429 432 #Subrutina de Remocion de la Interferencia
430 433 for ich in range(num_channel):
431 434 #Se ordena los espectros segun su potencia (menor a mayor)
432 435 power = jspectra[ich,mask_prof,:]
433 436 power = power[:,hei_interf]
434 437 power = power.sum(axis = 0)
435 438 psort = power.ravel().argsort()
436 439
437 440 #Se estima la interferencia promedio en los Espectros de Potencia empleando
438 441 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
439 442
440 443 if noise_exist:
441 444 # tmp_noise = jnoise[ich] / num_prof
442 445 tmp_noise = jnoise[ich]
443 446 junkspc_interf = junkspc_interf - tmp_noise
444 447 #junkspc_interf[:,comp_mask_prof] = 0
445 448
446 449 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
447 450 jspc_interf = jspc_interf.transpose()
448 451 #Calculando el espectro de interferencia promedio
449 452 noiseid = numpy.where(jspc_interf <= tmp_noise/ math.sqrt(num_incoh))
450 453 noiseid = noiseid[0]
451 454 cnoiseid = noiseid.size
452 455 interfid = numpy.where(jspc_interf > tmp_noise/ math.sqrt(num_incoh))
453 456 interfid = interfid[0]
454 457 cinterfid = interfid.size
455 458
456 459 if (cnoiseid > 0): jspc_interf[noiseid] = 0
457 460
458 461 #Expandiendo los perfiles a limpiar
459 462 if (cinterfid > 0):
460 463 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
461 464 new_interfid = numpy.asarray(new_interfid)
462 465 new_interfid = {x for x in new_interfid}
463 466 new_interfid = numpy.array(list(new_interfid))
464 467 new_cinterfid = new_interfid.size
465 468 else: new_cinterfid = 0
466 469
467 470 for ip in range(new_cinterfid):
468 471 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
469 472 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
470 473
471 474
472 475 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
473 476
474 477 #Removiendo la interferencia del punto de mayor interferencia
475 478 ListAux = jspc_interf[mask_prof].tolist()
476 479 maxid = ListAux.index(max(ListAux))
477 480
478 481
479 482 if cinterfid > 0:
480 483 for ip in range(cinterfid*(interf == 2) - 1):
481 484 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/math.sqrt(num_incoh))).nonzero()
482 485 cind = len(ind)
483 486
484 487 if (cind > 0):
485 488 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/math.sqrt(num_incoh))
486 489
487 490 ind = numpy.array([-2,-1,1,2])
488 491 xx = numpy.zeros([4,4])
489 492
490 493 for id1 in range(4):
491 494 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
492 495
493 496 xx_inv = numpy.linalg.inv(xx)
494 497 xx = xx_inv[:,0]
495 498 ind = (ind + maxid + num_mask_prof)%num_mask_prof
496 499 yy = jspectra[ich,mask_prof[ind],:]
497 500 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
498 501
499 502
500 503 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/math.sqrt(num_incoh))).nonzero()
501 504 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/math.sqrt(num_incoh))
502 505
503 506 #Remocion de Interferencia en el Cross Spectra
504 507 if jcspectra == None: return jspectra, jcspectra
505 508 num_pairs = jcspectra.size/(num_prof*num_hei)
506 509 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
507 510
508 511 for ip in range(num_pairs):
509 512
510 513 #-------------------------------------------
511 514
512 515 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
513 516 cspower = cspower[:,hei_interf]
514 517 cspower = cspower.sum(axis = 0)
515 518
516 519 cspsort = cspower.ravel().argsort()
517 520 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
518 521 junkcspc_interf = junkcspc_interf.transpose()
519 522 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
520 523
521 524 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
522 525
523 526 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
524 527 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
525 528 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
526 529
527 530 for iprof in range(num_prof):
528 531 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
529 532 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
530 533
531 534 #Removiendo la Interferencia
532 535 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
533 536
534 537 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
535 538 maxid = ListAux.index(max(ListAux))
536 539
537 540 ind = numpy.array([-2,-1,1,2])
538 541 xx = numpy.zeros([4,4])
539 542
540 543 for id1 in range(4):
541 544 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
542 545
543 546 xx_inv = numpy.linalg.inv(xx)
544 547 xx = xx_inv[:,0]
545 548
546 549 ind = (ind + maxid + num_mask_prof)%num_mask_prof
547 550 yy = jcspectra[ip,mask_prof[ind],:]
548 551 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
549 552
550 553 #Guardar Resultados
551 554 self.dataOut.data_spc = jspectra
552 555 self.dataOut.data_cspc = jcspectra
553 556
554 557 return 1
555 558
556 559 def setRadarFrequency(self, frequency=None):
557 560 if frequency != None:
558 561 self.dataOut.frequency = frequency
559 562
560 563 return 1
561 564
562 565 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
563 566 #validacion de rango
564 567 if minHei == None:
565 568 minHei = self.dataOut.heightList[0]
566 569
567 570 if maxHei == None:
568 571 maxHei = self.dataOut.heightList[-1]
569 572
570 573 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
571 574 print 'minHei: %.2f is out of the heights range'%(minHei)
572 575 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
573 576 minHei = self.dataOut.heightList[0]
574 577
575 578 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
576 579 print 'maxHei: %.2f is out of the heights range'%(maxHei)
577 580 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
578 581 maxHei = self.dataOut.heightList[-1]
579 582
580 583 # validacion de velocidades
581 584 velrange = self.dataOut.getVelRange(1)
582 585
583 586 if minVel == None:
584 587 minVel = velrange[0]
585 588
586 589 if maxVel == None:
587 590 maxVel = velrange[-1]
588 591
589 592 if (minVel < velrange[0]) or (minVel > maxVel):
590 593 print 'minVel: %.2f is out of the velocity range'%(minVel)
591 594 print 'minVel is setting to %.2f'%(velrange[0])
592 595 minVel = velrange[0]
593 596
594 597 if (maxVel > velrange[-1]) or (maxVel < minVel):
595 598 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
596 599 print 'maxVel is setting to %.2f'%(velrange[-1])
597 600 maxVel = velrange[-1]
598 601
599 602 # seleccion de indices para rango
600 603 minIndex = 0
601 604 maxIndex = 0
602 605 heights = self.dataOut.heightList
603 606
604 607 inda = numpy.where(heights >= minHei)
605 608 indb = numpy.where(heights <= maxHei)
606 609
607 610 try:
608 611 minIndex = inda[0][0]
609 612 except:
610 613 minIndex = 0
611 614
612 615 try:
613 616 maxIndex = indb[0][-1]
614 617 except:
615 618 maxIndex = len(heights)
616 619
617 620 if (minIndex < 0) or (minIndex > maxIndex):
618 621 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
619 622
620 623 if (maxIndex >= self.dataOut.nHeights):
621 624 maxIndex = self.dataOut.nHeights-1
622 625
623 626 # seleccion de indices para velocidades
624 627 indminvel = numpy.where(velrange >= minVel)
625 628 indmaxvel = numpy.where(velrange <= maxVel)
626 629 try:
627 630 minIndexVel = indminvel[0][0]
628 631 except:
629 632 minIndexVel = 0
630 633
631 634 try:
632 635 maxIndexVel = indmaxvel[0][-1]
633 636 except:
634 637 maxIndexVel = len(velrange)
635 638
636 639 #seleccion del espectro
637 640 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
638 641 #estimacion de ruido
639 642 noise = numpy.zeros(self.dataOut.nChannels)
640 643
641 644 for channel in range(self.dataOut.nChannels):
642 645 daux = data_spc[channel,:,:]
643 646 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
644 647
645 648 self.dataOut.noise_estimation = noise.copy()
646 649
647 650 return 1
648 651
649 652 class IncohInt(Operation):
650 653
651 654
652 655 __profIndex = 0
653 656 __withOverapping = False
654 657
655 658 __byTime = False
656 659 __initime = None
657 660 __lastdatatime = None
658 661 __integrationtime = None
659 662
660 663 __buffer_spc = None
661 664 __buffer_cspc = None
662 665 __buffer_dc = None
663 666
664 667 __dataReady = False
665 668
666 669 __timeInterval = None
667 670
668 671 n = None
669 672
670 673
671 674
672 675 def __init__(self):
673 676
674 677 Operation.__init__(self)
675 678 # self.isConfig = False
676 679
677 680 def setup(self, n=None, timeInterval=None, overlapping=False):
678 681 """
679 682 Set the parameters of the integration class.
680 683
681 684 Inputs:
682 685
683 686 n : Number of coherent integrations
684 687 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
685 688 overlapping :
686 689
687 690 """
688 691
689 692 self.__initime = None
690 693 self.__lastdatatime = 0
691 694 self.__buffer_spc = None
692 695 self.__buffer_cspc = None
693 696 self.__buffer_dc = None
694 697 self.__dataReady = False
695 698
696 699
697 700 if n == None and timeInterval == None:
698 701 raise ValueError, "n or timeInterval should be specified ..."
699 702
700 703 if n != None:
701 704 self.n = n
702 705 self.__byTime = False
703 706 else:
704 707 self.__integrationtime = timeInterval #if (type(timeInterval)!=integer) -> change this line
705 708 self.n = 9999
706 709 self.__byTime = True
707 710
708 711 if overlapping:
709 712 self.__withOverapping = True
710 713 else:
711 714 self.__withOverapping = False
712 715 self.__buffer_spc = 0
713 716 self.__buffer_cspc = 0
714 717 self.__buffer_dc = 0
715 718
716 719 self.__profIndex = 0
717 720
718 721 def putData(self, data_spc, data_cspc, data_dc):
719 722
720 723 """
721 724 Add a profile to the __buffer_spc and increase in one the __profileIndex
722 725
723 726 """
724 727
725 728 if not self.__withOverapping:
726 729 self.__buffer_spc += data_spc
727 730
728 731 if data_cspc == None:
729 732 self.__buffer_cspc = None
730 733 else:
731 734 self.__buffer_cspc += data_cspc
732 735
733 736 if data_dc == None:
734 737 self.__buffer_dc = None
735 738 else:
736 739 self.__buffer_dc += data_dc
737 740
738 741 self.__profIndex += 1
739 742 return
740 743
741 744 #Overlapping data
742 745 nChannels, nFFTPoints, nHeis = data_spc.shape
743 746 data_spc = numpy.reshape(data_spc, (1, nChannels, nFFTPoints, nHeis))
744 747 if data_cspc != None:
745 748 data_cspc = numpy.reshape(data_cspc, (1, -1, nFFTPoints, nHeis))
746 749 if data_dc != None:
747 750 data_dc = numpy.reshape(data_dc, (1, -1, nHeis))
748 751
749 752 #If the buffer is empty then it takes the data value
750 753 if self.__buffer_spc == None:
751 754 self.__buffer_spc = data_spc
752 755
753 756 if data_cspc == None:
754 757 self.__buffer_cspc = None
755 758 else:
756 759 self.__buffer_cspc += data_cspc
757 760
758 761 if data_dc == None:
759 762 self.__buffer_dc = None
760 763 else:
761 764 self.__buffer_dc += data_dc
762 765
763 766 self.__profIndex += 1
764 767 return
765 768
766 769 #If the buffer length is lower than n then stakcing the data value
767 770 if self.__profIndex < self.n:
768 771 self.__buffer_spc = numpy.vstack((self.__buffer_spc, data_spc))
769 772
770 773 if data_cspc != None:
771 774 self.__buffer_cspc = numpy.vstack((self.__buffer_cspc, data_cspc))
772 775
773 776 if data_dc != None:
774 777 self.__buffer_dc = numpy.vstack((self.__buffer_dc, data_dc))
775 778
776 779 self.__profIndex += 1
777 780 return
778 781
779 782 #If the buffer length is equal to n then replacing the last buffer value with the data value
780 783 self.__buffer_spc = numpy.roll(self.__buffer_spc, -1, axis=0)
781 784 self.__buffer_spc[self.n-1] = data_spc
782 785
783 786 if data_cspc != None:
784 787 self.__buffer_cspc = numpy.roll(self.__buffer_cspc, -1, axis=0)
785 788 self.__buffer_cspc[self.n-1] = data_cspc
786 789
787 790 if data_dc != None:
788 791 self.__buffer_dc = numpy.roll(self.__buffer_dc, -1, axis=0)
789 792 self.__buffer_dc[self.n-1] = data_dc
790 793
791 794 self.__profIndex = self.n
792 795 return
793 796
794 797
795 798 def pushData(self):
796 799 """
797 800 Return the sum of the last profiles and the profiles used in the sum.
798 801
799 802 Affected:
800 803
801 804 self.__profileIndex
802 805
803 806 """
804 807 data_spc = None
805 808 data_cspc = None
806 809 data_dc = None
807 810
808 811 if not self.__withOverapping:
809 812 data_spc = self.__buffer_spc
810 813 data_cspc = self.__buffer_cspc
811 814 data_dc = self.__buffer_dc
812 815
813 816 n = self.__profIndex
814 817
815 818 self.__buffer_spc = 0
816 819 self.__buffer_cspc = 0
817 820 self.__buffer_dc = 0
818 821 self.__profIndex = 0
819 822
820 823 return data_spc, data_cspc, data_dc, n
821 824
822 825 #Integration with Overlapping
823 826 data_spc = numpy.sum(self.__buffer_spc, axis=0)
824 827
825 828 if self.__buffer_cspc != None:
826 829 data_cspc = numpy.sum(self.__buffer_cspc, axis=0)
827 830
828 831 if self.__buffer_dc != None:
829 832 data_dc = numpy.sum(self.__buffer_dc, axis=0)
830 833
831 834 n = self.__profIndex
832 835
833 836 return data_spc, data_cspc, data_dc, n
834 837
835 838 def byProfiles(self, *args):
836 839
837 840 self.__dataReady = False
838 841 avgdata_spc = None
839 842 avgdata_cspc = None
840 843 avgdata_dc = None
841 844 # n = None
842 845
843 846 self.putData(*args)
844 847
845 848 if self.__profIndex == self.n:
846 849
847 850 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
848 851 self.__dataReady = True
849 852
850 853 return avgdata_spc, avgdata_cspc, avgdata_dc
851 854
852 855 def byTime(self, datatime, *args):
853 856
854 857 self.__dataReady = False
855 858 avgdata_spc = None
856 859 avgdata_cspc = None
857 860 avgdata_dc = None
858 861 n = None
859 862
860 863 self.putData(*args)
861 864
862 865 if (datatime - self.__initime) >= self.__integrationtime:
863 866 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
864 867 self.n = n
865 868 self.__dataReady = True
866 869
867 870 return avgdata_spc, avgdata_cspc, avgdata_dc
868 871
869 872 def integrate(self, datatime, *args):
870 873
871 874 if self.__initime == None:
872 875 self.__initime = datatime
873 876
874 877 if self.__byTime:
875 878 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
876 879 else:
877 880 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
878 881
879 882 self.__lastdatatime = datatime
880 883
881 884 if avgdata_spc == None:
882 885 return None, None, None, None
883 886
884 887 avgdatatime = self.__initime
885 888 try:
886 889 self.__timeInterval = (self.__lastdatatime - self.__initime)/(self.n - 1)
887 890 except:
888 891 self.__timeInterval = self.__lastdatatime - self.__initime
889 892
890 893 deltatime = datatime -self.__lastdatatime
891 894
892 895 if not self.__withOverapping:
893 896 self.__initime = datatime
894 897 else:
895 898 self.__initime += deltatime
896 899
897 900 return avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc
898 901
899 902 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
900 903
901 904 if n==1:
902 905 dataOut.flagNoData = False
903 906 return
904 907
905 908 if not self.isConfig:
906 909 self.setup(n, timeInterval, overlapping)
907 910 self.isConfig = True
908 911
909 912 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
910 913 dataOut.data_spc,
911 914 dataOut.data_cspc,
912 915 dataOut.data_dc)
913 916
914 917 # dataOut.timeInterval *= n
915 918 dataOut.flagNoData = True
916 919
917 920 if self.__dataReady:
918 921
919 922 dataOut.data_spc = avgdata_spc
920 923 dataOut.data_cspc = avgdata_cspc
921 924 dataOut.data_dc = avgdata_dc
922 925
923 926 dataOut.nIncohInt *= self.n
924 927 dataOut.utctime = avgdatatime
925 928 #dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt * dataOut.nIncohInt * dataOut.nFFTPoints
926 929 dataOut.timeInterval = self.__timeInterval*self.n
927 930 dataOut.flagNoData = False
@@ -1,751 +1,754
1 1 import numpy
2 2
3 3 from jroproc_base import ProcessingUnit, Operation
4 4 from model.data.jrodata import Voltage
5 5
6 6 class VoltageProc(ProcessingUnit):
7 7
8 8
9 9 def __init__(self):
10 10
11 11 ProcessingUnit.__init__(self)
12 12
13 13 # self.objectDict = {}
14 14 self.dataOut = Voltage()
15 15 self.flip = 1
16 16
17 17 def run(self):
18 18 if self.dataIn.type == 'AMISR':
19 19 self.__updateObjFromAmisrInput()
20 20
21 21 if self.dataIn.type == 'Voltage':
22 22 self.dataOut.copy(self.dataIn)
23 23
24 24 # self.dataOut.copy(self.dataIn)
25 25
26 26 def __updateObjFromAmisrInput(self):
27 27
28 28 self.dataOut.timeZone = self.dataIn.timeZone
29 29 self.dataOut.dstFlag = self.dataIn.dstFlag
30 30 self.dataOut.errorCount = self.dataIn.errorCount
31 31 self.dataOut.useLocalTime = self.dataIn.useLocalTime
32 32
33 33 self.dataOut.flagNoData = self.dataIn.flagNoData
34 34 self.dataOut.data = self.dataIn.data
35 35 self.dataOut.utctime = self.dataIn.utctime
36 36 self.dataOut.channelList = self.dataIn.channelList
37 37 self.dataOut.timeInterval = self.dataIn.timeInterval
38 38 self.dataOut.heightList = self.dataIn.heightList
39 39 self.dataOut.nProfiles = self.dataIn.nProfiles
40 40
41 41 self.dataOut.nCohInt = self.dataIn.nCohInt
42 42 self.dataOut.ippSeconds = self.dataIn.ippSeconds
43 43 self.dataOut.frequency = self.dataIn.frequency
44
45 self.dataOut.azimuth = self.dataIn.azimuth
46 self.dataOut.zenith = self.dataIn.zenith
44 47 #
45 48 # pass#
46 49 #
47 50 # def init(self):
48 51 #
49 52 #
50 53 # if self.dataIn.type == 'AMISR':
51 54 # self.__updateObjFromAmisrInput()
52 55 #
53 56 # if self.dataIn.type == 'Voltage':
54 57 # self.dataOut.copy(self.dataIn)
55 58 # # No necesita copiar en cada init() los atributos de dataIn
56 59 # # la copia deberia hacerse por cada nuevo bloque de datos
57 60
58 61 def selectChannels(self, channelList):
59 62
60 63 channelIndexList = []
61 64
62 65 for channel in channelList:
63 66 index = self.dataOut.channelList.index(channel)
64 67 channelIndexList.append(index)
65 68
66 69 self.selectChannelsByIndex(channelIndexList)
67 70
68 71 def selectChannelsByIndex(self, channelIndexList):
69 72 """
70 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
71 74
72 75 Input:
73 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
74 77
75 78 Affected:
76 79 self.dataOut.data
77 80 self.dataOut.channelIndexList
78 81 self.dataOut.nChannels
79 82 self.dataOut.m_ProcessingHeader.totalSpectra
80 83 self.dataOut.systemHeaderObj.numChannels
81 84 self.dataOut.m_ProcessingHeader.blockSize
82 85
83 86 Return:
84 87 None
85 88 """
86 89
87 90 for channelIndex in channelIndexList:
88 91 if channelIndex not in self.dataOut.channelIndexList:
89 92 print channelIndexList
90 93 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
91 94
92 95 # nChannels = len(channelIndexList)
93 96
94 97 data = self.dataOut.data[channelIndexList,:]
95 98
96 99 self.dataOut.data = data
97 100 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
98 101 # self.dataOut.nChannels = nChannels
99 102
100 103 return 1
101 104
102 105 def selectHeights(self, minHei=None, maxHei=None):
103 106 """
104 107 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
105 108 minHei <= height <= maxHei
106 109
107 110 Input:
108 111 minHei : valor minimo de altura a considerar
109 112 maxHei : valor maximo de altura a considerar
110 113
111 114 Affected:
112 115 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
113 116
114 117 Return:
115 118 1 si el metodo se ejecuto con exito caso contrario devuelve 0
116 119 """
117 120
118 121 if minHei == None:
119 122 minHei = self.dataOut.heightList[0]
120 123
121 124 if maxHei == None:
122 125 maxHei = self.dataOut.heightList[-1]
123 126
124 127 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
125 128 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
126 129
127 130
128 131 if (maxHei > self.dataOut.heightList[-1]):
129 132 maxHei = self.dataOut.heightList[-1]
130 133 # raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
131 134
132 135 minIndex = 0
133 136 maxIndex = 0
134 137 heights = self.dataOut.heightList
135 138
136 139 inda = numpy.where(heights >= minHei)
137 140 indb = numpy.where(heights <= maxHei)
138 141
139 142 try:
140 143 minIndex = inda[0][0]
141 144 except:
142 145 minIndex = 0
143 146
144 147 try:
145 148 maxIndex = indb[0][-1]
146 149 except:
147 150 maxIndex = len(heights)
148 151
149 152 self.selectHeightsByIndex(minIndex, maxIndex)
150 153
151 154 return 1
152 155
153 156
154 157 def selectHeightsByIndex(self, minIndex, maxIndex):
155 158 """
156 159 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
157 160 minIndex <= index <= maxIndex
158 161
159 162 Input:
160 163 minIndex : valor de indice minimo de altura a considerar
161 164 maxIndex : valor de indice maximo de altura a considerar
162 165
163 166 Affected:
164 167 self.dataOut.data
165 168 self.dataOut.heightList
166 169
167 170 Return:
168 171 1 si el metodo se ejecuto con exito caso contrario devuelve 0
169 172 """
170 173
171 174 if (minIndex < 0) or (minIndex > maxIndex):
172 175 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
173 176
174 177 if (maxIndex >= self.dataOut.nHeights):
175 178 maxIndex = self.dataOut.nHeights-1
176 179 # raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
177 180
178 181 # nHeights = maxIndex - minIndex + 1
179 182
180 183 #voltage
181 184 data = self.dataOut.data[:,minIndex:maxIndex+1]
182 185
183 186 # firstHeight = self.dataOut.heightList[minIndex]
184 187
185 188 self.dataOut.data = data
186 189 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
187 190
188 191 return 1
189 192
190 193
191 194 def filterByHeights(self, window, axis=1):
192 195 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
193 196
194 197 if window == None:
195 198 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
196 199
197 200 newdelta = deltaHeight * window
198 201 r = self.dataOut.data.shape[axis] % window
199 202 if axis == 1:
200 203 buffer = self.dataOut.data[:,0:self.dataOut.data.shape[axis]-r]
201 204 buffer = buffer.reshape(self.dataOut.data.shape[0],self.dataOut.data.shape[axis]/window,window)
202 205 buffer = numpy.sum(buffer,axis+1)
203 206
204 207 elif axis == 2:
205 208 buffer = self.dataOut.data[:, :, 0:self.dataOut.data.shape[axis]-r]
206 209 buffer = buffer.reshape(self.dataOut.data.shape[0],self.dataOut.data.shape[1],self.dataOut.data.shape[axis]/window,window)
207 210 buffer = numpy.sum(buffer,axis+1)
208 211
209 212 else:
210 213 raise ValueError, "axis value should be 1 or 2, the input value %d is not valid" % (axis)
211 214
212 215 self.dataOut.data = buffer.copy()
213 216 self.dataOut.heightList = numpy.arange(self.dataOut.heightList[0],newdelta*(self.dataOut.nHeights-r)/window,newdelta)
214 217 self.dataOut.windowOfFilter = window
215 218
216 219 return 1
217 220
218 221 def deFlip(self):
219 222 self.dataOut.data *= self.flip
220 223 self.flip *= -1.
221 224
222 225 def setRadarFrequency(self, frequency=None):
223 226 if frequency != None:
224 227 self.dataOut.frequency = frequency
225 228
226 229 return 1
227 230
228 231 class CohInt(Operation):
229 232
230 233 isConfig = False
231 234
232 235 __profIndex = 0
233 236 __withOverapping = False
234 237
235 238 __byTime = False
236 239 __initime = None
237 240 __lastdatatime = None
238 241 __integrationtime = None
239 242
240 243 __buffer = None
241 244
242 245 __dataReady = False
243 246
244 247 n = None
245 248
246 249
247 250 def __init__(self):
248 251
249 252 Operation.__init__(self)
250 253
251 254 # self.isConfig = False
252 255
253 256 def setup(self, n=None, timeInterval=None, overlapping=False, byblock=False):
254 257 """
255 258 Set the parameters of the integration class.
256 259
257 260 Inputs:
258 261
259 262 n : Number of coherent integrations
260 263 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
261 264 overlapping :
262 265
263 266 """
264 267
265 268 self.__initime = None
266 269 self.__lastdatatime = 0
267 270 self.__buffer = None
268 271 self.__dataReady = False
269 272 self.byblock = byblock
270 273
271 274 if n == None and timeInterval == None:
272 275 raise ValueError, "n or timeInterval should be specified ..."
273 276
274 277 if n != None:
275 278 self.n = n
276 279 self.__byTime = False
277 280 else:
278 281 self.__integrationtime = timeInterval * 60. #if (type(timeInterval)!=integer) -> change this line
279 282 self.n = 9999
280 283 self.__byTime = True
281 284
282 285 if overlapping:
283 286 self.__withOverapping = True
284 287 self.__buffer = None
285 288 else:
286 289 self.__withOverapping = False
287 290 self.__buffer = 0
288 291
289 292 self.__profIndex = 0
290 293
291 294 def putData(self, data):
292 295
293 296 """
294 297 Add a profile to the __buffer and increase in one the __profileIndex
295 298
296 299 """
297 300
298 301 if not self.__withOverapping:
299 302 self.__buffer += data.copy()
300 303 self.__profIndex += 1
301 304 return
302 305
303 306 #Overlapping data
304 307 nChannels, nHeis = data.shape
305 308 data = numpy.reshape(data, (1, nChannels, nHeis))
306 309
307 310 #If the buffer is empty then it takes the data value
308 311 if self.__buffer == None:
309 312 self.__buffer = data
310 313 self.__profIndex += 1
311 314 return
312 315
313 316 #If the buffer length is lower than n then stakcing the data value
314 317 if self.__profIndex < self.n:
315 318 self.__buffer = numpy.vstack((self.__buffer, data))
316 319 self.__profIndex += 1
317 320 return
318 321
319 322 #If the buffer length is equal to n then replacing the last buffer value with the data value
320 323 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
321 324 self.__buffer[self.n-1] = data
322 325 self.__profIndex = self.n
323 326 return
324 327
325 328
326 329 def pushData(self):
327 330 """
328 331 Return the sum of the last profiles and the profiles used in the sum.
329 332
330 333 Affected:
331 334
332 335 self.__profileIndex
333 336
334 337 """
335 338
336 339 if not self.__withOverapping:
337 340 data = self.__buffer
338 341 n = self.__profIndex
339 342
340 343 self.__buffer = 0
341 344 self.__profIndex = 0
342 345
343 346 return data, n
344 347
345 348 #Integration with Overlapping
346 349 data = numpy.sum(self.__buffer, axis=0)
347 350 n = self.__profIndex
348 351
349 352 return data, n
350 353
351 354 def byProfiles(self, data):
352 355
353 356 self.__dataReady = False
354 357 avgdata = None
355 358 # n = None
356 359
357 360 self.putData(data)
358 361
359 362 if self.__profIndex == self.n:
360 363
361 364 avgdata, n = self.pushData()
362 365 self.__dataReady = True
363 366
364 367 return avgdata
365 368
366 369 def byTime(self, data, datatime):
367 370
368 371 self.__dataReady = False
369 372 avgdata = None
370 373 n = None
371 374
372 375 self.putData(data)
373 376
374 377 if (datatime - self.__initime) >= self.__integrationtime:
375 378 avgdata, n = self.pushData()
376 379 self.n = n
377 380 self.__dataReady = True
378 381
379 382 return avgdata
380 383
381 384 def integrate(self, data, datatime=None):
382 385
383 386 if self.__initime == None:
384 387 self.__initime = datatime
385 388
386 389 if self.__byTime:
387 390 avgdata = self.byTime(data, datatime)
388 391 else:
389 392 avgdata = self.byProfiles(data)
390 393
391 394
392 395 self.__lastdatatime = datatime
393 396
394 397 if avgdata == None:
395 398 return None, None
396 399
397 400 avgdatatime = self.__initime
398 401
399 402 deltatime = datatime -self.__lastdatatime
400 403
401 404 if not self.__withOverapping:
402 405 self.__initime = datatime
403 406 else:
404 407 self.__initime += deltatime
405 408
406 409 return avgdata, avgdatatime
407 410
408 411 def integrateByBlock(self, dataOut):
409 412 times = int(dataOut.data.shape[1]/self.n)
410 413 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
411 414
412 415 id_min = 0
413 416 id_max = self.n
414 417
415 418 for i in range(times):
416 419 junk = dataOut.data[:,id_min:id_max,:]
417 420 avgdata[:,i,:] = junk.sum(axis=1)
418 421 id_min += self.n
419 422 id_max += self.n
420 423
421 424 timeInterval = dataOut.ippSeconds*self.n
422 425 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
423 426 self.__dataReady = True
424 427 return avgdata, avgdatatime
425 428
426 429 def run(self, dataOut, **kwargs):
427 430
428 431 if not self.isConfig:
429 432 self.setup(**kwargs)
430 433 self.isConfig = True
431 434
432 435 if self.byblock:
433 436 avgdata, avgdatatime = self.integrateByBlock(dataOut)
434 437 else:
435 438 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
436 439
437 440 # dataOut.timeInterval *= n
438 441 dataOut.flagNoData = True
439 442
440 443 if self.__dataReady:
441 444 dataOut.data = avgdata
442 445 dataOut.nCohInt *= self.n
443 446 dataOut.utctime = avgdatatime
444 447 dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
445 448 dataOut.flagNoData = False
446 449
447 450 class Decoder(Operation):
448 451
449 452 isConfig = False
450 453 __profIndex = 0
451 454
452 455 code = None
453 456
454 457 nCode = None
455 458 nBaud = None
456 459
457 460
458 461 def __init__(self):
459 462
460 463 Operation.__init__(self)
461 464
462 465 self.times = None
463 466 self.osamp = None
464 467 self.__setValues = False
465 468 # self.isConfig = False
466 469
467 470 def setup(self, code, shape, times, osamp):
468 471
469 472 self.__profIndex = 0
470 473
471 474 self.code = code
472 475
473 476 self.nCode = len(code)
474 477 self.nBaud = len(code[0])
475 478
476 479 if times != None:
477 480 self.times = times
478 481
479 482 if ((osamp != None) and (osamp >1)):
480 483 self.osamp = osamp
481 484 self.code = numpy.repeat(code, repeats=self.osamp,axis=1)
482 485 self.nBaud = self.nBaud*self.osamp
483 486
484 487 if len(shape) == 2:
485 488 self.__nChannels, self.__nHeis = shape
486 489
487 490 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
488 491
489 492 __codeBuffer[:,0:self.nBaud] = self.code
490 493
491 494 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
492 495
493 496 self.ndatadec = self.__nHeis - self.nBaud + 1
494 497
495 498 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
496 499 else:
497 500 self.__nChannels, self.__nProfiles, self.__nHeis = shape
498 501
499 502 self.ndatadec = self.__nHeis - self.nBaud + 1
500 503
501 504 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
502 505
503 506
504 507
505 508 def convolutionInFreq(self, data):
506 509
507 510 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
508 511
509 512 fft_data = numpy.fft.fft(data, axis=1)
510 513
511 514 conv = fft_data*fft_code
512 515
513 516 data = numpy.fft.ifft(conv,axis=1)
514 517
515 518 datadec = data[:,:-self.nBaud+1]
516 519
517 520 return datadec
518 521
519 522 def convolutionInFreqOpt(self, data):
520 523
521 524 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
522 525
523 526 data = cfunctions.decoder(fft_code, data)
524 527
525 528 datadec = data[:,:-self.nBaud+1]
526 529
527 530 return datadec
528 531
529 532 def convolutionInTime(self, data):
530 533
531 534 code = self.code[self.__profIndex]
532 535
533 536 for i in range(self.__nChannels):
534 537 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='valid')
535 538
536 539 return self.datadecTime
537 540
538 541 def convolutionByBlockInTime(self, data):
539 542 junk = numpy.lib.stride_tricks.as_strided(self.code, (self.times, self.code.size), (0, self.code.itemsize))
540 543 junk = junk.flatten()
541 544 code_block = numpy.reshape(junk, (self.nCode*self.times,self.nBaud))
542 545
543 546 for i in range(self.__nChannels):
544 547 for j in range(self.__nProfiles):
545 548 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='valid')
546 549
547 550 return self.datadecTime
548 551
549 552 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, times=None, osamp=None):
550 553
551 554 if code == None:
552 555 code = dataOut.code
553 556 else:
554 557 code = numpy.array(code).reshape(nCode,nBaud)
555 558
556 559
557 560
558 561 if not self.isConfig:
559 562
560 563 self.setup(code, dataOut.data.shape, times, osamp)
561 564
562 565 dataOut.code = code
563 566 dataOut.nCode = nCode
564 567 dataOut.nBaud = nBaud
565 568 dataOut.radarControllerHeaderObj.code = code
566 569 dataOut.radarControllerHeaderObj.nCode = nCode
567 570 dataOut.radarControllerHeaderObj.nBaud = nBaud
568 571
569 572 self.isConfig = True
570 573
571 574 if mode == 0:
572 575 datadec = self.convolutionInTime(dataOut.data)
573 576
574 577 if mode == 1:
575 578 datadec = self.convolutionInFreq(dataOut.data)
576 579
577 580 if mode == 2:
578 581 datadec = self.convolutionInFreqOpt(dataOut.data)
579 582
580 583 if mode == 3:
581 584 datadec = self.convolutionByBlockInTime(dataOut.data)
582 585
583 586 if not(self.__setValues):
584 587 dataOut.code = self.code
585 588 dataOut.nCode = self.nCode
586 589 dataOut.nBaud = self.nBaud
587 590 dataOut.radarControllerHeaderObj.code = self.code
588 591 dataOut.radarControllerHeaderObj.nCode = self.nCode
589 592 dataOut.radarControllerHeaderObj.nBaud = self.nBaud
590 593 self.__setValues = True
591 594
592 595 dataOut.data = datadec
593 596
594 597 dataOut.heightList = dataOut.heightList[0:self.ndatadec]
595 598
596 599 dataOut.flagDecodeData = True #asumo q la data no esta decodificada
597 600
598 601 if self.__profIndex == self.nCode-1:
599 602 self.__profIndex = 0
600 603 return 1
601 604
602 605 self.__profIndex += 1
603 606
604 607 return 1
605 608 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
606 609
607 610
608 611 class ProfileConcat(Operation):
609 612
610 613 isConfig = False
611 614 buffer = None
612 615
613 616 def __init__(self):
614 617
615 618 Operation.__init__(self)
616 619 self.profileIndex = 0
617 620
618 621 def reset(self):
619 622 self.buffer = numpy.zeros_like(self.buffer)
620 623 self.start_index = 0
621 624 self.times = 1
622 625
623 626 def setup(self, data, m, n=1):
624 627 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
625 628 self.profiles = data.shape[1]
626 629 self.start_index = 0
627 630 self.times = 1
628 631
629 632 def concat(self, data):
630 633
631 634 self.buffer[:,self.start_index:self.profiles*self.times] = data.copy()
632 635 self.start_index = self.start_index + self.profiles
633 636
634 637 def run(self, dataOut, m):
635 638
636 639 dataOut.flagNoData = True
637 640
638 641 if not self.isConfig:
639 642 self.setup(dataOut.data, m, 1)
640 643 self.isConfig = True
641 644
642 645 self.concat(dataOut.data)
643 646 self.times += 1
644 647 if self.times > m:
645 648 dataOut.data = self.buffer
646 649 self.reset()
647 650 dataOut.flagNoData = False
648 651 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
649 652 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
650 653 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * 5
651 654 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
652 655
653 656 class ProfileSelector(Operation):
654 657
655 658 profileIndex = None
656 659 # Tamanho total de los perfiles
657 660 nProfiles = None
658 661
659 662 def __init__(self):
660 663
661 664 Operation.__init__(self)
662 665 self.profileIndex = 0
663 666
664 667 def incIndex(self):
665 668 self.profileIndex += 1
666 669
667 670 if self.profileIndex >= self.nProfiles:
668 671 self.profileIndex = 0
669 672
670 673 def isProfileInRange(self, minIndex, maxIndex):
671 674
672 675 if self.profileIndex < minIndex:
673 676 return False
674 677
675 678 if self.profileIndex > maxIndex:
676 679 return False
677 680
678 681 return True
679 682
680 683 def isProfileInList(self, profileList):
681 684
682 685 if self.profileIndex not in profileList:
683 686 return False
684 687
685 688 return True
686 689
687 690 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False):
688 691
689 692 dataOut.flagNoData = True
690 693 self.nProfiles = dataOut.nProfiles
691 694
692 695 if byblock:
693 696
694 697 if profileList != None:
695 698 dataOut.data = dataOut.data[:,profileList,:]
696 699 pass
697 700 else:
698 701 pmin = profileRangeList[0]
699 702 pmax = profileRangeList[1]
700 703 dataOut.data = dataOut.data[:,pmin:pmax+1,:]
701 704 dataOut.flagNoData = False
702 705 self.profileIndex = 0
703 706 return 1
704 707
705 708 if profileList != None:
706 709 if self.isProfileInList(profileList):
707 710 dataOut.flagNoData = False
708 711
709 712 self.incIndex()
710 713 return 1
711 714
712 715
713 716 elif profileRangeList != None:
714 717 minIndex = profileRangeList[0]
715 718 maxIndex = profileRangeList[1]
716 719 if self.isProfileInRange(minIndex, maxIndex):
717 720 dataOut.flagNoData = False
718 721
719 722 self.incIndex()
720 723 return 1
721 724 elif beam != None: #beam is only for AMISR data
722 725 if self.isProfileInList(dataOut.beamRangeDict[beam]):
723 726 dataOut.flagNoData = False
724 727
725 728 self.incIndex()
726 729 return 1
727 730
728 731 else:
729 732 raise ValueError, "ProfileSelector needs profileList or profileRangeList"
730 733
731 734 return 0
732 735
733 736
734 737
735 738 class Reshaper(Operation):
736 739 def __init__(self):
737 740 Operation.__init__(self)
738 741 self.updateNewHeights = False
739 742
740 743 def run(self, dataOut, shape):
741 744 shape_tuple = tuple(shape)
742 745 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
743 746 dataOut.flagNoData = False
744 747
745 748 if not(self.updateNewHeights):
746 749 old_nheights = dataOut.nHeights
747 750 new_nheights = dataOut.data.shape[2]
748 751 factor = new_nheights / old_nheights
749 752 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
750 753 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * factor
751 754 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight) No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now