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