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