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