##// END OF EJS Templates
schain
RSS
cambios para amisr ISR getNoise
joabAM -
r1468:1fb12846b034
parent child
Show More
Side by Side Unified
Context file:
r1468:1fb12846b034 Loading diff...:
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
@@ -1,961 +1,962
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13 from itertools import combinations
14 14
15 15
16 16 class SpectraPlot(Plot):
17 17 '''
18 18 Plot for Spectra data
19 19 '''
20 20
21 21 CODE = 'spc'
22 22 colormap = 'jet'
23 23 plot_type = 'pcolor'
24 24 buffering = False
25 25 channelList = []
26 26
27 27 def setup(self):
28 28
29 29 self.nplots = len(self.data.channels)
30 30 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
31 31 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
32 32 self.height = 2.6 * self.nrows
33 33
34 34 self.cb_label = 'dB'
35 35 if self.showprofile:
36 36 self.width = 4 * self.ncols
37 37 else:
38 38 self.width = 3.5 * self.ncols
39 39 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
40 40 self.ylabel = 'Range [km]'
41 41
42 42
43 43 def update_list(self,dataOut):
44 44 if len(self.channelList) == 0:
45 45 self.channelList = dataOut.channelList
46 46
47 47 def update(self, dataOut):
48 48
49 49 self.update_list(dataOut)
50 50 data = {}
51 51 meta = {}
52 52 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
53 53 data['spc'] = spc
54 54 data['rti'] = dataOut.getPower()
55 55 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
56 56 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
57 57 if self.CODE == 'spc_moments':
58 58 data['moments'] = dataOut.moments
59 59
60 60 return data, meta
61 61
62 62 def plot(self):
63 63 if self.xaxis == "frequency":
64 64 x = self.data.xrange[0]
65 65 self.xlabel = "Frequency (kHz)"
66 66 elif self.xaxis == "time":
67 67 x = self.data.xrange[1]
68 68 self.xlabel = "Time (ms)"
69 69 else:
70 70 x = self.data.xrange[2]
71 71 self.xlabel = "Velocity (m/s)"
72 72
73 73 if self.CODE == 'spc_moments':
74 74 x = self.data.xrange[2]
75 75 self.xlabel = "Velocity (m/s)"
76 76
77 77 self.titles = []
78 78 y = self.data.yrange
79 79 self.y = y
80 80
81 81 data = self.data[-1]
82 82 z = data['spc']
83 83
84 84 for n, ax in enumerate(self.axes):
85 85 noise = data['noise'][n]
86 86 if self.CODE == 'spc_moments':
87 87 mean = data['moments'][n, 1]
88 88 if ax.firsttime:
89 89 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
90 90 self.xmin = self.xmin if self.xmin else -self.xmax
91 91 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
92 92 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
93 93 ax.plt = ax.pcolormesh(x, y, z[n].T,
94 94 vmin=self.zmin,
95 95 vmax=self.zmax,
96 96 cmap=plt.get_cmap(self.colormap)
97 97 )
98 98
99 99 if self.showprofile:
100 100 ax.plt_profile = self.pf_axes[n].plot(
101 101 data['rti'][n], y)[0]
102 102 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
103 103 color="k", linestyle="dashed", lw=1)[0]
104 104 if self.CODE == 'spc_moments':
105 105 ax.plt_mean = ax.plot(mean, y, color='k')[0]
106 106 else:
107 107 ax.plt.set_array(z[n].T.ravel())
108 108 if self.showprofile:
109 109 ax.plt_profile.set_data(data['rti'][n], y)
110 110 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
111 111 if self.CODE == 'spc_moments':
112 112 ax.plt_mean.set_data(mean, y)
113 113 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
114 114
115 115
116 116 class CrossSpectraPlot(Plot):
117 117
118 118 CODE = 'cspc'
119 119 colormap = 'jet'
120 120 plot_type = 'pcolor'
121 121 zmin_coh = None
122 122 zmax_coh = None
123 123 zmin_phase = None
124 124 zmax_phase = None
125 125 realChannels = None
126 126 crossPairs = None
127 127
128 128 def setup(self):
129 129
130 130 self.ncols = 4
131 131 self.nplots = len(self.data.pairs) * 2
132 132 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
133 133 self.width = 3.1 * self.ncols
134 134 self.height = 2.6 * self.nrows
135 135 self.ylabel = 'Range [km]'
136 136 self.showprofile = False
137 137 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
138 138
139 139 def update(self, dataOut):
140 140
141 141 data = {}
142 142 meta = {}
143 143
144 144 spc = dataOut.data_spc
145 145 cspc = dataOut.data_cspc
146 146 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
147 147 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
148 148 meta['pairs'] = rawPairs
149 149
150 150 if self.crossPairs == None:
151 151 self.crossPairs = dataOut.pairsList
152 152
153 153 tmp = []
154 154
155 155 for n, pair in enumerate(meta['pairs']):
156 156
157 157 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
158 158 coh = numpy.abs(out)
159 159 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
160 160 tmp.append(coh)
161 161 tmp.append(phase)
162 162
163 163 data['cspc'] = numpy.array(tmp)
164 164
165 165 return data, meta
166 166
167 167 def plot(self):
168 168
169 169 if self.xaxis == "frequency":
170 170 x = self.data.xrange[0]
171 171 self.xlabel = "Frequency (kHz)"
172 172 elif self.xaxis == "time":
173 173 x = self.data.xrange[1]
174 174 self.xlabel = "Time (ms)"
175 175 else:
176 176 x = self.data.xrange[2]
177 177 self.xlabel = "Velocity (m/s)"
178 178
179 179 self.titles = []
180 180
181 181 y = self.data.yrange
182 182 self.y = y
183 183
184 184 data = self.data[-1]
185 185 cspc = data['cspc']
186 186
187 187 for n in range(len(self.data.pairs)):
188 188
189 189 pair = self.crossPairs[n]
190 190
191 191 coh = cspc[n*2]
192 192 phase = cspc[n*2+1]
193 193 ax = self.axes[2 * n]
194 194
195 195 if ax.firsttime:
196 196 ax.plt = ax.pcolormesh(x, y, coh.T,
197 197 vmin=self.zmin_coh,
198 198 vmax=self.zmax_coh,
199 199 cmap=plt.get_cmap(self.colormap_coh)
200 200 )
201 201 else:
202 202 ax.plt.set_array(coh.T.ravel())
203 203 self.titles.append(
204 204 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
205 205
206 206 ax = self.axes[2 * n + 1]
207 207 if ax.firsttime:
208 208 ax.plt = ax.pcolormesh(x, y, phase.T,
209 209 vmin=-180,
210 210 vmax=180,
211 211 cmap=plt.get_cmap(self.colormap_phase)
212 212 )
213 213 else:
214 214 ax.plt.set_array(phase.T.ravel())
215 215
216 216 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
217 217
218 218
219 219 class RTIPlot(Plot):
220 220 '''
221 221 Plot for RTI data
222 222 '''
223 223
224 224 CODE = 'rti'
225 225 colormap = 'jet'
226 226 plot_type = 'pcolorbuffer'
227 227 titles = None
228 228 channelList = []
229 229
230 230 def setup(self):
231 231 self.xaxis = 'time'
232 232 self.ncols = 1
233 233 #print("dataChannels ",self.data.channels)
234 234 self.nrows = len(self.data.channels)
235 235 self.nplots = len(self.data.channels)
236 236 self.ylabel = 'Range [km]'
237 237 self.xlabel = 'Time'
238 238 self.cb_label = 'dB'
239 239 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
240 240 self.titles = ['{} Channel {}'.format(
241 241 self.CODE.upper(), x) for x in range(self.nplots)]
242 242
243 243 def update_list(self,dataOut):
244 244
245 245 self.channelList = dataOut.channelList
246 246
247 247
248 248 def update(self, dataOut):
249 249 if len(self.channelList) == 0:
250 250 self.update_list(dataOut)
251 251 data = {}
252 252 meta = {}
253 253 data['rti'] = dataOut.getPower()
254 254 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
255 255 return data, meta
256 256
257 257 def plot(self):
258 258
259 259 self.x = self.data.times
260 260 self.y = self.data.yrange
261 261 self.z = self.data[self.CODE]
262 262 self.z = numpy.array(self.z, dtype=float)
263 263 self.z = numpy.ma.masked_invalid(self.z)
264 264
265 265 try:
266 266 if self.channelList != None:
267 267 self.titles = ['{} Channel {}'.format(
268 268 self.CODE.upper(), x) for x in self.channelList]
269 269 except:
270 270 if self.channelList.any() != None:
271 271 self.titles = ['{} Channel {}'.format(
272 272 self.CODE.upper(), x) for x in self.channelList]
273 273 if self.decimation is None:
274 274 x, y, z = self.fill_gaps(self.x, self.y, self.z)
275 275 else:
276 276 x, y, z = self.fill_gaps(*self.decimate())
277 277 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
278 278 for n, ax in enumerate(self.axes):
279 279 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
280 280 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
281 281 data = self.data[-1]
282 282 if ax.firsttime:
283 283 ax.plt = ax.pcolormesh(x, y, z[n].T,
284 284 vmin=self.zmin,
285 285 vmax=self.zmax,
286 286 cmap=plt.get_cmap(self.colormap)
287 287 )
288 288 if self.showprofile:
289 289 ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0]
290 290
291 291 if "noise" in self.data:
292 292 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
293 293 color="k", linestyle="dashed", lw=1)[0]
294 294 else:
295 295 ax.collections.remove(ax.collections[0])
296 296 ax.plt = ax.pcolormesh(x, y, z[n].T,
297 297 vmin=self.zmin,
298 298 vmax=self.zmax,
299 299 cmap=plt.get_cmap(self.colormap)
300 300 )
301 301 if self.showprofile:
302 302 ax.plot_profile.set_data(data[self.CODE][n], self.y)
303 303 if "noise" in self.data:
304 304 ax.plot_noise.set_data(numpy.repeat(
305 305 data['noise'][n], len(self.y)), self.y)
306 306
307 307
308 308 class CoherencePlot(RTIPlot):
309 309 '''
310 310 Plot for Coherence data
311 311 '''
312 312
313 313 CODE = 'coh'
314 314
315 315 def setup(self):
316 316 self.xaxis = 'time'
317 317 self.ncols = 1
318 318 self.nrows = len(self.data.pairs)
319 319 self.nplots = len(self.data.pairs)
320 320 self.ylabel = 'Range [km]'
321 321 self.xlabel = 'Time'
322 322 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
323 323 if self.CODE == 'coh':
324 324 self.cb_label = ''
325 325 self.titles = [
326 326 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
327 327 else:
328 328 self.cb_label = 'Degrees'
329 329 self.titles = [
330 330 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
331 331
332 332 def update(self, dataOut):
333 333 self.update_list(dataOut)
334 334 data = {}
335 335 meta = {}
336 336 data['coh'] = dataOut.getCoherence()
337 337 meta['pairs'] = dataOut.pairsList
338 338
339 339
340 340 return data, meta
341 341
342 342 class PhasePlot(CoherencePlot):
343 343 '''
344 344 Plot for Phase map data
345 345 '''
346 346
347 347 CODE = 'phase'
348 348 colormap = 'seismic'
349 349
350 350 def update(self, dataOut):
351 351
352 352 data = {}
353 353 meta = {}
354 354 data['phase'] = dataOut.getCoherence(phase=True)
355 355 meta['pairs'] = dataOut.pairsList
356 356
357 357 return data, meta
358 358
359 359 class NoisePlot(Plot):
360 360 '''
361 361 Plot for noise
362 362 '''
363 363
364 364 CODE = 'noise'
365 365 plot_type = 'scatterbuffer'
366 366
367 367 def setup(self):
368 368 self.xaxis = 'time'
369 369 self.ncols = 1
370 370 self.nrows = 1
371 371 self.nplots = 1
372 372 self.ylabel = 'Intensity [dB]'
373 373 self.xlabel = 'Time'
374 374 self.titles = ['Noise']
375 375 self.colorbar = False
376 376 self.plots_adjust.update({'right': 0.85 })
377 377
378 378 def update(self, dataOut):
379 379
380 380 data = {}
381 381 meta = {}
382 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
382 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
383 data['noise'] = noise
383 384 meta['yrange'] = numpy.array([])
384 385
385 386 return data, meta
386 387
387 388 def plot(self):
388 389
389 390 x = self.data.times
390 391 xmin = self.data.min_time
391 392 xmax = xmin + self.xrange * 60 * 60
392 393 Y = self.data['noise']
393 394
394 395 if self.axes[0].firsttime:
395 self.ymin = numpy.nanmin(Y) - 5
396 self.ymax = numpy.nanmax(Y) + 5
396 if self.ymin is None: self.ymin = numpy.nanmin(Y) - 5
397 if self.ymax is None: self.ymax = numpy.nanmax(Y) + 5
397 398 for ch in self.data.channels:
398 399 y = Y[ch]
399 400 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
400 401 plt.legend(bbox_to_anchor=(1.18, 1.0))
401 402 else:
402 403 for ch in self.data.channels:
403 404 y = Y[ch]
404 405 self.axes[0].lines[ch].set_data(x, y)
405 406
406 407
407 408 class PowerProfilePlot(Plot):
408 409
409 410 CODE = 'pow_profile'
410 411 plot_type = 'scatter'
411 412
412 413 def setup(self):
413 414
414 415 self.ncols = 1
415 416 self.nrows = 1
416 417 self.nplots = 1
417 418 self.height = 4
418 419 self.width = 3
419 420 self.ylabel = 'Range [km]'
420 421 self.xlabel = 'Intensity [dB]'
421 422 self.titles = ['Power Profile']
422 423 self.colorbar = False
423 424
424 425 def update(self, dataOut):
425 426
426 427 data = {}
427 428 meta = {}
428 429 data[self.CODE] = dataOut.getPower()
429 430
430 431 return data, meta
431 432
432 433 def plot(self):
433 434
434 435 y = self.data.yrange
435 436 self.y = y
436 437
437 438 x = self.data[-1][self.CODE]
438 439
439 440 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
440 441 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
441 442
442 443 if self.axes[0].firsttime:
443 444 for ch in self.data.channels:
444 445 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
445 446 plt.legend()
446 447 else:
447 448 for ch in self.data.channels:
448 449 self.axes[0].lines[ch].set_data(x[ch], y)
449 450
450 451
451 452 class SpectraCutPlot(Plot):
452 453
453 454 CODE = 'spc_cut'
454 455 plot_type = 'scatter'
455 456 buffering = False
456 457 heights = []
457 458 channelList = []
458 459 maintitle = "Spectra Cuts"
459 460
460 461 def setup(self):
461 462
462 463 self.nplots = len(self.data.channels)
463 464 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
464 465 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
465 466 self.width = 3.6 * self.ncols + 1.5
466 467 self.height = 3 * self.nrows
467 468 self.ylabel = 'Power [dB]'
468 469 self.colorbar = False
469 470 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
470 471 if self.selectedHeight:
471 472 self.maintitle = "Spectra Cut for %d km " %(int(self.selectedHeight))
472 473
473 474 def update(self, dataOut):
474 475 if len(self.channelList) == 0:
475 476 self.channelList = dataOut.channelList
476 477
477 478 self.heights = dataOut.heightList
478 479 if self.selectedHeight:
479 480 index_list = numpy.where(self.heights >= self.selectedHeight)
480 481 self.height_index = index_list[0]
481 482 self.height_index = self.height_index[0]
482 483 #print(self.height_index)
483 484 data = {}
484 485 meta = {}
485 486 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
486 487 if self.selectedHeight:
487 488 data['spc'] = spc[:,:,self.height_index]
488 489 else:
489 490 data['spc'] = spc
490 491 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
491 492
492 493 return data, meta
493 494
494 495 def plot(self):
495 496 if self.xaxis == "frequency":
496 497 x = self.data.xrange[0][1:]
497 498 self.xlabel = "Frequency (kHz)"
498 499 elif self.xaxis == "time":
499 500 x = self.data.xrange[1]
500 501 self.xlabel = "Time (ms)"
501 502 else:
502 503 x = self.data.xrange[2]
503 504 self.xlabel = "Velocity (m/s)"
504 505
505 506 self.titles = []
506 507
507 508 y = self.data.yrange
508 509 z = self.data[-1]['spc']
509 510 #print(z.shape)
510 511 if self.height_index:
511 512 index = numpy.array(self.height_index)
512 513 else:
513 514 index = numpy.arange(0, len(y), int((len(y))/9))
514 515
515 516 for n, ax in enumerate(self.axes):
516 517 if ax.firsttime:
517 518 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
518 519 self.xmin = self.xmin if self.xmin else -self.xmax
519 520 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
520 521 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
521 522 if self.selectedHeight:
522 523 ax.plt = ax.plot(x, z[n,:])
523 524
524 525 else:
525 526 ax.plt = ax.plot(x, z[n, :, index].T)
526 527 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
527 528 self.figures[0].legend(ax.plt, labels, loc='center right')
528 529 else:
529 530 for i, line in enumerate(ax.plt):
530 531 if self.selectedHeight:
531 532 line.set_data(x, z[n, :])
532 533 else:
533 534 line.set_data(x, z[n, :, index[i]])
534 535 self.titles.append('CH {}'.format(self.channelList[n]))
535 536 plt.suptitle(self.maintitle)
536 537
537 538 class BeaconPhase(Plot):
538 539
539 540 __isConfig = None
540 541 __nsubplots = None
541 542
542 543 PREFIX = 'beacon_phase'
543 544
544 545 def __init__(self):
545 546 Plot.__init__(self)
546 547 self.timerange = 24*60*60
547 548 self.isConfig = False
548 549 self.__nsubplots = 1
549 550 self.counter_imagwr = 0
550 551 self.WIDTH = 800
551 552 self.HEIGHT = 400
552 553 self.WIDTHPROF = 120
553 554 self.HEIGHTPROF = 0
554 555 self.xdata = None
555 556 self.ydata = None
556 557
557 558 self.PLOT_CODE = BEACON_CODE
558 559
559 560 self.FTP_WEI = None
560 561 self.EXP_CODE = None
561 562 self.SUB_EXP_CODE = None
562 563 self.PLOT_POS = None
563 564
564 565 self.filename_phase = None
565 566
566 567 self.figfile = None
567 568
568 569 self.xmin = None
569 570 self.xmax = None
570 571
571 572 def getSubplots(self):
572 573
573 574 ncol = 1
574 575 nrow = 1
575 576
576 577 return nrow, ncol
577 578
578 579 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
579 580
580 581 self.__showprofile = showprofile
581 582 self.nplots = nplots
582 583
583 584 ncolspan = 7
584 585 colspan = 6
585 586 self.__nsubplots = 2
586 587
587 588 self.createFigure(id = id,
588 589 wintitle = wintitle,
589 590 widthplot = self.WIDTH+self.WIDTHPROF,
590 591 heightplot = self.HEIGHT+self.HEIGHTPROF,
591 592 show=show)
592 593
593 594 nrow, ncol = self.getSubplots()
594 595
595 596 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
596 597
597 598 def save_phase(self, filename_phase):
598 599 f = open(filename_phase,'w+')
599 600 f.write('\n\n')
600 601 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
601 602 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
602 603 f.close()
603 604
604 605 def save_data(self, filename_phase, data, data_datetime):
605 606 f=open(filename_phase,'a')
606 607 timetuple_data = data_datetime.timetuple()
607 608 day = str(timetuple_data.tm_mday)
608 609 month = str(timetuple_data.tm_mon)
609 610 year = str(timetuple_data.tm_year)
610 611 hour = str(timetuple_data.tm_hour)
611 612 minute = str(timetuple_data.tm_min)
612 613 second = str(timetuple_data.tm_sec)
613 614 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
614 615 f.close()
615 616
616 617 def plot(self):
617 618 log.warning('TODO: Not yet implemented...')
618 619
619 620 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
620 621 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
621 622 timerange=None,
622 623 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
623 624 server=None, folder=None, username=None, password=None,
624 625 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
625 626
626 627 if dataOut.flagNoData:
627 628 return dataOut
628 629
629 630 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
630 631 return
631 632
632 633 if pairsList == None:
633 634 pairsIndexList = dataOut.pairsIndexList[:10]
634 635 else:
635 636 pairsIndexList = []
636 637 for pair in pairsList:
637 638 if pair not in dataOut.pairsList:
638 639 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
639 640 pairsIndexList.append(dataOut.pairsList.index(pair))
640 641
641 642 if pairsIndexList == []:
642 643 return
643 644
644 645 # if len(pairsIndexList) > 4:
645 646 # pairsIndexList = pairsIndexList[0:4]
646 647
647 648 hmin_index = None
648 649 hmax_index = None
649 650
650 651 if hmin != None and hmax != None:
651 652 indexes = numpy.arange(dataOut.nHeights)
652 653 hmin_list = indexes[dataOut.heightList >= hmin]
653 654 hmax_list = indexes[dataOut.heightList <= hmax]
654 655
655 656 if hmin_list.any():
656 657 hmin_index = hmin_list[0]
657 658
658 659 if hmax_list.any():
659 660 hmax_index = hmax_list[-1]+1
660 661
661 662 x = dataOut.getTimeRange()
662 663
663 664 thisDatetime = dataOut.datatime
664 665
665 666 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
666 667 xlabel = "Local Time"
667 668 ylabel = "Phase (degrees)"
668 669
669 670 update_figfile = False
670 671
671 672 nplots = len(pairsIndexList)
672 673 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
673 674 phase_beacon = numpy.zeros(len(pairsIndexList))
674 675 for i in range(nplots):
675 676 pair = dataOut.pairsList[pairsIndexList[i]]
676 677 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
677 678 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
678 679 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
679 680 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
680 681 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
681 682
682 683 if dataOut.beacon_heiIndexList:
683 684 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
684 685 else:
685 686 phase_beacon[i] = numpy.average(phase)
686 687
687 688 if not self.isConfig:
688 689
689 690 nplots = len(pairsIndexList)
690 691
691 692 self.setup(id=id,
692 693 nplots=nplots,
693 694 wintitle=wintitle,
694 695 showprofile=showprofile,
695 696 show=show)
696 697
697 698 if timerange != None:
698 699 self.timerange = timerange
699 700
700 701 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
701 702
702 703 if ymin == None: ymin = 0
703 704 if ymax == None: ymax = 360
704 705
705 706 self.FTP_WEI = ftp_wei
706 707 self.EXP_CODE = exp_code
707 708 self.SUB_EXP_CODE = sub_exp_code
708 709 self.PLOT_POS = plot_pos
709 710
710 711 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
711 712 self.isConfig = True
712 713 self.figfile = figfile
713 714 self.xdata = numpy.array([])
714 715 self.ydata = numpy.array([])
715 716
716 717 update_figfile = True
717 718
718 719 #open file beacon phase
719 720 path = '%s%03d' %(self.PREFIX, self.id)
720 721 beacon_file = os.path.join(path,'%s.txt'%self.name)
721 722 self.filename_phase = os.path.join(figpath,beacon_file)
722 723 #self.save_phase(self.filename_phase)
723 724
724 725
725 726 #store data beacon phase
726 727 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
727 728
728 729 self.setWinTitle(title)
729 730
730 731
731 732 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
732 733
733 734 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
734 735
735 736 axes = self.axesList[0]
736 737
737 738 self.xdata = numpy.hstack((self.xdata, x[0:1]))
738 739
739 740 if len(self.ydata)==0:
740 741 self.ydata = phase_beacon.reshape(-1,1)
741 742 else:
742 743 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
743 744
744 745
745 746 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
746 747 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
747 748 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
748 749 XAxisAsTime=True, grid='both'
749 750 )
750 751
751 752 self.draw()
752 753
753 754 if dataOut.ltctime >= self.xmax:
754 755 self.counter_imagwr = wr_period
755 756 self.isConfig = False
756 757 update_figfile = True
757 758
758 759 self.save(figpath=figpath,
759 760 figfile=figfile,
760 761 save=save,
761 762 ftp=ftp,
762 763 wr_period=wr_period,
763 764 thisDatetime=thisDatetime,
764 765 update_figfile=update_figfile)
765 766
766 767 return dataOut
767 768
768 769 class NoiselessSpectraPlot(Plot):
769 770 '''
770 771 Plot for Spectra data, subtracting
771 772 the noise in all channels, using for
772 773 amisr-14 data
773 774 '''
774 775
775 776 CODE = 'noiseless_spc'
776 777 colormap = 'nipy_spectral'
777 778 plot_type = 'pcolor'
778 779 buffering = False
779 780 channelList = []
780 781
781 782 def setup(self):
782 783
783 784 self.nplots = len(self.data.channels)
784 785 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
785 786 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
786 787 self.height = 2.6 * self.nrows
787 788
788 789 self.cb_label = 'dB'
789 790 if self.showprofile:
790 791 self.width = 4 * self.ncols
791 792 else:
792 793 self.width = 3.5 * self.ncols
793 794 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
794 795 self.ylabel = 'Range [km]'
795 796
796 797
797 798 def update_list(self,dataOut):
798 799 if len(self.channelList) == 0:
799 800 self.channelList = dataOut.channelList
800 801
801 802 def update(self, dataOut):
802 803
803 804 self.update_list(dataOut)
804 805 data = {}
805 806 meta = {}
806 807 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
807 808 (nch, nff, nh) = dataOut.data_spc.shape
808 809 n1 = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
809 810 noise = numpy.repeat(n1,nff, axis=1).reshape((nch,nff,nh))
810 811 #print(noise.shape, "noise", noise)
811 812
812 813 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - noise
813 814
814 815 data['spc'] = spc
815 816 data['rti'] = dataOut.getPower() - n1
816 817
817 818 data['noise'] = n0
818 819 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
819 820
820 821 return data, meta
821 822
822 823 def plot(self):
823 824 if self.xaxis == "frequency":
824 825 x = self.data.xrange[0]
825 826 self.xlabel = "Frequency (kHz)"
826 827 elif self.xaxis == "time":
827 828 x = self.data.xrange[1]
828 829 self.xlabel = "Time (ms)"
829 830 else:
830 831 x = self.data.xrange[2]
831 832 self.xlabel = "Velocity (m/s)"
832 833
833 834 self.titles = []
834 835 y = self.data.yrange
835 836 self.y = y
836 837
837 838 data = self.data[-1]
838 839 z = data['spc']
839 840
840 841 for n, ax in enumerate(self.axes):
841 842 noise = data['noise'][n]
842 843
843 844 if ax.firsttime:
844 845 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
845 846 self.xmin = self.xmin if self.xmin else -self.xmax
846 847 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
847 848 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
848 849 ax.plt = ax.pcolormesh(x, y, z[n].T,
849 850 vmin=self.zmin,
850 851 vmax=self.zmax,
851 852 cmap=plt.get_cmap(self.colormap)
852 853 )
853 854
854 855 if self.showprofile:
855 856 ax.plt_profile = self.pf_axes[n].plot(
856 857 data['rti'][n], y)[0]
857 858 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
858 859 color="k", linestyle="dashed", lw=1)[0]
859 860
860 861 else:
861 862 ax.plt.set_array(z[n].T.ravel())
862 863 if self.showprofile:
863 864 ax.plt_profile.set_data(data['rti'][n], y)
864 865 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
865 866
866 867 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
867 868
868 869
869 870 class NoiselessRTIPlot(Plot):
870 871 '''
871 872 Plot for RTI data
872 873 '''
873 874
874 875 CODE = 'noiseless_rti'
875 876 colormap = 'jet'
876 877 plot_type = 'pcolorbuffer'
877 878 titles = None
878 879 channelList = []
879 880
880 881 def setup(self):
881 882 self.xaxis = 'time'
882 883 self.ncols = 1
883 884 #print("dataChannels ",self.data.channels)
884 885 self.nrows = len(self.data.channels)
885 886 self.nplots = len(self.data.channels)
886 887 self.ylabel = 'Range [km]'
887 888 self.xlabel = 'Time'
888 889 self.cb_label = 'dB'
889 890 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
890 891 self.titles = ['{} Channel {}'.format(
891 892 self.CODE.upper(), x) for x in range(self.nplots)]
892 893
893 894 def update_list(self,dataOut):
894 895
895 896 self.channelList = dataOut.channelList
896 897
897 898
898 899 def update(self, dataOut):
899 900 if len(self.channelList) == 0:
900 901 self.update_list(dataOut)
901 902 data = {}
902 903 meta = {}
903 904
904 905 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
905 906 (nch, nff, nh) = dataOut.data_spc.shape
906 907 noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
907 908
908 909
909 910 data['noiseless_rti'] = dataOut.getPower() - noise
910 911 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
911 912 return data, meta
912 913
913 914 def plot(self):
914 915
915 916 self.x = self.data.times
916 917 self.y = self.data.yrange
917 918 self.z = self.data['noiseless_rti']
918 919 self.z = numpy.array(self.z, dtype=float)
919 920 self.z = numpy.ma.masked_invalid(self.z)
920 921
921 922 try:
922 923 if self.channelList != None:
923 924 self.titles = ['{} Channel {}'.format(
924 925 self.CODE.upper(), x) for x in self.channelList]
925 926 except:
926 927 if self.channelList.any() != None:
927 928 self.titles = ['{} Channel {}'.format(
928 929 self.CODE.upper(), x) for x in self.channelList]
929 930 if self.decimation is None:
930 931 x, y, z = self.fill_gaps(self.x, self.y, self.z)
931 932 else:
932 933 x, y, z = self.fill_gaps(*self.decimate())
933 934 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
934 935 for n, ax in enumerate(self.axes):
935 936 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
936 937 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
937 938 data = self.data[-1]
938 939 if ax.firsttime:
939 940 ax.plt = ax.pcolormesh(x, y, z[n].T,
940 941 vmin=self.zmin,
941 942 vmax=self.zmax,
942 943 cmap=plt.get_cmap(self.colormap)
943 944 )
944 945 if self.showprofile:
945 946 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
946 947
947 948 if "noise" in self.data:
948 949 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
949 950 color="k", linestyle="dashed", lw=1)[0]
950 951 else:
951 952 ax.collections.remove(ax.collections[0])
952 953 ax.plt = ax.pcolormesh(x, y, z[n].T,
953 954 vmin=self.zmin,
954 955 vmax=self.zmax,
955 956 cmap=plt.get_cmap(self.colormap)
956 957 )
957 958 if self.showprofile:
958 959 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
959 960 if "noise" in self.data:
960 961 ax.plot_noise.set_data(numpy.repeat(
961 962 data['noise'][n], len(self.y)), self.y)
Show file before | Show file after | Hide comments
@@ -1,1688 +1,1814
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 import math
16 16
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 18 from schainpy.model.data.jrodata import Spectra
19 19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 20 from schainpy.utils import log
21 21
22 22 from scipy.optimize import curve_fit
23 23
24 24 class SpectraProc(ProcessingUnit):
25 25
26 26 def __init__(self):
27 27
28 28 ProcessingUnit.__init__(self)
29 29
30 30 self.buffer = None
31 31 self.firstdatatime = None
32 32 self.profIndex = 0
33 33 self.dataOut = Spectra()
34 34 self.id_min = None
35 35 self.id_max = None
36 36 self.setupReq = False #Agregar a todas las unidades de proc
37 37
38 38 def __updateSpecFromVoltage(self):
39 39
40 40 self.dataOut.timeZone = self.dataIn.timeZone
41 41 self.dataOut.dstFlag = self.dataIn.dstFlag
42 42 self.dataOut.errorCount = self.dataIn.errorCount
43 43 self.dataOut.useLocalTime = self.dataIn.useLocalTime
44 44 try:
45 45 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
46 46 except:
47 47 pass
48 48 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
49 49 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
50 50 self.dataOut.channelList = self.dataIn.channelList
51 51 self.dataOut.heightList = self.dataIn.heightList
52 52 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
53 53 self.dataOut.nProfiles = self.dataOut.nFFTPoints
54 54 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
55 55 self.dataOut.utctime = self.firstdatatime
56 56 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
57 57 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
58 58 self.dataOut.flagShiftFFT = False
59 59 self.dataOut.nCohInt = self.dataIn.nCohInt
60 60 self.dataOut.nIncohInt = 1
61 61 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
62 62 self.dataOut.frequency = self.dataIn.frequency
63 63 self.dataOut.realtime = self.dataIn.realtime
64 64 self.dataOut.azimuth = self.dataIn.azimuth
65 65 self.dataOut.zenith = self.dataIn.zenith
66 66 self.dataOut.codeList = self.dataIn.codeList
67 67 self.dataOut.azimuthList = self.dataIn.azimuthList
68 68 self.dataOut.elevationList = self.dataIn.elevationList
69 69
70 70
71 71 def __getFft(self):
72 72 """
73 73 Convierte valores de Voltaje a Spectra
74 74
75 75 Affected:
76 76 self.dataOut.data_spc
77 77 self.dataOut.data_cspc
78 78 self.dataOut.data_dc
79 79 self.dataOut.heightList
80 80 self.profIndex
81 81 self.buffer
82 82 self.dataOut.flagNoData
83 83 """
84 84 fft_volt = numpy.fft.fft(
85 85 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
86 86 fft_volt = fft_volt.astype(numpy.dtype('complex'))
87 87 dc = fft_volt[:, 0, :]
88 88
89 89 # calculo de self-spectra
90 90 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
91 91 spc = fft_volt * numpy.conjugate(fft_volt)
92 92 spc = spc.real
93 93
94 94 blocksize = 0
95 95 blocksize += dc.size
96 96 blocksize += spc.size
97 97
98 98 cspc = None
99 99 pairIndex = 0
100 100 if self.dataOut.pairsList != None:
101 101 # calculo de cross-spectra
102 102 cspc = numpy.zeros(
103 103 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
104 104 for pair in self.dataOut.pairsList:
105 105 if pair[0] not in self.dataOut.channelList:
106 106 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
107 107 str(pair), str(self.dataOut.channelList)))
108 108 if pair[1] not in self.dataOut.channelList:
109 109 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
110 110 str(pair), str(self.dataOut.channelList)))
111 111
112 112 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
113 113 numpy.conjugate(fft_volt[pair[1], :, :])
114 114 pairIndex += 1
115 115 blocksize += cspc.size
116 116
117 117 self.dataOut.data_spc = spc
118 118 self.dataOut.data_cspc = cspc
119 119 self.dataOut.data_dc = dc
120 120 self.dataOut.blockSize = blocksize
121 121 self.dataOut.flagShiftFFT = False
122 122
123 123 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
124 124
125 125 if self.dataIn.type == "Spectra":
126 126
127 127 try:
128 128 self.dataOut.copy(self.dataIn)
129 129
130 130 except Exception as e:
131 131 print(e)
132 132
133 133 if shift_fft:
134 134 #desplaza a la derecha en el eje 2 determinadas posiciones
135 135 shift = int(self.dataOut.nFFTPoints/2)
136 136 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
137 137
138 138 if self.dataOut.data_cspc is not None:
139 139 #desplaza a la derecha en el eje 2 determinadas posiciones
140 140 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
141 141 if pairsList:
142 142 self.__selectPairs(pairsList)
143 143
144 144
145 145 elif self.dataIn.type == "Voltage":
146 146
147 147 self.dataOut.flagNoData = True
148 148
149 149 if nFFTPoints == None:
150 150 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
151 151
152 152 if nProfiles == None:
153 153 nProfiles = nFFTPoints
154 154
155 155 if ippFactor == None:
156 156 self.dataOut.ippFactor = 1
157 157
158 158 self.dataOut.nFFTPoints = nFFTPoints
159 159
160 160 if self.buffer is None:
161 161 self.buffer = numpy.zeros((self.dataIn.nChannels,
162 162 nProfiles,
163 163 self.dataIn.nHeights),
164 164 dtype='complex')
165 165
166 166 if self.dataIn.flagDataAsBlock:
167 167 nVoltProfiles = self.dataIn.data.shape[1]
168 168
169 169 if nVoltProfiles == nProfiles:
170 170 self.buffer = self.dataIn.data.copy()
171 171 self.profIndex = nVoltProfiles
172 172
173 173 elif nVoltProfiles < nProfiles:
174 174
175 175 if self.profIndex == 0:
176 176 self.id_min = 0
177 177 self.id_max = nVoltProfiles
178 178
179 179 self.buffer[:, self.id_min:self.id_max,
180 180 :] = self.dataIn.data
181 181 self.profIndex += nVoltProfiles
182 182 self.id_min += nVoltProfiles
183 183 self.id_max += nVoltProfiles
184 184 else:
185 185 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
186 186 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
187 187 self.dataOut.flagNoData = True
188 188 else:
189 189 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
190 190 self.profIndex += 1
191 191
192 192 if self.firstdatatime == None:
193 193 self.firstdatatime = self.dataIn.utctime
194 194
195 195 if self.profIndex == nProfiles:
196 196 self.__updateSpecFromVoltage()
197 197 if pairsList == None:
198 198 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
199 199 else:
200 200 self.dataOut.pairsList = pairsList
201 201 self.__getFft()
202 202 self.dataOut.flagNoData = False
203 203 self.firstdatatime = None
204 204 self.profIndex = 0
205 self.dataOut.noise_estimation = None
205 206 else:
206 207 raise ValueError("The type of input object '%s' is not valid".format(
207 208 self.dataIn.type))
208 209
209 210 def __selectPairs(self, pairsList):
210 211
211 212 if not pairsList:
212 213 return
213 214
214 215 pairs = []
215 216 pairsIndex = []
216 217
217 218 for pair in pairsList:
218 219 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
219 220 continue
220 221 pairs.append(pair)
221 222 pairsIndex.append(pairs.index(pair))
222 223
223 224 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
224 225 self.dataOut.pairsList = pairs
225 226
226 227 return
227 228
228 229 def selectFFTs(self, minFFT, maxFFT ):
229 230 """
230 231 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
231 232 minFFT<= FFT <= maxFFT
232 233 """
233 234
234 235 if (minFFT > maxFFT):
235 236 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
236 237
237 238 if (minFFT < self.dataOut.getFreqRange()[0]):
238 239 minFFT = self.dataOut.getFreqRange()[0]
239 240
240 241 if (maxFFT > self.dataOut.getFreqRange()[-1]):
241 242 maxFFT = self.dataOut.getFreqRange()[-1]
242 243
243 244 minIndex = 0
244 245 maxIndex = 0
245 246 FFTs = self.dataOut.getFreqRange()
246 247
247 248 inda = numpy.where(FFTs >= minFFT)
248 249 indb = numpy.where(FFTs <= maxFFT)
249 250
250 251 try:
251 252 minIndex = inda[0][0]
252 253 except:
253 254 minIndex = 0
254 255
255 256 try:
256 257 maxIndex = indb[0][-1]
257 258 except:
258 259 maxIndex = len(FFTs)
259 260
260 261 self.selectFFTsByIndex(minIndex, maxIndex)
261 262
262 263 return 1
263 264
264 265 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
265 266 newheis = numpy.where(
266 267 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
267 268
268 269 if hei_ref != None:
269 270 newheis = numpy.where(self.dataOut.heightList > hei_ref)
270 271
271 272 minIndex = min(newheis[0])
272 273 maxIndex = max(newheis[0])
273 274 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
274 275 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
275 276
276 277 # determina indices
277 278 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
278 279 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
279 280 avg_dB = 10 * \
280 281 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
281 282 beacon_dB = numpy.sort(avg_dB)[-nheis:]
282 283 beacon_heiIndexList = []
283 284 for val in avg_dB.tolist():
284 285 if val >= beacon_dB[0]:
285 286 beacon_heiIndexList.append(avg_dB.tolist().index(val))
286 287
287 288 #data_spc = data_spc[:,:,beacon_heiIndexList]
288 289 data_cspc = None
289 290 if self.dataOut.data_cspc is not None:
290 291 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
291 292 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
292 293
293 294 data_dc = None
294 295 if self.dataOut.data_dc is not None:
295 296 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
296 297 #data_dc = data_dc[:,beacon_heiIndexList]
297 298
298 299 self.dataOut.data_spc = data_spc
299 300 self.dataOut.data_cspc = data_cspc
300 301 self.dataOut.data_dc = data_dc
301 302 self.dataOut.heightList = heightList
302 303 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
303 304
304 305 return 1
305 306
306 307 def selectFFTsByIndex(self, minIndex, maxIndex):
307 308 """
308 309
309 310 """
310 311
311 312 if (minIndex < 0) or (minIndex > maxIndex):
312 313 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
313 314
314 315 if (maxIndex >= self.dataOut.nProfiles):
315 316 maxIndex = self.dataOut.nProfiles-1
316 317
317 318 #Spectra
318 319 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
319 320
320 321 data_cspc = None
321 322 if self.dataOut.data_cspc is not None:
322 323 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
323 324
324 325 data_dc = None
325 326 if self.dataOut.data_dc is not None:
326 327 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
327 328
328 329 self.dataOut.data_spc = data_spc
329 330 self.dataOut.data_cspc = data_cspc
330 331 self.dataOut.data_dc = data_dc
331 332
332 333 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
333 334 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
334 335 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
335 336
336 337 return 1
337 338
338 339 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
339 340 # validacion de rango
340 341 if minHei == None:
341 342 minHei = self.dataOut.heightList[0]
342 343
343 344 if maxHei == None:
344 345 maxHei = self.dataOut.heightList[-1]
345 346
346 347 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
347 348 print('minHei: %.2f is out of the heights range' % (minHei))
348 349 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
349 350 minHei = self.dataOut.heightList[0]
350 351
351 352 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
352 353 print('maxHei: %.2f is out of the heights range' % (maxHei))
353 354 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
354 355 maxHei = self.dataOut.heightList[-1]
355 356
356 357 # validacion de velocidades
357 358 velrange = self.dataOut.getVelRange(1)
358 359
359 360 if minVel == None:
360 361 minVel = velrange[0]
361 362
362 363 if maxVel == None:
363 364 maxVel = velrange[-1]
364 365
365 366 if (minVel < velrange[0]) or (minVel > maxVel):
366 367 print('minVel: %.2f is out of the velocity range' % (minVel))
367 368 print('minVel is setting to %.2f' % (velrange[0]))
368 369 minVel = velrange[0]
369 370
370 371 if (maxVel > velrange[-1]) or (maxVel < minVel):
371 372 print('maxVel: %.2f is out of the velocity range' % (maxVel))
372 373 print('maxVel is setting to %.2f' % (velrange[-1]))
373 374 maxVel = velrange[-1]
374 375
375 376 # seleccion de indices para rango
376 377 minIndex = 0
377 378 maxIndex = 0
378 379 heights = self.dataOut.heightList
379 380
380 381 inda = numpy.where(heights >= minHei)
381 382 indb = numpy.where(heights <= maxHei)
382 383
383 384 try:
384 385 minIndex = inda[0][0]
385 386 except:
386 387 minIndex = 0
387 388
388 389 try:
389 390 maxIndex = indb[0][-1]
390 391 except:
391 392 maxIndex = len(heights)
392 393
393 394 if (minIndex < 0) or (minIndex > maxIndex):
394 395 raise ValueError("some value in (%d,%d) is not valid" % (
395 396 minIndex, maxIndex))
396 397
397 398 if (maxIndex >= self.dataOut.nHeights):
398 399 maxIndex = self.dataOut.nHeights - 1
399 400
400 401 # seleccion de indices para velocidades
401 402 indminvel = numpy.where(velrange >= minVel)
402 403 indmaxvel = numpy.where(velrange <= maxVel)
403 404 try:
404 405 minIndexVel = indminvel[0][0]
405 406 except:
406 407 minIndexVel = 0
407 408
408 409 try:
409 410 maxIndexVel = indmaxvel[0][-1]
410 411 except:
411 412 maxIndexVel = len(velrange)
412 413
413 414 # seleccion del espectro
414 415 data_spc = self.dataOut.data_spc[:,
415 416 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
416 417 # estimacion de ruido
417 418 noise = numpy.zeros(self.dataOut.nChannels)
418 419
419 420 for channel in range(self.dataOut.nChannels):
420 421 daux = data_spc[channel, :, :]
421 422 sortdata = numpy.sort(daux, axis=None)
422 423 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
423 424
424 425 self.dataOut.noise_estimation = noise.copy()
425 426
426 427 return 1
427 428
428 429 class removeDC(Operation):
429 430
430 431 def run(self, dataOut, mode=2):
431 432 self.dataOut = dataOut
432 433 jspectra = self.dataOut.data_spc
433 434 jcspectra = self.dataOut.data_cspc
434 435
435 436 num_chan = jspectra.shape[0]
436 437 num_hei = jspectra.shape[2]
437 438
438 439 if jcspectra is not None:
439 440 jcspectraExist = True
440 441 num_pairs = jcspectra.shape[0]
441 442 else:
442 443 jcspectraExist = False
443 444
444 445 freq_dc = int(jspectra.shape[1] / 2)
445 446 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
446 447 ind_vel = ind_vel.astype(int)
447 448
448 449 if ind_vel[0] < 0:
449 450 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
450 451
451 452 if mode == 1:
452 453 jspectra[:, freq_dc, :] = (
453 454 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
454 455
455 456 if jcspectraExist:
456 457 jcspectra[:, freq_dc, :] = (
457 458 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
458 459
459 460 if mode == 2:
460 461
461 462 vel = numpy.array([-2, -1, 1, 2])
462 463 xx = numpy.zeros([4, 4])
463 464
464 465 for fil in range(4):
465 466 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
466 467
467 468 xx_inv = numpy.linalg.inv(xx)
468 469 xx_aux = xx_inv[0, :]
469 470
470 471 for ich in range(num_chan):
471 472 yy = jspectra[ich, ind_vel, :]
472 473 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
473 474
474 475 junkid = jspectra[ich, freq_dc, :] <= 0
475 476 cjunkid = sum(junkid)
476 477
477 478 if cjunkid.any():
478 479 jspectra[ich, freq_dc, junkid.nonzero()] = (
479 480 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
480 481
481 482 if jcspectraExist:
482 483 for ip in range(num_pairs):
483 484 yy = jcspectra[ip, ind_vel, :]
484 485 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
485 486
486 487 self.dataOut.data_spc = jspectra
487 488 self.dataOut.data_cspc = jcspectra
488 489
489 490 return self.dataOut
490 491
492 class getNoise(Operation):
493 def __init__(self):
494
495 Operation.__init__(self)
496
497 def run(self, dataOut, minHei=None, maxHei=None, minVel=None, maxVel=None, minFreq= None, maxFreq=None,):
498 self.dataOut = dataOut.copy()
499 print("1: ",dataOut.noise_estimation, dataOut.normFactor)
500
501 if minHei == None:
502 minHei = self.dataOut.heightList[0]
503
504 if maxHei == None:
505 maxHei = self.dataOut.heightList[-1]
506
507 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
508 print('minHei: %.2f is out of the heights range' % (minHei))
509 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
510 minHei = self.dataOut.heightList[0]
511
512 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
513 print('maxHei: %.2f is out of the heights range' % (maxHei))
514 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
515 maxHei = self.dataOut.heightList[-1]
516
517
518 #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
519 minIndexFFT = 0
520 maxIndexFFT = 0
521 # validacion de velocidades
522 indminPoint = None
523 indmaxPoint = None
524
525 if minVel == None and maxVel == None:
526
527 freqrange = self.dataOut.getFreqRange(1)
528
529 if minFreq == None:
530 minFreq = freqrange[0]
531
532 if maxFreq == None:
533 maxFreq = freqrange[-1]
534
535 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
536 print('minFreq: %.2f is out of the frequency range' % (minFreq))
537 print('minFreq is setting to %.2f' % (freqrange[0]))
538 minFreq = freqrange[0]
539
540 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
541 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
542 print('maxFreq is setting to %.2f' % (freqrange[-1]))
543 maxFreq = freqrange[-1]
544
545 indminPoint = numpy.where(freqrange >= minFreq)
546 indmaxPoint = numpy.where(freqrange <= maxFreq)
547
548 else:
549 velrange = self.dataOut.getVelRange(1)
550
551 if minVel == None:
552 minVel = velrange[0]
553
554 if maxVel == None:
555 maxVel = velrange[-1]
556
557 if (minVel < velrange[0]) or (minVel > maxVel):
558 print('minVel: %.2f is out of the velocity range' % (minVel))
559 print('minVel is setting to %.2f' % (velrange[0]))
560 minVel = velrange[0]
561
562 if (maxVel > velrange[-1]) or (maxVel < minVel):
563 print('maxVel: %.2f is out of the velocity range' % (maxVel))
564 print('maxVel is setting to %.2f' % (velrange[-1]))
565 maxVel = velrange[-1]
566
567 indminPoint = numpy.where(velrange >= minVel)
568 indmaxPoint = numpy.where(velrange <= maxVel)
569
570
571 # seleccion de indices para rango
572 minIndex = 0
573 maxIndex = 0
574 heights = self.dataOut.heightList
575
576 inda = numpy.where(heights >= minHei)
577 indb = numpy.where(heights <= maxHei)
578
579 try:
580 minIndex = inda[0][0]
581 except:
582 minIndex = 0
583
584 try:
585 maxIndex = indb[0][-1]
586 except:
587 maxIndex = len(heights)
588
589 if (minIndex < 0) or (minIndex > maxIndex):
590 raise ValueError("some value in (%d,%d) is not valid" % (
591 minIndex, maxIndex))
592
593 if (maxIndex >= self.dataOut.nHeights):
594 maxIndex = self.dataOut.nHeights - 1
595 #############################################################3
596 # seleccion de indices para velocidades
597
598 try:
599 minIndexFFT = indminPoint[0][0]
600 except:
601 minIndexFFT = 0
602
603 try:
604 maxIndexFFT = indmaxPoint[0][-1]
605 except:
606 maxIndexFFT = len( self.dataOut.getFreqRange(1))
607
608 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
609 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
610
611 self.dataOut.noise_estimation = noise.copy()
612 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
613 return self.dataOut
614
615
616
491 617 # import matplotlib.pyplot as plt
492 618
493 619 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
494 620 z = (x - a1) / a2
495 621 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
496 622 return y
497 623
498 624
499 625 class CleanRayleigh(Operation):
500 626
501 627 def __init__(self):
502 628
503 629 Operation.__init__(self)
504 630 self.i=0
505 631 self.isConfig = False
506 632 self.__dataReady = False
507 633 self.__profIndex = 0
508 634 self.byTime = False
509 635 self.byProfiles = False
510 636
511 637 self.bloques = None
512 638 self.bloque0 = None
513 639
514 640 self.index = 0
515 641
516 642 self.buffer = 0
517 643 self.buffer2 = 0
518 644 self.buffer3 = 0
519 645
520 646
521 647 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
522 648
523 649 self.nChannels = dataOut.nChannels
524 650 self.nProf = dataOut.nProfiles
525 651 self.nPairs = dataOut.data_cspc.shape[0]
526 652 self.pairsArray = numpy.array(dataOut.pairsList)
527 653 self.spectra = dataOut.data_spc
528 654 self.cspectra = dataOut.data_cspc
529 655 self.heights = dataOut.heightList #alturas totales
530 656 self.nHeights = len(self.heights)
531 657 self.min_hei = min_hei
532 658 self.max_hei = max_hei
533 659 if (self.min_hei == None):
534 660 self.min_hei = 0
535 661 if (self.max_hei == None):
536 662 self.max_hei = dataOut.heightList[-1]
537 663 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
538 664 self.heightsClean = self.heights[self.hval] #alturas filtradas
539 665 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
540 666 self.nHeightsClean = len(self.heightsClean)
541 667 self.channels = dataOut.channelList
542 668 self.nChan = len(self.channels)
543 669 self.nIncohInt = dataOut.nIncohInt
544 670 self.__initime = dataOut.utctime
545 671 self.maxAltInd = self.hval[-1]+1
546 672 self.minAltInd = self.hval[0]
547 673
548 674 self.crosspairs = dataOut.pairsList
549 675 self.nPairs = len(self.crosspairs)
550 676 self.normFactor = dataOut.normFactor
551 677 self.nFFTPoints = dataOut.nFFTPoints
552 678 self.ippSeconds = dataOut.ippSeconds
553 679 self.currentTime = self.__initime
554 680 self.pairsArray = numpy.array(dataOut.pairsList)
555 681 self.factor_stdv = factor_stdv
556 682
557 683 if n != None :
558 684 self.byProfiles = True
559 685 self.nIntProfiles = n
560 686 else:
561 687 self.__integrationtime = timeInterval
562 688
563 689 self.__dataReady = False
564 690 self.isConfig = True
565 691
566 692
567 693
568 694 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
569 695
570 696 if not self.isConfig :
571 697
572 698 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
573 699
574 700 tini=dataOut.utctime
575 701
576 702 if self.byProfiles:
577 703 if self.__profIndex == self.nIntProfiles:
578 704 self.__dataReady = True
579 705 else:
580 706 if (tini - self.__initime) >= self.__integrationtime:
581 707
582 708 self.__dataReady = True
583 709 self.__initime = tini
584 710
585 711 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
586 712
587 713 if self.__dataReady:
588 714
589 715 self.__profIndex = 0
590 716 jspc = self.buffer
591 717 jcspc = self.buffer2
592 718 #jnoise = self.buffer3
593 719 self.buffer = dataOut.data_spc
594 720 self.buffer2 = dataOut.data_cspc
595 721 #self.buffer3 = dataOut.noise
596 722 self.currentTime = dataOut.utctime
597 723 if numpy.any(jspc) :
598 724 #print( jspc.shape, jcspc.shape)
599 725 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
600 726 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
601 727 self.__dataReady = False
602 728 #print( jspc.shape, jcspc.shape)
603 729 dataOut.flagNoData = False
604 730 else:
605 731 dataOut.flagNoData = True
606 732 self.__dataReady = False
607 733 return dataOut
608 734 else:
609 735 #print( len(self.buffer))
610 736 if numpy.any(self.buffer):
611 737 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
612 738 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
613 739 self.buffer3 += dataOut.data_dc
614 740 else:
615 741 self.buffer = dataOut.data_spc
616 742 self.buffer2 = dataOut.data_cspc
617 743 self.buffer3 = dataOut.data_dc
618 744 #print self.index, self.fint
619 745 #print self.buffer2.shape
620 746 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
621 747 self.__profIndex += 1
622 748 return dataOut ## NOTE: REV
623 749
624 750
625 751 #index = tini.tm_hour*12+tini.tm_min/5
626 752 '''REVISAR'''
627 753 # jspc = jspc/self.nFFTPoints/self.normFactor
628 754 # jcspc = jcspc/self.nFFTPoints/self.normFactor
629 755
630 756
631 757
632 758 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
633 759 dataOut.data_spc = tmp_spectra
634 760 dataOut.data_cspc = tmp_cspectra
635 761
636 762 #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
637 763
638 764 dataOut.data_dc = self.buffer3
639 765 dataOut.nIncohInt *= self.nIntProfiles
640 766 dataOut.utctime = self.currentTime #tiempo promediado
641 767 #print("Time: ",time.localtime(dataOut.utctime))
642 768 # dataOut.data_spc = sat_spectra
643 769 # dataOut.data_cspc = sat_cspectra
644 770 self.buffer = 0
645 771 self.buffer2 = 0
646 772 self.buffer3 = 0
647 773
648 774 return dataOut
649 775
650 776 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
651 777 #print("OP cleanRayleigh")
652 778 #import matplotlib.pyplot as plt
653 779 #for k in range(149):
654 780 #channelsProcssd = []
655 781 #channelA_ok = False
656 782 #rfunc = cspectra.copy() #self.bloques
657 783 rfunc = spectra.copy()
658 784 #rfunc = cspectra
659 785 #val_spc = spectra*0.0 #self.bloque0*0.0
660 786 #val_cspc = cspectra*0.0 #self.bloques*0.0
661 787 #in_sat_spectra = spectra.copy() #self.bloque0
662 788 #in_sat_cspectra = cspectra.copy() #self.bloques
663 789
664 790
665 791 ###ONLY FOR TEST:
666 792 raxs = math.ceil(math.sqrt(self.nPairs))
667 793 caxs = math.ceil(self.nPairs/raxs)
668 794 if self.nPairs <4:
669 795 raxs = 2
670 796 caxs = 2
671 797 #print(raxs, caxs)
672 798 fft_rev = 14 #nFFT to plot
673 799 hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
674 800 hei_rev = hei_rev[0]
675 801 #print(hei_rev)
676 802
677 803 #print numpy.absolute(rfunc[:,0,0,14])
678 804
679 805 gauss_fit, covariance = None, None
680 806 for ih in range(self.minAltInd,self.maxAltInd):
681 807 for ifreq in range(self.nFFTPoints):
682 808 '''
683 809 ###ONLY FOR TEST:
684 810 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
685 811 fig, axs = plt.subplots(raxs, caxs)
686 812 fig2, axs2 = plt.subplots(raxs, caxs)
687 813 col_ax = 0
688 814 row_ax = 0
689 815 '''
690 816 #print(self.nPairs)
691 817 for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
692 818 # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
693 819 # continue
694 820 # if not self.crosspairs[ii][0] in channelsProcssd:
695 821 # channelA_ok = True
696 822 #print("pair: ",self.crosspairs[ii])
697 823 '''
698 824 ###ONLY FOR TEST:
699 825 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
700 826 col_ax = 0
701 827 row_ax += 1
702 828 '''
703 829 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
704 830 #print(func2clean.shape)
705 831 val = (numpy.isfinite(func2clean)==True).nonzero()
706 832
707 833 if len(val)>0: #limitador
708 834 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
709 835 if min_val <= -40 :
710 836 min_val = -40
711 837 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
712 838 if max_val >= 200 :
713 839 max_val = 200
714 840 #print min_val, max_val
715 841 step = 1
716 842 #print("Getting bins and the histogram")
717 843 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
718 844 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
719 845 #print(len(y_dist),len(binstep[:-1]))
720 846 #print(row_ax,col_ax, " ..")
721 847 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
722 848 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
723 849 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
724 850 parg = [numpy.amax(y_dist),mean,sigma]
725 851
726 852 newY = None
727 853
728 854 try :
729 855 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
730 856 mode = gauss_fit[1]
731 857 stdv = gauss_fit[2]
732 858 #print(" FIT OK",gauss_fit)
733 859 '''
734 860 ###ONLY FOR TEST:
735 861 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
736 862 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
737 863 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
738 864 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
739 865 axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
740 866 '''
741 867 except:
742 868 mode = mean
743 869 stdv = sigma
744 870 #print("FIT FAIL")
745 871 #continue
746 872
747 873
748 874 #print(mode,stdv)
749 875 #Removing echoes greater than mode + std_factor*stdv
750 876 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
751 877 #noval tiene los indices que se van a remover
752 878 #print("Chan ",ii," novals: ",len(noval[0]))
753 879 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
754 880 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
755 881 #print(novall)
756 882 #print(" ",self.pairsArray[ii])
757 883 #cross_pairs = self.pairsArray[ii]
758 884 #Getting coherent echoes which are removed.
759 885 # if len(novall[0]) > 0:
760 886 #
761 887 # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
762 888 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
763 889 # val_cspc[novall[0],ii,ifreq,ih] = 1
764 890 #print("OUT NOVALL 1")
765 891 try:
766 892 pair = (self.channels[ii],self.channels[ii + 1])
767 893 except:
768 894 pair = (99,99)
769 895 #print("par ", pair)
770 896 if ( pair in self.crosspairs):
771 897 q = self.crosspairs.index(pair)
772 898 #print("está aqui: ", q, (ii,ii + 1))
773 899 new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
774 900 cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
775 901
776 902 #if channelA_ok:
777 903 #chA = self.channels.index(cross_pairs[0])
778 904 new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
779 905 spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
780 906 #channelA_ok = False
781 907
782 908 # chB = self.channels.index(cross_pairs[1])
783 909 # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
784 910 # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
785 911 #
786 912 # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
787 913 # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
788 914 '''
789 915 ###ONLY FOR TEST:
790 916 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
791 917 func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
792 918 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
793 919 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
794 920 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
795 921 axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
796 922 '''
797 923 '''
798 924 ###ONLY FOR TEST:
799 925 col_ax += 1 #contador de ploteo columnas
800 926 ##print(col_ax)
801 927 ###ONLY FOR TEST:
802 928 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
803 929 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
804 930 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
805 931 fig.suptitle(title)
806 932 fig2.suptitle(title2)
807 933 plt.show()
808 934 '''
809 935 ##################################################################################################
810 936
811 937 #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
812 938 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
813 939 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
814 940 for ih in range(self.nHeights):
815 941 for ifreq in range(self.nFFTPoints):
816 942 for ich in range(self.nChan):
817 943 tmp = spectra[:,ich,ifreq,ih]
818 944 valid = (numpy.isfinite(tmp[:])==True).nonzero()
819 945
820 946 if len(valid[0]) >0 :
821 947 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
822 948
823 949 for icr in range(self.nPairs):
824 950 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
825 951 valid = (numpy.isfinite(tmp)==True).nonzero()
826 952 if len(valid[0]) > 0:
827 953 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
828 954
829 955 return out_spectra, out_cspectra
830 956
831 957 def REM_ISOLATED_POINTS(self,array,rth):
832 958 # import matplotlib.pyplot as plt
833 959 if rth == None :
834 960 rth = 4
835 961 #print("REM ISO")
836 962 num_prof = len(array[0,:,0])
837 963 num_hei = len(array[0,0,:])
838 964 n2d = len(array[:,0,0])
839 965
840 966 for ii in range(n2d) :
841 967 #print ii,n2d
842 968 tmp = array[ii,:,:]
843 969 #print tmp.shape, array[ii,101,:],array[ii,102,:]
844 970
845 971 # fig = plt.figure(figsize=(6,5))
846 972 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
847 973 # ax = fig.add_axes([left, bottom, width, height])
848 974 # x = range(num_prof)
849 975 # y = range(num_hei)
850 976 # cp = ax.contour(y,x,tmp)
851 977 # ax.clabel(cp, inline=True,fontsize=10)
852 978 # plt.show()
853 979
854 980 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
855 981 tmp = numpy.reshape(tmp,num_prof*num_hei)
856 982 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
857 983 indxs2 = (tmp > 0).nonzero()
858 984
859 985 indxs1 = (indxs1[0])
860 986 indxs2 = indxs2[0]
861 987 #indxs1 = numpy.array(indxs1[0])
862 988 #indxs2 = numpy.array(indxs2[0])
863 989 indxs = None
864 990 #print indxs1 , indxs2
865 991 for iv in range(len(indxs2)):
866 992 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
867 993 #print len(indxs2), indv
868 994 if len(indv[0]) > 0 :
869 995 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
870 996 # print indxs
871 997 indxs = indxs[1:]
872 998 #print(indxs, len(indxs))
873 999 if len(indxs) < 4 :
874 1000 array[ii,:,:] = 0.
875 1001 return
876 1002
877 1003 xpos = numpy.mod(indxs ,num_hei)
878 1004 ypos = (indxs / num_hei)
879 1005 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
880 1006 #print sx
881 1007 xpos = xpos[sx]
882 1008 ypos = ypos[sx]
883 1009
884 1010 # *********************************** Cleaning isolated points **********************************
885 1011 ic = 0
886 1012 while True :
887 1013 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
888 1014 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
889 1015 #plt.plot(r)
890 1016 #plt.show()
891 1017 no_coh1 = (numpy.isfinite(r)==True).nonzero()
892 1018 no_coh2 = (r <= rth).nonzero()
893 1019 #print r, no_coh1, no_coh2
894 1020 no_coh1 = numpy.array(no_coh1[0])
895 1021 no_coh2 = numpy.array(no_coh2[0])
896 1022 no_coh = None
897 1023 #print valid1 , valid2
898 1024 for iv in range(len(no_coh2)):
899 1025 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
900 1026 if len(indv[0]) > 0 :
901 1027 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
902 1028 no_coh = no_coh[1:]
903 1029 #print len(no_coh), no_coh
904 1030 if len(no_coh) < 4 :
905 1031 #print xpos[ic], ypos[ic], ic
906 1032 # plt.plot(r)
907 1033 # plt.show()
908 1034 xpos[ic] = numpy.nan
909 1035 ypos[ic] = numpy.nan
910 1036
911 1037 ic = ic + 1
912 1038 if (ic == len(indxs)) :
913 1039 break
914 1040 #print( xpos, ypos)
915 1041
916 1042 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
917 1043 #print indxs[0]
918 1044 if len(indxs[0]) < 4 :
919 1045 array[ii,:,:] = 0.
920 1046 return
921 1047
922 1048 xpos = xpos[indxs[0]]
923 1049 ypos = ypos[indxs[0]]
924 1050 for i in range(0,len(ypos)):
925 1051 ypos[i]=int(ypos[i])
926 1052 junk = tmp
927 1053 tmp = junk*0.0
928 1054
929 1055 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
930 1056 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
931 1057
932 1058 #print array.shape
933 1059 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
934 1060 #print tmp.shape
935 1061
936 1062 # fig = plt.figure(figsize=(6,5))
937 1063 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
938 1064 # ax = fig.add_axes([left, bottom, width, height])
939 1065 # x = range(num_prof)
940 1066 # y = range(num_hei)
941 1067 # cp = ax.contour(y,x,array[ii,:,:])
942 1068 # ax.clabel(cp, inline=True,fontsize=10)
943 1069 # plt.show()
944 1070 return array
945 1071
946 1072
947 1073 class IntegrationFaradaySpectra(Operation):
948 1074
949 1075 __profIndex = 0
950 1076 __withOverapping = False
951 1077
952 1078 __byTime = False
953 1079 __initime = None
954 1080 __lastdatatime = None
955 1081 __integrationtime = None
956 1082
957 1083 __buffer_spc = None
958 1084 __buffer_cspc = None
959 1085 __buffer_dc = None
960 1086
961 1087 __dataReady = False
962 1088
963 1089 __timeInterval = None
964 1090
965 1091 n = None
966 1092
967 1093 def __init__(self):
968 1094
969 1095 Operation.__init__(self)
970 1096
971 1097 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None):
972 1098 """
973 1099 Set the parameters of the integration class.
974 1100
975 1101 Inputs:
976 1102
977 1103 n : Number of coherent integrations
978 1104 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
979 1105 overlapping :
980 1106
981 1107 """
982 1108
983 1109 self.__initime = None
984 1110 self.__lastdatatime = 0
985 1111
986 1112 self.__buffer_spc = []
987 1113 self.__buffer_cspc = []
988 1114 self.__buffer_dc = 0
989 1115
990 1116 self.__profIndex = 0
991 1117 self.__dataReady = False
992 1118 self.__byTime = False
993 1119
994 1120 #self.ByLags = dataOut.ByLags ###REDEFINIR
995 1121 self.ByLags = False
996 1122
997 1123 if DPL != None:
998 1124 self.DPL=DPL
999 1125 else:
1000 1126 #self.DPL=dataOut.DPL ###REDEFINIR
1001 1127 self.DPL=0
1002 1128
1003 1129 if n is None and timeInterval is None:
1004 1130 raise ValueError("n or timeInterval should be specified ...")
1005 1131
1006 1132 if n is not None:
1007 1133 self.n = int(n)
1008 1134 else:
1009 1135
1010 1136 self.__integrationtime = int(timeInterval)
1011 1137 self.n = None
1012 1138 self.__byTime = True
1013 1139
1014 1140 def putData(self, data_spc, data_cspc, data_dc):
1015 1141 """
1016 1142 Add a profile to the __buffer_spc and increase in one the __profileIndex
1017 1143
1018 1144 """
1019 1145
1020 1146 self.__buffer_spc.append(data_spc)
1021 1147
1022 1148 if data_cspc is None:
1023 1149 self.__buffer_cspc = None
1024 1150 else:
1025 1151 self.__buffer_cspc.append(data_cspc)
1026 1152
1027 1153 if data_dc is None:
1028 1154 self.__buffer_dc = None
1029 1155 else:
1030 1156 self.__buffer_dc += data_dc
1031 1157
1032 1158 self.__profIndex += 1
1033 1159
1034 1160 return
1035 1161
1036 1162 def hildebrand_sekhon_Integration(self,data,navg):
1037 1163
1038 1164 sortdata = numpy.sort(data, axis=None)
1039 1165 sortID=data.argsort()
1040 1166 lenOfData = len(sortdata)
1041 1167 nums_min = lenOfData*0.75
1042 1168 if nums_min <= 5:
1043 1169 nums_min = 5
1044 1170 sump = 0.
1045 1171 sumq = 0.
1046 1172 j = 0
1047 1173 cont = 1
1048 1174 while((cont == 1)and(j < lenOfData)):
1049 1175 sump += sortdata[j]
1050 1176 sumq += sortdata[j]**2
1051 1177 if j > nums_min:
1052 1178 rtest = float(j)/(j-1) + 1.0/navg
1053 1179 if ((sumq*j) > (rtest*sump**2)):
1054 1180 j = j - 1
1055 1181 sump = sump - sortdata[j]
1056 1182 sumq = sumq - sortdata[j]**2
1057 1183 cont = 0
1058 1184 j += 1
1059 1185 #lnoise = sump / j
1060 1186
1061 1187 return j,sortID
1062 1188
1063 1189 def pushData(self):
1064 1190 """
1065 1191 Return the sum of the last profiles and the profiles used in the sum.
1066 1192
1067 1193 Affected:
1068 1194
1069 1195 self.__profileIndex
1070 1196
1071 1197 """
1072 1198 bufferH=None
1073 1199 buffer=None
1074 1200 buffer1=None
1075 1201 buffer_cspc=None
1076 1202 self.__buffer_spc=numpy.array(self.__buffer_spc)
1077 1203 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1078 1204 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1079 1205 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1080 1206 for k in range(7,self.nHeights):
1081 1207 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1082 1208 outliers_IDs_cspc=[]
1083 1209 cspc_outliers_exist=False
1084 1210 for i in range(self.nChannels):#dataOut.nChannels):
1085 1211
1086 1212 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1087 1213 indexes=[]
1088 1214 #sortIDs=[]
1089 1215 outliers_IDs=[]
1090 1216
1091 1217 for j in range(self.nProfiles):
1092 1218 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1093 1219 # continue
1094 1220 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1095 1221 # continue
1096 1222 buffer=buffer1[:,j]
1097 1223 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1098 1224
1099 1225 indexes.append(index)
1100 1226 #sortIDs.append(sortID)
1101 1227 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1102 1228
1103 1229 outliers_IDs=numpy.array(outliers_IDs)
1104 1230 outliers_IDs=outliers_IDs.ravel()
1105 1231 outliers_IDs=numpy.unique(outliers_IDs)
1106 1232 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1107 1233 indexes=numpy.array(indexes)
1108 1234 indexmin=numpy.min(indexes)
1109 1235
1110 1236 if indexmin != buffer1.shape[0]:
1111 1237 cspc_outliers_exist=True
1112 1238 ###sortdata=numpy.sort(buffer1,axis=0)
1113 1239 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
1114 1240 lt=outliers_IDs
1115 1241 avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1116 1242
1117 1243 for p in list(outliers_IDs):
1118 1244 buffer1[p,:]=avg
1119 1245
1120 1246 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1121 1247 ###cspc IDs
1122 1248 #indexmin_cspc+=indexmin_cspc
1123 1249 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1124 1250
1125 1251 #if not breakFlag:
1126 1252 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1127 1253 if cspc_outliers_exist:
1128 1254 #sortdata=numpy.sort(buffer_cspc,axis=0)
1129 1255 #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0)
1130 1256 lt=outliers_IDs_cspc
1131 1257
1132 1258 avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1133 1259 for p in list(outliers_IDs_cspc):
1134 1260 buffer_cspc[p,:]=avg
1135 1261
1136 1262 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1137 1263 #else:
1138 1264 #break
1139 1265
1140 1266
1141 1267
1142 1268
1143 1269 buffer=None
1144 1270 bufferH=None
1145 1271 buffer1=None
1146 1272 buffer_cspc=None
1147 1273
1148 1274 #print("cpsc",self.__buffer_cspc[:,0,0,0,0])
1149 1275 #print(self.__profIndex)
1150 1276 #exit()
1151 1277
1152 1278 buffer=None
1153 1279 #print(self.__buffer_spc[:,1,3,20,0])
1154 1280 #print(self.__buffer_spc[:,1,5,37,0])
1155 1281 data_spc = numpy.sum(self.__buffer_spc,axis=0)
1156 1282 data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1157 1283
1158 1284 #print(numpy.shape(data_spc))
1159 1285 #data_spc[1,4,20,0]=numpy.nan
1160 1286
1161 1287 #data_cspc = self.__buffer_cspc
1162 1288 data_dc = self.__buffer_dc
1163 1289 n = self.__profIndex
1164 1290
1165 1291 self.__buffer_spc = []
1166 1292 self.__buffer_cspc = []
1167 1293 self.__buffer_dc = 0
1168 1294 self.__profIndex = 0
1169 1295
1170 1296 return data_spc, data_cspc, data_dc, n
1171 1297
1172 1298 def byProfiles(self, *args):
1173 1299
1174 1300 self.__dataReady = False
1175 1301 avgdata_spc = None
1176 1302 avgdata_cspc = None
1177 1303 avgdata_dc = None
1178 1304
1179 1305 self.putData(*args)
1180 1306
1181 1307 if self.__profIndex == self.n:
1182 1308
1183 1309 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1184 1310 self.n = n
1185 1311 self.__dataReady = True
1186 1312
1187 1313 return avgdata_spc, avgdata_cspc, avgdata_dc
1188 1314
1189 1315 def byTime(self, datatime, *args):
1190 1316
1191 1317 self.__dataReady = False
1192 1318 avgdata_spc = None
1193 1319 avgdata_cspc = None
1194 1320 avgdata_dc = None
1195 1321
1196 1322 self.putData(*args)
1197 1323
1198 1324 if (datatime - self.__initime) >= self.__integrationtime:
1199 1325 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1200 1326 self.n = n
1201 1327 self.__dataReady = True
1202 1328
1203 1329 return avgdata_spc, avgdata_cspc, avgdata_dc
1204 1330
1205 1331 def integrate(self, datatime, *args):
1206 1332
1207 1333 if self.__profIndex == 0:
1208 1334 self.__initime = datatime
1209 1335
1210 1336 if self.__byTime:
1211 1337 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1212 1338 datatime, *args)
1213 1339 else:
1214 1340 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1215 1341
1216 1342 if not self.__dataReady:
1217 1343 return None, None, None, None
1218 1344
1219 1345 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1220 1346
1221 1347 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False):
1222 1348 if n == 1:
1223 1349 return dataOut
1224 1350
1225 1351 dataOut.flagNoData = True
1226 1352
1227 1353 if not self.isConfig:
1228 1354 self.setup(dataOut, n, timeInterval, overlapping,DPL )
1229 1355 self.isConfig = True
1230 1356
1231 1357 if not self.ByLags:
1232 1358 self.nProfiles=dataOut.nProfiles
1233 1359 self.nChannels=dataOut.nChannels
1234 1360 self.nHeights=dataOut.nHeights
1235 1361 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1236 1362 dataOut.data_spc,
1237 1363 dataOut.data_cspc,
1238 1364 dataOut.data_dc)
1239 1365 else:
1240 1366 self.nProfiles=dataOut.nProfiles
1241 1367 self.nChannels=dataOut.nChannels
1242 1368 self.nHeights=dataOut.nHeights
1243 1369 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1244 1370 dataOut.dataLag_spc,
1245 1371 dataOut.dataLag_cspc,
1246 1372 dataOut.dataLag_dc)
1247 1373
1248 1374 if self.__dataReady:
1249 1375
1250 1376 if not self.ByLags:
1251 1377
1252 1378 dataOut.data_spc = numpy.squeeze(avgdata_spc)
1253 1379 dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1254 1380 dataOut.data_dc = avgdata_dc
1255 1381 else:
1256 1382 dataOut.dataLag_spc = avgdata_spc
1257 1383 dataOut.dataLag_cspc = avgdata_cspc
1258 1384 dataOut.dataLag_dc = avgdata_dc
1259 1385
1260 1386 dataOut.data_spc=dataOut.dataLag_spc[:,:,:,dataOut.LagPlot]
1261 1387 dataOut.data_cspc=dataOut.dataLag_cspc[:,:,:,dataOut.LagPlot]
1262 1388 dataOut.data_dc=dataOut.dataLag_dc[:,:,dataOut.LagPlot]
1263 1389
1264 1390
1265 1391 dataOut.nIncohInt *= self.n
1266 1392 dataOut.utctime = avgdatatime
1267 1393 dataOut.flagNoData = False
1268 1394
1269 1395 return dataOut
1270 1396
1271 1397 class removeInterference(Operation):
1272 1398
1273 1399 def removeInterference2(self):
1274 1400
1275 1401 cspc = self.dataOut.data_cspc
1276 1402 spc = self.dataOut.data_spc
1277 1403 Heights = numpy.arange(cspc.shape[2])
1278 1404 realCspc = numpy.abs(cspc)
1279 1405
1280 1406 for i in range(cspc.shape[0]):
1281 1407 LinePower= numpy.sum(realCspc[i], axis=0)
1282 1408 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1283 1409 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1284 1410 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1285 1411 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1286 1412 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1287 1413
1288 1414
1289 1415 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1290 1416 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1291 1417 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1292 1418 cspc[i,InterferenceRange,:] = numpy.NaN
1293 1419
1294 1420 self.dataOut.data_cspc = cspc
1295 1421
1296 1422 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1297 1423
1298 1424 jspectra = self.dataOut.data_spc
1299 1425 jcspectra = self.dataOut.data_cspc
1300 1426 jnoise = self.dataOut.getNoise()
1301 1427 num_incoh = self.dataOut.nIncohInt
1302 1428
1303 1429 num_channel = jspectra.shape[0]
1304 1430 num_prof = jspectra.shape[1]
1305 1431 num_hei = jspectra.shape[2]
1306 1432
1307 1433 # hei_interf
1308 1434 if hei_interf is None:
1309 1435 count_hei = int(num_hei / 2)
1310 1436 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1311 1437 hei_interf = numpy.asarray(hei_interf)[0]
1312 1438 # nhei_interf
1313 1439 if (nhei_interf == None):
1314 1440 nhei_interf = 5
1315 1441 if (nhei_interf < 1):
1316 1442 nhei_interf = 1
1317 1443 if (nhei_interf > count_hei):
1318 1444 nhei_interf = count_hei
1319 1445 if (offhei_interf == None):
1320 1446 offhei_interf = 0
1321 1447
1322 1448 ind_hei = list(range(num_hei))
1323 1449 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1324 1450 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1325 1451 mask_prof = numpy.asarray(list(range(num_prof)))
1326 1452 num_mask_prof = mask_prof.size
1327 1453 comp_mask_prof = [0, num_prof / 2]
1328 1454
1329 1455 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1330 1456 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1331 1457 jnoise = numpy.nan
1332 1458 noise_exist = jnoise[0] < numpy.Inf
1333 1459
1334 1460 # Subrutina de Remocion de la Interferencia
1335 1461 for ich in range(num_channel):
1336 1462 # Se ordena los espectros segun su potencia (menor a mayor)
1337 1463 power = jspectra[ich, mask_prof, :]
1338 1464 power = power[:, hei_interf]
1339 1465 power = power.sum(axis=0)
1340 1466 psort = power.ravel().argsort()
1341 1467
1342 1468 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1343 1469 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1344 1470 offhei_interf, nhei_interf + offhei_interf))]]]
1345 1471
1346 1472 if noise_exist:
1347 1473 # tmp_noise = jnoise[ich] / num_prof
1348 1474 tmp_noise = jnoise[ich]
1349 1475 junkspc_interf = junkspc_interf - tmp_noise
1350 1476 #junkspc_interf[:,comp_mask_prof] = 0
1351 1477
1352 1478 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1353 1479 jspc_interf = jspc_interf.transpose()
1354 1480 # Calculando el espectro de interferencia promedio
1355 1481 noiseid = numpy.where(
1356 1482 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1357 1483 noiseid = noiseid[0]
1358 1484 cnoiseid = noiseid.size
1359 1485 interfid = numpy.where(
1360 1486 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1361 1487 interfid = interfid[0]
1362 1488 cinterfid = interfid.size
1363 1489
1364 1490 if (cnoiseid > 0):
1365 1491 jspc_interf[noiseid] = 0
1366 1492
1367 1493 # Expandiendo los perfiles a limpiar
1368 1494 if (cinterfid > 0):
1369 1495 new_interfid = (
1370 1496 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1371 1497 new_interfid = numpy.asarray(new_interfid)
1372 1498 new_interfid = {x for x in new_interfid}
1373 1499 new_interfid = numpy.array(list(new_interfid))
1374 1500 new_cinterfid = new_interfid.size
1375 1501 else:
1376 1502 new_cinterfid = 0
1377 1503
1378 1504 for ip in range(new_cinterfid):
1379 1505 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1380 1506 jspc_interf[new_interfid[ip]
1381 1507 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1382 1508
1383 1509 jspectra[ich, :, ind_hei] = jspectra[ich, :,
1384 1510 ind_hei] - jspc_interf # Corregir indices
1385 1511
1386 1512 # Removiendo la interferencia del punto de mayor interferencia
1387 1513 ListAux = jspc_interf[mask_prof].tolist()
1388 1514 maxid = ListAux.index(max(ListAux))
1389 1515
1390 1516 if cinterfid > 0:
1391 1517 for ip in range(cinterfid * (interf == 2) - 1):
1392 1518 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1393 1519 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1394 1520 cind = len(ind)
1395 1521
1396 1522 if (cind > 0):
1397 1523 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1398 1524 (1 + (numpy.random.uniform(cind) - 0.5) /
1399 1525 numpy.sqrt(num_incoh))
1400 1526
1401 1527 ind = numpy.array([-2, -1, 1, 2])
1402 1528 xx = numpy.zeros([4, 4])
1403 1529
1404 1530 for id1 in range(4):
1405 1531 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1406 1532
1407 1533 xx_inv = numpy.linalg.inv(xx)
1408 1534 xx = xx_inv[:, 0]
1409 1535 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1410 1536 yy = jspectra[ich, mask_prof[ind], :]
1411 1537 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
1412 1538 yy.transpose(), xx)
1413 1539
1414 1540 indAux = (jspectra[ich, :, :] < tmp_noise *
1415 1541 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1416 1542 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1417 1543 (1 - 1 / numpy.sqrt(num_incoh))
1418 1544
1419 1545 # Remocion de Interferencia en el Cross Spectra
1420 1546 if jcspectra is None:
1421 1547 return jspectra, jcspectra
1422 1548 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1423 1549 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1424 1550
1425 1551 for ip in range(num_pairs):
1426 1552
1427 1553 #-------------------------------------------
1428 1554
1429 1555 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1430 1556 cspower = cspower[:, hei_interf]
1431 1557 cspower = cspower.sum(axis=0)
1432 1558
1433 1559 cspsort = cspower.ravel().argsort()
1434 1560 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1435 1561 offhei_interf, nhei_interf + offhei_interf))]]]
1436 1562 junkcspc_interf = junkcspc_interf.transpose()
1437 1563 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1438 1564
1439 1565 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1440 1566
1441 1567 median_real = int(numpy.median(numpy.real(
1442 1568 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1443 1569 median_imag = int(numpy.median(numpy.imag(
1444 1570 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1445 1571 comp_mask_prof = [int(e) for e in comp_mask_prof]
1446 1572 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1447 1573 median_real, median_imag)
1448 1574
1449 1575 for iprof in range(num_prof):
1450 1576 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1451 1577 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1452 1578
1453 1579 # Removiendo la Interferencia
1454 1580 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1455 1581 :, ind_hei] - jcspc_interf
1456 1582
1457 1583 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1458 1584 maxid = ListAux.index(max(ListAux))
1459 1585
1460 1586 ind = numpy.array([-2, -1, 1, 2])
1461 1587 xx = numpy.zeros([4, 4])
1462 1588
1463 1589 for id1 in range(4):
1464 1590 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1465 1591
1466 1592 xx_inv = numpy.linalg.inv(xx)
1467 1593 xx = xx_inv[:, 0]
1468 1594
1469 1595 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1470 1596 yy = jcspectra[ip, mask_prof[ind], :]
1471 1597 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1472 1598
1473 1599 # Guardar Resultados
1474 1600 self.dataOut.data_spc = jspectra
1475 1601 self.dataOut.data_cspc = jcspectra
1476 1602
1477 1603 return 1
1478 1604
1479 1605 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
1480 1606
1481 1607 self.dataOut = dataOut
1482 1608
1483 1609 if mode == 1:
1484 1610 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1485 1611 elif mode == 2:
1486 1612 self.removeInterference2()
1487 1613
1488 1614 return self.dataOut
1489 1615
1490 1616
1491 1617 class IncohInt(Operation):
1492 1618
1493 1619 __profIndex = 0
1494 1620 __withOverapping = False
1495 1621
1496 1622 __byTime = False
1497 1623 __initime = None
1498 1624 __lastdatatime = None
1499 1625 __integrationtime = None
1500 1626
1501 1627 __buffer_spc = None
1502 1628 __buffer_cspc = None
1503 1629 __buffer_dc = None
1504 1630
1505 1631 __dataReady = False
1506 1632
1507 1633 __timeInterval = None
1508 1634
1509 1635 n = None
1510 1636
1511 1637 def __init__(self):
1512 1638
1513 1639 Operation.__init__(self)
1514 1640
1515 1641 def setup(self, n=None, timeInterval=None, overlapping=False):
1516 1642 """
1517 1643 Set the parameters of the integration class.
1518 1644
1519 1645 Inputs:
1520 1646
1521 1647 n : Number of coherent integrations
1522 1648 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1523 1649 overlapping :
1524 1650
1525 1651 """
1526 1652
1527 1653 self.__initime = None
1528 1654 self.__lastdatatime = 0
1529 1655
1530 1656 self.__buffer_spc = 0
1531 1657 self.__buffer_cspc = 0
1532 1658 self.__buffer_dc = 0
1533 1659
1534 1660 self.__profIndex = 0
1535 1661 self.__dataReady = False
1536 1662 self.__byTime = False
1537 1663
1538 1664 if n is None and timeInterval is None:
1539 1665 raise ValueError("n or timeInterval should be specified ...")
1540 1666
1541 1667 if n is not None:
1542 1668 self.n = int(n)
1543 1669 else:
1544 1670
1545 1671 self.__integrationtime = int(timeInterval)
1546 1672 self.n = None
1547 1673 self.__byTime = True
1548 1674
1549 1675 def putData(self, data_spc, data_cspc, data_dc):
1550 1676 """
1551 1677 Add a profile to the __buffer_spc and increase in one the __profileIndex
1552 1678
1553 1679 """
1554 1680
1555 1681 self.__buffer_spc += data_spc
1556 1682
1557 1683 if data_cspc is None:
1558 1684 self.__buffer_cspc = None
1559 1685 else:
1560 1686 self.__buffer_cspc += data_cspc
1561 1687
1562 1688 if data_dc is None:
1563 1689 self.__buffer_dc = None
1564 1690 else:
1565 1691 self.__buffer_dc += data_dc
1566 1692
1567 1693 self.__profIndex += 1
1568 1694
1569 1695 return
1570 1696
1571 1697 def pushData(self):
1572 1698 """
1573 1699 Return the sum of the last profiles and the profiles used in the sum.
1574 1700
1575 1701 Affected:
1576 1702
1577 1703 self.__profileIndex
1578 1704
1579 1705 """
1580 1706
1581 1707 data_spc = self.__buffer_spc
1582 1708 data_cspc = self.__buffer_cspc
1583 1709 data_dc = self.__buffer_dc
1584 1710 n = self.__profIndex
1585 1711
1586 1712 self.__buffer_spc = 0
1587 1713 self.__buffer_cspc = 0
1588 1714 self.__buffer_dc = 0
1589 1715 self.__profIndex = 0
1590 1716
1591 1717 return data_spc, data_cspc, data_dc, n
1592 1718
1593 1719 def byProfiles(self, *args):
1594 1720
1595 1721 self.__dataReady = False
1596 1722 avgdata_spc = None
1597 1723 avgdata_cspc = None
1598 1724 avgdata_dc = None
1599 1725
1600 1726 self.putData(*args)
1601 1727
1602 1728 if self.__profIndex == self.n:
1603 1729
1604 1730 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1605 1731 self.n = n
1606 1732 self.__dataReady = True
1607 1733
1608 1734 return avgdata_spc, avgdata_cspc, avgdata_dc
1609 1735
1610 1736 def byTime(self, datatime, *args):
1611 1737
1612 1738 self.__dataReady = False
1613 1739 avgdata_spc = None
1614 1740 avgdata_cspc = None
1615 1741 avgdata_dc = None
1616 1742
1617 1743 self.putData(*args)
1618 1744
1619 1745 if (datatime - self.__initime) >= self.__integrationtime:
1620 1746 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1621 1747 self.n = n
1622 1748 self.__dataReady = True
1623 1749
1624 1750 return avgdata_spc, avgdata_cspc, avgdata_dc
1625 1751
1626 1752 def integrate(self, datatime, *args):
1627 1753
1628 1754 if self.__profIndex == 0:
1629 1755 self.__initime = datatime
1630 1756
1631 1757 if self.__byTime:
1632 1758 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1633 1759 datatime, *args)
1634 1760 else:
1635 1761 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1636 1762
1637 1763 if not self.__dataReady:
1638 1764 return None, None, None, None
1639 1765
1640 1766 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1641 1767
1642 1768 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1643 1769 if n == 1:
1644 1770 return dataOut
1645 1771
1646 1772 dataOut.flagNoData = True
1647 1773
1648 1774 if not self.isConfig:
1649 1775 self.setup(n, timeInterval, overlapping)
1650 1776 self.isConfig = True
1651 1777
1652 1778 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1653 1779 dataOut.data_spc,
1654 1780 dataOut.data_cspc,
1655 1781 dataOut.data_dc)
1656 1782
1657 1783 if self.__dataReady:
1658 1784
1659 1785 dataOut.data_spc = avgdata_spc
1660 1786 dataOut.data_cspc = avgdata_cspc
1661 1787 dataOut.data_dc = avgdata_dc
1662 1788 dataOut.nIncohInt *= self.n
1663 1789 dataOut.utctime = avgdatatime
1664 1790 dataOut.flagNoData = False
1665 1791
1666 1792 return dataOut
1667 1793
1668 1794 class dopplerFlip(Operation):
1669 1795
1670 1796 def run(self, dataOut):
1671 1797 # arreglo 1: (num_chan, num_profiles, num_heights)
1672 1798 self.dataOut = dataOut
1673 1799 # JULIA-oblicua, indice 2
1674 1800 # arreglo 2: (num_profiles, num_heights)
1675 1801 jspectra = self.dataOut.data_spc[2]
1676 1802 jspectra_tmp = numpy.zeros(jspectra.shape)
1677 1803 num_profiles = jspectra.shape[0]
1678 1804 freq_dc = int(num_profiles / 2)
1679 1805 # Flip con for
1680 1806 for j in range(num_profiles):
1681 1807 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1682 1808 # Intercambio perfil de DC con perfil inmediato anterior
1683 1809 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1684 1810 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1685 1811 # canal modificado es re-escrito en el arreglo de canales
1686 1812 self.dataOut.data_spc[2] = jspectra_tmp
1687 1813
1688 1814 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now