##// END OF EJS Templates
schain
RSS
Valley Densities Proc
rflores -
r1707:21f41e756eb1
parent child
Show More
Side by Side Unified
Context file:
r1707:21f41e756eb1 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,1345 +1,1428
1 1
2 2 import os
3 3 import time
4 4 import math
5 5 import datetime
6 6 import numpy
7 7
8 8 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator #YONG
9 9
10 10 from .jroplot_spectra import RTIPlot, NoisePlot
11 11
12 12 from schainpy.utils import log
13 13 from .plotting_codes import *
14 14
15 15 from schainpy.model.graphics.jroplot_base import Plot, plt
16 16
17 17 import matplotlib.pyplot as plt
18 18 import matplotlib.colors as colors
19 19 from matplotlib.ticker import MultipleLocator, LogLocator, NullFormatter
20 20
21 21 class RTIDPPlot(RTIPlot):
22 22 '''
23 23 Written by R. Flores
24 24 '''
25 25 '''Plot for RTI Double Pulse Experiment Using Cross Products Analysis
26 26 '''
27 27
28 28 CODE = 'RTIDP'
29 29 colormap = 'jet'
30 30 plot_name = 'RTI'
31 31 plot_type = 'pcolorbuffer'
32 32
33 33 def setup(self):
34 34 self.xaxis = 'time'
35 35 self.ncols = 1
36 36 self.nrows = 3
37 37 self.nplots = self.nrows
38 38
39 39 self.ylabel = 'Range [km]'
40 40 self.xlabel = 'Time (LT)'
41 41
42 42 self.cb_label = 'Intensity (dB)'
43 43
44 44 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
45 45
46 46 self.titles = ['{} Channel {}'.format(
47 47 self.plot_name.upper(), '0x1'),'{} Channel {}'.format(
48 48 self.plot_name.upper(), '0'),'{} Channel {}'.format(
49 49 self.plot_name.upper(), '1')]
50 50
51 51 def update(self, dataOut):
52 52
53 53 data = {}
54 54 meta = {}
55 55 data['rti'] = dataOut.data_for_RTI_DP
56 56 data['NDP'] = dataOut.NDP
57 57
58 58 return data, meta
59 59
60 60 def plot(self):
61 61
62 62 NDP = self.data['NDP'][-1]
63 63 self.x = self.data.times
64 64 self.y = self.data.yrange[0:NDP]
65 65 self.z = self.data['rti']
66 66 self.z = numpy.ma.masked_invalid(self.z)
67 67
68 68 if self.decimation is None:
69 69 x, y, z = self.fill_gaps(self.x, self.y, self.z)
70 70 else:
71 71 x, y, z = self.fill_gaps(*self.decimate())
72 72
73 73 for n, ax in enumerate(self.axes):
74 74
75 75 self.zmax = self.zmax if self.zmax is not None else numpy.max(
76 76 self.z[1][0,12:40])
77 77 self.zmin = self.zmin if self.zmin is not None else numpy.min(
78 78 self.z[1][0,12:40])
79 79
80 80 if ax.firsttime:
81 81
82 82 if self.zlimits is not None:
83 83 self.zmin, self.zmax = self.zlimits[n]
84 84
85 85 ax.plt = ax.pcolormesh(x, y, z[n].T,
86 86 vmin=self.zmin,
87 87 vmax=self.zmax,
88 88 cmap=plt.get_cmap(self.colormap)
89 89 )
90 90 else:
91 91 #if self.zlimits is not None:
92 92 #self.zmin, self.zmax = self.zlimits[n]
93 93 ax.collections.remove(ax.collections[0])
94 94 ax.plt = ax.pcolormesh(x, y, z[n].T,
95 95 vmin=self.zmin,
96 96 vmax=self.zmax,
97 97 cmap=plt.get_cmap(self.colormap)
98 98 )
99 99
100 100
101 101 class RTILPPlot(RTIPlot):
102 102 '''
103 103 Written by R. Flores
104 104 '''
105 105 '''
106 106 Plot for RTI Long Pulse Using Cross Products Analysis
107 107 '''
108 108
109 109 CODE = 'RTILP'
110 110 colormap = 'jet'
111 111 plot_name = 'RTI LP'
112 112 plot_type = 'pcolorbuffer'
113 113
114 114 def setup(self):
115 115 self.xaxis = 'time'
116 116 self.ncols = 1
117 117 self.nrows = 2
118 118 self.nplots = self.nrows
119 119
120 120 self.ylabel = 'Range [km]'
121 121 self.xlabel = 'Time (LT)'
122 122
123 123 self.cb_label = 'Intensity (dB)'
124 124
125 125 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
126 126
127 127
128 128 self.titles = ['{} Channel {}'.format(
129 129 self.plot_name.upper(), '0'),'{} Channel {}'.format(
130 130 self.plot_name.upper(), '1'),'{} Channel {}'.format(
131 131 self.plot_name.upper(), '2'),'{} Channel {}'.format(
132 132 self.plot_name.upper(), '3')]
133 133
134 134
135 135 def update(self, dataOut):
136 136
137 137 data = {}
138 138 meta = {}
139 139 data['rti'] = dataOut.data_for_RTI_LP
140 140 data['NRANGE'] = dataOut.NRANGE
141 141
142 142 return data, meta
143 143
144 144 def plot(self):
145 145
146 146 NRANGE = self.data['NRANGE'][-1]
147 147 self.x = self.data.times
148 148 self.y = self.data.yrange[0:NRANGE]
149 149
150 150 self.z = self.data['rti']
151 151
152 152 self.z = numpy.ma.masked_invalid(self.z)
153 153
154 154 if self.decimation is None:
155 155 x, y, z = self.fill_gaps(self.x, self.y, self.z)
156 156 else:
157 157 x, y, z = self.fill_gaps(*self.decimate())
158 158
159 159 for n, ax in enumerate(self.axes):
160 160
161 161 self.zmax = self.zmax if self.zmax is not None else numpy.max(
162 162 self.z[1][0,12:40])
163 163 self.zmin = self.zmin if self.zmin is not None else numpy.min(
164 164 self.z[1][0,12:40])
165 165
166 166 if ax.firsttime:
167 167
168 168 if self.zlimits is not None:
169 169 self.zmin, self.zmax = self.zlimits[n]
170 170
171 171
172 172 ax.plt = ax.pcolormesh(x, y, z[n].T,
173 173 vmin=self.zmin,
174 174 vmax=self.zmax,
175 175 cmap=plt.get_cmap(self.colormap)
176 176 )
177 177
178 178 else:
179 179 if self.zlimits is not None:
180 180 self.zmin, self.zmax = self.zlimits[n]
181 181 ax.collections.remove(ax.collections[0])
182 182 ax.plt = ax.pcolormesh(x, y, z[n].T,
183 183 vmin=self.zmin,
184 184 vmax=self.zmax,
185 185 cmap=plt.get_cmap(self.colormap)
186 186 )
187 187
188 188
189 189 class DenRTIPlot(RTIPlot):
190 190 '''
191 191 Written by R. Flores
192 192 '''
193 193 '''
194 194 Plot for Den
195 195 '''
196 196
197 197 CODE = 'denrti'
198 198 colormap = 'jet'
199 199
200 200 def setup(self):
201 201 self.xaxis = 'time'
202 202 self.ncols = 1
203 203 self.nrows = self.data.shape(self.CODE)[0]
204 204 self.nplots = self.nrows
205 205
206 206 self.ylabel = 'Range [km]'
207 207 self.xlabel = 'Time (LT)'
208 208
209 209 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
210 210
211 211 if self.CODE == 'denrti':
212 212 self.cb_label = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
213 213
214 214 self.titles = ['Electron Density RTI']
215 215
216 216 def update(self, dataOut):
217 217
218 218 data = {}
219 219 meta = {}
220 220
221 221 data['denrti'] = dataOut.DensityFinal*1.e-6 #To Plot in cm^-3
222 222
223 223 return data, meta
224 224
225 225 def plot(self):
226 226
227 227 self.x = self.data.times
228 228 self.y = self.data.yrange
229 229
230 230 self.z = self.data[self.CODE]
231 231
232 232 self.z = numpy.ma.masked_invalid(self.z)
233 233
234 234 if self.decimation is None:
235 235 x, y, z = self.fill_gaps(self.x, self.y, self.z)
236 236 else:
237 237 x, y, z = self.fill_gaps(*self.decimate())
238 238
239 239 for n, ax in enumerate(self.axes):
240 240
241 241 self.zmax = self.zmax if self.zmax is not None else numpy.max(
242 242 self.z[n])
243 243 self.zmin = self.zmin if self.zmin is not None else numpy.min(
244 244 self.z[n])
245 245
246 246 if ax.firsttime:
247 247
248 248 if self.zlimits is not None:
249 249 self.zmin, self.zmax = self.zlimits[n]
250 250 if numpy.log10(self.zmin)<0:
251 251 self.zmin=1
252 252 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
253 253 vmin=self.zmin,
254 254 vmax=self.zmax,
255 255 cmap=self.cmaps[n],
256 256 norm=colors.LogNorm()
257 257 )
258 258
259 259 else:
260 260 if self.zlimits is not None:
261 261 self.zmin, self.zmax = self.zlimits[n]
262 262 ax.collections.remove(ax.collections[0])
263 263 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
264 264 vmin=self.zmin,
265 265 vmax=self.zmax,
266 266 cmap=self.cmaps[n],
267 267 norm=colors.LogNorm()
268 268 )
269 269
270 270
271 271 class ETempRTIPlot(RTIPlot):
272 272 '''
273 273 Written by R. Flores
274 274 '''
275 275 '''
276 276 Plot for Electron Temperature
277 277 '''
278 278
279 279 CODE = 'ETemp'
280 280 colormap = 'jet'
281 281
282 282 def setup(self):
283 283 self.xaxis = 'time'
284 284 self.ncols = 1
285 285 self.nrows = self.data.shape(self.CODE)[0]
286 286 self.nplots = self.nrows
287 287
288 288 self.ylabel = 'Range [km]'
289 289 self.xlabel = 'Time (LT)'
290 290 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
291 291 if self.CODE == 'ETemp':
292 292 self.cb_label = 'Electron Temperature (K)'
293 293 self.titles = ['Electron Temperature RTI']
294 294 if self.CODE == 'ITemp':
295 295 self.cb_label = 'Ion Temperature (K)'
296 296 self.titles = ['Ion Temperature RTI']
297 297 if self.CODE == 'HeFracLP':
298 298 self.cb_label ='He+ Fraction'
299 299 self.titles = ['He+ Fraction RTI']
300 300 self.zmax=0.16
301 301 if self.CODE == 'HFracLP':
302 302 self.cb_label ='H+ Fraction'
303 303 self.titles = ['H+ Fraction RTI']
304 304
305 305 def update(self, dataOut):
306 306
307 307 data = {}
308 308 meta = {}
309 309
310 310 data['ETemp'] = dataOut.ElecTempFinal
311 311
312 312 return data, meta
313 313
314 314 def plot(self):
315 315
316 316 self.x = self.data.times
317 317 self.y = self.data.yrange
318 318 self.z = self.data[self.CODE]
319 319
320 320 self.z = numpy.ma.masked_invalid(self.z)
321 321
322 322 if self.decimation is None:
323 323 x, y, z = self.fill_gaps(self.x, self.y, self.z)
324 324 else:
325 325 x, y, z = self.fill_gaps(*self.decimate())
326 326
327 327 for n, ax in enumerate(self.axes):
328 328
329 329 self.zmax = self.zmax if self.zmax is not None else numpy.max(
330 330 self.z[n])
331 331 self.zmin = self.zmin if self.zmin is not None else numpy.min(
332 332 self.z[n])
333 333
334 334 if ax.firsttime:
335 335
336 336 if self.zlimits is not None:
337 337 self.zmin, self.zmax = self.zlimits[n]
338 338
339 339 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
340 340 vmin=self.zmin,
341 341 vmax=self.zmax,
342 342 cmap=self.cmaps[n]
343 343 )
344 344 #plt.tight_layout()
345 345
346 346 else:
347 347 if self.zlimits is not None:
348 348 self.zmin, self.zmax = self.zlimits[n]
349 349 ax.collections.remove(ax.collections[0])
350 350 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
351 351 vmin=self.zmin,
352 352 vmax=self.zmax,
353 353 cmap=self.cmaps[n]
354 354 )
355 355
356 356
357 357 class ITempRTIPlot(ETempRTIPlot):
358 358 '''
359 359 Written by R. Flores
360 360 '''
361 361 '''
362 362 Plot for Ion Temperature
363 363 '''
364 364
365 365 CODE = 'ITemp'
366 366 colormap = 'jet'
367 367 plot_name = 'Ion Temperature'
368 368
369 369 def update(self, dataOut):
370 370
371 371 data = {}
372 372 meta = {}
373 373
374 374 data['ITemp'] = dataOut.IonTempFinal
375 375
376 376 return data, meta
377 377
378 378
379 379 class HFracRTIPlot(ETempRTIPlot):
380 380 '''
381 381 Written by R. Flores
382 382 '''
383 383 '''
384 384 Plot for H+ LP
385 385 '''
386 386
387 387 CODE = 'HFracLP'
388 388 colormap = 'jet'
389 389 plot_name = 'H+ Frac'
390 390
391 391 def update(self, dataOut):
392 392
393 393 data = {}
394 394 meta = {}
395 395 data['HFracLP'] = dataOut.PhyFinal
396 396
397 397 return data, meta
398 398
399 399
400 400 class HeFracRTIPlot(ETempRTIPlot):
401 401 '''
402 402 Written by R. Flores
403 403 '''
404 404 '''
405 405 Plot for He+ LP
406 406 '''
407 407
408 408 CODE = 'HeFracLP'
409 409 colormap = 'jet'
410 410 plot_name = 'He+ Frac'
411 411
412 412 def update(self, dataOut):
413 413
414 414 data = {}
415 415 meta = {}
416 416 data['HeFracLP'] = dataOut.PheFinal
417 417
418 418 return data, meta
419 419
420 420
421 421 class TempsDPPlot(Plot):
422 422 '''
423 423 Written by R. Flores
424 424 '''
425 425 '''
426 426 Plot for Electron - Ion Temperatures
427 427 '''
428 428
429 429 CODE = 'tempsDP'
430 430 #plot_name = 'Temperatures'
431 431 plot_type = 'scatterbuffer'
432 432
433 433 def setup(self):
434 434
435 435 self.ncols = 1
436 436 self.nrows = 1
437 437 self.nplots = 1
438 438 self.ylabel = 'Range [km]'
439 439 self.xlabel = 'Temperature (K)'
440 440 self.titles = ['Electron/Ion Temperatures']
441 441 self.width = 3.5
442 442 self.height = 5.5
443 443 self.colorbar = False
444 444 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
445 445
446 446 def update(self, dataOut):
447 447 data = {}
448 448 meta = {}
449 449
450 450 data['Te'] = dataOut.te2
451 451 data['Ti'] = dataOut.ti2
452 452 data['Te_error'] = dataOut.ete2
453 453 data['Ti_error'] = dataOut.eti2
454 454
455 455 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
456 456
457 457 return data, meta
458 458
459 459 def plot(self):
460 460
461 461 y = self.data.yrange
462 462
463 463 self.xmin = -100
464 464 self.xmax = 5000
465 465
466 466 ax = self.axes[0]
467 467
468 468 data = self.data[-1]
469 469
470 470 Te = data['Te']
471 471 Ti = data['Ti']
472 472 errTe = data['Te_error']
473 473 errTi = data['Ti_error']
474 474
475 475 if ax.firsttime:
476 476 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
477 477 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
478 478 plt.legend(loc='lower right')
479 479 self.ystep_given = 50
480 480 ax.yaxis.set_minor_locator(MultipleLocator(15))
481 481 ax.grid(which='minor')
482 482
483 483 else:
484 484 self.clear_figures()
485 485 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
486 486 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
487 487 plt.legend(loc='lower right')
488 488 ax.yaxis.set_minor_locator(MultipleLocator(15))
489 489
490 490
491 491 class TempsHPPlot(Plot):
492 492 '''
493 493 Written by R. Flores
494 494 '''
495 495 '''
496 496 Plot for Temperatures Hybrid Experiment
497 497 '''
498 498
499 499 CODE = 'temps_LP'
500 500 #plot_name = 'Temperatures'
501 501 plot_type = 'scatterbuffer'
502 502
503 503
504 504 def setup(self):
505 505
506 506 self.ncols = 1
507 507 self.nrows = 1
508 508 self.nplots = 1
509 509 self.ylabel = 'Range [km]'
510 510 self.xlabel = 'Temperature (K)'
511 511 self.titles = ['Electron/Ion Temperatures']
512 512 self.width = 3.5
513 513 self.height = 6.5
514 514 self.colorbar = False
515 515 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
516 516
517 517 def update(self, dataOut):
518 518 data = {}
519 519 meta = {}
520 520
521 521
522 522 data['Te'] = numpy.concatenate((dataOut.te2[:dataOut.cut],dataOut.te[dataOut.cut:]))
523 523 data['Ti'] = numpy.concatenate((dataOut.ti2[:dataOut.cut],dataOut.ti[dataOut.cut:]))
524 524 data['Te_error'] = numpy.concatenate((dataOut.ete2[:dataOut.cut],dataOut.ete[dataOut.cut:]))
525 525 data['Ti_error'] = numpy.concatenate((dataOut.eti2[:dataOut.cut],dataOut.eti[dataOut.cut:]))
526 526
527 527 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
528 528
529 529 return data, meta
530 530
531 531 def plot(self):
532 532
533 533
534 534 self.y = self.data.yrange
535 535 self.xmin = -100
536 536 self.xmax = 4500
537 537 ax = self.axes[0]
538 538
539 539 data = self.data[-1]
540 540
541 541 Te = data['Te']
542 542 Ti = data['Ti']
543 543 errTe = data['Te_error']
544 544 errTi = data['Ti_error']
545 545
546 546 if ax.firsttime:
547 547
548 548 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
549 549 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
550 550 plt.legend(loc='lower right')
551 551 self.ystep_given = 200
552 552 ax.yaxis.set_minor_locator(MultipleLocator(15))
553 553 ax.grid(which='minor')
554 554
555 555 else:
556 556 self.clear_figures()
557 557 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
558 558 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
559 559 plt.legend(loc='lower right')
560 560 ax.yaxis.set_minor_locator(MultipleLocator(15))
561 561 ax.grid(which='minor')
562 562
563 563
564 564 class FracsHPPlot(Plot):
565 565 '''
566 566 Written by R. Flores
567 567 '''
568 568 '''
569 569 Plot for Composition LP
570 570 '''
571 571
572 572 CODE = 'fracs_LP'
573 573 plot_type = 'scatterbuffer'
574 574
575 575
576 576 def setup(self):
577 577
578 578 self.ncols = 1
579 579 self.nrows = 1
580 580 self.nplots = 1
581 581 self.ylabel = 'Range [km]'
582 582 self.xlabel = 'Frac'
583 583 self.titles = ['Composition']
584 584 self.width = 3.5
585 585 self.height = 6.5
586 586 self.colorbar = False
587 587 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
588 588
589 589 def update(self, dataOut):
590 590 data = {}
591 591 meta = {}
592 592
593 593 #aux_nan=numpy.zeros(dataOut.cut,'float32')
594 594 #aux_nan[:]=numpy.nan
595 595 #data['ph'] = numpy.concatenate((aux_nan,dataOut.ph[dataOut.cut:]))
596 596 #data['eph'] = numpy.concatenate((aux_nan,dataOut.eph[dataOut.cut:]))
597 597
598 598 data['ph'] = dataOut.ph[dataOut.cut:]
599 599 data['eph'] = dataOut.eph[dataOut.cut:]
600 600 data['phe'] = dataOut.phe[dataOut.cut:]
601 601 data['ephe'] = dataOut.ephe[dataOut.cut:]
602 602
603 603 data['cut'] = dataOut.cut
604 604
605 605 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
606 606
607 607
608 608 return data, meta
609 609
610 610 def plot(self):
611 611
612 612 data = self.data[-1]
613 613
614 614 ph = data['ph']
615 615 eph = data['eph']
616 616 phe = data['phe']
617 617 ephe = data['ephe']
618 618 cut = data['cut']
619 619 self.y = self.data.yrange
620 620
621 621 self.xmin = 0
622 622 self.xmax = 1
623 623 ax = self.axes[0]
624 624
625 625 if ax.firsttime:
626 626
627 627 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
628 628 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
629 629 plt.legend(loc='lower right')
630 630 self.xstep_given = 0.2
631 631 self.ystep_given = 200
632 632 ax.yaxis.set_minor_locator(MultipleLocator(15))
633 633 ax.grid(which='minor')
634 634
635 635 else:
636 636 self.clear_figures()
637 637 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
638 638 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
639 639 plt.legend(loc='lower right')
640 640 ax.yaxis.set_minor_locator(MultipleLocator(15))
641 641 ax.grid(which='minor')
642 642
643 643 class EDensityPlot(Plot):
644 644 '''
645 645 Written by R. Flores
646 646 '''
647 647 '''
648 648 Plot for electron density
649 649 '''
650 650
651 651 CODE = 'den'
652 652 #plot_name = 'Electron Density'
653 653 plot_type = 'scatterbuffer'
654 654
655 655 def setup(self):
656 656
657 657 self.ncols = 1
658 658 self.nrows = 1
659 659 self.nplots = 1
660 660 self.ylabel = 'Range [km]'
661 661 self.xlabel = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
662 662 self.titles = ['Electron Density']
663 663 self.width = 3.5
664 664 self.height = 5.5
665 665 self.colorbar = False
666 666 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
667 667
668 668 def update(self, dataOut):
669 669 data = {}
670 670 meta = {}
671 671
672 672 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
673 673 data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS]
674 674 data['den_error'] = dataOut.sdp2[:dataOut.NSHTS]
675 675 #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS]
676 676 #print(numpy.shape(data['den_power']))
677 677 #print(numpy.shape(data['den_Faraday']))
678 678 #print(numpy.shape(data['den_error']))
679 679
680 680 data['NSHTS'] = dataOut.NSHTS
681 681
682 682 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
683 683
684 684 return data, meta
685 685
686 686 def plot(self):
687 687
688 688 y = self.data.yrange
689 689
690 690 #self.xmin = 1e3
691 691 #self.xmax = 1e7
692 692
693 693 ax = self.axes[0]
694 694
695 695 data = self.data[-1]
696 696
697 697 DenPow = data['den_power']
698 698 DenFar = data['den_Faraday']
699 699 errDenPow = data['den_error']
700 700 #errFaraday = data['err_Faraday']
701 701
702 702 NSHTS = data['NSHTS']
703 703
704 704 if self.CODE == 'denLP':
705 705 DenPowLP = data['den_LP']
706 706 errDenPowLP = data['den_LP_error']
707 707 cut = data['cut']
708 708
709 709 if ax.firsttime:
710 710 self.autoxticks=False
711 711 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
712 712 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2,linestyle='-')
713 713 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
714 714 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2,linestyle='-')
715 715
716 716 if self.CODE=='denLP':
717 717 ax.errorbar(DenPowLP[cut:], y[cut:], xerr=errDenPowLP[cut:], fmt='r^-',elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
718 718
719 719 plt.legend(loc='upper left',fontsize=8.5)
720 720 #plt.legend(loc='lower left',fontsize=8.5)
721 721 ax.set_xscale("log", nonposx='clip')
722 722 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
723 723 self.ystep_given=100
724 724 if self.CODE=='denLP':
725 725 self.ystep_given=200
726 726 ax.set_yticks(grid_y_ticks,minor=True)
727 727 locmaj = LogLocator(base=10,numticks=12)
728 728 ax.xaxis.set_major_locator(locmaj)
729 729 locmin = LogLocator(base=10.0,subs=(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12)
730 730 ax.xaxis.set_minor_locator(locmin)
731 731 ax.xaxis.set_minor_formatter(NullFormatter())
732 732 ax.grid(which='minor')
733 733
734 734 else:
735 735 dataBefore = self.data[-2]
736 736 DenPowBefore = dataBefore['den_power']
737 737 self.clear_figures()
738 738 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
739 739 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2,linestyle='-')
740 740 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
741 741 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2,linestyle='-')
742 742 ax.errorbar(DenPowBefore, y[:NSHTS], elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed")
743 743
744 744 if self.CODE=='denLP':
745 745 ax.errorbar(DenPowLP[cut:], y[cut:], fmt='r^-', xerr=errDenPowLP[cut:],elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
746 746
747 747 ax.set_xscale("log", nonposx='clip')
748 748 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
749 749 ax.set_yticks(grid_y_ticks,minor=True)
750 750 locmaj = LogLocator(base=10,numticks=12)
751 751 ax.xaxis.set_major_locator(locmaj)
752 752 locmin = LogLocator(base=10.0,subs=(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12)
753 753 ax.xaxis.set_minor_locator(locmin)
754 754 ax.xaxis.set_minor_formatter(NullFormatter())
755 755 ax.grid(which='minor')
756 756 plt.legend(loc='upper left',fontsize=8.5)
757 757 #plt.legend(loc='lower left',fontsize=8.5)
758 758
759 class RelativeDenPlot(Plot):
760 '''
761 Written by R. Flores
762 '''
763 '''
764 Plot for electron density
765 '''
766
767 CODE = 'den'
768 #plot_name = 'Electron Density'
769 plot_type = 'scatterbuffer'
770
771 def setup(self):
772
773 self.ncols = 1
774 self.nrows = 1
775 self.nplots = 1
776 self.ylabel = 'Range [km]'
777 self.xlabel = r'$\mathrm{N_e}$ Relative Electron Density ($\mathrm{1/cm^3}$)'
778 self.titles = ['Electron Density']
779 self.width = 3.5
780 self.height = 5.5
781 self.colorbar = False
782 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
783
784 def update(self, dataOut):
785 data = {}
786 meta = {}
787
788 data['den_power'] = dataOut.ph2
789 data['den_error'] = dataOut.sdp2
790
791 meta['yrange'] = dataOut.heightList
792
793 return data, meta
794
795 def plot(self):
796
797 y = self.data.yrange
798
799 ax = self.axes[0]
800
801 data = self.data[-1]
802
803 DenPow = data['den_power']
804 errDenPow = data['den_error']
805
806 if ax.firsttime:
807 self.autoxticks=False
808 ax.errorbar(DenPow, y, fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2,linestyle='-')
809
810 plt.legend(loc='upper left',fontsize=8.5)
811 #plt.legend(loc='lower left',fontsize=8.5)
812 ax.set_xscale("log", nonposx='clip')
813 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
814 self.ystep_given=100
815 ax.set_yticks(grid_y_ticks,minor=True)
816 locmaj = LogLocator(base=10,numticks=12)
817 ax.xaxis.set_major_locator(locmaj)
818 locmin = LogLocator(base=10.0,subs=(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12)
819 ax.xaxis.set_minor_locator(locmin)
820 ax.xaxis.set_minor_formatter(NullFormatter())
821 ax.grid(which='minor')
822
823 else:
824 dataBefore = self.data[-2]
825 DenPowBefore = dataBefore['den_power']
826 self.clear_figures()
827 ax.errorbar(DenPow, y, fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2,linestyle='-')
828 ax.errorbar(DenPowBefore, y, elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed")
829
830 ax.set_xscale("log", nonposx='clip')
831 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
832 ax.set_yticks(grid_y_ticks,minor=True)
833 locmaj = LogLocator(base=10,numticks=12)
834 ax.xaxis.set_major_locator(locmaj)
835 locmin = LogLocator(base=10.0,subs=(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12)
836 ax.xaxis.set_minor_locator(locmin)
837 ax.xaxis.set_minor_formatter(NullFormatter())
838 ax.grid(which='minor')
839 plt.legend(loc='upper left',fontsize=8.5)
840 #plt.legend(loc='lower left',fontsize=8.5)
841
759 842 class FaradayAnglePlot(Plot):
760 843 '''
761 844 Written by R. Flores
762 845 '''
763 846 '''
764 847 Plot for electron density
765 848 '''
766 849
767 850 CODE = 'angle'
768 851 plot_name = 'Faraday Angle'
769 852 plot_type = 'scatterbuffer'
770 853
771 854 def setup(self):
772 855
773 856 self.ncols = 1
774 857 self.nrows = 1
775 858 self.nplots = 1
776 859 self.ylabel = 'Range [km]'
777 860 self.xlabel = 'Faraday Angle (º)'
778 861 self.titles = ['Electron Density']
779 862 self.width = 3.5
780 863 self.height = 5.5
781 864 self.colorbar = False
782 865 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
783 866
784 867 def update(self, dataOut):
785 868 data = {}
786 869 meta = {}
787 870
788 871 data['angle'] = numpy.degrees(dataOut.phi)
789 872 #'''
790 873 #print(dataOut.phi_uwrp)
791 874 #print(data['angle'])
792 875 #exit(1)
793 876 #'''
794 877 data['dphi'] = dataOut.dphi_uc*10
795 878 #print(dataOut.dphi)
796 879
797 880 #data['NSHTS'] = dataOut.NSHTS
798 881
799 882 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
800 883
801 884 return data, meta
802 885
803 886 def plot(self):
804 887
805 888 data = self.data[-1]
806 889 self.x = data[self.CODE]
807 890 dphi = data['dphi']
808 891 self.y = self.data.yrange
809 892 self.xmin = -360#-180
810 893 self.xmax = 360#180
811 894 ax = self.axes[0]
812 895
813 896 if ax.firsttime:
814 897 self.autoxticks=False
815 898 #if self.CODE=='den':
816 899 ax.plot(self.x, self.y,marker='o',color='g',linewidth=1.0,markersize=2)
817 900 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
818 901
819 902 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
820 903 self.ystep_given=100
821 904 if self.CODE=='denLP':
822 905 self.ystep_given=200
823 906 ax.set_yticks(grid_y_ticks,minor=True)
824 907 ax.grid(which='minor')
825 908 #plt.tight_layout()
826 909 else:
827 910
828 911 self.clear_figures()
829 912 #if self.CODE=='den':
830 913 #print(numpy.shape(self.x))
831 914 ax.plot(self.x, self.y, marker='o',color='g',linewidth=1.0, markersize=2)
832 915 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
833 916
834 917 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
835 918 ax.set_yticks(grid_y_ticks,minor=True)
836 919 ax.grid(which='minor')
837 920
838 921 class EDensityHPPlot(EDensityPlot):
839 922 '''
840 923 Written by R. Flores
841 924 '''
842 925 '''
843 926 Plot for Electron Density Hybrid Experiment
844 927 '''
845 928
846 929 CODE = 'denLP'
847 930 plot_name = 'Electron Density'
848 931 plot_type = 'scatterbuffer'
849 932
850 933 def update(self, dataOut):
851 934 data = {}
852 935 meta = {}
853 936
854 937 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
855 938 data['den_Faraday']=dataOut.dphi[:dataOut.NSHTS]
856 939 data['den_error']=dataOut.sdp2[:dataOut.NSHTS]
857 940 data['den_LP']=dataOut.ne[:dataOut.NACF]
858 941 data['den_LP_error']=dataOut.ene[:dataOut.NACF]*dataOut.ne[:dataOut.NACF]*0.434
859 942 #self.ene=10**dataOut.ene[:dataOut.NACF]
860 943 data['NSHTS']=dataOut.NSHTS
861 944 data['cut']=dataOut.cut
862 945
863 946 return data, meta
864 947
865 948
866 949 class ACFsPlot(Plot):
867 950 '''
868 951 Written by R. Flores
869 952 '''
870 953 '''
871 954 Plot for ACFs Double Pulse Experiment
872 955 '''
873 956
874 957 CODE = 'acfs'
875 958 #plot_name = 'ACF'
876 959 plot_type = 'scatterbuffer'
877 960
878 961
879 962 def setup(self):
880 963 self.ncols = 1
881 964 self.nrows = 1
882 965 self.nplots = 1
883 966 self.ylabel = 'Range [km]'
884 967 self.xlabel = 'Lag (ms)'
885 968 self.titles = ['ACFs']
886 969 self.width = 3.5
887 970 self.height = 5.5
888 971 self.colorbar = False
889 972 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
890 973
891 974 def update(self, dataOut):
892 975 data = {}
893 976 meta = {}
894 977
895 978 data['ACFs'] = dataOut.acfs_to_plot
896 979 data['ACFs_error'] = dataOut.acfs_error_to_plot
897 980 data['lags'] = dataOut.lags_to_plot
898 981 data['Lag_contaminated_1'] = dataOut.x_igcej_to_plot
899 982 data['Lag_contaminated_2'] = dataOut.x_ibad_to_plot
900 983 data['Height_contaminated_1'] = dataOut.y_igcej_to_plot
901 984 data['Height_contaminated_2'] = dataOut.y_ibad_to_plot
902 985
903 986 meta['yrange'] = numpy.array([])
904 987 #meta['NSHTS'] = dataOut.NSHTS
905 988 #meta['DPL'] = dataOut.DPL
906 989 data['NSHTS'] = dataOut.NSHTS #This is metadata
907 990 data['DPL'] = dataOut.DPL #This is metadata
908 991
909 992 return data, meta
910 993
911 994 def plot(self):
912 995
913 996 data = self.data[-1]
914 997 #NSHTS = self.meta['NSHTS']
915 998 #DPL = self.meta['DPL']
916 999 NSHTS = data['NSHTS'] #This is metadata
917 1000 DPL = data['DPL'] #This is metadata
918 1001
919 1002 lags = data['lags']
920 1003 ACFs = data['ACFs']
921 1004 errACFs = data['ACFs_error']
922 1005 BadLag1 = data['Lag_contaminated_1']
923 1006 BadLag2 = data['Lag_contaminated_2']
924 1007 BadHei1 = data['Height_contaminated_1']
925 1008 BadHei2 = data['Height_contaminated_2']
926 1009
927 1010 self.xmin = 0.0
928 1011 self.xmax = 2.0
929 1012 self.y = ACFs
930 1013
931 1014 ax = self.axes[0]
932 1015
933 1016 if ax.firsttime:
934 1017
935 1018 for i in range(NSHTS):
936 1019 x_aux = numpy.isfinite(lags[i,:])
937 1020 y_aux = numpy.isfinite(ACFs[i,:])
938 1021 yerr_aux = numpy.isfinite(errACFs[i,:])
939 1022 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
940 1023 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
941 1024 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
942 1025 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
943 1026 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
944 1027 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',marker='o',linewidth=1.0,markersize=2)
945 1028 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
946 1029 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
947 1030
948 1031 self.xstep_given = (self.xmax-self.xmin)/(DPL-1)
949 1032 self.ystep_given = 50
950 1033 ax.yaxis.set_minor_locator(MultipleLocator(15))
951 1034 ax.grid(which='minor')
952 1035
953 1036 else:
954 1037 self.clear_figures()
955 1038 for i in range(NSHTS):
956 1039 x_aux = numpy.isfinite(lags[i,:])
957 1040 y_aux = numpy.isfinite(ACFs[i,:])
958 1041 yerr_aux = numpy.isfinite(errACFs[i,:])
959 1042 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
960 1043 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
961 1044 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
962 1045 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
963 1046 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
964 1047 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],linewidth=1.0,markersize=2,color='b',marker='o')
965 1048 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
966 1049 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
967 1050 ax.yaxis.set_minor_locator(MultipleLocator(15))
968 1051
969 1052 class ACFsLPPlot(Plot):
970 1053 '''
971 1054 Written by R. Flores
972 1055 '''
973 1056 '''
974 1057 Plot for ACFs Double Pulse Experiment
975 1058 '''
976 1059
977 1060 CODE = 'acfs_LP'
978 1061 #plot_name = 'ACF'
979 1062 plot_type = 'scatterbuffer'
980 1063
981 1064
982 1065 def setup(self):
983 1066 self.ncols = 1
984 1067 self.nrows = 1
985 1068 self.nplots = 1
986 1069 self.ylabel = 'Range [km]'
987 1070 self.xlabel = 'Lag (ms)'
988 1071 self.titles = ['ACFs']
989 1072 self.width = 3.5
990 1073 self.height = 5.5
991 1074 self.colorbar = False
992 1075 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
993 1076
994 1077 def update(self, dataOut):
995 1078 data = {}
996 1079 meta = {}
997 1080
998 1081 aux=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
999 1082 errors=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
1000 1083 lags_LP_to_plot=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
1001 1084
1002 1085 for i in range(dataOut.NACF):
1003 1086 for j in range(dataOut.IBITS):
1004 1087 if numpy.abs(dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0])<1.0:
1005 1088 aux[i,j]=dataOut.output_LP_integrated.real[j,i,0]/dataOut.output_LP_integrated.real[0,i,0]
1006 1089 aux[i,j]=max(min(aux[i,j],1.0),-1.0)*dataOut.DH+dataOut.heightList[i]
1007 1090 lags_LP_to_plot[i,j]=dataOut.lags_LP[j]
1008 1091 errors[i,j]=dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0]*dataOut.DH
1009 1092 else:
1010 1093 aux[i,j]=numpy.nan
1011 1094 lags_LP_to_plot[i,j]=numpy.nan
1012 1095 errors[i,j]=numpy.nan
1013 1096
1014 1097 data['ACFs'] = aux
1015 1098 data['ACFs_error'] = errors
1016 1099 data['lags'] = lags_LP_to_plot
1017 1100
1018 1101 meta['yrange'] = numpy.array([])
1019 1102 #meta['NACF'] = dataOut.NACF
1020 1103 #meta['NLAG'] = dataOut.NLAG
1021 1104 data['NACF'] = dataOut.NACF #This is metadata
1022 1105 data['NLAG'] = dataOut.NLAG #This is metadata
1023 1106
1024 1107 return data, meta
1025 1108
1026 1109 def plot(self):
1027 1110
1028 1111 data = self.data[-1]
1029 1112 #NACF = self.meta['NACF']
1030 1113 #NLAG = self.meta['NLAG']
1031 1114 NACF = data['NACF'] #This is metadata
1032 1115 NLAG = data['NLAG'] #This is metadata
1033 1116
1034 1117 lags = data['lags']
1035 1118 ACFs = data['ACFs']
1036 1119 errACFs = data['ACFs_error']
1037 1120
1038 1121 self.xmin = 0.0
1039 1122 self.xmax = 1.5
1040 1123
1041 1124 self.y = ACFs
1042 1125
1043 1126 ax = self.axes[0]
1044 1127
1045 1128 if ax.firsttime:
1046 1129
1047 1130 for i in range(NACF):
1048 1131 x_aux = numpy.isfinite(lags[i,:])
1049 1132 y_aux = numpy.isfinite(ACFs[i,:])
1050 1133 yerr_aux = numpy.isfinite(errACFs[i,:])
1051 1134
1052 1135 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1053 1136 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1054 1137
1055 1138 #self.xstep_given = (self.xmax-self.xmin)/(self.data.NLAG-1)
1056 1139 self.xstep_given=0.3
1057 1140 self.ystep_given = 200
1058 1141 ax.yaxis.set_minor_locator(MultipleLocator(15))
1059 1142 ax.grid(which='minor')
1060 1143
1061 1144 else:
1062 1145 self.clear_figures()
1063 1146
1064 1147 for i in range(NACF):
1065 1148 x_aux = numpy.isfinite(lags[i,:])
1066 1149 y_aux = numpy.isfinite(ACFs[i,:])
1067 1150 yerr_aux = numpy.isfinite(errACFs[i,:])
1068 1151
1069 1152 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1070 1153 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1071 1154
1072 1155 ax.yaxis.set_minor_locator(MultipleLocator(15))
1073 1156
1074 1157
1075 1158 class CrossProductsPlot(Plot):
1076 1159 '''
1077 1160 Written by R. Flores
1078 1161 '''
1079 1162 '''
1080 1163 Plot for cross products
1081 1164 '''
1082 1165
1083 1166 CODE = 'crossprod'
1084 1167 plot_name = 'Cross Products'
1085 1168 plot_type = 'scatterbuffer'
1086 1169
1087 1170 def setup(self):
1088 1171
1089 1172 self.ncols = 3
1090 1173 self.nrows = 1
1091 1174 self.nplots = 3
1092 1175 self.ylabel = 'Range [km]'
1093 1176 self.titles = []
1094 1177 self.width = 3.5*self.nplots
1095 1178 self.height = 5.5
1096 1179 self.colorbar = False
1097 1180 self.plots_adjust.update({'wspace':.3, 'left': 0.12, 'right': 0.92, 'bottom': 0.1})
1098 1181
1099 1182
1100 1183 def update(self, dataOut):
1101 1184
1102 1185 data = {}
1103 1186 meta = {}
1104 1187
1105 1188 data['crossprod'] = dataOut.crossprods
1106 1189 data['NDP'] = dataOut.NDP
1107 1190
1108 1191 return data, meta
1109 1192
1110 1193 def plot(self):
1111 1194
1112 1195 NDP = self.data['NDP'][-1]
1113 1196 x = self.data['crossprod'][:,-1,:,:,:,:]
1114 1197 y = self.data.yrange[0:NDP]
1115 1198
1116 1199 for n, ax in enumerate(self.axes):
1117 1200
1118 1201 self.xmin=numpy.min(numpy.concatenate((x[n][0,20:30,0,0],x[n][1,20:30,0,0],x[n][2,20:30,0,0],x[n][3,20:30,0,0])))
1119 1202 self.xmax=numpy.max(numpy.concatenate((x[n][0,20:30,0,0],x[n][1,20:30,0,0],x[n][2,20:30,0,0],x[n][3,20:30,0,0])))
1120 1203
1121 1204 if ax.firsttime:
1122 1205
1123 1206 self.autoxticks=False
1124 1207 if n==0:
1125 1208 label1='kax'
1126 1209 label2='kay'
1127 1210 label3='kbx'
1128 1211 label4='kby'
1129 1212 self.xlimits=[(self.xmin,self.xmax)]
1130 1213 elif n==1:
1131 1214 label1='kax2'
1132 1215 label2='kay2'
1133 1216 label3='kbx2'
1134 1217 label4='kby2'
1135 1218 self.xlimits.append((self.xmin,self.xmax))
1136 1219 elif n==2:
1137 1220 label1='kaxay'
1138 1221 label2='kbxby'
1139 1222 label3='kaxbx'
1140 1223 label4='kaxby'
1141 1224 self.xlimits.append((self.xmin,self.xmax))
1142 1225
1143 1226 ax.plotline1 = ax.plot(x[n][0,:,0,0], y, color='r',linewidth=2.0, label=label1)
1144 1227 ax.plotline2 = ax.plot(x[n][1,:,0,0], y, color='k',linewidth=2.0, label=label2)
1145 1228 ax.plotline3 = ax.plot(x[n][2,:,0,0], y, color='b',linewidth=2.0, label=label3)
1146 1229 ax.plotline4 = ax.plot(x[n][3,:,0,0], y, color='m',linewidth=2.0, label=label4)
1147 1230 ax.legend(loc='upper right')
1148 1231 ax.set_xlim(self.xmin, self.xmax)
1149 1232 self.titles.append('{}'.format(self.plot_name.upper()))
1150 1233
1151 1234 else:
1152 1235
1153 1236 if n==0:
1154 1237 self.xlimits=[(self.xmin,self.xmax)]
1155 1238 else:
1156 1239 self.xlimits.append((self.xmin,self.xmax))
1157 1240
1158 1241 ax.set_xlim(self.xmin, self.xmax)
1159 1242
1160 1243 ax.plotline1[0].set_data(x[n][0,:,0,0],y)
1161 1244 ax.plotline2[0].set_data(x[n][1,:,0,0],y)
1162 1245 ax.plotline3[0].set_data(x[n][2,:,0,0],y)
1163 1246 ax.plotline4[0].set_data(x[n][3,:,0,0],y)
1164 1247 self.titles.append('{}'.format(self.plot_name.upper()))
1165 1248
1166 1249
1167 1250 class CrossProductsLPPlot(Plot):
1168 1251 '''
1169 1252 Written by R. Flores
1170 1253 '''
1171 1254 '''
1172 1255 Plot for cross products LP
1173 1256 '''
1174 1257
1175 1258 CODE = 'crossprodslp'
1176 1259 plot_name = 'Cross Products LP'
1177 1260 plot_type = 'scatterbuffer'
1178 1261
1179 1262
1180 1263 def setup(self):
1181 1264
1182 1265 self.ncols = 2
1183 1266 self.nrows = 1
1184 1267 self.nplots = 2
1185 1268 self.ylabel = 'Range [km]'
1186 1269 self.xlabel = 'dB'
1187 1270 self.width = 3.5*self.nplots
1188 1271 self.height = 5.5
1189 1272 self.colorbar = False
1190 1273 self.titles = []
1191 1274 self.plots_adjust.update({'wspace': .8 ,'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1192 1275
1193 1276 def update(self, dataOut):
1194 1277 data = {}
1195 1278 meta = {}
1196 1279
1197 1280 data['crossprodslp'] = 10*numpy.log10(numpy.abs(dataOut.output_LP))
1198 1281
1199 1282 data['NRANGE'] = dataOut.NRANGE #This is metadata
1200 1283 data['NLAG'] = dataOut.NLAG #This is metadata
1201 1284
1202 1285 return data, meta
1203 1286
1204 1287 def plot(self):
1205 1288
1206 1289 NRANGE = self.data['NRANGE'][-1]
1207 1290 NLAG = self.data['NLAG'][-1]
1208 1291
1209 1292 x = self.data[self.CODE][:,-1,:,:]
1210 1293 self.y = self.data.yrange[0:NRANGE]
1211 1294
1212 1295 label_array=numpy.array(['lag '+ str(x) for x in range(NLAG)])
1213 1296 color_array=['r','k','g','b','c','m','y','orange','steelblue','purple','peru','darksalmon','grey','limegreen','olive','midnightblue']
1214 1297
1215 1298
1216 1299 for n, ax in enumerate(self.axes):
1217 1300
1218 1301 self.xmin=28#30
1219 1302 self.xmax=70#70
1220 1303 #self.xmin=numpy.min(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1221 1304 #self.xmax=numpy.max(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1222 1305
1223 1306 if ax.firsttime:
1224 1307
1225 1308 self.autoxticks=False
1226 1309 if n == 0:
1227 1310 self.plotline_array=numpy.zeros((2,NLAG),dtype=object)
1228 1311
1229 1312 for i in range(NLAG):
1230 1313 self.plotline_array[n,i], = ax.plot(x[i,:,n], self.y, color=color_array[i],linewidth=1.0, label=label_array[i])
1231 1314
1232 1315 ax.legend(loc='upper right')
1233 1316 ax.set_xlim(self.xmin, self.xmax)
1234 1317 if n==0:
1235 1318 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1236 1319 if n==1:
1237 1320 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1238 1321 else:
1239 1322 for i in range(NLAG):
1240 1323 self.plotline_array[n,i].set_data(x[i,:,n],self.y)
1241 1324
1242 1325 if n==0:
1243 1326 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1244 1327 if n==1:
1245 1328 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1246 1329
1247 1330
1248 1331 class NoiseDPPlot(NoisePlot):
1249 1332 '''
1250 1333 Written by R. Flores
1251 1334 '''
1252 1335 '''
1253 1336 Plot for noise Double Pulse
1254 1337 '''
1255 1338
1256 1339 CODE = 'noise'
1257 1340 #plot_name = 'Noise'
1258 1341 #plot_type = 'scatterbuffer'
1259 1342
1260 1343 def update(self, dataOut):
1261 1344
1262 1345 data = {}
1263 1346 meta = {}
1264 1347 data['noise'] = 10*numpy.log10(dataOut.noise_final)
1265 1348
1266 1349 return data, meta
1267 1350
1268 1351
1269 1352 class XmitWaveformPlot(Plot):
1270 1353 '''
1271 1354 Written by R. Flores
1272 1355 '''
1273 1356 '''
1274 1357 Plot for xmit waveform
1275 1358 '''
1276 1359
1277 1360 CODE = 'xmit'
1278 1361 plot_name = 'Xmit Waveform'
1279 1362 plot_type = 'scatterbuffer'
1280 1363
1281 1364
1282 1365 def setup(self):
1283 1366
1284 1367 self.ncols = 1
1285 1368 self.nrows = 1
1286 1369 self.nplots = 1
1287 1370 self.ylabel = ''
1288 1371 self.xlabel = 'Number of Lag'
1289 1372 self.width = 5.5
1290 1373 self.height = 3.5
1291 1374 self.colorbar = False
1292 1375 self.plots_adjust.update({'right': 0.85 })
1293 1376 self.titles = [self.plot_name]
1294 1377 #self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1295 1378
1296 1379 #if not self.titles:
1297 1380 #self.titles = self.data.parameters \
1298 1381 #if self.data.parameters else ['{}'.format(self.plot_name.upper())]
1299 1382
1300 1383 def update(self, dataOut):
1301 1384
1302 1385 data = {}
1303 1386 meta = {}
1304 1387
1305 1388 y_1=numpy.arctan2(dataOut.output_LP[:,0,2].imag,dataOut.output_LP[:,0,2].real)* 180 / (numpy.pi*10)
1306 1389 y_2=numpy.abs(dataOut.output_LP[:,0,2])
1307 1390 norm=numpy.max(y_2)
1308 1391 norm=max(norm,0.1)
1309 1392 y_2=y_2/norm
1310 1393
1311 1394 meta['yrange'] = numpy.array([])
1312 1395
1313 1396 data['xmit'] = numpy.vstack((y_1,y_2))
1314 1397 data['NLAG'] = dataOut.NLAG
1315 1398
1316 1399 return data, meta
1317 1400
1318 1401 def plot(self):
1319 1402
1320 1403 data = self.data[-1]
1321 1404 NLAG = data['NLAG']
1322 1405 x = numpy.arange(0,NLAG,1,'float32')
1323 1406 y = data['xmit']
1324 1407
1325 1408 self.xmin = 0
1326 1409 self.xmax = NLAG-1
1327 1410 self.ymin = -1.0
1328 1411 self.ymax = 1.0
1329 1412 ax = self.axes[0]
1330 1413
1331 1414 if ax.firsttime:
1332 1415 ax.plotline0=ax.plot(x,y[0,:],color='blue')
1333 1416 ax.plotline1=ax.plot(x,y[1,:],color='red')
1334 1417 secax=ax.secondary_xaxis(location=0.5)
1335 1418 secax.xaxis.tick_bottom()
1336 1419 secax.tick_params( labelleft=False, labeltop=False,
1337 1420 labelright=False, labelbottom=False)
1338 1421
1339 1422 self.xstep_given = 3
1340 1423 self.ystep_given = .25
1341 1424 secax.set_xticks(numpy.linspace(self.xmin, self.xmax, 6)) #only works on matplotlib.version>3.2
1342 1425
1343 1426 else:
1344 1427 ax.plotline0[0].set_data(x,y[0,:])
1345 1428 ax.plotline1[0].set_data(x,y[1,:])
Show file before | Show file after | Hide comments
@@ -1,626 +1,636
1 1 import os
2 2 import time
3 3 import datetime
4 4
5 5 import numpy
6 6 import h5py
7 7
8 8 import schainpy.admin
9 9 from schainpy.model.data.jrodata import *
10 10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 11 from schainpy.model.io.jroIO_base import *
12 12 from schainpy.utils import log
13 13
14 14
15 15 class HDFReader(Reader, ProcessingUnit):
16 16 """Processing unit to read HDF5 format files
17 17
18 18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 20 attributes.
21 21 It is possible to read any HDF5 file by given the structure in the `description`
22 22 parameter, also you can add extra values to metadata with the parameter `extras`.
23 23
24 24 Parameters:
25 25 -----------
26 26 path : str
27 27 Path where files are located.
28 28 startDate : date
29 29 Start date of the files
30 30 endDate : list
31 31 End date of the files
32 32 startTime : time
33 33 Start time of the files
34 34 endTime : time
35 35 End time of the files
36 36 description : dict, optional
37 37 Dictionary with the description of the HDF5 file
38 38 extras : dict, optional
39 39 Dictionary with extra metadata to be be added to `dataOut`
40 40
41 41 Examples
42 42 --------
43 43
44 44 desc = {
45 45 'Data': {
46 46 'data_output': ['u', 'v', 'w'],
47 47 'utctime': 'timestamps',
48 48 } ,
49 49 'Metadata': {
50 50 'heightList': 'heights'
51 51 }
52 52 }
53 53
54 54 desc = {
55 55 'Data': {
56 56 'data_output': 'winds',
57 57 'utctime': 'timestamps'
58 58 },
59 59 'Metadata': {
60 60 'heightList': 'heights'
61 61 }
62 62 }
63 63
64 64 extras = {
65 65 'timeZone': 300
66 66 }
67 67
68 68 reader = project.addReadUnit(
69 69 name='HDFReader',
70 70 path='/path/to/files',
71 71 startDate='2019/01/01',
72 72 endDate='2019/01/31',
73 73 startTime='00:00:00',
74 74 endTime='23:59:59',
75 75 # description=json.dumps(desc),
76 76 # extras=json.dumps(extras),
77 77 )
78 78
79 79 """
80 80
81 81 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
82 82
83 83 def __init__(self):
84 84 ProcessingUnit.__init__(self)
85 85 self.dataOut = Parameters()
86 86 self.ext = ".hdf5"
87 87 self.optchar = "D"
88 88 self.meta = {}
89 89 self.data = {}
90 90 self.open_file = h5py.File
91 91 self.open_mode = 'r'
92 92 self.description = {}
93 93 self.extras = {}
94 94 self.filefmt = "*%Y%j***"
95 95 self.folderfmt = "*%Y%j"
96 96 self.utcoffset = 0
97 97
98 98 def setup(self, **kwargs):
99 99
100 100 self.set_kwargs(**kwargs)
101 101 if not self.ext.startswith('.'):
102 102 self.ext = '.{}'.format(self.ext)
103 103
104 104 if self.online:
105 105 log.log("Searching files in online mode...", self.name)
106 106
107 107 for nTries in range(self.nTries):
108 108 fullpath = self.searchFilesOnLine(self.path, self.startDate,
109 109 self.endDate, self.expLabel, self.ext, self.walk,
110 110 self.filefmt, self.folderfmt)
111 111 try:
112 112 fullpath = next(fullpath)
113 113 except:
114 114 fullpath = None
115 115
116 116 if fullpath:
117 117 break
118 118
119 119 log.warning(
120 120 'Waiting {} sec for a valid file in {}: try {} ...'.format(
121 121 self.delay, self.path, nTries + 1),
122 122 self.name)
123 123 time.sleep(self.delay)
124 124
125 125 if not(fullpath):
126 126 raise schainpy.admin.SchainError(
127 127 'There isn\'t any valid file in {}'.format(self.path))
128 128
129 129 pathname, filename = os.path.split(fullpath)
130 130 self.year = int(filename[1:5])
131 131 self.doy = int(filename[5:8])
132 132 self.set = int(filename[8:11]) - 1
133 133 else:
134 134 log.log("Searching files in {}".format(self.path), self.name)
135 135 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
136 136 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
137 137
138 138 self.setNextFile()
139 139
140 140 return
141 141
142 142 def readFirstHeader(self):
143 143 '''Read metadata and data'''
144 144
145 145 self.__readMetadata()
146 146 self.__readData()
147 147 self.__setBlockList()
148 148
149 149 if 'type' in self.meta:
150 150 self.dataOut = eval(self.meta['type'])()
151 151
152 152 for attr in self.meta:
153 153 setattr(self.dataOut, attr, self.meta[attr])
154 154
155 155 self.blockIndex = 0
156 156
157 157 return
158 158
159 159 def __setBlockList(self):
160 160 '''
161 161 Selects the data within the times defined
162 162
163 163 self.fp
164 164 self.startTime
165 165 self.endTime
166 166 self.blockList
167 167 self.blocksPerFile
168 168
169 169 '''
170 170
171 171 startTime = self.startTime
172 172 endTime = self.endTime
173 173 thisUtcTime = self.data['utctime'] + self.utcoffset
174 174 self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
175 175 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
176 176
177 177 thisDate = thisDatetime.date()
178 178 thisTime = thisDatetime.time()
179 179
180 180 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
181 181 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
182 182
183 183 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
184 184
185 185 self.blockList = ind
186 186 self.blocksPerFile = len(ind)
187 187 return
188 188
189 189 def __readMetadata(self):
190 190 '''
191 191 Reads Metadata
192 192 '''
193 193
194 194 meta = {}
195 195
196 196 if self.description:
197 197 for key, value in self.description['Metadata'].items():
198 198 meta[key] = self.fp[value][()]
199 199 else:
200 200 grp = self.fp['Metadata']
201 201 for name in grp:
202 202 meta[name] = grp[name][()]
203 203
204 204 if self.extras:
205 205 for key, value in self.extras.items():
206 206 meta[key] = value
207 207 self.meta = meta
208 208
209 209 return
210 210
211 211 def __readData(self):
212 212
213 213 data = {}
214 214
215 215 if self.description:
216 216 for key, value in self.description['Data'].items():
217 217 if isinstance(value, str):
218 218 if isinstance(self.fp[value], h5py.Dataset):
219 219 data[key] = self.fp[value][()]
220 220 elif isinstance(self.fp[value], h5py.Group):
221 221 array = []
222 222 for ch in self.fp[value]:
223 223 array.append(self.fp[value][ch][()])
224 224 data[key] = numpy.array(array)
225 225 elif isinstance(value, list):
226 226 array = []
227 227 for ch in value:
228 228 array.append(self.fp[ch][()])
229 229 data[key] = numpy.array(array)
230 230 else:
231 231 grp = self.fp['Data']
232 232 for name in grp:
233 233 if isinstance(grp[name], h5py.Dataset):
234 234 array = grp[name][()]
235 235 elif isinstance(grp[name], h5py.Group):
236 236 array = []
237 237 for ch in grp[name]:
238 238 array.append(grp[name][ch][()])
239 239 array = numpy.array(array)
240 240 else:
241 241 log.warning('Unknown type: {}'.format(name))
242 242
243 243 if name in self.description:
244 244 key = self.description[name]
245 245 else:
246 246 key = name
247 247 data[key] = array
248 248
249 249 self.data = data
250 250 return
251 251
252 252 def getData(self):
253 253
254 254 for attr in self.data:
255 255 if self.data[attr].ndim == 1:
256 256 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
257 257 else:
258 258 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
259 259
260 260 self.dataOut.flagNoData = False
261 261 self.blockIndex += 1
262 262
263 263 log.log("Block No. {}/{} -> {}".format(
264 264 self.blockIndex,
265 265 self.blocksPerFile,
266 266 self.dataOut.datatime.ctime()), self.name)
267 267
268 268 return
269 269
270 270 def run(self, **kwargs):
271 271
272 272 if not(self.isConfig):
273 273 self.setup(**kwargs)
274 274 self.isConfig = True
275 275
276 276 if self.blockIndex == self.blocksPerFile:
277 277 self.setNextFile()
278 278
279 279 self.getData()
280 280
281 281 return
282 282
283 283 @MPDecorator
284 284 class HDFWriter(Operation):
285 285 """Operation to write HDF5 files.
286 286
287 287 The HDF5 file contains by default two groups Data and Metadata where
288 288 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
289 289 parameters, data attributes are normaly time dependent where the metadata
290 290 are not.
291 291 It is possible to customize the structure of the HDF5 file with the
292 292 optional description parameter see the examples.
293 293
294 294 Parameters:
295 295 -----------
296 296 path : str
297 297 Path where files will be saved.
298 298 blocksPerFile : int
299 299 Number of blocks per file
300 300 metadataList : list
301 301 List of the dataOut attributes that will be saved as metadata
302 302 dataList : int
303 303 List of the dataOut attributes that will be saved as data
304 304 setType : bool
305 305 If True the name of the files corresponds to the timestamp of the data
306 306 description : dict, optional
307 307 Dictionary with the desired description of the HDF5 file
308 308
309 309 Examples
310 310 --------
311 311
312 312 desc = {
313 313 'data_output': {'winds': ['z', 'w', 'v']},
314 314 'utctime': 'timestamps',
315 315 'heightList': 'heights'
316 316 }
317 317 desc = {
318 318 'data_output': ['z', 'w', 'v'],
319 319 'utctime': 'timestamps',
320 320 'heightList': 'heights'
321 321 }
322 322 desc = {
323 323 'Data': {
324 324 'data_output': 'winds',
325 325 'utctime': 'timestamps'
326 326 },
327 327 'Metadata': {
328 328 'heightList': 'heights'
329 329 }
330 330 }
331 331
332 332 writer = proc_unit.addOperation(name='HDFWriter')
333 333 writer.addParameter(name='path', value='/path/to/file')
334 334 writer.addParameter(name='blocksPerFile', value='32')
335 335 writer.addParameter(name='metadataList', value='heightList,timeZone')
336 336 writer.addParameter(name='dataList',value='data_output,utctime')
337 337 # writer.addParameter(name='description',value=json.dumps(desc))
338 338
339 339 """
340 340
341 341 ext = ".hdf5"
342 342 optchar = "D"
343 343 filename = None
344 344 path = None
345 345 setFile = None
346 346 fp = None
347 347 firsttime = True
348 348 #Configurations
349 349 blocksPerFile = None
350 350 blockIndex = None
351 351 dataOut = None
352 352 #Data Arrays
353 353 dataList = None
354 354 metadataList = None
355 355 currentDay = None
356 356 lastTime = None
357 357
358 358 def __init__(self):
359 359
360 360 Operation.__init__(self)
361 361 return
362 362
363 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None):
363 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None, uniqueChannel=False):
364 364 self.path = path
365 365 self.blocksPerFile = blocksPerFile
366 366 self.metadataList = metadataList
367 367 self.dataList = [s.strip() for s in dataList]
368 368 self.setType = setType
369 369 self.description = description
370 self.uniqueChannel = uniqueChannel
370 371
371 372 if self.metadataList is None:
372 373 self.metadataList = self.dataOut.metadata_list
373 374
374 375 tableList = []
375 376 dsList = []
376 377
377 378 for i in range(len(self.dataList)):
378 379 dsDict = {}
379 380 if hasattr(self.dataOut, self.dataList[i]):
380 381 dataAux = getattr(self.dataOut, self.dataList[i])
381 382 dsDict['variable'] = self.dataList[i]
382 383 else:
383 384 log.warning('Attribute {} not found in dataOut', self.name)
384 385 continue
385 386
386 387 if dataAux is None:
387 388 continue
388 389 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
389 390 dsDict['nDim'] = 0
390 391 else:
392 if uniqueChannel: #Creates extra dimension to avoid the creation of multiple channels
393 dataAux = numpy.expand_dims(dataAux, axis=0)
394
391 395 dsDict['nDim'] = len(dataAux.shape)
392 396 dsDict['shape'] = dataAux.shape
393 397 dsDict['dsNumber'] = dataAux.shape[0]
394 398 dsDict['dtype'] = dataAux.dtype
395 399
396 400 dsList.append(dsDict)
397 401
398 402 self.dsList = dsList
399 403 self.currentDay = self.dataOut.datatime.date()
400 404
401 405 def timeFlag(self):
402 406 currentTime = self.dataOut.utctime
403 407 timeTuple = time.localtime(currentTime)
404 408 dataDay = timeTuple.tm_yday
405 409
406 410 if self.lastTime is None:
407 411 self.lastTime = currentTime
408 412 self.currentDay = dataDay
409 413 return False
410 414
411 415 timeDiff = currentTime - self.lastTime
412 416
413 417 #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
414 418 if dataDay != self.currentDay:
415 419 self.currentDay = dataDay
416 420 return True
417 421 elif timeDiff > 3*60*60:
418 422 self.lastTime = currentTime
419 423 return True
420 424 else:
421 425 self.lastTime = currentTime
422 426 return False
423 427
424 428 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
425 dataList=[], setType=None, description={}):
429 dataList=[], setType=None, description={}, uniqueChannel= False):
426 430
427 431 self.dataOut = dataOut
428 432 if not(self.isConfig):
429 433 self.setup(path=path, blocksPerFile=blocksPerFile,
430 434 metadataList=metadataList, dataList=dataList,
431 setType=setType, description=description)
435 setType=setType, description=description, uniqueChannel=uniqueChannel)
432 436
433 437 self.isConfig = True
434 438 self.setNextFile()
435 439
436 440 self.putData()
441
437 442 return
438 443
439 444 def setNextFile(self):
440 445
441 446 ext = self.ext
442 447 path = self.path
443 448 setFile = self.setFile
444 449
445 450 timeTuple = time.localtime(self.dataOut.utctime)
446 451 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
447 452 fullpath = os.path.join(path, subfolder)
448 453
449 454 if os.path.exists(fullpath):
450 455 filesList = os.listdir(fullpath)
451 456 filesList = [k for k in filesList if k.startswith(self.optchar)]
452 457 if len( filesList ) > 0:
453 458 filesList = sorted(filesList, key=str.lower)
454 459 filen = filesList[-1]
455 460 # el filename debera tener el siguiente formato
456 461 # 0 1234 567 89A BCDE (hex)
457 462 # x YYYY DDD SSS .ext
458 463 if isNumber(filen[8:11]):
459 464 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
460 465 else:
461 466 setFile = -1
462 467 else:
463 468 setFile = -1 #inicializo mi contador de seteo
464 469 else:
465 470 os.makedirs(fullpath)
466 471 setFile = -1 #inicializo mi contador de seteo
467 472
468 473 if self.setType is None:
469 474 setFile += 1
470 475 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
471 476 timeTuple.tm_year,
472 477 timeTuple.tm_yday,
473 478 setFile,
474 479 ext )
475 480 else:
476 481 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
477 482 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
478 483 timeTuple.tm_year,
479 484 timeTuple.tm_yday,
480 485 setFile,
481 486 ext )
482 487
483 488 self.filename = os.path.join( path, subfolder, file )
484 489
485 490 #Setting HDF5 File
486 491 self.fp = h5py.File(self.filename, 'w')
487 492 #write metadata
488 493 self.writeMetadata(self.fp)
489 494 #Write data
490 495 self.writeData(self.fp)
491 496
492 497 def getLabel(self, name, x=None):
493
498 #print("x: ", x)
494 499 if x is None:
495 500 if 'Data' in self.description:
496 501 data = self.description['Data']
497 502 if 'Metadata' in self.description:
498 503 data.update(self.description['Metadata'])
499 504 else:
500 505 data = self.description
501 506 if name in data:
502 507 if isinstance(data[name], str):
503 508 return data[name]
504 509 elif isinstance(data[name], list):
505 510 return None
506 511 elif isinstance(data[name], dict):
507 512 for key, value in data[name].items():
508 513 return key
509 514 return name
510 515 else:
511 516 if 'Metadata' in self.description:
512 517 meta = self.description['Metadata']
513 518 else:
514 519 meta = self.description
515 520 if name in meta:
516 521 if isinstance(meta[name], list):
517 522 return meta[name][x]
518 523 elif isinstance(meta[name], dict):
519 524 for key, value in meta[name].items():
520 525 return value[x]
521 526 if 'cspc' in name:
522 527 return 'pair{:02d}'.format(x)
523 528 else:
524 529 return 'channel{:02d}'.format(x)
525 530
526 531 def writeMetadata(self, fp):
527 532
528 533 if self.description:
529 534 if 'Metadata' in self.description:
530 535 grp = fp.create_group('Metadata')
531 536 else:
532 537 grp = fp
533 538 else:
534 539 grp = fp.create_group('Metadata')
535 540
536 541 for i in range(len(self.metadataList)):
537 542 if not hasattr(self.dataOut, self.metadataList[i]):
538 543 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
539 544 continue
540 545 value = getattr(self.dataOut, self.metadataList[i])
541 546 if isinstance(value, bool):
542 547 if value is True:
543 548 value = 1
544 549 else:
545 550 value = 0
546 551 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
547 552 return
548 553
549 554 def writeData(self, fp):
550 555
551 556 if self.description:
552 557 if 'Data' in self.description:
553 558 grp = fp.create_group('Data')
554 559 else:
555 560 grp = fp
556 561 else:
557 562 grp = fp.create_group('Data')
558 563
559 564 dtsets = []
560 565 data = []
561
566 #print("self.dsList: ", self.dsList)
562 567 for dsInfo in self.dsList:
563 568 if dsInfo['nDim'] == 0:
564 569 ds = grp.create_dataset(
565 570 self.getLabel(dsInfo['variable']),
566 571 (self.blocksPerFile, ),
567 572 chunks=True,
568 573 dtype=numpy.float64)
569 574 dtsets.append(ds)
570 575 data.append((dsInfo['variable'], -1))
571 576 else:
572 577 label = self.getLabel(dsInfo['variable'])
573 578 if label is not None:
574 579 sgrp = grp.create_group(label)
575 580 else:
576 581 sgrp = grp
577 582 for i in range(dsInfo['dsNumber']):
578 583 ds = sgrp.create_dataset(
579 584 self.getLabel(dsInfo['variable'], i),
580 585 (self.blocksPerFile, ) + dsInfo['shape'][1:],
581 586 chunks=True,
582 587 dtype=dsInfo['dtype'])
583 588 dtsets.append(ds)
584 589 data.append((dsInfo['variable'], i))
590
591 if self.uniqueChannel: #Deletes extra dimension created to avoid the creation of multiple channels
592 dataAux = getattr(self.dataOut, dsInfo['variable'])
593 dataAux = dataAux[0]
594
585 595 fp.flush()
586 596
587 597 log.log('Creating file: {}'.format(fp.filename), self.name)
588 598
589 599 self.ds = dtsets
590 600 self.data = data
591 601 self.firsttime = True
592 602 self.blockIndex = 0
593 603 return
594 604
595 605 def putData(self):
596 606
597 607 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
598 608 self.closeFile()
599 609 self.setNextFile()
600 610
601 611 for i, ds in enumerate(self.ds):
602 612 attr, ch = self.data[i]
603 613 if ch == -1:
604 614 ds[self.blockIndex] = getattr(self.dataOut, attr)
605 615 else:
606 616 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
607 617
608 618 self.fp.flush()
609 619 self.blockIndex += 1
610 620 log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
611 621
612 622 return
613 623
614 624 def closeFile(self):
615 625
616 626 if self.blockIndex != self.blocksPerFile:
617 627 for ds in self.ds:
618 628 ds.resize(self.blockIndex, axis=0)
619 629
620 630 if self.fp:
621 631 self.fp.flush()
622 632 self.fp.close()
623 633
624 634 def close(self):
625 635
626 636 self.closeFile()
Show file before | Show file after | Hide comments
@@ -1,982 +1,1049
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 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15
16 16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 17 from schainpy.model.data.jrodata import Spectra
18 18 from schainpy.model.data.jrodata import hildebrand_sekhon
19 19 from schainpy.utils import log
20 20
21 21
22 22 class SpectraProc(ProcessingUnit):
23 23
24 24 def __init__(self):
25 25
26 26 ProcessingUnit.__init__(self)
27 27
28 28 self.buffer = None
29 29 self.firstdatatime = None
30 30 self.profIndex = 0
31 31 self.dataOut = Spectra()
32 32 self.id_min = None
33 33 self.id_max = None
34 34 self.setupReq = False #Agregar a todas las unidades de proc
35 35
36 36 def __updateSpecFromVoltage(self):
37 37
38 38 self.dataOut.timeZone = self.dataIn.timeZone
39 39 self.dataOut.dstFlag = self.dataIn.dstFlag
40 40 self.dataOut.errorCount = self.dataIn.errorCount
41 41 self.dataOut.useLocalTime = self.dataIn.useLocalTime
42 42 try:
43 43 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
44 44 except:
45 45 pass
46 46 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
47 47 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
48 48 self.dataOut.channelList = self.dataIn.channelList
49 49 self.dataOut.heightList = self.dataIn.heightList
50 50 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
51 51 self.dataOut.nProfiles = self.dataOut.nFFTPoints
52 52 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
53 53 self.dataOut.utctime = self.firstdatatime
54 54 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
55 55 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
56 56 self.dataOut.flagShiftFFT = False
57 57 self.dataOut.nCohInt = self.dataIn.nCohInt
58 58 self.dataOut.nIncohInt = 1
59 59 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
60 60 self.dataOut.frequency = self.dataIn.frequency
61 61 self.dataOut.realtime = self.dataIn.realtime
62 62 self.dataOut.azimuth = self.dataIn.azimuth
63 63 self.dataOut.zenith = self.dataIn.zenith
64 64 self.dataOut.beam.codeList = self.dataIn.beam.codeList
65 65 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
66 66 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
67 67 self.dataOut.runNextUnit = self.dataIn.runNextUnit
68 68 try:
69 69 self.dataOut.step = self.dataIn.step
70 70 except:
71 71 pass
72 72
73 73 def __getFft(self):
74 74 """
75 75 Convierte valores de Voltaje a Spectra
76 76
77 77 Affected:
78 78 self.dataOut.data_spc
79 79 self.dataOut.data_cspc
80 80 self.dataOut.data_dc
81 81 self.dataOut.heightList
82 82 self.profIndex
83 83 self.buffer
84 84 self.dataOut.flagNoData
85 85 """
86 86 fft_volt = numpy.fft.fft(
87 87 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
88 88 fft_volt = fft_volt.astype(numpy.dtype('complex'))
89 89 dc = fft_volt[:, 0, :]
90 90
91 91 # calculo de self-spectra
92 92 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
93 93 spc = fft_volt * numpy.conjugate(fft_volt)
94 94 spc = spc.real
95 95
96 96 blocksize = 0
97 97 blocksize += dc.size
98 98 blocksize += spc.size
99 99
100 100 cspc = None
101 101 pairIndex = 0
102 102 if self.dataOut.pairsList != None:
103 103 # calculo de cross-spectra
104 104 cspc = numpy.zeros(
105 105 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
106 106 for pair in self.dataOut.pairsList:
107 107 if pair[0] not in self.dataOut.channelList:
108 108 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
109 109 str(pair), str(self.dataOut.channelList)))
110 110 if pair[1] not in self.dataOut.channelList:
111 111 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
112 112 str(pair), str(self.dataOut.channelList)))
113 113
114 114 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
115 115 numpy.conjugate(fft_volt[pair[1], :, :])
116 116 pairIndex += 1
117 117 blocksize += cspc.size
118 118
119 119 self.dataOut.data_spc = spc
120 120 self.dataOut.data_cspc = cspc
121 121 self.dataOut.data_dc = dc
122 122 self.dataOut.blockSize = blocksize
123 123 self.dataOut.flagShiftFFT = False
124 124
125 125 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
126 126
127 127 self.dataIn.runNextUnit = runNextUnit
128 128 if self.dataIn.type == "Spectra":
129 129
130 130 self.dataOut.copy(self.dataIn)
131 131 if shift_fft:
132 132 #desplaza a la derecha en el eje 2 determinadas posiciones
133 133 shift = int(self.dataOut.nFFTPoints/2)
134 134 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
135 135
136 136 if self.dataOut.data_cspc is not None:
137 137 #desplaza a la derecha en el eje 2 determinadas posiciones
138 138 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
139 139 if pairsList:
140 140 self.__selectPairs(pairsList)
141 141
142 142 elif self.dataIn.type == "Voltage":
143 143
144 144 self.dataOut.flagNoData = True
145 145
146 146 if nFFTPoints == None:
147 147 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
148 148
149 149 if nProfiles == None:
150 150 nProfiles = nFFTPoints
151 151 #print(self.dataOut.ipp)
152 152 #exit(1)
153 153 if ippFactor == None:
154 154 self.dataOut.ippFactor = 1
155 155 #if ippFactor is not None:
156 156 #self.dataOut.ippFactor = ippFactor
157 157 #print(ippFactor)
158 158 #print(self.dataOut.ippFactor)
159 159 #exit(1)
160 160
161 161 self.dataOut.nFFTPoints = nFFTPoints
162 162
163 163 if self.buffer is None:
164 164 self.buffer = numpy.zeros((self.dataIn.nChannels,
165 165 nProfiles,
166 166 self.dataIn.nHeights),
167 167 dtype='complex')
168 168
169 169 if self.dataIn.flagDataAsBlock:
170 170 nVoltProfiles = self.dataIn.data.shape[1]
171 171
172 172 if nVoltProfiles == nProfiles:
173 173 self.buffer = self.dataIn.data.copy()
174 174 self.profIndex = nVoltProfiles
175 175
176 176 elif nVoltProfiles < nProfiles:
177 177
178 178 if self.profIndex == 0:
179 179 self.id_min = 0
180 180 self.id_max = nVoltProfiles
181 181 #print(self.id_min)
182 182 #print(self.id_max)
183 183 #print(numpy.shape(self.buffer))
184 184 self.buffer[:, self.id_min:self.id_max,
185 185 :] = self.dataIn.data
186 186 self.profIndex += nVoltProfiles
187 187 self.id_min += nVoltProfiles
188 188 self.id_max += nVoltProfiles
189 189 else:
190 190 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
191 191 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
192 192 self.dataOut.flagNoData = True
193 193 else:
194 194 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
195 195 self.profIndex += 1
196 196
197 197 if self.firstdatatime == None:
198 198 self.firstdatatime = self.dataIn.utctime
199 199
200 200 if self.profIndex == nProfiles:
201 201 self.__updateSpecFromVoltage()
202 202 if pairsList == None:
203 203 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
204 204 else:
205 205 self.dataOut.pairsList = pairsList
206 206 self.__getFft()
207 207 self.dataOut.flagNoData = False
208 208 self.firstdatatime = None
209 209 self.profIndex = 0
210 210 else:
211 211 raise ValueError("The type of input object '%s' is not valid".format(
212 212 self.dataIn.type))
213 213
214 214
215 215 def __selectPairs(self, pairsList):
216 216
217 217 if not pairsList:
218 218 return
219 219
220 220 pairs = []
221 221 pairsIndex = []
222 222
223 223 for pair in pairsList:
224 224 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
225 225 continue
226 226 pairs.append(pair)
227 227 pairsIndex.append(pairs.index(pair))
228 228
229 229 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
230 230 self.dataOut.pairsList = pairs
231 231
232 232 return
233 233
234 234 def selectFFTs(self, minFFT, maxFFT ):
235 235 """
236 236 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
237 237 minFFT<= FFT <= maxFFT
238 238 """
239 239
240 240 if (minFFT > maxFFT):
241 241 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
242 242
243 243 if (minFFT < self.dataOut.getFreqRange()[0]):
244 244 minFFT = self.dataOut.getFreqRange()[0]
245 245
246 246 if (maxFFT > self.dataOut.getFreqRange()[-1]):
247 247 maxFFT = self.dataOut.getFreqRange()[-1]
248 248
249 249 minIndex = 0
250 250 maxIndex = 0
251 251 FFTs = self.dataOut.getFreqRange()
252 252
253 253 inda = numpy.where(FFTs >= minFFT)
254 254 indb = numpy.where(FFTs <= maxFFT)
255 255
256 256 try:
257 257 minIndex = inda[0][0]
258 258 except:
259 259 minIndex = 0
260 260
261 261 try:
262 262 maxIndex = indb[0][-1]
263 263 except:
264 264 maxIndex = len(FFTs)
265 265
266 266 self.selectFFTsByIndex(minIndex, maxIndex)
267 267
268 268 return 1
269 269
270 270 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
271 271 newheis = numpy.where(
272 272 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
273 273
274 274 if hei_ref != None:
275 275 newheis = numpy.where(self.dataOut.heightList > hei_ref)
276 276
277 277 minIndex = min(newheis[0])
278 278 maxIndex = max(newheis[0])
279 279 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
280 280 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
281 281
282 282 # determina indices
283 283 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
284 284 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
285 285 avg_dB = 10 * \
286 286 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
287 287 beacon_dB = numpy.sort(avg_dB)[-nheis:]
288 288 beacon_heiIndexList = []
289 289 for val in avg_dB.tolist():
290 290 if val >= beacon_dB[0]:
291 291 beacon_heiIndexList.append(avg_dB.tolist().index(val))
292 292
293 293 #data_spc = data_spc[:,:,beacon_heiIndexList]
294 294 data_cspc = None
295 295 if self.dataOut.data_cspc is not None:
296 296 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
297 297 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
298 298
299 299 data_dc = None
300 300 if self.dataOut.data_dc is not None:
301 301 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
302 302 #data_dc = data_dc[:,beacon_heiIndexList]
303 303
304 304 self.dataOut.data_spc = data_spc
305 305 self.dataOut.data_cspc = data_cspc
306 306 self.dataOut.data_dc = data_dc
307 307 self.dataOut.heightList = heightList
308 308 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
309 309
310 310 return 1
311 311
312 312 def selectFFTsByIndex(self, minIndex, maxIndex):
313 313 """
314 314
315 315 """
316 316
317 317 if (minIndex < 0) or (minIndex > maxIndex):
318 318 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
319 319
320 320 if (maxIndex >= self.dataOut.nProfiles):
321 321 maxIndex = self.dataOut.nProfiles-1
322 322
323 323 #Spectra
324 324 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
325 325
326 326 data_cspc = None
327 327 if self.dataOut.data_cspc is not None:
328 328 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
329 329
330 330 data_dc = None
331 331 if self.dataOut.data_dc is not None:
332 332 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
333 333
334 334 self.dataOut.data_spc = data_spc
335 335 self.dataOut.data_cspc = data_cspc
336 336 self.dataOut.data_dc = data_dc
337 337
338 338 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
339 339 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
340 340 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
341 341
342 342 return 1
343 343
344 344 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
345 345 # validacion de rango
346 346 print("NOISeeee")
347 347 if minHei == None:
348 348 minHei = self.dataOut.heightList[0]
349 349
350 350 if maxHei == None:
351 351 maxHei = self.dataOut.heightList[-1]
352 352
353 353 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
354 354 print('minHei: %.2f is out of the heights range' % (minHei))
355 355 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
356 356 minHei = self.dataOut.heightList[0]
357 357
358 358 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
359 359 print('maxHei: %.2f is out of the heights range' % (maxHei))
360 360 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
361 361 maxHei = self.dataOut.heightList[-1]
362 362
363 363 # validacion de velocidades
364 364 velrange = self.dataOut.getVelRange(1)
365 365
366 366 if minVel == None:
367 367 minVel = velrange[0]
368 368
369 369 if maxVel == None:
370 370 maxVel = velrange[-1]
371 371
372 372 if (minVel < velrange[0]) or (minVel > maxVel):
373 373 print('minVel: %.2f is out of the velocity range' % (minVel))
374 374 print('minVel is setting to %.2f' % (velrange[0]))
375 375 minVel = velrange[0]
376 376
377 377 if (maxVel > velrange[-1]) or (maxVel < minVel):
378 378 print('maxVel: %.2f is out of the velocity range' % (maxVel))
379 379 print('maxVel is setting to %.2f' % (velrange[-1]))
380 380 maxVel = velrange[-1]
381 381
382 382 # seleccion de indices para rango
383 383 minIndex = 0
384 384 maxIndex = 0
385 385 heights = self.dataOut.heightList
386 386
387 387 inda = numpy.where(heights >= minHei)
388 388 indb = numpy.where(heights <= maxHei)
389 389
390 390 try:
391 391 minIndex = inda[0][0]
392 392 except:
393 393 minIndex = 0
394 394
395 395 try:
396 396 maxIndex = indb[0][-1]
397 397 except:
398 398 maxIndex = len(heights)
399 399
400 400 if (minIndex < 0) or (minIndex > maxIndex):
401 401 raise ValueError("some value in (%d,%d) is not valid" % (
402 402 minIndex, maxIndex))
403 403
404 404 if (maxIndex >= self.dataOut.nHeights):
405 405 maxIndex = self.dataOut.nHeights - 1
406 406
407 407 # seleccion de indices para velocidades
408 408 indminvel = numpy.where(velrange >= minVel)
409 409 indmaxvel = numpy.where(velrange <= maxVel)
410 410 try:
411 411 minIndexVel = indminvel[0][0]
412 412 except:
413 413 minIndexVel = 0
414 414
415 415 try:
416 416 maxIndexVel = indmaxvel[0][-1]
417 417 except:
418 418 maxIndexVel = len(velrange)
419 419
420 420 # seleccion del espectro
421 421 data_spc = self.dataOut.data_spc[:,
422 422 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
423 423 # estimacion de ruido
424 424 noise = numpy.zeros(self.dataOut.nChannels)
425 425
426 426 for channel in range(self.dataOut.nChannels):
427 427 daux = data_spc[channel, :, :]
428 428 sortdata = numpy.sort(daux, axis=None)
429 429 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
430 430
431 431 self.dataOut.noise_estimation = noise.copy()
432 432
433 433 return 1
434 434
435 435 class GetSNR(Operation):
436 436 '''
437 437 Written by R. Flores
438 438 '''
439 439 """Operation to get SNR.
440 440
441 441 Parameters:
442 442 -----------
443 443
444 444 Example
445 445 --------
446 446
447 447 op = proc_unit.addOperation(name='GetSNR', optype='other')
448 448
449 449 """
450 450
451 451 def __init__(self, **kwargs):
452 452
453 453 Operation.__init__(self, **kwargs)
454 454
455 455
456 456 def run(self,dataOut):
457 457
458 458 #noise = dataOut.getNoise()
459 459 noise = dataOut.getNoise(ymin_index=-10) #Región superior donde solo debería de haber ruido
460 460 #print("Noise: ", noise)
461 461 #print("Noise_dB: ", 10*numpy.log10(noise/dataOut.normFactor))
462 462 #print("Heights: ", dataOut.heightList)
463 463 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.normFactor)
464 464 ################dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints) #Before 12Jan2023
465 465 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None])/(noise[:,None])
466 466 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
467 467 dataOut.snl = numpy.log10(dataOut.data_snr)
468 468 #print("snl: ", dataOut.snl)
469 469 #exit(1)
470 470 #print(dataOut.heightList[-11])
471 471 #print(numpy.shape(dataOut.heightList))
472 472 #print(dataOut.data_snr)
473 473 #print(dataOut.data_snr[0,-11])
474 474 #exit(1)
475 475 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr)
476 476 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.1, numpy.nan, dataOut.data_snr)
477 477 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.0, numpy.nan, dataOut.data_snr)
478 478 #dataOut.data_snr = numpy.where(dataOut.data_snr<.05, numpy.nan, dataOut.data_snr)
479 479 #dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
480 480 dataOut.snl = numpy.where(dataOut.snl<-1, numpy.nan, dataOut.snl)
481 481 '''
482 482 import matplotlib.pyplot as plt
483 483 #plt.plot(10*numpy.log10(dataOut.data_snr[0]),dataOut.heightList)
484 484 plt.plot(dataOut.data_snr[0],dataOut.heightList)#,marker='*')
485 485 plt.xlim(-1,10)
486 486 plt.axvline(1,color='k')
487 487 plt.axvline(.1,color='k',linestyle='--')
488 488 plt.grid()
489 489 plt.show()
490 490 '''
491 491 #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr)
492 492 #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0)
493 493 #print(dataOut.data_snr.shape)
494 494 #exit(1)
495 495 #print("Before: ", dataOut.data_snr[0])
496 496
497 497
498 498 return dataOut
499 499
500 500 class removeDC(Operation):
501 501
502 502 def run(self, dataOut, mode=2):
503 503 self.dataOut = dataOut
504 504 jspectra = self.dataOut.data_spc
505 505 jcspectra = self.dataOut.data_cspc
506 506
507 507 num_chan = jspectra.shape[0]
508 508 num_hei = jspectra.shape[2]
509 509
510 510 if jcspectra is not None:
511 511 jcspectraExist = True
512 512 num_pairs = jcspectra.shape[0]
513 513 else:
514 514 jcspectraExist = False
515 515
516 516 freq_dc = int(jspectra.shape[1] / 2)
517 517 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
518 518 ind_vel = ind_vel.astype(int)
519 519
520 520 if ind_vel[0] < 0:
521 521 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
522 522
523 523 if mode == 1:
524 524 jspectra[:, freq_dc, :] = (
525 525 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
526 526
527 527 if jcspectraExist:
528 528 jcspectra[:, freq_dc, :] = (
529 529 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
530 530
531 531 if mode == 2:
532 532
533 533 vel = numpy.array([-2, -1, 1, 2])
534 534 xx = numpy.zeros([4, 4])
535 535
536 536 for fil in range(4):
537 537 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
538 538
539 539 xx_inv = numpy.linalg.inv(xx)
540 540 xx_aux = xx_inv[0, :]
541 541
542 542 for ich in range(num_chan):
543 543 yy = jspectra[ich, ind_vel, :]
544 544 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
545 545
546 546 junkid = jspectra[ich, freq_dc, :] <= 0
547 547 cjunkid = sum(junkid)
548 548
549 549 if cjunkid.any():
550 550 jspectra[ich, freq_dc, junkid.nonzero()] = (
551 551 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
552 552
553 553 if jcspectraExist:
554 554 for ip in range(num_pairs):
555 555 yy = jcspectra[ip, ind_vel, :]
556 556 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
557 557
558 558 self.dataOut.data_spc = jspectra
559 559 self.dataOut.data_cspc = jcspectra
560 560
561 561 return self.dataOut
562 562
563 563 class removeInterference(Operation):
564 564
565 565 def removeInterference2(self):
566 566
567 567 cspc = self.dataOut.data_cspc
568 568 spc = self.dataOut.data_spc
569 569 Heights = numpy.arange(cspc.shape[2])
570 570 realCspc = numpy.abs(cspc)
571 571
572 572 for i in range(cspc.shape[0]):
573 573 LinePower= numpy.sum(realCspc[i], axis=0)
574 574 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
575 575 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
576 576 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
577 577 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
578 578 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
579 579
580 580
581 581 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
582 582 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
583 583 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
584 584 cspc[i,InterferenceRange,:] = numpy.NaN
585 585
586 586 self.dataOut.data_cspc = cspc
587 587
588 588 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
589 589
590 590 jspectra = self.dataOut.data_spc
591 591 jcspectra = self.dataOut.data_cspc
592 592 jnoise = self.dataOut.getNoise()
593 593 num_incoh = self.dataOut.nIncohInt
594 594
595 595 num_channel = jspectra.shape[0]
596 596 num_prof = jspectra.shape[1]
597 597 num_hei = jspectra.shape[2]
598 598
599 599 # hei_interf
600 600 if hei_interf is None:
601 601 count_hei = int(num_hei / 2)
602 602 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
603 603 hei_interf = numpy.asarray(hei_interf)[0]
604 604 # nhei_interf
605 605 if (nhei_interf == None):
606 606 nhei_interf = 5
607 607 if (nhei_interf < 1):
608 608 nhei_interf = 1
609 609 if (nhei_interf > count_hei):
610 610 nhei_interf = count_hei
611 611 if (offhei_interf == None):
612 612 offhei_interf = 0
613 613
614 614 ind_hei = list(range(num_hei))
615 615 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
616 616 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
617 617 mask_prof = numpy.asarray(list(range(num_prof)))
618 618 num_mask_prof = mask_prof.size
619 619 comp_mask_prof = [0, num_prof / 2]
620 620
621 621 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
622 622 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
623 623 jnoise = numpy.nan
624 624 noise_exist = jnoise[0] < numpy.Inf
625 625
626 626 # Subrutina de Remocion de la Interferencia
627 627 for ich in range(num_channel):
628 628 # Se ordena los espectros segun su potencia (menor a mayor)
629 629 power = jspectra[ich, mask_prof, :]
630 630 power = power[:, hei_interf]
631 631 power = power.sum(axis=0)
632 632 psort = power.ravel().argsort()
633 633
634 634 # Se estima la interferencia promedio en los Espectros de Potencia empleando
635 635 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
636 636 offhei_interf, nhei_interf + offhei_interf))]]]
637 637
638 638 if noise_exist:
639 639 # tmp_noise = jnoise[ich] / num_prof
640 640 tmp_noise = jnoise[ich]
641 641 junkspc_interf = junkspc_interf - tmp_noise
642 642 #junkspc_interf[:,comp_mask_prof] = 0
643 643
644 644 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
645 645 jspc_interf = jspc_interf.transpose()
646 646 # Calculando el espectro de interferencia promedio
647 647 noiseid = numpy.where(
648 648 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
649 649 noiseid = noiseid[0]
650 650 cnoiseid = noiseid.size
651 651 interfid = numpy.where(
652 652 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
653 653 interfid = interfid[0]
654 654 cinterfid = interfid.size
655 655
656 656 if (cnoiseid > 0):
657 657 jspc_interf[noiseid] = 0
658 658
659 659 # Expandiendo los perfiles a limpiar
660 660 if (cinterfid > 0):
661 661 new_interfid = (
662 662 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
663 663 new_interfid = numpy.asarray(new_interfid)
664 664 new_interfid = {x for x in new_interfid}
665 665 new_interfid = numpy.array(list(new_interfid))
666 666 new_cinterfid = new_interfid.size
667 667 else:
668 668 new_cinterfid = 0
669 669
670 670 for ip in range(new_cinterfid):
671 671 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
672 672 jspc_interf[new_interfid[ip]
673 673 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
674 674
675 675 jspectra[ich, :, ind_hei] = jspectra[ich, :,
676 676 ind_hei] - jspc_interf # Corregir indices
677 677
678 678 # Removiendo la interferencia del punto de mayor interferencia
679 679 ListAux = jspc_interf[mask_prof].tolist()
680 680 maxid = ListAux.index(max(ListAux))
681 681
682 682 if cinterfid > 0:
683 683 for ip in range(cinterfid * (interf == 2) - 1):
684 684 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
685 685 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
686 686 cind = len(ind)
687 687
688 688 if (cind > 0):
689 689 jspectra[ich, interfid[ip], ind] = tmp_noise * \
690 690 (1 + (numpy.random.uniform(cind) - 0.5) /
691 691 numpy.sqrt(num_incoh))
692 692
693 693 ind = numpy.array([-2, -1, 1, 2])
694 694 xx = numpy.zeros([4, 4])
695 695
696 696 for id1 in range(4):
697 697 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
698 698
699 699 xx_inv = numpy.linalg.inv(xx)
700 700 xx = xx_inv[:, 0]
701 701 ind = (ind + maxid + num_mask_prof) % num_mask_prof
702 702 yy = jspectra[ich, mask_prof[ind], :]
703 703 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
704 704 yy.transpose(), xx)
705 705
706 706 indAux = (jspectra[ich, :, :] < tmp_noise *
707 707 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
708 708 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
709 709 (1 - 1 / numpy.sqrt(num_incoh))
710 710
711 711 # Remocion de Interferencia en el Cross Spectra
712 712 if jcspectra is None:
713 713 return jspectra, jcspectra
714 714 num_pairs = int(jcspectra.size / (num_prof * num_hei))
715 715 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
716 716
717 717 for ip in range(num_pairs):
718 718
719 719 #-------------------------------------------
720 720
721 721 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
722 722 cspower = cspower[:, hei_interf]
723 723 cspower = cspower.sum(axis=0)
724 724
725 725 cspsort = cspower.ravel().argsort()
726 726 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
727 727 offhei_interf, nhei_interf + offhei_interf))]]]
728 728 junkcspc_interf = junkcspc_interf.transpose()
729 729 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
730 730
731 731 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
732 732
733 733 median_real = int(numpy.median(numpy.real(
734 734 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
735 735 median_imag = int(numpy.median(numpy.imag(
736 736 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
737 737 comp_mask_prof = [int(e) for e in comp_mask_prof]
738 738 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
739 739 median_real, median_imag)
740 740
741 741 for iprof in range(num_prof):
742 742 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
743 743 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
744 744
745 745 # Removiendo la Interferencia
746 746 jcspectra[ip, :, ind_hei] = jcspectra[ip,
747 747 :, ind_hei] - jcspc_interf
748 748
749 749 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
750 750 maxid = ListAux.index(max(ListAux))
751 751
752 752 ind = numpy.array([-2, -1, 1, 2])
753 753 xx = numpy.zeros([4, 4])
754 754
755 755 for id1 in range(4):
756 756 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
757 757
758 758 xx_inv = numpy.linalg.inv(xx)
759 759 xx = xx_inv[:, 0]
760 760
761 761 ind = (ind + maxid + num_mask_prof) % num_mask_prof
762 762 yy = jcspectra[ip, mask_prof[ind], :]
763 763 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
764 764
765 765 # Guardar Resultados
766 766 self.dataOut.data_spc = jspectra
767 767 self.dataOut.data_cspc = jcspectra
768 768
769 769 return 1
770 770
771 771 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
772 772
773 773 self.dataOut = dataOut
774 774
775 775 if mode == 1:
776 776 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
777 777 elif mode == 2:
778 778 self.removeInterference2()
779 779
780 780 return self.dataOut
781 781
782 class removeInterferenceAtFreq(Operation):
783 '''
784 Written by R. Flores
785 '''
786 """Operation to remove interfernce at a known frequency(s).
787
788 Parameters:
789 -----------
790 None
791
792 Example
793 --------
794
795 op = proc_unit.addOperation(name='removeInterferenceAtFreq')
796
797 """
798
799 def __init__(self):
800
801 Operation.__init__(self)
802
803 def run(self, dataOut, freq = None, freqList = None):
804
805 VelRange = dataOut.getVelRange()
806 #print("VelRange: ", VelRange)
807
808 freq_ids = []
809
810 if freq is not None:
811 #print("freq")
812 #if freq < 0:
813 inda = numpy.where(VelRange >= freq)
814 minIndex = inda[0][0]
815 #print(numpy.shape(dataOut.dataLag_spc))
816 dataOut.data_spc[:,minIndex,:] = numpy.nan
817
818 #inda = numpy.where(VelRange >= ymin_noise)
819 #indb = numpy.where(VelRange <= ymax_noise)
820
821 #minIndex = inda[0][0]
822 #maxIndex = indb[0][-1]
823
824 elif freqList is not None:
825 #print("freqList")
826 for freq in freqList:
827 #if freq < 0:
828 inda = numpy.where(VelRange >= freq)
829 minIndex = inda[0][0]
830 #print(numpy.shape(dataOut.dataLag_spc))
831 if freq > 0:
832 #dataOut.data_spc[:,minIndex-1,:] = numpy.nan
833 freq_ids.append(minIndex-1)
834 else:
835 #dataOut.data_spc[:,minIndex,:] = numpy.nan
836 freq_ids.append(minIndex)
837 else:
838 raise ValueError("freq or freqList should be specified ...")
839
840 #freq_ids = numpy.array(freq_ids).flatten()
841
842 avg = numpy.mean(dataOut.data_spc[:,[t for t in range(dataOut.data_spc.shape[0]) if t not in freq_ids],:],axis=1)
843
844 for p in list(freq_ids):
845 dataOut.data_spc[:,p,:] = avg#numpy.nan
846
847
848 return dataOut
782 849
783 850 class IncohInt(Operation):
784 851
785 852 __profIndex = 0
786 853 __withOverapping = False
787 854
788 855 __byTime = False
789 856 __initime = None
790 857 __lastdatatime = None
791 858 __integrationtime = None
792 859
793 860 __buffer_spc = None
794 861 __buffer_cspc = None
795 862 __buffer_dc = None
796 863
797 864 __dataReady = False
798 865
799 866 __timeInterval = None
800 867
801 868 n = None
802 869
803 870 def __init__(self):
804 871
805 872 Operation.__init__(self)
806 873
807 874 def setup(self, n=None, timeInterval=None, overlapping=False):
808 875 """
809 876 Set the parameters of the integration class.
810 877
811 878 Inputs:
812 879
813 880 n : Number of coherent integrations
814 881 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
815 882 overlapping :
816 883
817 884 """
818 885
819 886 self.__initime = None
820 887 self.__lastdatatime = 0
821 888
822 889 self.__buffer_spc = 0
823 890 self.__buffer_cspc = 0
824 891 self.__buffer_dc = 0
825 892
826 893 self.__profIndex = 0
827 894 self.__dataReady = False
828 895 self.__byTime = False
829 896
830 897 if n is None and timeInterval is None:
831 898 raise ValueError("n or timeInterval should be specified ...")
832 899
833 900 if n is not None:
834 901 self.n = int(n)
835 902 else:
836 903
837 904 self.__integrationtime = int(timeInterval)
838 905 self.n = None
839 906 self.__byTime = True
840 907
841 908 def putData(self, data_spc, data_cspc, data_dc):
842 909 """
843 910 Add a profile to the __buffer_spc and increase in one the __profileIndex
844 911
845 912 """
846 913
847 914 self.__buffer_spc += data_spc
848 915
849 916 if data_cspc is None:
850 917 self.__buffer_cspc = None
851 918 else:
852 919 self.__buffer_cspc += data_cspc
853 920
854 921 if data_dc is None:
855 922 self.__buffer_dc = None
856 923 else:
857 924 self.__buffer_dc += data_dc
858 925
859 926 self.__profIndex += 1
860 927
861 928 return
862 929
863 930 def pushData(self):
864 931 """
865 932 Return the sum of the last profiles and the profiles used in the sum.
866 933
867 934 Affected:
868 935
869 936 self.__profileIndex
870 937
871 938 """
872 939
873 940 data_spc = self.__buffer_spc
874 941 data_cspc = self.__buffer_cspc
875 942 data_dc = self.__buffer_dc
876 943 n = self.__profIndex
877 944
878 945 self.__buffer_spc = 0
879 946 self.__buffer_cspc = 0
880 947 self.__buffer_dc = 0
881 948 self.__profIndex = 0
882 949
883 950 return data_spc, data_cspc, data_dc, n
884 951
885 952 def byProfiles(self, *args):
886 953
887 954 self.__dataReady = False
888 955 avgdata_spc = None
889 956 avgdata_cspc = None
890 957 avgdata_dc = None
891 958
892 959 self.putData(*args)
893 960
894 961 if self.__profIndex == self.n:
895 962
896 963 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
897 964 self.n = n
898 965 self.__dataReady = True
899 966
900 967 return avgdata_spc, avgdata_cspc, avgdata_dc
901 968
902 969 def byTime(self, datatime, *args):
903 970
904 971 self.__dataReady = False
905 972 avgdata_spc = None
906 973 avgdata_cspc = None
907 974 avgdata_dc = None
908 975
909 976 self.putData(*args)
910 977
911 978 if (datatime - self.__initime) >= self.__integrationtime:
912 979 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
913 980 self.n = n
914 981 self.__dataReady = True
915 982
916 983 return avgdata_spc, avgdata_cspc, avgdata_dc
917 984
918 985 def integrate(self, datatime, *args):
919 986
920 987 if self.__profIndex == 0:
921 988 self.__initime = datatime
922 989
923 990 if self.__byTime:
924 991 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
925 992 datatime, *args)
926 993 else:
927 994 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
928 995
929 996 if not self.__dataReady:
930 997 return None, None, None, None
931 998
932 999 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
933 1000
934 1001 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
935 1002 if n == 1:
936 1003 return dataOut
937 1004 print("JERE")
938 1005 dataOut.flagNoData = True
939 1006
940 1007 if not self.isConfig:
941 1008 self.setup(n, timeInterval, overlapping)
942 1009 self.isConfig = True
943 1010
944 1011 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
945 1012 dataOut.data_spc,
946 1013 dataOut.data_cspc,
947 1014 dataOut.data_dc)
948 1015
949 1016 if self.__dataReady:
950 1017
951 1018 dataOut.data_spc = avgdata_spc
952 1019 print(numpy.sum(dataOut.data_spc))
953 1020 exit(1)
954 1021 dataOut.data_cspc = avgdata_cspc
955 1022 dataOut.data_dc = avgdata_dc
956 1023 dataOut.nIncohInt *= self.n
957 1024 dataOut.utctime = avgdatatime
958 1025 dataOut.flagNoData = False
959 1026
960 1027 return dataOut
961 1028
962 1029 class dopplerFlip(Operation):
963 1030
964 1031 def run(self, dataOut, chann = None):
965 1032 # arreglo 1: (num_chan, num_profiles, num_heights)
966 1033 self.dataOut = dataOut
967 1034 # JULIA-oblicua, indice 2
968 1035 # arreglo 2: (num_profiles, num_heights)
969 1036 jspectra = self.dataOut.data_spc[chann]
970 1037 jspectra_tmp = numpy.zeros(jspectra.shape)
971 1038 num_profiles = jspectra.shape[0]
972 1039 freq_dc = int(num_profiles / 2)
973 1040 # Flip con for
974 1041 for j in range(num_profiles):
975 1042 jspectra_tmp[num_profiles-j-1]= jspectra[j]
976 1043 # Intercambio perfil de DC con perfil inmediato anterior
977 1044 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
978 1045 jspectra_tmp[freq_dc]= jspectra[freq_dc]
979 1046 # canal modificado es re-escrito en el arreglo de canales
980 1047 self.dataOut.data_spc[chann] = jspectra_tmp
981 1048
982 1049 return self.dataOut
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
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now