##// END OF EJS Templates
schain
RSS
Update2 for EW-Drifts
Percy Condor -
r1383:3e971ac8dea1
parent child
Show More
Side by Side Unified
Context file:
r1383:3e971ac8dea1 Loading diff...:

The requested changes are too big and content was truncated. Show full diff

Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
@@ -1,702 +1,703
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13
14 14
15 15 class SpectraPlot(Plot):
16 16 '''
17 17 Plot for Spectra data
18 18 '''
19 19
20 CODE = 'spc'
20 CODE = 'spc_moments'
21 21 colormap = 'jet'
22 22 plot_type = 'pcolor'
23 23 buffering = False
24 24
25 25 def setup(self):
26 26 self.nplots = len(self.data.channels)
27 27 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
28 28 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
29 29 self.height = 2.6 * self.nrows
30 30 self.cb_label = 'dB'
31 31 if self.showprofile:
32 32 self.width = 4 * self.ncols
33 33 else:
34 34 self.width = 3.5 * self.ncols
35 35 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
36 36 self.ylabel = 'Range [km]'
37 37
38 38 def update(self, dataOut):
39 39
40 40 data = {}
41 41 meta = {}
42 42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 43 data['spc'] = spc
44 44 data['rti'] = dataOut.getPower()
45 45 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 46 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
47 47 if self.CODE == 'spc_moments':
48 48 data['moments'] = dataOut.moments
49 49
50 50 return data, meta
51 51
52 52 def plot(self):
53 53 if self.xaxis == "frequency":
54 54 x = self.data.xrange[0]
55 55 self.xlabel = "Frequency (kHz)"
56 56 elif self.xaxis == "time":
57 57 x = self.data.xrange[1]
58 58 self.xlabel = "Time (ms)"
59 59 else:
60 60 x = self.data.xrange[2]
61 61 self.xlabel = "Velocity (m/s)"
62 62
63 63 if self.CODE == 'spc_moments':
64 64 x = self.data.xrange[2]
65 65 self.xlabel = "Velocity (m/s)"
66 66
67 67 self.titles = []
68 68
69 69 y = self.data.yrange
70 70 self.y = y
71 71
72 72 data = self.data[-1]
73 73 z = data['spc']
74
74 #self.CODE = 'spc_moments'
75 75 for n, ax in enumerate(self.axes):
76 76 noise = data['noise'][n]
77 print(n,self.CODE)
77 78 if self.CODE == 'spc_moments':
78 79 mean = data['moments'][n, 1]
79 80 if ax.firsttime:
80 81 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
81 82 self.xmin = self.xmin if self.xmin else -self.xmax
82 83 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
83 84 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
84 85 ax.plt = ax.pcolormesh(x, y, z[n].T,
85 86 vmin=self.zmin,
86 87 vmax=self.zmax,
87 88 cmap=plt.get_cmap(self.colormap)
88 89 )
89 90
90 91 if self.showprofile:
91 92 ax.plt_profile = self.pf_axes[n].plot(
92 93 data['rti'][n], y)[0]
93 94 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
94 95 color="k", linestyle="dashed", lw=1)[0]
95 96 if self.CODE == 'spc_moments':
96 97 ax.plt_mean = ax.plot(mean, y, color='k')[0]
97 98 else:
98 99 ax.plt.set_array(z[n].T.ravel())
99 100 if self.showprofile:
100 101 ax.plt_profile.set_data(data['rti'][n], y)
101 102 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
102 103 if self.CODE == 'spc_moments':
103 104 ax.plt_mean.set_data(mean, y)
104 105 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
105 106
106 107
107 108 class CrossSpectraPlot(Plot):
108 109
109 110 CODE = 'cspc'
110 111 colormap = 'jet'
111 112 plot_type = 'pcolor'
112 113 zmin_coh = None
113 114 zmax_coh = None
114 115 zmin_phase = None
115 116 zmax_phase = None
116 117
117 118 def setup(self):
118 119
119 120 self.ncols = 4
120 121 self.nplots = len(self.data.pairs) * 2
121 122 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
122 123 self.width = 3.1 * self.ncols
123 124 self.height = 2.6 * self.nrows
124 125 self.ylabel = 'Range [km]'
125 126 self.showprofile = False
126 127 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
127 128
128 129 def update(self, dataOut):
129 130
130 131 data = {}
131 132 meta = {}
132 133
133 134 spc = dataOut.data_spc
134 135 cspc = dataOut.data_cspc
135 136 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
136 137 meta['pairs'] = dataOut.pairsList
137 138
138 139 tmp = []
139 140
140 141 for n, pair in enumerate(meta['pairs']):
141 142 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
142 143 coh = numpy.abs(out)
143 144 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
144 145 tmp.append(coh)
145 146 tmp.append(phase)
146 147
147 148 data['cspc'] = numpy.array(tmp)
148 149
149 150 return data, meta
150 151
151 152 def plot(self):
152 153
153 154 if self.xaxis == "frequency":
154 155 x = self.data.xrange[0]
155 156 self.xlabel = "Frequency (kHz)"
156 157 elif self.xaxis == "time":
157 158 x = self.data.xrange[1]
158 159 self.xlabel = "Time (ms)"
159 160 else:
160 161 x = self.data.xrange[2]
161 162 self.xlabel = "Velocity (m/s)"
162 163
163 164 self.titles = []
164 165
165 166 y = self.data.yrange
166 167 self.y = y
167 168
168 169 data = self.data[-1]
169 170 cspc = data['cspc']
170 171
171 172 for n in range(len(self.data.pairs)):
172 173 pair = self.data.pairs[n]
173 174 coh = cspc[n*2]
174 175 phase = cspc[n*2+1]
175 176 ax = self.axes[2 * n]
176 177 if ax.firsttime:
177 178 ax.plt = ax.pcolormesh(x, y, coh.T,
178 179 vmin=0,
179 180 vmax=1,
180 181 cmap=plt.get_cmap(self.colormap_coh)
181 182 )
182 183 else:
183 184 ax.plt.set_array(coh.T.ravel())
184 185 self.titles.append(
185 186 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
186 187
187 188 ax = self.axes[2 * n + 1]
188 189 if ax.firsttime:
189 190 ax.plt = ax.pcolormesh(x, y, phase.T,
190 191 vmin=-180,
191 192 vmax=180,
192 193 cmap=plt.get_cmap(self.colormap_phase)
193 194 )
194 195 else:
195 196 ax.plt.set_array(phase.T.ravel())
196 197 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
197 198
198 199
199 200 class RTIPlot(Plot):
200 201 '''
201 202 Plot for RTI data
202 203 '''
203 204
204 205 CODE = 'rti'
205 206 colormap = 'jet'
206 207 plot_type = 'pcolorbuffer'
207 208
208 209 def setup(self):
209 210 self.xaxis = 'time'
210 211 self.ncols = 1
211 212 self.nrows = len(self.data.channels)
212 213 self.nplots = len(self.data.channels)
213 214 self.ylabel = 'Range [km]'
214 215 self.xlabel = 'Time'
215 216 self.cb_label = 'dB'
216 217 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
217 218 self.titles = ['{} Channel {}'.format(
218 219 self.CODE.upper(), x) for x in range(self.nrows)]
219 220
220 221 def update(self, dataOut):
221 222
222 223 data = {}
223 224 meta = {}
224 225 data['rti'] = dataOut.getPower()
225 226 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
226 227
227 228 return data, meta
228 229
229 230 def plot(self):
230 231 self.x = self.data.times
231 232 self.y = self.data.yrange
232 233 self.z = self.data[self.CODE]
233 234 self.z = numpy.ma.masked_invalid(self.z)
234 235
235 236 if self.decimation is None:
236 237 x, y, z = self.fill_gaps(self.x, self.y, self.z)
237 238 else:
238 239 x, y, z = self.fill_gaps(*self.decimate())
239 240
240 241 for n, ax in enumerate(self.axes):
241 242 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
242 243 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
243 244 data = self.data[-1]
244 245 if ax.firsttime:
245 246 ax.plt = ax.pcolormesh(x, y, z[n].T,
246 247 vmin=self.zmin,
247 248 vmax=self.zmax,
248 249 cmap=plt.get_cmap(self.colormap)
249 250 )
250 251 if self.showprofile:
251 252 ax.plot_profile = self.pf_axes[n].plot(
252 253 data['rti'][n], self.y)[0]
253 254 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
254 255 color="k", linestyle="dashed", lw=1)[0]
255 256 else:
256 257 ax.collections.remove(ax.collections[0])
257 258 ax.plt = ax.pcolormesh(x, y, z[n].T,
258 259 vmin=self.zmin,
259 260 vmax=self.zmax,
260 261 cmap=plt.get_cmap(self.colormap)
261 262 )
262 263 if self.showprofile:
263 264 ax.plot_profile.set_data(data['rti'][n], self.y)
264 265 ax.plot_noise.set_data(numpy.repeat(
265 266 data['noise'][n], len(self.y)), self.y)
266 267
267 268
268 269 class CoherencePlot(RTIPlot):
269 270 '''
270 271 Plot for Coherence data
271 272 '''
272 273
273 274 CODE = 'coh'
274 275
275 276 def setup(self):
276 277 self.xaxis = 'time'
277 278 self.ncols = 1
278 279 self.nrows = len(self.data.pairs)
279 280 self.nplots = len(self.data.pairs)
280 281 self.ylabel = 'Range [km]'
281 282 self.xlabel = 'Time'
282 283 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
283 284 if self.CODE == 'coh':
284 285 self.cb_label = ''
285 286 self.titles = [
286 287 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
287 288 else:
288 289 self.cb_label = 'Degrees'
289 290 self.titles = [
290 291 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
291 292
292 293 def update(self, dataOut):
293 294
294 295 data = {}
295 296 meta = {}
296 297 data['coh'] = dataOut.getCoherence()
297 298 meta['pairs'] = dataOut.pairsList
298 299
299 300 return data, meta
300 301
301 302 class PhasePlot(CoherencePlot):
302 303 '''
303 304 Plot for Phase map data
304 305 '''
305 306
306 307 CODE = 'phase'
307 308 colormap = 'seismic'
308 309
309 310 def update(self, dataOut):
310 311
311 312 data = {}
312 313 meta = {}
313 314 data['phase'] = dataOut.getCoherence(phase=True)
314 315 meta['pairs'] = dataOut.pairsList
315 316
316 317 return data, meta
317 318
318 319 class NoisePlot(Plot):
319 320 '''
320 321 Plot for noise
321 322 '''
322 323
323 324 CODE = 'noise'
324 325 plot_type = 'scatterbuffer'
325 326
326 327 def setup(self):
327 328 self.xaxis = 'time'
328 329 self.ncols = 1
329 330 self.nrows = 1
330 331 self.nplots = 1
331 332 self.ylabel = 'Intensity [dB]'
332 333 self.xlabel = 'Time'
333 334 self.titles = ['Noise']
334 335 self.colorbar = False
335 336 self.plots_adjust.update({'right': 0.85 })
336 337
337 338 def update(self, dataOut):
338 339
339 340 data = {}
340 341 meta = {}
341 342 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
342 343 meta['yrange'] = numpy.array([])
343 344
344 345 return data, meta
345 346
346 347 def plot(self):
347 348
348 349 x = self.data.times
349 350 xmin = self.data.min_time
350 351 xmax = xmin + self.xrange * 60 * 60
351 352 Y = self.data['noise']
352 353
353 354 if self.axes[0].firsttime:
354 355 self.ymin = numpy.nanmin(Y) - 5
355 356 self.ymax = numpy.nanmax(Y) + 5
356 357 for ch in self.data.channels:
357 358 y = Y[ch]
358 359 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
359 360 plt.legend(bbox_to_anchor=(1.18, 1.0))
360 361 else:
361 362 for ch in self.data.channels:
362 363 y = Y[ch]
363 364 self.axes[0].lines[ch].set_data(x, y)
364 365
365 366
366 367 class PowerProfilePlot(Plot):
367 368
368 369 CODE = 'pow_profile'
369 370 plot_type = 'scatter'
370 371
371 372 def setup(self):
372 373
373 374 self.ncols = 1
374 375 self.nrows = 1
375 376 self.nplots = 1
376 377 self.height = 4
377 378 self.width = 3
378 379 self.ylabel = 'Range [km]'
379 380 self.xlabel = 'Intensity [dB]'
380 381 self.titles = ['Power Profile']
381 382 self.colorbar = False
382 383
383 384 def update(self, dataOut):
384 385
385 386 data = {}
386 387 meta = {}
387 388 data[self.CODE] = dataOut.getPower()
388 389
389 390 return data, meta
390 391
391 392 def plot(self):
392 393
393 394 y = self.data.yrange
394 395 self.y = y
395 396
396 397 x = self.data[-1][self.CODE]
397 398
398 399 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
399 400 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
400 401
401 402 if self.axes[0].firsttime:
402 403 for ch in self.data.channels:
403 404 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
404 405 plt.legend()
405 406 else:
406 407 for ch in self.data.channels:
407 408 self.axes[0].lines[ch].set_data(x[ch], y)
408 409
409 410
410 411 class SpectraCutPlot(Plot):
411 412
412 413 CODE = 'spc_cut'
413 414 plot_type = 'scatter'
414 415 buffering = False
415 416
416 417 def setup(self):
417 418
418 419 self.nplots = len(self.data.channels)
419 420 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
420 421 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
421 422 self.width = 3.4 * self.ncols + 1.5
422 423 self.height = 3 * self.nrows
423 424 self.ylabel = 'Power [dB]'
424 425 self.colorbar = False
425 426 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
426 427
427 428 def update(self, dataOut):
428 429
429 430 data = {}
430 431 meta = {}
431 432 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
432 433 data['spc'] = spc
433 434 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
434 435
435 436 return data, meta
436 437
437 438 def plot(self):
438 439 if self.xaxis == "frequency":
439 440 x = self.data.xrange[0][1:]
440 441 self.xlabel = "Frequency (kHz)"
441 442 elif self.xaxis == "time":
442 443 x = self.data.xrange[1]
443 444 self.xlabel = "Time (ms)"
444 445 else:
445 446 x = self.data.xrange[2]
446 447 self.xlabel = "Velocity (m/s)"
447 448
448 449 self.titles = []
449 450
450 451 y = self.data.yrange
451 452 z = self.data[-1]['spc']
452 453
453 454 if self.height_index:
454 455 index = numpy.array(self.height_index)
455 456 else:
456 457 index = numpy.arange(0, len(y), int((len(y))/9))
457 458
458 459 for n, ax in enumerate(self.axes):
459 460 if ax.firsttime:
460 461 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
461 462 self.xmin = self.xmin if self.xmin else -self.xmax
462 463 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
463 464 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
464 465 ax.plt = ax.plot(x, z[n, :, index].T)
465 466 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
466 467 self.figures[0].legend(ax.plt, labels, loc='center right')
467 468 else:
468 469 for i, line in enumerate(ax.plt):
469 470 line.set_data(x, z[n, :, index[i]])
470 471 self.titles.append('CH {}'.format(n))
471 472
472 473
473 474 class BeaconPhase(Plot):
474 475
475 476 __isConfig = None
476 477 __nsubplots = None
477 478
478 479 PREFIX = 'beacon_phase'
479 480
480 481 def __init__(self):
481 482 Plot.__init__(self)
482 483 self.timerange = 24*60*60
483 484 self.isConfig = False
484 485 self.__nsubplots = 1
485 486 self.counter_imagwr = 0
486 487 self.WIDTH = 800
487 488 self.HEIGHT = 400
488 489 self.WIDTHPROF = 120
489 490 self.HEIGHTPROF = 0
490 491 self.xdata = None
491 492 self.ydata = None
492 493
493 494 self.PLOT_CODE = BEACON_CODE
494 495
495 496 self.FTP_WEI = None
496 497 self.EXP_CODE = None
497 498 self.SUB_EXP_CODE = None
498 499 self.PLOT_POS = None
499 500
500 501 self.filename_phase = None
501 502
502 503 self.figfile = None
503 504
504 505 self.xmin = None
505 506 self.xmax = None
506 507
507 508 def getSubplots(self):
508 509
509 510 ncol = 1
510 511 nrow = 1
511 512
512 513 return nrow, ncol
513 514
514 515 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
515 516
516 517 self.__showprofile = showprofile
517 518 self.nplots = nplots
518 519
519 520 ncolspan = 7
520 521 colspan = 6
521 522 self.__nsubplots = 2
522 523
523 524 self.createFigure(id = id,
524 525 wintitle = wintitle,
525 526 widthplot = self.WIDTH+self.WIDTHPROF,
526 527 heightplot = self.HEIGHT+self.HEIGHTPROF,
527 528 show=show)
528 529
529 530 nrow, ncol = self.getSubplots()
530 531
531 532 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
532 533
533 534 def save_phase(self, filename_phase):
534 535 f = open(filename_phase,'w+')
535 536 f.write('\n\n')
536 537 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
537 538 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
538 539 f.close()
539 540
540 541 def save_data(self, filename_phase, data, data_datetime):
541 542 f=open(filename_phase,'a')
542 543 timetuple_data = data_datetime.timetuple()
543 544 day = str(timetuple_data.tm_mday)
544 545 month = str(timetuple_data.tm_mon)
545 546 year = str(timetuple_data.tm_year)
546 547 hour = str(timetuple_data.tm_hour)
547 548 minute = str(timetuple_data.tm_min)
548 549 second = str(timetuple_data.tm_sec)
549 550 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
550 551 f.close()
551 552
552 553 def plot(self):
553 554 log.warning('TODO: Not yet implemented...')
554 555
555 556 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
556 557 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
557 558 timerange=None,
558 559 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
559 560 server=None, folder=None, username=None, password=None,
560 561 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
561 562
562 563 if dataOut.flagNoData:
563 564 return dataOut
564 565
565 566 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
566 567 return
567 568
568 569 if pairsList == None:
569 570 pairsIndexList = dataOut.pairsIndexList[:10]
570 571 else:
571 572 pairsIndexList = []
572 573 for pair in pairsList:
573 574 if pair not in dataOut.pairsList:
574 575 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
575 576 pairsIndexList.append(dataOut.pairsList.index(pair))
576 577
577 578 if pairsIndexList == []:
578 579 return
579 580
580 581 # if len(pairsIndexList) > 4:
581 582 # pairsIndexList = pairsIndexList[0:4]
582 583
583 584 hmin_index = None
584 585 hmax_index = None
585 586
586 587 if hmin != None and hmax != None:
587 588 indexes = numpy.arange(dataOut.nHeights)
588 589 hmin_list = indexes[dataOut.heightList >= hmin]
589 590 hmax_list = indexes[dataOut.heightList <= hmax]
590 591
591 592 if hmin_list.any():
592 593 hmin_index = hmin_list[0]
593 594
594 595 if hmax_list.any():
595 596 hmax_index = hmax_list[-1]+1
596 597
597 598 x = dataOut.getTimeRange()
598 599
599 600 thisDatetime = dataOut.datatime
600 601
601 602 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
602 603 xlabel = "Local Time"
603 604 ylabel = "Phase (degrees)"
604 605
605 606 update_figfile = False
606 607
607 608 nplots = len(pairsIndexList)
608 609 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
609 610 phase_beacon = numpy.zeros(len(pairsIndexList))
610 611 for i in range(nplots):
611 612 pair = dataOut.pairsList[pairsIndexList[i]]
612 613 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
613 614 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
614 615 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
615 616 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
616 617 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
617 618
618 619 if dataOut.beacon_heiIndexList:
619 620 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
620 621 else:
621 622 phase_beacon[i] = numpy.average(phase)
622 623
623 624 if not self.isConfig:
624 625
625 626 nplots = len(pairsIndexList)
626 627
627 628 self.setup(id=id,
628 629 nplots=nplots,
629 630 wintitle=wintitle,
630 631 showprofile=showprofile,
631 632 show=show)
632 633
633 634 if timerange != None:
634 635 self.timerange = timerange
635 636
636 637 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
637 638
638 639 if ymin == None: ymin = 0
639 640 if ymax == None: ymax = 360
640 641
641 642 self.FTP_WEI = ftp_wei
642 643 self.EXP_CODE = exp_code
643 644 self.SUB_EXP_CODE = sub_exp_code
644 645 self.PLOT_POS = plot_pos
645 646
646 647 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
647 648 self.isConfig = True
648 649 self.figfile = figfile
649 650 self.xdata = numpy.array([])
650 651 self.ydata = numpy.array([])
651 652
652 653 update_figfile = True
653 654
654 655 #open file beacon phase
655 656 path = '%s%03d' %(self.PREFIX, self.id)
656 657 beacon_file = os.path.join(path,'%s.txt'%self.name)
657 658 self.filename_phase = os.path.join(figpath,beacon_file)
658 659 #self.save_phase(self.filename_phase)
659 660
660 661
661 662 #store data beacon phase
662 663 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
663 664
664 665 self.setWinTitle(title)
665 666
666 667
667 668 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
668 669
669 670 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
670 671
671 672 axes = self.axesList[0]
672 673
673 674 self.xdata = numpy.hstack((self.xdata, x[0:1]))
674 675
675 676 if len(self.ydata)==0:
676 677 self.ydata = phase_beacon.reshape(-1,1)
677 678 else:
678 679 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
679 680
680 681
681 682 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
682 683 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
683 684 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
684 685 XAxisAsTime=True, grid='both'
685 686 )
686 687
687 688 self.draw()
688 689
689 690 if dataOut.ltctime >= self.xmax:
690 691 self.counter_imagwr = wr_period
691 692 self.isConfig = False
692 693 update_figfile = True
693 694
694 695 self.save(figpath=figpath,
695 696 figfile=figfile,
696 697 save=save,
697 698 ftp=ftp,
698 699 wr_period=wr_period,
699 700 thisDatetime=thisDatetime,
700 701 update_figfile=update_figfile)
701 702
702 703 return dataOut
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
@@ -1,1625 +1,1631
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62 62 for channel in channelList:
63 63 if channel not in self.dataOut.channelList:
64 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
65 65
66 66 index = self.dataOut.channelList.index(channel)
67 67 channelIndexList.append(index)
68 68 self.selectChannelsByIndex(channelIndexList)
69 69 return self.dataOut
70 70
71 71 def selectChannelsByIndex(self, channelIndexList):
72 72 """
73 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
74 74
75 75 Input:
76 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
77 77
78 78 Affected:
79 79 self.dataOut.data
80 80 self.dataOut.channelIndexList
81 81 self.dataOut.nChannels
82 82 self.dataOut.m_ProcessingHeader.totalSpectra
83 83 self.dataOut.systemHeaderObj.numChannels
84 84 self.dataOut.m_ProcessingHeader.blockSize
85 85
86 86 Return:
87 87 None
88 88 """
89 89
90 90 for channelIndex in channelIndexList:
91 91 if channelIndex not in self.dataOut.channelIndexList:
92 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
93 93
94 94 if self.dataOut.type == 'Voltage':
95 95 if self.dataOut.flagDataAsBlock:
96 96 """
97 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
98 98 """
99 99 data = self.dataOut.data[channelIndexList,:,:]
100 100 else:
101 101 data = self.dataOut.data[channelIndexList,:]
102 102
103 103 self.dataOut.data = data
104 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
105 105 self.dataOut.channelList = range(len(channelIndexList))
106 106
107 107 elif self.dataOut.type == 'Spectra':
108 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
109 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
110 110
111 111 self.dataOut.data_spc = data_spc
112 112 self.dataOut.data_dc = data_dc
113 113
114 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 self.dataOut.channelList = range(len(channelIndexList))
116 116 self.__selectPairsByChannel(channelIndexList)
117 117
118 118 return 1
119 119
120 120 def __selectPairsByChannel(self, channelList=None):
121 121
122 122 if channelList == None:
123 123 return
124 124
125 125 pairsIndexListSelected = []
126 126 for pairIndex in self.dataOut.pairsIndexList:
127 127 # First pair
128 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
129 129 continue
130 130 # Second pair
131 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
132 132 continue
133 133
134 134 pairsIndexListSelected.append(pairIndex)
135 135
136 136 if not pairsIndexListSelected:
137 137 self.dataOut.data_cspc = None
138 138 self.dataOut.pairsList = []
139 139 return
140 140
141 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
142 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
143 143 for i in pairsIndexListSelected]
144 144
145 145 return
146 146
147 147 class selectHeights(Operation):
148 148
149 149 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
150 150 """
151 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
152 152 minHei <= height <= maxHei
153 153
154 154 Input:
155 155 minHei : valor minimo de altura a considerar
156 156 maxHei : valor maximo de altura a considerar
157 157
158 158 Affected:
159 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
160 160
161 161 Return:
162 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
163 163 """
164 164
165 165 self.dataOut = dataOut
166 166
167 if minHei and maxHei:
167 #if minHei and maxHei:
168 if 1:
169 if minHei == None:
170 minHei = self.dataOut.heightList[0]
171
172 if maxHei == None:
173 maxHei = self.dataOut.heightList[-1]
168 174
169 175 if (minHei < self.dataOut.heightList[0]):
170 176 minHei = self.dataOut.heightList[0]
171 177
172 178 if (maxHei > self.dataOut.heightList[-1]):
173 179 maxHei = self.dataOut.heightList[-1]
174 180
175 181 minIndex = 0
176 182 maxIndex = 0
177 183 heights = self.dataOut.heightList
178 184
179 185 inda = numpy.where(heights >= minHei)
180 186 indb = numpy.where(heights <= maxHei)
181 187
182 188 try:
183 189 minIndex = inda[0][0]
184 190 except:
185 191 minIndex = 0
186 192
187 193 try:
188 194 maxIndex = indb[0][-1]
189 195 except:
190 196 maxIndex = len(heights)
191 197
192 198 self.selectHeightsByIndex(minIndex, maxIndex)
193 199
194 200 return self.dataOut
195 201
196 202 def selectHeightsByIndex(self, minIndex, maxIndex):
197 203 """
198 204 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
199 205 minIndex <= index <= maxIndex
200 206
201 207 Input:
202 208 minIndex : valor de indice minimo de altura a considerar
203 209 maxIndex : valor de indice maximo de altura a considerar
204 210
205 211 Affected:
206 212 self.dataOut.data
207 213 self.dataOut.heightList
208 214
209 215 Return:
210 216 1 si el metodo se ejecuto con exito caso contrario devuelve 0
211 217 """
212 218
213 219 if self.dataOut.type == 'Voltage':
214 220 if (minIndex < 0) or (minIndex > maxIndex):
215 221 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
216 222
217 223 if (maxIndex >= self.dataOut.nHeights):
218 224 maxIndex = self.dataOut.nHeights
219 225
220 226 #voltage
221 227 if self.dataOut.flagDataAsBlock:
222 228 """
223 229 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
224 230 """
225 231 data = self.dataOut.data[:,:, minIndex:maxIndex]
226 232 else:
227 233 data = self.dataOut.data[:, minIndex:maxIndex]
228 234
229 235 # firstHeight = self.dataOut.heightList[minIndex]
230 236
231 237 self.dataOut.data = data
232 238 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
233 239
234 240 if self.dataOut.nHeights <= 1:
235 241 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
236 242 elif self.dataOut.type == 'Spectra':
237 243 if (minIndex < 0) or (minIndex > maxIndex):
238 244 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
239 245 minIndex, maxIndex))
240 246
241 247 if (maxIndex >= self.dataOut.nHeights):
242 248 maxIndex = self.dataOut.nHeights - 1
243 249
244 250 # Spectra
245 251 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
246 252
247 253 data_cspc = None
248 254 if self.dataOut.data_cspc is not None:
249 255 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
250 256
251 257 data_dc = None
252 258 if self.dataOut.data_dc is not None:
253 259 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
254 260
255 261 self.dataOut.data_spc = data_spc
256 262 self.dataOut.data_cspc = data_cspc
257 263 self.dataOut.data_dc = data_dc
258 264
259 265 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
260 266
261 267 return 1
262 268
263 269
264 270 class filterByHeights(Operation):
265 271
266 272 def run(self, dataOut, window):
267 273
268 274 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
269 275
270 276 if window == None:
271 277 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
272 278
273 279 newdelta = deltaHeight * window
274 280 r = dataOut.nHeights % window
275 281 newheights = (dataOut.nHeights-r)/window
276 282
277 283 if newheights <= 1:
278 284 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
279 285
280 286 if dataOut.flagDataAsBlock:
281 287 """
282 288 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
283 289 """
284 290 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
285 291 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
286 292 buffer = numpy.sum(buffer,3)
287 293
288 294 else:
289 295 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
290 296 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
291 297 buffer = numpy.sum(buffer,2)
292 298
293 299 dataOut.data = buffer
294 300 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
295 301 dataOut.windowOfFilter = window
296 302
297 303 return dataOut
298 304
299 305
300 306 class setH0(Operation):
301 307
302 308 def run(self, dataOut, h0, deltaHeight = None):
303 309
304 310 if not deltaHeight:
305 311 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
306 312
307 313 nHeights = dataOut.nHeights
308 314
309 315 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
310 316
311 317 dataOut.heightList = newHeiRange
312 318
313 319 return dataOut
314 320
315 321
316 322 class deFlip(Operation):
317 323
318 324 def run(self, dataOut, channelList = []):
319 325
320 326 data = dataOut.data.copy()
321 327
322 328 if dataOut.flagDataAsBlock:
323 329 flip = self.flip
324 330 profileList = list(range(dataOut.nProfiles))
325 331
326 332 if not channelList:
327 333 for thisProfile in profileList:
328 334 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
329 335 flip *= -1.0
330 336 else:
331 337 for thisChannel in channelList:
332 338 if thisChannel not in dataOut.channelList:
333 339 continue
334 340
335 341 for thisProfile in profileList:
336 342 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
337 343 flip *= -1.0
338 344
339 345 self.flip = flip
340 346
341 347 else:
342 348 if not channelList:
343 349 data[:,:] = data[:,:]*self.flip
344 350 else:
345 351 for thisChannel in channelList:
346 352 if thisChannel not in dataOut.channelList:
347 353 continue
348 354
349 355 data[thisChannel,:] = data[thisChannel,:]*self.flip
350 356
351 357 self.flip *= -1.
352 358
353 359 dataOut.data = data
354 360
355 361 return dataOut
356 362
357 363
358 364 class setAttribute(Operation):
359 365 '''
360 366 Set an arbitrary attribute(s) to dataOut
361 367 '''
362 368
363 369 def __init__(self):
364 370
365 371 Operation.__init__(self)
366 372 self._ready = False
367 373
368 374 def run(self, dataOut, **kwargs):
369 375
370 376 for key, value in kwargs.items():
371 377 setattr(dataOut, key, value)
372 378
373 379 return dataOut
374 380
375 381
376 382 @MPDecorator
377 383 class printAttribute(Operation):
378 384 '''
379 385 Print an arbitrary attribute of dataOut
380 386 '''
381 387
382 388 def __init__(self):
383 389
384 390 Operation.__init__(self)
385 391
386 392 def run(self, dataOut, attributes):
387 393
388 394 if isinstance(attributes, str):
389 395 attributes = [attributes]
390 396 for attr in attributes:
391 397 if hasattr(dataOut, attr):
392 398 log.log(getattr(dataOut, attr), attr)
393 399
394 400
395 401 class interpolateHeights(Operation):
396 402
397 403 def run(self, dataOut, topLim, botLim):
398 404 #69 al 72 para julia
399 405 #82-84 para meteoros
400 406 if len(numpy.shape(dataOut.data))==2:
401 407 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
402 408 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
403 409 #dataOut.data[:,botLim:limSup+1] = sampInterp
404 410 dataOut.data[:,botLim:topLim+1] = sampInterp
405 411 else:
406 412 nHeights = dataOut.data.shape[2]
407 413 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
408 414 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
409 415 f = interpolate.interp1d(x, y, axis = 2)
410 416 xnew = numpy.arange(botLim,topLim+1)
411 417 ynew = f(xnew)
412 418 dataOut.data[:,:,botLim:topLim+1] = ynew
413 419
414 420 return dataOut
415 421
416 422
417 423 class CohInt(Operation):
418 424
419 425 isConfig = False
420 426 __profIndex = 0
421 427 __byTime = False
422 428 __initime = None
423 429 __lastdatatime = None
424 430 __integrationtime = None
425 431 __buffer = None
426 432 __bufferStride = []
427 433 __dataReady = False
428 434 __profIndexStride = 0
429 435 __dataToPutStride = False
430 436 n = None
431 437
432 438 def __init__(self, **kwargs):
433 439
434 440 Operation.__init__(self, **kwargs)
435 441
436 442 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
437 443 """
438 444 Set the parameters of the integration class.
439 445
440 446 Inputs:
441 447
442 448 n : Number of coherent integrations
443 449 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
444 450 overlapping :
445 451 """
446 452
447 453 self.__initime = None
448 454 self.__lastdatatime = 0
449 455 self.__buffer = None
450 456 self.__dataReady = False
451 457 self.byblock = byblock
452 458 self.stride = stride
453 459
454 460 if n == None and timeInterval == None:
455 461 raise ValueError("n or timeInterval should be specified ...")
456 462
457 463 if n != None:
458 464 self.n = n
459 465 self.__byTime = False
460 466 else:
461 467 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
462 468 self.n = 9999
463 469 self.__byTime = True
464 470
465 471 if overlapping:
466 472 self.__withOverlapping = True
467 473 self.__buffer = None
468 474 else:
469 475 self.__withOverlapping = False
470 476 self.__buffer = 0
471 477
472 478 self.__profIndex = 0
473 479
474 480 def putData(self, data):
475 481
476 482 """
477 483 Add a profile to the __buffer and increase in one the __profileIndex
478 484
479 485 """
480 486
481 487 if not self.__withOverlapping:
482 488 self.__buffer += data.copy()
483 489 self.__profIndex += 1
484 490 return
485 491
486 492 #Overlapping data
487 493 nChannels, nHeis = data.shape
488 494 data = numpy.reshape(data, (1, nChannels, nHeis))
489 495
490 496 #If the buffer is empty then it takes the data value
491 497 if self.__buffer is None:
492 498 self.__buffer = data
493 499 self.__profIndex += 1
494 500 return
495 501
496 502 #If the buffer length is lower than n then stakcing the data value
497 503 if self.__profIndex < self.n:
498 504 self.__buffer = numpy.vstack((self.__buffer, data))
499 505 self.__profIndex += 1
500 506 return
501 507
502 508 #If the buffer length is equal to n then replacing the last buffer value with the data value
503 509 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
504 510 self.__buffer[self.n-1] = data
505 511 self.__profIndex = self.n
506 512 return
507 513
508 514
509 515 def pushData(self):
510 516 """
511 517 Return the sum of the last profiles and the profiles used in the sum.
512 518
513 519 Affected:
514 520
515 521 self.__profileIndex
516 522
517 523 """
518 524
519 525 if not self.__withOverlapping:
520 526 data = self.__buffer
521 527 n = self.__profIndex
522 528
523 529 self.__buffer = 0
524 530 self.__profIndex = 0
525 531
526 532 return data, n
527 533
528 534 #Integration with Overlapping
529 535 data = numpy.sum(self.__buffer, axis=0)
530 536 # print data
531 537 # raise
532 538 n = self.__profIndex
533 539
534 540 return data, n
535 541
536 542 def byProfiles(self, data):
537 543
538 544 self.__dataReady = False
539 545 avgdata = None
540 546 # n = None
541 547 # print data
542 548 # raise
543 549 self.putData(data)
544 550
545 551 if self.__profIndex == self.n:
546 552 avgdata, n = self.pushData()
547 553 self.__dataReady = True
548 554
549 555 return avgdata
550 556
551 557 def byTime(self, data, datatime):
552 558
553 559 self.__dataReady = False
554 560 avgdata = None
555 561 n = None
556 562
557 563 self.putData(data)
558 564
559 565 if (datatime - self.__initime) >= self.__integrationtime:
560 566 avgdata, n = self.pushData()
561 567 self.n = n
562 568 self.__dataReady = True
563 569
564 570 return avgdata
565 571
566 572 def integrateByStride(self, data, datatime):
567 573 # print data
568 574 if self.__profIndex == 0:
569 575 self.__buffer = [[data.copy(), datatime]]
570 576 else:
571 577 self.__buffer.append([data.copy(),datatime])
572 578 self.__profIndex += 1
573 579 self.__dataReady = False
574 580
575 581 if self.__profIndex == self.n * self.stride :
576 582 self.__dataToPutStride = True
577 583 self.__profIndexStride = 0
578 584 self.__profIndex = 0
579 585 self.__bufferStride = []
580 586 for i in range(self.stride):
581 587 current = self.__buffer[i::self.stride]
582 588 data = numpy.sum([t[0] for t in current], axis=0)
583 589 avgdatatime = numpy.average([t[1] for t in current])
584 590 # print data
585 591 self.__bufferStride.append((data, avgdatatime))
586 592
587 593 if self.__dataToPutStride:
588 594 self.__dataReady = True
589 595 self.__profIndexStride += 1
590 596 if self.__profIndexStride == self.stride:
591 597 self.__dataToPutStride = False
592 598 # print self.__bufferStride[self.__profIndexStride - 1]
593 599 # raise
594 600 return self.__bufferStride[self.__profIndexStride - 1]
595 601
596 602
597 603 return None, None
598 604
599 605 def integrate(self, data, datatime=None):
600 606
601 607 if self.__initime == None:
602 608 self.__initime = datatime
603 609
604 610 if self.__byTime:
605 611 avgdata = self.byTime(data, datatime)
606 612 else:
607 613 avgdata = self.byProfiles(data)
608 614
609 615
610 616 self.__lastdatatime = datatime
611 617
612 618 if avgdata is None:
613 619 return None, None
614 620
615 621 avgdatatime = self.__initime
616 622
617 623 deltatime = datatime - self.__lastdatatime
618 624
619 625 if not self.__withOverlapping:
620 626 self.__initime = datatime
621 627 else:
622 628 self.__initime += deltatime
623 629
624 630 return avgdata, avgdatatime
625 631
626 632 def integrateByBlock(self, dataOut):
627 633
628 634 times = int(dataOut.data.shape[1]/self.n)
629 635 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
630 636
631 637 id_min = 0
632 638 id_max = self.n
633 639
634 640 for i in range(times):
635 641 junk = dataOut.data[:,id_min:id_max,:]
636 642 avgdata[:,i,:] = junk.sum(axis=1)
637 643 id_min += self.n
638 644 id_max += self.n
639 645
640 646 timeInterval = dataOut.ippSeconds*self.n
641 647 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
642 648 self.__dataReady = True
643 649 return avgdata, avgdatatime
644 650
645 651 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
646 652
647 653 if not self.isConfig:
648 654 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
649 655 self.isConfig = True
650 656
651 657 if dataOut.flagDataAsBlock:
652 658 """
653 659 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
654 660 """
655 661 avgdata, avgdatatime = self.integrateByBlock(dataOut)
656 662 dataOut.nProfiles /= self.n
657 663 else:
658 664 if stride is None:
659 665 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
660 666 else:
661 667 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
662 668
663 669
664 670 # dataOut.timeInterval *= n
665 671 dataOut.flagNoData = True
666 672
667 673 if self.__dataReady:
668 674 dataOut.data = avgdata
669 675 if not dataOut.flagCohInt:
670 676 dataOut.nCohInt *= self.n
671 677 dataOut.flagCohInt = True
672 678 dataOut.utctime = avgdatatime
673 679 # print avgdata, avgdatatime
674 680 # raise
675 681 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
676 682 dataOut.flagNoData = False
677 683 return dataOut
678 684
679 685 class Decoder(Operation):
680 686
681 687 isConfig = False
682 688 __profIndex = 0
683 689
684 690 code = None
685 691
686 692 nCode = None
687 693 nBaud = None
688 694
689 695 def __init__(self, **kwargs):
690 696
691 697 Operation.__init__(self, **kwargs)
692 698
693 699 self.times = None
694 700 self.osamp = None
695 701 # self.__setValues = False
696 702 self.isConfig = False
697 703 self.setupReq = False
698 704 def setup(self, code, osamp, dataOut):
699 705
700 706 self.__profIndex = 0
701 707
702 708 self.code = code
703 709
704 710 self.nCode = len(code)
705 711 self.nBaud = len(code[0])
706 712
707 713 if (osamp != None) and (osamp >1):
708 714 self.osamp = osamp
709 715 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
710 716 self.nBaud = self.nBaud*self.osamp
711 717
712 718 self.__nChannels = dataOut.nChannels
713 719 self.__nProfiles = dataOut.nProfiles
714 720 self.__nHeis = dataOut.nHeights
715 721
716 722 if self.__nHeis < self.nBaud:
717 723 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
718 724
719 725 #Frequency
720 726 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
721 727
722 728 __codeBuffer[:,0:self.nBaud] = self.code
723 729
724 730 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
725 731
726 732 if dataOut.flagDataAsBlock:
727 733
728 734 self.ndatadec = self.__nHeis #- self.nBaud + 1
729 735
730 736 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
731 737
732 738 else:
733 739
734 740 #Time
735 741 self.ndatadec = self.__nHeis #- self.nBaud + 1
736 742
737 743 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
738 744
739 745 def __convolutionInFreq(self, data):
740 746
741 747 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
742 748
743 749 fft_data = numpy.fft.fft(data, axis=1)
744 750
745 751 conv = fft_data*fft_code
746 752
747 753 data = numpy.fft.ifft(conv,axis=1)
748 754
749 755 return data
750 756
751 757 def __convolutionInFreqOpt(self, data):
752 758
753 759 raise NotImplementedError
754 760
755 761 def __convolutionInTime(self, data):
756 762
757 763 code = self.code[self.__profIndex]
758 764 for i in range(self.__nChannels):
759 765 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
760 766
761 767 return self.datadecTime
762 768
763 769 def __convolutionByBlockInTime(self, data):
764 770
765 771 repetitions = int(self.__nProfiles / self.nCode)
766 772 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
767 773 junk = junk.flatten()
768 774 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
769 775 profilesList = range(self.__nProfiles)
770 776
771 777 for i in range(self.__nChannels):
772 778 for j in profilesList:
773 779 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
774 780 return self.datadecTime
775 781
776 782 def __convolutionByBlockInFreq(self, data):
777 783
778 784 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
779 785
780 786
781 787 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
782 788
783 789 fft_data = numpy.fft.fft(data, axis=2)
784 790
785 791 conv = fft_data*fft_code
786 792
787 793 data = numpy.fft.ifft(conv,axis=2)
788 794
789 795 return data
790 796
791 797
792 798 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
793 799
794 800 if dataOut.flagDecodeData:
795 801 print("This data is already decoded, recoding again ...")
796 802
797 803 if not self.isConfig:
798 804
799 805 if code is None:
800 806 if dataOut.code is None:
801 807 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
802 808
803 809 code = dataOut.code
804 810 else:
805 811 code = numpy.array(code).reshape(nCode,nBaud)
806 812 self.setup(code, osamp, dataOut)
807 813
808 814 self.isConfig = True
809 815
810 816 if mode == 3:
811 817 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
812 818
813 819 if times != None:
814 820 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
815 821
816 822 if self.code is None:
817 823 print("Fail decoding: Code is not defined.")
818 824 return
819 825
820 826 self.__nProfiles = dataOut.nProfiles
821 827 datadec = None
822 828
823 829 if mode == 3:
824 830 mode = 0
825 831
826 832 if dataOut.flagDataAsBlock:
827 833 """
828 834 Decoding when data have been read as block,
829 835 """
830 836
831 837 if mode == 0:
832 838 datadec = self.__convolutionByBlockInTime(dataOut.data)
833 839 if mode == 1:
834 840 datadec = self.__convolutionByBlockInFreq(dataOut.data)
835 841 else:
836 842 """
837 843 Decoding when data have been read profile by profile
838 844 """
839 845 if mode == 0:
840 846 datadec = self.__convolutionInTime(dataOut.data)
841 847
842 848 if mode == 1:
843 849 datadec = self.__convolutionInFreq(dataOut.data)
844 850
845 851 if mode == 2:
846 852 datadec = self.__convolutionInFreqOpt(dataOut.data)
847 853
848 854 if datadec is None:
849 855 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
850 856
851 857 dataOut.code = self.code
852 858 dataOut.nCode = self.nCode
853 859 dataOut.nBaud = self.nBaud
854 860
855 861 dataOut.data = datadec
856 862
857 863 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
858 864
859 865 dataOut.flagDecodeData = True #asumo q la data esta decodificada
860 866
861 867 if self.__profIndex == self.nCode-1:
862 868 self.__profIndex = 0
863 869 return dataOut
864 870
865 871 self.__profIndex += 1
866 872
867 873 return dataOut
868 874 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
869 875
870 876
871 877 class ProfileConcat(Operation):
872 878
873 879 isConfig = False
874 880 buffer = None
875 881
876 882 def __init__(self, **kwargs):
877 883
878 884 Operation.__init__(self, **kwargs)
879 885 self.profileIndex = 0
880 886
881 887 def reset(self):
882 888 self.buffer = numpy.zeros_like(self.buffer)
883 889 self.start_index = 0
884 890 self.times = 1
885 891
886 892 def setup(self, data, m, n=1):
887 893 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
888 894 self.nHeights = data.shape[1]#.nHeights
889 895 self.start_index = 0
890 896 self.times = 1
891 897
892 898 def concat(self, data):
893 899
894 900 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
895 901 self.start_index = self.start_index + self.nHeights
896 902
897 903 def run(self, dataOut, m):
898 904 dataOut.flagNoData = True
899 905
900 906 if not self.isConfig:
901 907 self.setup(dataOut.data, m, 1)
902 908 self.isConfig = True
903 909
904 910 if dataOut.flagDataAsBlock:
905 911 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
906 912
907 913 else:
908 914 self.concat(dataOut.data)
909 915 self.times += 1
910 916 if self.times > m:
911 917 dataOut.data = self.buffer
912 918 self.reset()
913 919 dataOut.flagNoData = False
914 920 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
915 921 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
916 922 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
917 923 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
918 924 dataOut.ippSeconds *= m
919 925 return dataOut
920 926
921 927 class ProfileSelector(Operation):
922 928
923 929 profileIndex = None
924 930 # Tamanho total de los perfiles
925 931 nProfiles = None
926 932
927 933 def __init__(self, **kwargs):
928 934
929 935 Operation.__init__(self, **kwargs)
930 936 self.profileIndex = 0
931 937
932 938 def incProfileIndex(self):
933 939
934 940 self.profileIndex += 1
935 941
936 942 if self.profileIndex >= self.nProfiles:
937 943 self.profileIndex = 0
938 944
939 945 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
940 946
941 947 if profileIndex < minIndex:
942 948 return False
943 949
944 950 if profileIndex > maxIndex:
945 951 return False
946 952
947 953 return True
948 954
949 955 def isThisProfileInList(self, profileIndex, profileList):
950 956
951 957 if profileIndex not in profileList:
952 958 return False
953 959
954 960 return True
955 961
956 962 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
957 963
958 964 """
959 965 ProfileSelector:
960 966
961 967 Inputs:
962 968 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
963 969
964 970 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
965 971
966 972 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
967 973
968 974 """
969 975
970 976 if rangeList is not None:
971 977 if type(rangeList[0]) not in (tuple, list):
972 978 rangeList = [rangeList]
973 979
974 980 dataOut.flagNoData = True
975 981
976 982 if dataOut.flagDataAsBlock:
977 983 """
978 984 data dimension = [nChannels, nProfiles, nHeis]
979 985 """
980 986 if profileList != None:
981 987 dataOut.data = dataOut.data[:,profileList,:]
982 988
983 989 if profileRangeList != None:
984 990 minIndex = profileRangeList[0]
985 991 maxIndex = profileRangeList[1]
986 992 profileList = list(range(minIndex, maxIndex+1))
987 993
988 994 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
989 995
990 996 if rangeList != None:
991 997
992 998 profileList = []
993 999
994 1000 for thisRange in rangeList:
995 1001 minIndex = thisRange[0]
996 1002 maxIndex = thisRange[1]
997 1003
998 1004 profileList.extend(list(range(minIndex, maxIndex+1)))
999 1005
1000 1006 dataOut.data = dataOut.data[:,profileList,:]
1001 1007
1002 1008 dataOut.nProfiles = len(profileList)
1003 1009 dataOut.profileIndex = dataOut.nProfiles - 1
1004 1010 dataOut.flagNoData = False
1005 1011
1006 1012 return dataOut
1007 1013
1008 1014 """
1009 1015 data dimension = [nChannels, nHeis]
1010 1016 """
1011 1017
1012 1018 if profileList != None:
1013 1019
1014 1020 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1015 1021
1016 1022 self.nProfiles = len(profileList)
1017 1023 dataOut.nProfiles = self.nProfiles
1018 1024 dataOut.profileIndex = self.profileIndex
1019 1025 dataOut.flagNoData = False
1020 1026
1021 1027 self.incProfileIndex()
1022 1028 return dataOut
1023 1029
1024 1030 if profileRangeList != None:
1025 1031
1026 1032 minIndex = profileRangeList[0]
1027 1033 maxIndex = profileRangeList[1]
1028 1034
1029 1035 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1030 1036
1031 1037 self.nProfiles = maxIndex - minIndex + 1
1032 1038 dataOut.nProfiles = self.nProfiles
1033 1039 dataOut.profileIndex = self.profileIndex
1034 1040 dataOut.flagNoData = False
1035 1041
1036 1042 self.incProfileIndex()
1037 1043 return dataOut
1038 1044
1039 1045 if rangeList != None:
1040 1046
1041 1047 nProfiles = 0
1042 1048
1043 1049 for thisRange in rangeList:
1044 1050 minIndex = thisRange[0]
1045 1051 maxIndex = thisRange[1]
1046 1052
1047 1053 nProfiles += maxIndex - minIndex + 1
1048 1054
1049 1055 for thisRange in rangeList:
1050 1056
1051 1057 minIndex = thisRange[0]
1052 1058 maxIndex = thisRange[1]
1053 1059
1054 1060 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1055 1061
1056 1062 self.nProfiles = nProfiles
1057 1063 dataOut.nProfiles = self.nProfiles
1058 1064 dataOut.profileIndex = self.profileIndex
1059 1065 dataOut.flagNoData = False
1060 1066
1061 1067 self.incProfileIndex()
1062 1068
1063 1069 break
1064 1070
1065 1071 return dataOut
1066 1072
1067 1073
1068 1074 if beam != None: #beam is only for AMISR data
1069 1075 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1070 1076 dataOut.flagNoData = False
1071 1077 dataOut.profileIndex = self.profileIndex
1072 1078
1073 1079 self.incProfileIndex()
1074 1080
1075 1081 return dataOut
1076 1082
1077 1083 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1078 1084
1079 1085
1080 1086 class Reshaper(Operation):
1081 1087
1082 1088 def __init__(self, **kwargs):
1083 1089
1084 1090 Operation.__init__(self, **kwargs)
1085 1091
1086 1092 self.__buffer = None
1087 1093 self.__nitems = 0
1088 1094
1089 1095 def __appendProfile(self, dataOut, nTxs):
1090 1096
1091 1097 if self.__buffer is None:
1092 1098 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1093 1099 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1094 1100
1095 1101 ini = dataOut.nHeights * self.__nitems
1096 1102 end = ini + dataOut.nHeights
1097 1103
1098 1104 self.__buffer[:, ini:end] = dataOut.data
1099 1105
1100 1106 self.__nitems += 1
1101 1107
1102 1108 return int(self.__nitems*nTxs)
1103 1109
1104 1110 def __getBuffer(self):
1105 1111
1106 1112 if self.__nitems == int(1./self.__nTxs):
1107 1113
1108 1114 self.__nitems = 0
1109 1115
1110 1116 return self.__buffer.copy()
1111 1117
1112 1118 return None
1113 1119
1114 1120 def __checkInputs(self, dataOut, shape, nTxs):
1115 1121
1116 1122 if shape is None and nTxs is None:
1117 1123 raise ValueError("Reshaper: shape of factor should be defined")
1118 1124
1119 1125 if nTxs:
1120 1126 if nTxs < 0:
1121 1127 raise ValueError("nTxs should be greater than 0")
1122 1128
1123 1129 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1124 1130 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1125 1131
1126 1132 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1127 1133
1128 1134 return shape, nTxs
1129 1135
1130 1136 if len(shape) != 2 and len(shape) != 3:
1131 1137 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))
1132 1138
1133 1139 if len(shape) == 2:
1134 1140 shape_tuple = [dataOut.nChannels]
1135 1141 shape_tuple.extend(shape)
1136 1142 else:
1137 1143 shape_tuple = list(shape)
1138 1144
1139 1145 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1140 1146
1141 1147 return shape_tuple, nTxs
1142 1148
1143 1149 def run(self, dataOut, shape=None, nTxs=None):
1144 1150
1145 1151 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1146 1152
1147 1153 dataOut.flagNoData = True
1148 1154 profileIndex = None
1149 1155
1150 1156 if dataOut.flagDataAsBlock:
1151 1157
1152 1158 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1153 1159 dataOut.flagNoData = False
1154 1160
1155 1161 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1156 1162
1157 1163 else:
1158 1164
1159 1165 if self.__nTxs < 1:
1160 1166
1161 1167 self.__appendProfile(dataOut, self.__nTxs)
1162 1168 new_data = self.__getBuffer()
1163 1169
1164 1170 if new_data is not None:
1165 1171 dataOut.data = new_data
1166 1172 dataOut.flagNoData = False
1167 1173
1168 1174 profileIndex = dataOut.profileIndex*nTxs
1169 1175
1170 1176 else:
1171 1177 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1172 1178
1173 1179 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1174 1180
1175 1181 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1176 1182
1177 1183 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1178 1184
1179 1185 dataOut.profileIndex = profileIndex
1180 1186
1181 1187 dataOut.ippSeconds /= self.__nTxs
1182 1188
1183 1189 return dataOut
1184 1190
1185 1191 class SplitProfiles(Operation):
1186 1192
1187 1193 def __init__(self, **kwargs):
1188 1194
1189 1195 Operation.__init__(self, **kwargs)
1190 1196
1191 1197 def run(self, dataOut, n):
1192 1198
1193 1199 dataOut.flagNoData = True
1194 1200 profileIndex = None
1195 1201
1196 1202 if dataOut.flagDataAsBlock:
1197 1203
1198 1204 #nchannels, nprofiles, nsamples
1199 1205 shape = dataOut.data.shape
1200 1206
1201 1207 if shape[2] % n != 0:
1202 1208 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1203 1209
1204 1210 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1205 1211
1206 1212 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1207 1213 dataOut.flagNoData = False
1208 1214
1209 1215 profileIndex = int(dataOut.nProfiles/n) - 1
1210 1216
1211 1217 else:
1212 1218
1213 1219 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1214 1220
1215 1221 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1216 1222
1217 1223 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1218 1224
1219 1225 dataOut.nProfiles = int(dataOut.nProfiles*n)
1220 1226
1221 1227 dataOut.profileIndex = profileIndex
1222 1228
1223 1229 dataOut.ippSeconds /= n
1224 1230
1225 1231 return dataOut
1226 1232
1227 1233 class CombineProfiles(Operation):
1228 1234 def __init__(self, **kwargs):
1229 1235
1230 1236 Operation.__init__(self, **kwargs)
1231 1237
1232 1238 self.__remData = None
1233 1239 self.__profileIndex = 0
1234 1240
1235 1241 def run(self, dataOut, n):
1236 1242
1237 1243 dataOut.flagNoData = True
1238 1244 profileIndex = None
1239 1245
1240 1246 if dataOut.flagDataAsBlock:
1241 1247
1242 1248 #nchannels, nprofiles, nsamples
1243 1249 shape = dataOut.data.shape
1244 1250 new_shape = shape[0], shape[1]/n, shape[2]*n
1245 1251
1246 1252 if shape[1] % n != 0:
1247 1253 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1248 1254
1249 1255 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1250 1256 dataOut.flagNoData = False
1251 1257
1252 1258 profileIndex = int(dataOut.nProfiles*n) - 1
1253 1259
1254 1260 else:
1255 1261
1256 1262 #nchannels, nsamples
1257 1263 if self.__remData is None:
1258 1264 newData = dataOut.data
1259 1265 else:
1260 1266 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1261 1267
1262 1268 self.__profileIndex += 1
1263 1269
1264 1270 if self.__profileIndex < n:
1265 1271 self.__remData = newData
1266 1272 #continue
1267 1273 return
1268 1274
1269 1275 self.__profileIndex = 0
1270 1276 self.__remData = None
1271 1277
1272 1278 dataOut.data = newData
1273 1279 dataOut.flagNoData = False
1274 1280
1275 1281 profileIndex = dataOut.profileIndex/n
1276 1282
1277 1283
1278 1284 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1279 1285
1280 1286 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1281 1287
1282 1288 dataOut.nProfiles = int(dataOut.nProfiles/n)
1283 1289
1284 1290 dataOut.profileIndex = profileIndex
1285 1291
1286 1292 dataOut.ippSeconds *= n
1287 1293
1288 1294 return dataOut
1289 1295
1290 1296 class PulsePairVoltage(Operation):
1291 1297 '''
1292 1298 Function PulsePair(Signal Power, Velocity)
1293 1299 The real component of Lag[0] provides Intensity Information
1294 1300 The imag component of Lag[1] Phase provides Velocity Information
1295 1301
1296 1302 Configuration Parameters:
1297 1303 nPRF = Number of Several PRF
1298 1304 theta = Degree Azimuth angel Boundaries
1299 1305
1300 1306 Input:
1301 1307 self.dataOut
1302 1308 lag[N]
1303 1309 Affected:
1304 1310 self.dataOut.spc
1305 1311 '''
1306 1312 isConfig = False
1307 1313 __profIndex = 0
1308 1314 __initime = None
1309 1315 __lastdatatime = None
1310 1316 __buffer = None
1311 1317 noise = None
1312 1318 __dataReady = False
1313 1319 n = None
1314 1320 __nch = 0
1315 1321 __nHeis = 0
1316 1322 removeDC = False
1317 1323 ipp = None
1318 1324 lambda_ = 0
1319 1325
1320 1326 def __init__(self,**kwargs):
1321 1327 Operation.__init__(self,**kwargs)
1322 1328
1323 1329 def setup(self, dataOut, n = None, removeDC=False):
1324 1330 '''
1325 1331 n= Numero de PRF's de entrada
1326 1332 '''
1327 1333 self.__initime = None
1328 1334 self.__lastdatatime = 0
1329 1335 self.__dataReady = False
1330 1336 self.__buffer = 0
1331 1337 self.__profIndex = 0
1332 1338 self.noise = None
1333 1339 self.__nch = dataOut.nChannels
1334 1340 self.__nHeis = dataOut.nHeights
1335 1341 self.removeDC = removeDC
1336 1342 self.lambda_ = 3.0e8/(9345.0e6)
1337 1343 self.ippSec = dataOut.ippSeconds
1338 1344 self.nCohInt = dataOut.nCohInt
1339 1345 print("IPPseconds",dataOut.ippSeconds)
1340 1346
1341 1347 print("ELVALOR DE n es:", n)
1342 1348 if n == None:
1343 1349 raise ValueError("n should be specified.")
1344 1350
1345 1351 if n != None:
1346 1352 if n<2:
1347 1353 raise ValueError("n should be greater than 2")
1348 1354
1349 1355 self.n = n
1350 1356 self.__nProf = n
1351 1357
1352 1358 self.__buffer = numpy.zeros((dataOut.nChannels,
1353 1359 n,
1354 1360 dataOut.nHeights),
1355 1361 dtype='complex')
1356 1362
1357 1363 def putData(self,data):
1358 1364 '''
1359 1365 Add a profile to he __buffer and increase in one the __profiel Index
1360 1366 '''
1361 1367 self.__buffer[:,self.__profIndex,:]= data
1362 1368 self.__profIndex += 1
1363 1369 return
1364 1370
1365 1371 def pushData(self,dataOut):
1366 1372 '''
1367 1373 Return the PULSEPAIR and the profiles used in the operation
1368 1374 Affected : self.__profileIndex
1369 1375 '''
1370 1376 #----------------- Remove DC-----------------------------------
1371 1377 if self.removeDC==True:
1372 1378 mean = numpy.mean(self.__buffer,1)
1373 1379 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1374 1380 dc= numpy.tile(tmp,[1,self.__nProf,1])
1375 1381 self.__buffer = self.__buffer - dc
1376 1382 #------------------Calculo de Potencia ------------------------
1377 1383 pair0 = self.__buffer*numpy.conj(self.__buffer)
1378 1384 pair0 = pair0.real
1379 1385 lag_0 = numpy.sum(pair0,1)
1380 1386 #------------------Calculo de Ruido x canal--------------------
1381 1387 self.noise = numpy.zeros(self.__nch)
1382 1388 for i in range(self.__nch):
1383 1389 daux = numpy.sort(pair0[i,:,:],axis= None)
1384 1390 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1385 1391
1386 1392 self.noise = self.noise.reshape(self.__nch,1)
1387 1393 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1388 1394 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1389 1395 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1390 1396 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1391 1397 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1392 1398 #-------------------- Power --------------------------------------------------
1393 1399 data_power = lag_0/(self.n*self.nCohInt)
1394 1400 #------------------ Senal ---------------------------------------------------
1395 1401 data_intensity = pair0 - noise_buffer
1396 1402 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1397 1403 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1398 1404 for i in range(self.__nch):
1399 1405 for j in range(self.__nHeis):
1400 1406 if data_intensity[i][j] < 0:
1401 1407 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1402 1408
1403 1409 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1404 1410 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1405 1411 lag_1 = numpy.sum(pair1,1)
1406 1412 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1407 1413 data_velocity = (self.lambda_/2.0)*data_freq
1408 1414
1409 1415 #---------------- Potencia promedio estimada de la Senal-----------
1410 1416 lag_0 = lag_0/self.n
1411 1417 S = lag_0-self.noise
1412 1418
1413 1419 #---------------- Frecuencia Doppler promedio ---------------------
1414 1420 lag_1 = lag_1/(self.n-1)
1415 1421 R1 = numpy.abs(lag_1)
1416 1422
1417 1423 #---------------- Calculo del SNR----------------------------------
1418 1424 data_snrPP = S/self.noise
1419 1425 for i in range(self.__nch):
1420 1426 for j in range(self.__nHeis):
1421 1427 if data_snrPP[i][j] < 1.e-20:
1422 1428 data_snrPP[i][j] = 1.e-20
1423 1429
1424 1430 #----------------- Calculo del ancho espectral ----------------------
1425 1431 L = S/R1
1426 1432 L = numpy.where(L<0,1,L)
1427 1433 L = numpy.log(L)
1428 1434 tmp = numpy.sqrt(numpy.absolute(L))
1429 1435 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1430 1436 n = self.__profIndex
1431 1437
1432 1438 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1433 1439 self.__profIndex = 0
1434 1440 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1435 1441
1436 1442
1437 1443 def pulsePairbyProfiles(self,dataOut):
1438 1444
1439 1445 self.__dataReady = False
1440 1446 data_power = None
1441 1447 data_intensity = None
1442 1448 data_velocity = None
1443 1449 data_specwidth = None
1444 1450 data_snrPP = None
1445 1451 self.putData(data=dataOut.data)
1446 1452 if self.__profIndex == self.n:
1447 1453 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1448 1454 self.__dataReady = True
1449 1455
1450 1456 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1451 1457
1452 1458
1453 1459 def pulsePairOp(self, dataOut, datatime= None):
1454 1460
1455 1461 if self.__initime == None:
1456 1462 self.__initime = datatime
1457 1463 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1458 1464 self.__lastdatatime = datatime
1459 1465
1460 1466 if data_power is None:
1461 1467 return None, None, None,None,None,None
1462 1468
1463 1469 avgdatatime = self.__initime
1464 1470 deltatime = datatime - self.__lastdatatime
1465 1471 self.__initime = datatime
1466 1472
1467 1473 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1468 1474
1469 1475 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1470 1476
1471 1477 if not self.isConfig:
1472 1478 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1473 1479 self.isConfig = True
1474 1480 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1475 1481 dataOut.flagNoData = True
1476 1482
1477 1483 if self.__dataReady:
1478 1484 dataOut.nCohInt *= self.n
1479 1485 dataOut.dataPP_POW = data_intensity # S
1480 1486 dataOut.dataPP_POWER = data_power # P
1481 1487 dataOut.dataPP_DOP = data_velocity
1482 1488 dataOut.dataPP_SNR = data_snrPP
1483 1489 dataOut.dataPP_WIDTH = data_specwidth
1484 1490 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1485 1491 dataOut.utctime = avgdatatime
1486 1492 dataOut.flagNoData = False
1487 1493 return dataOut
1488 1494
1489 1495
1490 1496
1491 1497 # import collections
1492 1498 # from scipy.stats import mode
1493 1499 #
1494 1500 # class Synchronize(Operation):
1495 1501 #
1496 1502 # isConfig = False
1497 1503 # __profIndex = 0
1498 1504 #
1499 1505 # def __init__(self, **kwargs):
1500 1506 #
1501 1507 # Operation.__init__(self, **kwargs)
1502 1508 # # self.isConfig = False
1503 1509 # self.__powBuffer = None
1504 1510 # self.__startIndex = 0
1505 1511 # self.__pulseFound = False
1506 1512 #
1507 1513 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1508 1514 #
1509 1515 # #Read data
1510 1516 #
1511 1517 # powerdB = dataOut.getPower(channel = channel)
1512 1518 # noisedB = dataOut.getNoise(channel = channel)[0]
1513 1519 #
1514 1520 # self.__powBuffer.extend(powerdB.flatten())
1515 1521 #
1516 1522 # dataArray = numpy.array(self.__powBuffer)
1517 1523 #
1518 1524 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1519 1525 #
1520 1526 # maxValue = numpy.nanmax(filteredPower)
1521 1527 #
1522 1528 # if maxValue < noisedB + 10:
1523 1529 # #No se encuentra ningun pulso de transmision
1524 1530 # return None
1525 1531 #
1526 1532 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1527 1533 #
1528 1534 # if len(maxValuesIndex) < 2:
1529 1535 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1530 1536 # return None
1531 1537 #
1532 1538 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1533 1539 #
1534 1540 # #Seleccionar solo valores con un espaciamiento de nSamples
1535 1541 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1536 1542 #
1537 1543 # if len(pulseIndex) < 2:
1538 1544 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1539 1545 # return None
1540 1546 #
1541 1547 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1542 1548 #
1543 1549 # #remover senales que se distancien menos de 10 unidades o muestras
1544 1550 # #(No deberian existir IPP menor a 10 unidades)
1545 1551 #
1546 1552 # realIndex = numpy.where(spacing > 10 )[0]
1547 1553 #
1548 1554 # if len(realIndex) < 2:
1549 1555 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1550 1556 # return None
1551 1557 #
1552 1558 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1553 1559 # realPulseIndex = pulseIndex[realIndex]
1554 1560 #
1555 1561 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1556 1562 #
1557 1563 # print "IPP = %d samples" %period
1558 1564 #
1559 1565 # self.__newNSamples = dataOut.nHeights #int(period)
1560 1566 # self.__startIndex = int(realPulseIndex[0])
1561 1567 #
1562 1568 # return 1
1563 1569 #
1564 1570 #
1565 1571 # def setup(self, nSamples, nChannels, buffer_size = 4):
1566 1572 #
1567 1573 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1568 1574 # maxlen = buffer_size*nSamples)
1569 1575 #
1570 1576 # bufferList = []
1571 1577 #
1572 1578 # for i in range(nChannels):
1573 1579 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1574 1580 # maxlen = buffer_size*nSamples)
1575 1581 #
1576 1582 # bufferList.append(bufferByChannel)
1577 1583 #
1578 1584 # self.__nSamples = nSamples
1579 1585 # self.__nChannels = nChannels
1580 1586 # self.__bufferList = bufferList
1581 1587 #
1582 1588 # def run(self, dataOut, channel = 0):
1583 1589 #
1584 1590 # if not self.isConfig:
1585 1591 # nSamples = dataOut.nHeights
1586 1592 # nChannels = dataOut.nChannels
1587 1593 # self.setup(nSamples, nChannels)
1588 1594 # self.isConfig = True
1589 1595 #
1590 1596 # #Append new data to internal buffer
1591 1597 # for thisChannel in range(self.__nChannels):
1592 1598 # bufferByChannel = self.__bufferList[thisChannel]
1593 1599 # bufferByChannel.extend(dataOut.data[thisChannel])
1594 1600 #
1595 1601 # if self.__pulseFound:
1596 1602 # self.__startIndex -= self.__nSamples
1597 1603 #
1598 1604 # #Finding Tx Pulse
1599 1605 # if not self.__pulseFound:
1600 1606 # indexFound = self.__findTxPulse(dataOut, channel)
1601 1607 #
1602 1608 # if indexFound == None:
1603 1609 # dataOut.flagNoData = True
1604 1610 # return
1605 1611 #
1606 1612 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1607 1613 # self.__pulseFound = True
1608 1614 # self.__startIndex = indexFound
1609 1615 #
1610 1616 # #If pulse was found ...
1611 1617 # for thisChannel in range(self.__nChannels):
1612 1618 # bufferByChannel = self.__bufferList[thisChannel]
1613 1619 # #print self.__startIndex
1614 1620 # x = numpy.array(bufferByChannel)
1615 1621 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1616 1622 #
1617 1623 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1618 1624 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1619 1625 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1620 1626 #
1621 1627 # dataOut.data = self.__arrayBuffer
1622 1628 #
1623 1629 # self.__startIndex += self.__newNSamples
1624 1630 #
1625 1631 # return
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
General Comments 0
You need to be logged in to leave comments. Login now