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