##// END OF EJS Templates
schain
RSS
update
rflores -
r1727:66f67e6e4b37
parent child
Show More
Side by Side Unified
Context file:
r1727:66f67e6e4b37 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,1348 +1,1349
1 1 # Copyright (c) 2012-2021 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 #import collections.abc
12 12
13 13 from schainpy.model.graphics.jroplot_base import Plot, plt, log
14 14
15 15 class SpectraPlot(Plot):
16 16 '''
17 17 Plot for Spectra data
18 18 '''
19 19
20 20 CODE = 'spc'
21 21 colormap = 'jet'
22 22 plot_type = 'pcolor'
23 23 buffering = False
24 24
25 25 def setup(self):
26 26
27 27 self.nplots = len(self.data.channels)
28 28 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
29 29 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
30 30 self.height = 2.6 * self.nrows
31 31 self.cb_label = 'dB'
32 32 if self.showprofile:
33 33 self.width = 4 * self.ncols
34 34 else:
35 35 self.width = 3.5 * self.ncols
36 36 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
37 37 self.ylabel = 'Range [km]'
38 38
39 39 def update(self, dataOut):
40 40
41 41 data = {}
42 42 meta = {}
43 43
44 44 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
45 45 #print("dataOut.normFactor: ", dataOut.normFactor)
46 46 #print("spc: ", dataOut.data_spc[0,0,0])
47 47 #spc = 10*numpy.log10(dataOut.data_spc)
48 48 #print("Spc: ",spc[0])
49 49 #exit(1)
50 50 data['spc'] = spc
51 51 data['rti'] = dataOut.getPower()
52 52 #print(data['rti'][0])
53 53 #exit(1)
54 54 #print("NormFactor: ",dataOut.normFactor)
55 55 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
56 56 if hasattr(dataOut, 'LagPlot'): #Double Pulse
57 57 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
58 58 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=46,ymax_index=max_hei_id)/dataOut.normFactor)
59 59 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=40,ymax_index=max_hei_id)/dataOut.normFactor)
60 60 data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)/dataOut.normFactor)
61 61 data['noise'][0] = 10*numpy.log10(dataOut.getNoise(ymin_index=53)[0]/dataOut.normFactor)
62 62 #data['noise'][1] = 22.035507
63 63 else:
64 64 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
65 65 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=26,ymax_index=44)/dataOut.normFactor)
66 66 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
67 67
68 68 if self.CODE == 'spc_moments':
69 69 data['moments'] = dataOut.moments
70 70 if self.CODE == 'gaussian_fit':
71 71 data['gaussfit'] = dataOut.DGauFitParams
72 72
73 73 return data, meta
74 74
75 75 def plot(self):
76 76
77 77 if self.xaxis == "frequency":
78 78 x = self.data.xrange[0]
79 79 self.xlabel = "Frequency (kHz)"
80 80 elif self.xaxis == "time":
81 81 x = self.data.xrange[1]
82 82 self.xlabel = "Time (ms)"
83 83 else:
84 84 x = self.data.xrange[2]
85 85 self.xlabel = "Velocity (m/s)"
86 86
87 87 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
88 88 x = self.data.xrange[2]
89 89 self.xlabel = "Velocity (m/s)"
90 90
91 91 self.titles = []
92 92
93 93 y = self.data.yrange
94 94 self.y = y
95 95
96 96 data = self.data[-1]
97 97 z = data['spc']
98 98
99 99 self.CODE2 = 'spc_oblique'
100 100
101 101 for n, ax in enumerate(self.axes):
102 102 noise = data['noise'][n]
103 103 if self.CODE == 'spc_moments':
104 104 mean = data['moments'][n, 1]
105 105 if self.CODE == 'gaussian_fit':
106 106 gau0 = data['gaussfit'][n][2,:,0]
107 107 gau1 = data['gaussfit'][n][2,:,1]
108 108 if ax.firsttime:
109 109 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
110 110 self.xmin = self.xmin if self.xmin else numpy.nanmin(x)#-self.xmax
111 111 #self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
112 112 #self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
113 113 if self.zlimits is not None:
114 114 self.zmin, self.zmax = self.zlimits[n]
115 115
116 116 ax.plt = ax.pcolormesh(x, y, z[n].T,
117 117 vmin=self.zmin,
118 118 vmax=self.zmax,
119 119 cmap=plt.get_cmap(self.colormap),
120 120 )
121 121
122 122 if self.showprofile:
123 123 ax.plt_profile = self.pf_axes[n].plot(
124 124 data['rti'][n], y)[0]
125 125 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
126 126 color="k", linestyle="dashed", lw=1)[0]
127 127 if self.CODE == 'spc_moments':
128 128 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
129 129 if self.CODE == 'gaussian_fit':
130 130 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
131 131 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
132 132 else:
133 133 if self.zlimits is not None:
134 134 self.zmin, self.zmax = self.zlimits[n]
135 135 ax.plt.set_array(z[n].T.ravel())
136 136 if self.showprofile:
137 137 ax.plt_profile.set_data(data['rti'][n], y)
138 138 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
139 139 if self.CODE == 'spc_moments':
140 140 ax.plt_mean.set_data(mean, y)
141 141 if self.CODE == 'gaussian_fit':
142 142 ax.plt_gau0.set_data(gau0, y)
143 143 ax.plt_gau1.set_data(gau1, y)
144 144 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
145 145
146 146 class SpectraObliquePlot(Plot):
147 147 '''
148 148 Plot for Spectra data
149 149 '''
150 150
151 151 CODE = 'spc_oblique'
152 152 colormap = 'jet'
153 153 plot_type = 'pcolor'
154 154
155 155 def setup(self):
156 156 self.xaxis = "oblique"
157 157 self.nplots = len(self.data.channels)
158 158 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
159 159 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
160 160 self.height = 2.6 * self.nrows
161 161 self.cb_label = 'dB'
162 162 if self.showprofile:
163 163 self.width = 4 * self.ncols
164 164 else:
165 165 self.width = 3.5 * self.ncols
166 166 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
167 167 self.ylabel = 'Range [km]'
168 168
169 169 def update(self, dataOut):
170 170
171 171 data = {}
172 172 meta = {}
173 173
174 174 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
175 175 data['spc'] = spc
176 176 data['rti'] = dataOut.getPower()
177 177 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
178 178 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
179 179 '''
180 180 data['shift1'] = dataOut.Oblique_params[0,-2,:]
181 181 data['shift2'] = dataOut.Oblique_params[0,-1,:]
182 182 data['shift1_error'] = dataOut.Oblique_param_errors[0,-2,:]
183 183 data['shift2_error'] = dataOut.Oblique_param_errors[0,-1,:]
184 184 '''
185 185 '''
186 186 data['shift1'] = dataOut.Oblique_params[0,1,:]
187 187 data['shift2'] = dataOut.Oblique_params[0,4,:]
188 188 data['shift1_error'] = dataOut.Oblique_param_errors[0,1,:]
189 189 data['shift2_error'] = dataOut.Oblique_param_errors[0,4,:]
190 190 '''
191 191 data['shift1'] = dataOut.Dop_EEJ_T1[0]
192 192 data['shift2'] = dataOut.Dop_EEJ_T2[0]
193 193 data['max_val_2'] = dataOut.Oblique_params[0,-1,:]
194 194 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
195 195 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
196 196
197 197 return data, meta
198 198
199 199 def plot(self):
200 200
201 201 if self.xaxis == "frequency":
202 202 x = self.data.xrange[0]
203 203 self.xlabel = "Frequency (kHz)"
204 204 elif self.xaxis == "time":
205 205 x = self.data.xrange[1]
206 206 self.xlabel = "Time (ms)"
207 207 else:
208 208 x = self.data.xrange[2]
209 209 self.xlabel = "Velocity (m/s)"
210 210
211 211 self.titles = []
212 212
213 213 y = self.data.yrange
214 214 self.y = y
215 215
216 216 data = self.data[-1]
217 217 z = data['spc']
218 218
219 219 for n, ax in enumerate(self.axes):
220 220 noise = self.data['noise'][n][-1]
221 221 shift1 = data['shift1']
222 222 #print(shift1)
223 223 shift2 = data['shift2']
224 224 max_val_2 = data['max_val_2']
225 225 err1 = data['shift1_error']
226 226 err2 = data['shift2_error']
227 227 if ax.firsttime:
228 228
229 229 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
230 230 self.xmin = self.xmin if self.xmin else -self.xmax
231 231 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
232 232 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
233 233 ax.plt = ax.pcolormesh(x, y, z[n].T,
234 234 vmin=self.zmin,
235 235 vmax=self.zmax,
236 236 cmap=plt.get_cmap(self.colormap)
237 237 )
238 238
239 239 if self.showprofile:
240 240 ax.plt_profile = self.pf_axes[n].plot(
241 241 self.data['rti'][n][-1], y)[0]
242 242 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
243 243 color="k", linestyle="dashed", lw=1)[0]
244 244
245 245 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
246 246 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
247 247 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
248 248
249 249 #print("plotter1: ", self.ploterr1,shift1)
250 250
251 251 else:
252 252 #print("else plotter1: ", self.ploterr1,shift1)
253 253 self.ploterr1.remove()
254 254 self.ploterr2.remove()
255 255 self.ploterr3.remove()
256 256 ax.plt.set_array(z[n].T.ravel())
257 257 if self.showprofile:
258 258 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
259 259 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
260 260 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
261 261 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
262 262 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
263 263
264 264 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
265 265
266 266
267 267 class CrossSpectraPlot(Plot):
268 268
269 269 CODE = 'cspc'
270 270 colormap = 'jet'
271 271 plot_type = 'pcolor'
272 272 zmin_coh = None
273 273 zmax_coh = None
274 274 zmin_phase = None
275 275 zmax_phase = None
276 276
277 277 def setup(self):
278 278
279 279 self.ncols = 4
280 280 self.nplots = len(self.data.pairs) * 2
281 281 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
282 282 self.width = 3.1 * self.ncols
283 283 self.height = 5 * self.nrows
284 284 self.ylabel = 'Range [km]'
285 285 self.showprofile = False
286 286 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
287 287
288 288 def update(self, dataOut):
289 289
290 290 data = {}
291 291 meta = {}
292 292
293 293 spc = dataOut.data_spc
294 294 cspc = dataOut.data_cspc
295 295 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
296 296 meta['pairs'] = dataOut.pairsList
297 297
298 298 tmp = []
299 299
300 300 for n, pair in enumerate(meta['pairs']):
301 301 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
302 302 coh = numpy.abs(out)
303 303 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
304 304 tmp.append(coh)
305 305 tmp.append(phase)
306 306
307 307 data['cspc'] = numpy.array(tmp)
308 308
309 309 return data, meta
310 310
311 311 def plot(self):
312 312
313 313 if self.xaxis == "frequency":
314 314 x = self.data.xrange[0]
315 315 self.xlabel = "Frequency (kHz)"
316 316 elif self.xaxis == "time":
317 317 x = self.data.xrange[1]
318 318 self.xlabel = "Time (ms)"
319 319 else:
320 320 x = self.data.xrange[2]
321 321 self.xlabel = "Velocity (m/s)"
322 322
323 323 self.titles = []
324 324
325 325 y = self.data.yrange
326 326 self.y = y
327 327
328 328 data = self.data[-1]
329 329 cspc = data['cspc']
330 330
331 331 for n in range(len(self.data.pairs)):
332 332 pair = self.data.pairs[n]
333 333 coh = cspc[n*2]
334 334 phase = cspc[n*2+1]
335 335 ax = self.axes[2 * n]
336 336 if ax.firsttime:
337 337 ax.plt = ax.pcolormesh(x, y, coh.T,
338 338 vmin=0,
339 339 vmax=1,
340 340 cmap=plt.get_cmap(self.colormap_coh)
341 341 )
342 342 else:
343 343 ax.plt.set_array(coh.T.ravel())
344 344 self.titles.append(
345 345 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
346 346
347 347 ax = self.axes[2 * n + 1]
348 348 if ax.firsttime:
349 349 ax.plt = ax.pcolormesh(x, y, phase.T,
350 350 vmin=-180,
351 351 vmax=180,
352 352 cmap=plt.get_cmap(self.colormap_phase)
353 353 )
354 354 else:
355 355 ax.plt.set_array(phase.T.ravel())
356 356 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
357 357
358 358
359 359 class CrossSpectra4Plot(Plot):
360 360
361 361 CODE = 'cspc'
362 362 colormap = 'jet'
363 363 plot_type = 'pcolor'
364 364 zmin_coh = None
365 365 zmax_coh = None
366 366 zmin_phase = None
367 367 zmax_phase = None
368 368
369 369 def setup(self):
370 370
371 371 self.ncols = 4
372 372 self.nrows = len(self.data.pairs)
373 373 self.nplots = self.nrows * 4
374 374 self.width = 3.1 * self.ncols
375 375 self.height = 5 * self.nrows
376 376 self.ylabel = 'Range [km]'
377 377 self.showprofile = False
378 378 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
379 379
380 380 def plot(self):
381 381
382 382 if self.xaxis == "frequency":
383 383 x = self.data.xrange[0]
384 384 self.xlabel = "Frequency (kHz)"
385 385 elif self.xaxis == "time":
386 386 x = self.data.xrange[1]
387 387 self.xlabel = "Time (ms)"
388 388 else:
389 389 x = self.data.xrange[2]
390 390 self.xlabel = "Velocity (m/s)"
391 391
392 392 self.titles = []
393 393
394 394
395 395 y = self.data.heights
396 396 self.y = y
397 397 nspc = self.data['spc']
398 398 #print(numpy.shape(self.data['spc']))
399 399 spc = self.data['cspc'][0]
400 400 #print(numpy.shape(nspc))
401 401 #exit()
402 402 #nspc[1,:,:] = numpy.flip(nspc[1,:,:],axis=0)
403 403 #print(numpy.shape(spc))
404 404 #exit()
405 405 cspc = self.data['cspc'][1]
406 406
407 407 #xflip=numpy.flip(x)
408 408 #print(numpy.shape(cspc))
409 409 #exit()
410 410
411 411 for n in range(self.nrows):
412 412 noise = self.data['noise'][:,-1]
413 413 pair = self.data.pairs[n]
414 414 #print(pair)
415 415 #exit()
416 416 ax = self.axes[4 * n]
417 417 if ax.firsttime:
418 418 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
419 419 self.xmin = self.xmin if self.xmin else -self.xmax
420 420 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
421 421 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
422 422 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
423 423 vmin=self.zmin,
424 424 vmax=self.zmax,
425 425 cmap=plt.get_cmap(self.colormap)
426 426 )
427 427 else:
428 428 #print(numpy.shape(nspc[pair[0]].T))
429 429 #exit()
430 430 ax.plt.set_array(nspc[pair[0]].T.ravel())
431 431 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
432 432
433 433 ax = self.axes[4 * n + 1]
434 434
435 435 if ax.firsttime:
436 436 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
437 437 vmin=self.zmin,
438 438 vmax=self.zmax,
439 439 cmap=plt.get_cmap(self.colormap)
440 440 )
441 441 else:
442 442
443 443 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
444 444 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
445 445
446 446 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
447 447 coh = numpy.abs(out)
448 448 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
449 449
450 450 ax = self.axes[4 * n + 2]
451 451 if ax.firsttime:
452 452 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
453 453 vmin=0,
454 454 vmax=1,
455 455 cmap=plt.get_cmap(self.colormap_coh)
456 456 )
457 457 else:
458 458 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
459 459 self.titles.append(
460 460 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
461 461
462 462 ax = self.axes[4 * n + 3]
463 463 if ax.firsttime:
464 464 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
465 465 vmin=-180,
466 466 vmax=180,
467 467 cmap=plt.get_cmap(self.colormap_phase)
468 468 )
469 469 else:
470 470 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
471 471 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
472 472
473 473
474 474 class CrossSpectra2Plot(Plot):
475 475
476 476 CODE = 'cspc'
477 477 colormap = 'jet'
478 478 plot_type = 'pcolor'
479 479 zmin_coh = None
480 480 zmax_coh = None
481 481 zmin_phase = None
482 482 zmax_phase = None
483 483
484 484 def setup(self):
485 485
486 486 self.ncols = 1
487 487 self.nrows = len(self.data.pairs)
488 488 self.nplots = self.nrows * 1
489 489 self.width = 3.1 * self.ncols
490 490 self.height = 5 * self.nrows
491 491 self.ylabel = 'Range [km]'
492 492 self.showprofile = False
493 493 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
494 494
495 495 def plot(self):
496 496
497 497 if self.xaxis == "frequency":
498 498 x = self.data.xrange[0]
499 499 self.xlabel = "Frequency (kHz)"
500 500 elif self.xaxis == "time":
501 501 x = self.data.xrange[1]
502 502 self.xlabel = "Time (ms)"
503 503 else:
504 504 x = self.data.xrange[2]
505 505 self.xlabel = "Velocity (m/s)"
506 506
507 507 self.titles = []
508 508
509 509
510 510 y = self.data.heights
511 511 self.y = y
512 512 #nspc = self.data['spc']
513 513 #print(numpy.shape(self.data['spc']))
514 514 #spc = self.data['cspc'][0]
515 515 #print(numpy.shape(spc))
516 516 #exit()
517 517 cspc = self.data['cspc'][1]
518 518 #print(numpy.shape(cspc))
519 519 #exit()
520 520
521 521 for n in range(self.nrows):
522 522 noise = self.data['noise'][:,-1]
523 523 pair = self.data.pairs[n]
524 524 #print(pair) #exit()
525 525
526 526
527 527
528 528 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
529 529
530 530 #print(out[:,53])
531 531 #exit()
532 532 cross = numpy.abs(out)
533 533 z = cross/self.data.nFactor
534 534 #print("here")
535 535 #print(dataOut.data_spc[0,0,0])
536 536 #exit()
537 537
538 538 cross = 10*numpy.log10(z)
539 539 #print(numpy.shape(cross))
540 540 #print(cross[0,:])
541 541 #print(self.data.nFactor)
542 542 #exit()
543 543 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
544 544
545 545 ax = self.axes[1 * n]
546 546 if ax.firsttime:
547 547 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
548 548 self.xmin = self.xmin if self.xmin else -self.xmax
549 549 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
550 550 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
551 551 ax.plt = ax.pcolormesh(x, y, cross.T,
552 552 vmin=self.zmin,
553 553 vmax=self.zmax,
554 554 cmap=plt.get_cmap(self.colormap)
555 555 )
556 556 else:
557 557 ax.plt.set_array(cross.T.ravel())
558 558 self.titles.append(
559 559 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
560 560
561 561
562 562 class CrossSpectra3Plot(Plot):
563 563
564 564 CODE = 'cspc'
565 565 colormap = 'jet'
566 566 plot_type = 'pcolor'
567 567 zmin_coh = None
568 568 zmax_coh = None
569 569 zmin_phase = None
570 570 zmax_phase = None
571 571
572 572 def setup(self):
573 573
574 574 self.ncols = 3
575 575 self.nrows = len(self.data.pairs)
576 576 self.nplots = self.nrows * 3
577 577 self.width = 3.1 * self.ncols
578 578 self.height = 5 * self.nrows
579 579 self.ylabel = 'Range [km]'
580 580 self.showprofile = False
581 581 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
582 582
583 583 def plot(self):
584 584
585 585 if self.xaxis == "frequency":
586 586 x = self.data.xrange[0]
587 587 self.xlabel = "Frequency (kHz)"
588 588 elif self.xaxis == "time":
589 589 x = self.data.xrange[1]
590 590 self.xlabel = "Time (ms)"
591 591 else:
592 592 x = self.data.xrange[2]
593 593 self.xlabel = "Velocity (m/s)"
594 594
595 595 self.titles = []
596 596
597 597
598 598 y = self.data.heights
599 599 self.y = y
600 600 #nspc = self.data['spc']
601 601 #print(numpy.shape(self.data['spc']))
602 602 #spc = self.data['cspc'][0]
603 603 #print(numpy.shape(spc))
604 604 #exit()
605 605 cspc = self.data['cspc'][1]
606 606 #print(numpy.shape(cspc))
607 607 #exit()
608 608
609 609 for n in range(self.nrows):
610 610 noise = self.data['noise'][:,-1]
611 611 pair = self.data.pairs[n]
612 612 #print(pair) #exit()
613 613
614 614
615 615
616 616 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
617 617
618 618 #print(out[:,53])
619 619 #exit()
620 620 cross = numpy.abs(out)
621 621 z = cross/self.data.nFactor
622 622 cross = 10*numpy.log10(z)
623 623
624 624 out_r= out.real/self.data.nFactor
625 625 #out_r = 10*numpy.log10(out_r)
626 626
627 627 out_i= out.imag/self.data.nFactor
628 628 #out_i = 10*numpy.log10(out_i)
629 629 #print(numpy.shape(cross))
630 630 #print(cross[0,:])
631 631 #print(self.data.nFactor)
632 632 #exit()
633 633 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
634 634
635 635 ax = self.axes[3 * n]
636 636 if ax.firsttime:
637 637 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
638 638 self.xmin = self.xmin if self.xmin else -self.xmax
639 639 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
640 640 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
641 641 ax.plt = ax.pcolormesh(x, y, cross.T,
642 642 vmin=self.zmin,
643 643 vmax=self.zmax,
644 644 cmap=plt.get_cmap(self.colormap)
645 645 )
646 646 else:
647 647 ax.plt.set_array(cross.T.ravel())
648 648 self.titles.append(
649 649 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
650 650
651 651 ax = self.axes[3 * n + 1]
652 652 if ax.firsttime:
653 653 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
654 654 self.xmin = self.xmin if self.xmin else -self.xmax
655 655 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
656 656 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
657 657 ax.plt = ax.pcolormesh(x, y, out_r.T,
658 658 vmin=-1.e6,
659 659 vmax=0,
660 660 cmap=plt.get_cmap(self.colormap)
661 661 )
662 662 else:
663 663 ax.plt.set_array(out_r.T.ravel())
664 664 self.titles.append(
665 665 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
666 666
667 667 ax = self.axes[3 * n + 2]
668 668
669 669
670 670 if ax.firsttime:
671 671 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
672 672 self.xmin = self.xmin if self.xmin else -self.xmax
673 673 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
674 674 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
675 675 ax.plt = ax.pcolormesh(x, y, out_i.T,
676 676 vmin=-1.e6,
677 677 vmax=1.e6,
678 678 cmap=plt.get_cmap(self.colormap)
679 679 )
680 680 else:
681 681 ax.plt.set_array(out_i.T.ravel())
682 682 self.titles.append(
683 683 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
684 684
685 685 class RTIPlot(Plot):
686 686 '''
687 687 Plot for RTI data
688 688 '''
689 689
690 690 CODE = 'rti'
691 691 colormap = 'jet'
692 692 plot_type = 'pcolorbuffer'
693 693
694 694 def setup(self):
695 695 self.xaxis = 'time'
696 696 self.ncols = 1
697 697 self.nrows = len(self.data.channels)
698 698 self.nplots = len(self.data.channels)
699 699 self.ylabel = 'Range [km]'
700 700 self.xlabel = 'Time'
701 701 self.cb_label = 'dB'
702 702 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
703 703 self.titles = ['{} Channel {}'.format(
704 704 self.CODE.upper(), x) for x in range(self.nrows)]
705 705
706 706 def update(self, dataOut):
707 707
708 708 data = {}
709 709 meta = {}
710 710 data['rti'] = dataOut.getPower()
711 711 #print(numpy.shape(data['rti']))
712 712
713 713 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
714 714
715 715 return data, meta
716 716
717 717 def plot(self):
718 718
719 719 self.x = self.data.times
720 720 self.y = self.data.yrange
721 721 self.z = self.data[self.CODE]
722 722 #print("Inside RTI: ", self.z)
723 723 self.z = numpy.ma.masked_invalid(self.z)
724 724
725 725 if self.decimation is None:
726 726 x, y, z = self.fill_gaps(self.x, self.y, self.z)
727 727 else:
728 728 x, y, z = self.fill_gaps(*self.decimate())
729 729 #print("self.z: ", self.z)
730 730 #exit(1)
731 731 '''
732 732 if not isinstance(self.zmin, collections.abc.Sequence):
733 733 if not self.zmin:
734 734 self.zmin = [numpy.min(self.z)]*len(self.axes)
735 735 else:
736 736 self.zmin = [self.zmin]*len(self.axes)
737 737
738 738 if not isinstance(self.zmax, collections.abc.Sequence):
739 739 if not self.zmax:
740 740 self.zmax = [numpy.max(self.z)]*len(self.axes)
741 741 else:
742 742 self.zmax = [self.zmax]*len(self.axes)
743 743 '''
744 744 for n, ax in enumerate(self.axes):
745 745
746 746 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
747 747 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
748 748
749 749 if ax.firsttime:
750 750 if self.zlimits is not None:
751 751 self.zmin, self.zmax = self.zlimits[n]
752 752 ax.plt = ax.pcolormesh(x, y, z[n].T,
753 753 vmin=self.zmin,
754 754 vmax=self.zmax,
755 755 cmap=plt.get_cmap(self.colormap)
756 756 )
757 757 if self.showprofile:
758 758 ax.plot_profile = self.pf_axes[n].plot(
759 759 self.data['rti'][n][-1], self.y)[0]
760 760 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
761 761 color="k", linestyle="dashed", lw=1)[0]
762 762 else:
763 763 if self.zlimits is not None:
764 764 self.zmin, self.zmax = self.zlimits[n]
765 765 ax.collections.remove(ax.collections[0])
766 766 ax.plt = ax.pcolormesh(x, y, z[n].T,
767 767 vmin=self.zmin,
768 768 vmax=self.zmax,
769 769 cmap=plt.get_cmap(self.colormap)
770 770 )
771 771 if self.showprofile:
772 772 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
773 773 ax.plot_noise.set_data(numpy.repeat(
774 774 self.data['noise'][n][-1], len(self.y)), self.y)
775 775
776 776
777 777 class SpectrogramPlot(Plot):
778 778 '''
779 779 Plot for Spectrogram data
780 780 '''
781 781
782 782 CODE = 'Spectrogram_Profile'
783 783 colormap = 'binary'
784 784 plot_type = 'pcolorbuffer'
785 785
786 786 def setup(self):
787 787 self.xaxis = 'time'
788 788 self.ncols = 1
789 789 self.nrows = len(self.data.channels)
790 790 self.nplots = len(self.data.channels)
791 791 self.xlabel = 'Time'
792 792 #self.cb_label = 'dB'
793 793 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
794 794 self.titles = []
795 795
796 796 #self.titles = ['{} Channel {} \n H = {} km ({} - {})'.format(
797 797 #self.CODE.upper(), x, self.data.heightList[self.data.hei], self.data.heightList[self.data.hei],self.data.heightList[self.data.hei]+(self.data.DH*self.data.nProfiles)) for x in range(self.nrows)]
798 798
799 799 self.titles = ['{} Channel {}'.format(
800 800 self.CODE.upper(), x) for x in range(self.nrows)]
801 801
802 802
803 803 def update(self, dataOut):
804 804 data = {}
805 805 meta = {}
806 806
807 807 maxHei = 1620#+12000
808 maxHei = 1180
808 809 indb = numpy.where(dataOut.heightList <= maxHei)
809 810 hei = indb[0][-1]
810 811 #print(dataOut.heightList)
811 812
812 813 factor = dataOut.nIncohInt
813 814 z = dataOut.data_spc[:,:,hei] / factor
814 815 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
815 816 #buffer = 10 * numpy.log10(z)
816 817
817 818 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
818 819
819 820
820 821 #self.hei = hei
821 822 #self.heightList = dataOut.heightList
822 823 #self.DH = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
823 824 #self.nProfiles = dataOut.nProfiles
824 825
825 826 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
826 827
827 828 data['hei'] = hei
828 829 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
829 830 data['nProfiles'] = dataOut.nProfiles
830 831 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
831 832 '''
832 833 import matplotlib.pyplot as plt
833 834 plt.plot(10 * numpy.log10(z[0,:]))
834 835 plt.show()
835 836
836 837 from time import sleep
837 838 sleep(10)
838 839 '''
839 840 return data, meta
840 841
841 842 def plot(self):
842 843
843 844 self.x = self.data.times
844 845 self.z = self.data[self.CODE]
845 846 self.y = self.data.xrange[0]
846 847
847 848 hei = self.data['hei'][-1]
848 849 DH = self.data['DH'][-1]
849 850 nProfiles = self.data['nProfiles'][-1]
850 851
851 852 self.ylabel = "Frequency (kHz)"
852 853
853 854 self.z = numpy.ma.masked_invalid(self.z)
854 855
855 856 if self.decimation is None:
856 857 x, y, z = self.fill_gaps(self.x, self.y, self.z)
857 858 else:
858 859 x, y, z = self.fill_gaps(*self.decimate())
859 860
860 861 for n, ax in enumerate(self.axes):
861 862 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
862 863 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
863 864 data = self.data[-1]
864 865 if ax.firsttime:
865 866 ax.plt = ax.pcolormesh(x, y, z[n].T,
866 867 vmin=self.zmin,
867 868 vmax=self.zmax,
868 869 cmap=plt.get_cmap(self.colormap)
869 870 )
870 871 else:
871 872 ax.collections.remove(ax.collections[0])
872 873 ax.plt = ax.pcolormesh(x, y, z[n].T,
873 874 vmin=self.zmin,
874 875 vmax=self.zmax,
875 876 cmap=plt.get_cmap(self.colormap)
876 877 )
877 878
878 879 #self.titles.append('Spectrogram')
879 880
880 881 #self.titles.append('{} Channel {} \n H = {} km ({} - {})'.format(
881 882 #self.CODE.upper(), x, y[hei], y[hei],y[hei]+(DH*nProfiles)))
882 883
883 884
884 885
885 886
886 887 class CoherencePlot(RTIPlot):
887 888 '''
888 889 Plot for Coherence data
889 890 '''
890 891
891 892 CODE = 'coh'
892 893
893 894 def setup(self):
894 895 self.xaxis = 'time'
895 896 self.ncols = 1
896 897 self.nrows = len(self.data.pairs)
897 898 self.nplots = len(self.data.pairs)
898 899 self.ylabel = 'Range [km]'
899 900 self.xlabel = 'Time'
900 901 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
901 902 if self.CODE == 'coh':
902 903 self.cb_label = ''
903 904 self.titles = [
904 905 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
905 906 else:
906 907 self.cb_label = 'Degrees'
907 908 self.titles = [
908 909 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
909 910
910 911 def update(self, dataOut):
911 912
912 913 data = {}
913 914 meta = {}
914 915 data['coh'] = dataOut.getCoherence()
915 916 meta['pairs'] = dataOut.pairsList
916 917
917 918 return data, meta
918 919
919 920 class PhasePlot(CoherencePlot):
920 921 '''
921 922 Plot for Phase map data
922 923 '''
923 924
924 925 CODE = 'phase'
925 926 colormap = 'seismic'
926 927
927 928 def update(self, dataOut):
928 929
929 930 data = {}
930 931 meta = {}
931 932 data['phase'] = dataOut.getCoherence(phase=True)
932 933 meta['pairs'] = dataOut.pairsList
933 934
934 935 return data, meta
935 936
936 937 class NoisePlot(Plot):
937 938 '''
938 939 Plot for noise
939 940 '''
940 941
941 942 CODE = 'noise'
942 943 plot_type = 'scatterbuffer'
943 944
944 945 def setup(self):
945 946 self.xaxis = 'time'
946 947 self.ncols = 1
947 948 self.nrows = 1
948 949 self.nplots = 1
949 950 self.ylabel = 'Intensity [dB]'
950 951 self.xlabel = 'Time'
951 952 self.titles = ['Noise']
952 953 self.colorbar = False
953 954 self.plots_adjust.update({'right': 0.85 })
954 955
955 956 def update(self, dataOut):
956 957
957 958 data = {}
958 959 meta = {}
959 960 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
960 961 meta['yrange'] = numpy.array([])
961 962
962 963 return data, meta
963 964
964 965 def plot(self):
965 966
966 967 x = self.data.times
967 968 xmin = self.data.min_time
968 969 xmax = xmin + self.xrange * 60 * 60
969 970 Y = self.data['noise']
970 971
971 972 if self.axes[0].firsttime:
972 973 self.ymin = numpy.nanmin(Y) - 5
973 974 self.ymax = numpy.nanmax(Y) + 5
974 975 for ch in self.data.channels:
975 976 y = Y[ch]
976 977 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
977 978 plt.legend(bbox_to_anchor=(1.18, 1.0))
978 979 else:
979 980 for ch in self.data.channels:
980 981 y = Y[ch]
981 982 self.axes[0].lines[ch].set_data(x, y)
982 983
983 984 self.ymin = numpy.nanmin(Y) - 5
984 985 self.ymax = numpy.nanmax(Y) + 10
985 986
986 987
987 988 class PowerProfilePlot(Plot):
988 989
989 990 CODE = 'pow_profile'
990 991 plot_type = 'scatter'
991 992
992 993 def setup(self):
993 994
994 995 self.ncols = 1
995 996 self.nrows = 1
996 997 self.nplots = 1
997 998 self.height = 4
998 999 self.width = 3
999 1000 self.ylabel = 'Range [km]'
1000 1001 self.xlabel = 'Intensity [dB]'
1001 1002 self.titles = ['Power Profile']
1002 1003 self.colorbar = False
1003 1004
1004 1005 def update(self, dataOut):
1005 1006
1006 1007 data = {}
1007 1008 meta = {}
1008 1009 data[self.CODE] = dataOut.getPower()
1009 1010
1010 1011 return data, meta
1011 1012
1012 1013 def plot(self):
1013 1014
1014 1015 y = self.data.yrange
1015 1016 self.y = y
1016 1017
1017 1018 x = self.data[-1][self.CODE]
1018 1019
1019 1020 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
1020 1021 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
1021 1022
1022 1023 if self.axes[0].firsttime:
1023 1024 for ch in self.data.channels:
1024 1025 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
1025 1026 plt.legend()
1026 1027 else:
1027 1028 for ch in self.data.channels:
1028 1029 self.axes[0].lines[ch].set_data(x[ch], y)
1029 1030
1030 1031
1031 1032 class SpectraCutPlot(Plot):
1032 1033
1033 1034 CODE = 'spc_cut'
1034 1035 plot_type = 'scatter'
1035 1036 buffering = False
1036 1037
1037 1038 def setup(self):
1038 1039
1039 1040 self.nplots = len(self.data.channels)
1040 1041 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
1041 1042 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
1042 1043 self.width = 3.4 * self.ncols + 1.5
1043 1044 self.height = 3 * self.nrows
1044 1045 self.ylabel = 'Power [dB]'
1045 1046 self.colorbar = False
1046 1047 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
1047 1048
1048 1049 def update(self, dataOut):
1049 1050
1050 1051 data = {}
1051 1052 meta = {}
1052 1053 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
1053 1054 data['spc'] = spc
1054 1055 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
1055 1056 if self.CODE == 'cut_gaussian_fit':
1056 1057 data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
1057 1058 data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
1058 1059 return data, meta
1059 1060
1060 1061 def plot(self):
1061 1062 if self.xaxis == "frequency":
1062 1063 x = self.data.xrange[0][1:]
1063 1064 self.xlabel = "Frequency (kHz)"
1064 1065 elif self.xaxis == "time":
1065 1066 x = self.data.xrange[1]
1066 1067 self.xlabel = "Time (ms)"
1067 1068 else:
1068 1069 x = self.data.xrange[2][:-1]
1069 1070 self.xlabel = "Velocity (m/s)"
1070 1071
1071 1072 if self.CODE == 'cut_gaussian_fit':
1072 1073 x = self.data.xrange[2][:-1]
1073 1074 self.xlabel = "Velocity (m/s)"
1074 1075
1075 1076 self.titles = []
1076 1077
1077 1078 y = self.data.yrange
1078 1079 data = self.data[-1]
1079 1080 z = data['spc']
1080 1081
1081 1082 if self.height_index:
1082 1083 index = numpy.array(self.height_index)
1083 1084 else:
1084 1085 index = numpy.arange(0, len(y), int((len(y))/9))
1085 1086
1086 1087 for n, ax in enumerate(self.axes):
1087 1088 if self.CODE == 'cut_gaussian_fit':
1088 1089 gau0 = data['gauss_fit0']
1089 1090 gau1 = data['gauss_fit1']
1090 1091 if ax.firsttime:
1091 1092 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1092 1093 self.xmin = self.xmin if self.xmin else -self.xmax
1093 1094 self.ymin = self.ymin if self.ymin else numpy.nanmin(z[:,:,index])
1094 1095 self.ymax = self.ymax if self.ymax else numpy.nanmax(z[:,:,index])
1095 1096 #print(self.ymax)
1096 1097 #print(z[n, :, index])
1097 1098 ax.plt = ax.plot(x, z[n, :, index].T, lw=0.25)
1098 1099 if self.CODE == 'cut_gaussian_fit':
1099 1100 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
1100 1101 for i, line in enumerate(ax.plt_gau0):
1101 1102 line.set_color(ax.plt[i].get_color())
1102 1103 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
1103 1104 for i, line in enumerate(ax.plt_gau1):
1104 1105 line.set_color(ax.plt[i].get_color())
1105 1106 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1106 1107 self.figures[0].legend(ax.plt, labels, loc='center right')
1107 1108 else:
1108 1109 for i, line in enumerate(ax.plt):
1109 1110 line.set_data(x, z[n, :, index[i]].T)
1110 1111 for i, line in enumerate(ax.plt_gau0):
1111 1112 line.set_data(x, gau0[n, :, index[i]].T)
1112 1113 line.set_color(ax.plt[i].get_color())
1113 1114 for i, line in enumerate(ax.plt_gau1):
1114 1115 line.set_data(x, gau1[n, :, index[i]].T)
1115 1116 line.set_color(ax.plt[i].get_color())
1116 1117 self.titles.append('CH {}'.format(n))
1117 1118
1118 1119
1119 1120 class BeaconPhase(Plot):
1120 1121
1121 1122 __isConfig = None
1122 1123 __nsubplots = None
1123 1124
1124 1125 PREFIX = 'beacon_phase'
1125 1126
1126 1127 def __init__(self):
1127 1128 Plot.__init__(self)
1128 1129 self.timerange = 24*60*60
1129 1130 self.isConfig = False
1130 1131 self.__nsubplots = 1
1131 1132 self.counter_imagwr = 0
1132 1133 self.WIDTH = 800
1133 1134 self.HEIGHT = 400
1134 1135 self.WIDTHPROF = 120
1135 1136 self.HEIGHTPROF = 0
1136 1137 self.xdata = None
1137 1138 self.ydata = None
1138 1139
1139 1140 self.PLOT_CODE = BEACON_CODE
1140 1141
1141 1142 self.FTP_WEI = None
1142 1143 self.EXP_CODE = None
1143 1144 self.SUB_EXP_CODE = None
1144 1145 self.PLOT_POS = None
1145 1146
1146 1147 self.filename_phase = None
1147 1148
1148 1149 self.figfile = None
1149 1150
1150 1151 self.xmin = None
1151 1152 self.xmax = None
1152 1153
1153 1154 def getSubplots(self):
1154 1155
1155 1156 ncol = 1
1156 1157 nrow = 1
1157 1158
1158 1159 return nrow, ncol
1159 1160
1160 1161 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1161 1162
1162 1163 self.__showprofile = showprofile
1163 1164 self.nplots = nplots
1164 1165
1165 1166 ncolspan = 7
1166 1167 colspan = 6
1167 1168 self.__nsubplots = 2
1168 1169
1169 1170 self.createFigure(id = id,
1170 1171 wintitle = wintitle,
1171 1172 widthplot = self.WIDTH+self.WIDTHPROF,
1172 1173 heightplot = self.HEIGHT+self.HEIGHTPROF,
1173 1174 show=show)
1174 1175
1175 1176 nrow, ncol = self.getSubplots()
1176 1177
1177 1178 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1178 1179
1179 1180 def save_phase(self, filename_phase):
1180 1181 f = open(filename_phase,'w+')
1181 1182 f.write('\n\n')
1182 1183 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1183 1184 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1184 1185 f.close()
1185 1186
1186 1187 def save_data(self, filename_phase, data, data_datetime):
1187 1188 f=open(filename_phase,'a')
1188 1189 timetuple_data = data_datetime.timetuple()
1189 1190 day = str(timetuple_data.tm_mday)
1190 1191 month = str(timetuple_data.tm_mon)
1191 1192 year = str(timetuple_data.tm_year)
1192 1193 hour = str(timetuple_data.tm_hour)
1193 1194 minute = str(timetuple_data.tm_min)
1194 1195 second = str(timetuple_data.tm_sec)
1195 1196 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1196 1197 f.close()
1197 1198
1198 1199 def plot(self):
1199 1200 log.warning('TODO: Not yet implemented...')
1200 1201
1201 1202 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1202 1203 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1203 1204 timerange=None,
1204 1205 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1205 1206 server=None, folder=None, username=None, password=None,
1206 1207 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1207 1208
1208 1209 if dataOut.flagNoData:
1209 1210 return dataOut
1210 1211
1211 1212 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1212 1213 return
1213 1214
1214 1215 if pairsList == None:
1215 1216 pairsIndexList = dataOut.pairsIndexList[:10]
1216 1217 else:
1217 1218 pairsIndexList = []
1218 1219 for pair in pairsList:
1219 1220 if pair not in dataOut.pairsList:
1220 1221 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1221 1222 pairsIndexList.append(dataOut.pairsList.index(pair))
1222 1223
1223 1224 if pairsIndexList == []:
1224 1225 return
1225 1226
1226 1227 # if len(pairsIndexList) > 4:
1227 1228 # pairsIndexList = pairsIndexList[0:4]
1228 1229
1229 1230 hmin_index = None
1230 1231 hmax_index = None
1231 1232
1232 1233 if hmin != None and hmax != None:
1233 1234 indexes = numpy.arange(dataOut.nHeights)
1234 1235 hmin_list = indexes[dataOut.heightList >= hmin]
1235 1236 hmax_list = indexes[dataOut.heightList <= hmax]
1236 1237
1237 1238 if hmin_list.any():
1238 1239 hmin_index = hmin_list[0]
1239 1240
1240 1241 if hmax_list.any():
1241 1242 hmax_index = hmax_list[-1]+1
1242 1243
1243 1244 x = dataOut.getTimeRange()
1244 1245
1245 1246 thisDatetime = dataOut.datatime
1246 1247
1247 1248 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1248 1249 xlabel = "Local Time"
1249 1250 ylabel = "Phase (degrees)"
1250 1251
1251 1252 update_figfile = False
1252 1253
1253 1254 nplots = len(pairsIndexList)
1254 1255 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1255 1256 phase_beacon = numpy.zeros(len(pairsIndexList))
1256 1257 for i in range(nplots):
1257 1258 pair = dataOut.pairsList[pairsIndexList[i]]
1258 1259 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1259 1260 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1260 1261 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1261 1262 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1262 1263 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1263 1264
1264 1265 if dataOut.beacon_heiIndexList:
1265 1266 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1266 1267 else:
1267 1268 phase_beacon[i] = numpy.average(phase)
1268 1269
1269 1270 if not self.isConfig:
1270 1271
1271 1272 nplots = len(pairsIndexList)
1272 1273
1273 1274 self.setup(id=id,
1274 1275 nplots=nplots,
1275 1276 wintitle=wintitle,
1276 1277 showprofile=showprofile,
1277 1278 show=show)
1278 1279
1279 1280 if timerange != None:
1280 1281 self.timerange = timerange
1281 1282
1282 1283 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1283 1284
1284 1285 if ymin == None: ymin = 0
1285 1286 if ymax == None: ymax = 360
1286 1287
1287 1288 self.FTP_WEI = ftp_wei
1288 1289 self.EXP_CODE = exp_code
1289 1290 self.SUB_EXP_CODE = sub_exp_code
1290 1291 self.PLOT_POS = plot_pos
1291 1292
1292 1293 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1293 1294 self.isConfig = True
1294 1295 self.figfile = figfile
1295 1296 self.xdata = numpy.array([])
1296 1297 self.ydata = numpy.array([])
1297 1298
1298 1299 update_figfile = True
1299 1300
1300 1301 #open file beacon phase
1301 1302 path = '%s%03d' %(self.PREFIX, self.id)
1302 1303 beacon_file = os.path.join(path,'%s.txt'%self.name)
1303 1304 self.filename_phase = os.path.join(figpath,beacon_file)
1304 1305 #self.save_phase(self.filename_phase)
1305 1306
1306 1307
1307 1308 #store data beacon phase
1308 1309 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1309 1310
1310 1311 self.setWinTitle(title)
1311 1312
1312 1313
1313 1314 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1314 1315
1315 1316 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1316 1317
1317 1318 axes = self.axesList[0]
1318 1319
1319 1320 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1320 1321
1321 1322 if len(self.ydata)==0:
1322 1323 self.ydata = phase_beacon.reshape(-1,1)
1323 1324 else:
1324 1325 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1325 1326
1326 1327
1327 1328 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1328 1329 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1329 1330 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1330 1331 XAxisAsTime=True, grid='both'
1331 1332 )
1332 1333
1333 1334 self.draw()
1334 1335
1335 1336 if dataOut.ltctime >= self.xmax:
1336 1337 self.counter_imagwr = wr_period
1337 1338 self.isConfig = False
1338 1339 update_figfile = True
1339 1340
1340 1341 self.save(figpath=figpath,
1341 1342 figfile=figfile,
1342 1343 save=save,
1343 1344 ftp=ftp,
1344 1345 wr_period=wr_period,
1345 1346 thisDatetime=thisDatetime,
1346 1347 update_figfile=update_figfile)
1347 1348
1348 1349 return dataOut
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,1049 +1,1049
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15
16 16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 17 from schainpy.model.data.jrodata import Spectra
18 18 from schainpy.model.data.jrodata import hildebrand_sekhon
19 19 from schainpy.utils import log
20 20
21 21
22 22 class SpectraProc(ProcessingUnit):
23 23
24 24 def __init__(self):
25 25
26 26 ProcessingUnit.__init__(self)
27 27
28 28 self.buffer = None
29 29 self.firstdatatime = None
30 30 self.profIndex = 0
31 31 self.dataOut = Spectra()
32 32 self.id_min = None
33 33 self.id_max = None
34 34 self.setupReq = False #Agregar a todas las unidades de proc
35 35
36 36 def __updateSpecFromVoltage(self):
37 37
38 38 self.dataOut.timeZone = self.dataIn.timeZone
39 39 self.dataOut.dstFlag = self.dataIn.dstFlag
40 40 self.dataOut.errorCount = self.dataIn.errorCount
41 41 self.dataOut.useLocalTime = self.dataIn.useLocalTime
42 42 try:
43 43 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
44 44 except:
45 45 pass
46 46 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
47 47 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
48 48 self.dataOut.channelList = self.dataIn.channelList
49 49 self.dataOut.heightList = self.dataIn.heightList
50 50 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
51 51 self.dataOut.nProfiles = self.dataOut.nFFTPoints
52 52 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
53 53 self.dataOut.utctime = self.firstdatatime
54 54 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
55 55 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
56 56 self.dataOut.flagShiftFFT = False
57 57 self.dataOut.nCohInt = self.dataIn.nCohInt
58 58 self.dataOut.nIncohInt = 1
59 59 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
60 60 self.dataOut.frequency = self.dataIn.frequency
61 61 self.dataOut.realtime = self.dataIn.realtime
62 62 self.dataOut.azimuth = self.dataIn.azimuth
63 63 self.dataOut.zenith = self.dataIn.zenith
64 64 self.dataOut.beam.codeList = self.dataIn.beam.codeList
65 65 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
66 66 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
67 67 self.dataOut.runNextUnit = self.dataIn.runNextUnit
68 68 try:
69 69 self.dataOut.step = self.dataIn.step
70 70 except:
71 71 pass
72 72
73 73 def __getFft(self):
74 74 """
75 75 Convierte valores de Voltaje a Spectra
76 76
77 77 Affected:
78 78 self.dataOut.data_spc
79 79 self.dataOut.data_cspc
80 80 self.dataOut.data_dc
81 81 self.dataOut.heightList
82 82 self.profIndex
83 83 self.buffer
84 84 self.dataOut.flagNoData
85 85 """
86 86 fft_volt = numpy.fft.fft(
87 87 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
88 88 fft_volt = fft_volt.astype(numpy.dtype('complex'))
89 89 dc = fft_volt[:, 0, :]
90 90
91 91 # calculo de self-spectra
92 92 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
93 93 spc = fft_volt * numpy.conjugate(fft_volt)
94 94 spc = spc.real
95 95
96 96 blocksize = 0
97 97 blocksize += dc.size
98 98 blocksize += spc.size
99 99
100 100 cspc = None
101 101 pairIndex = 0
102 102 if self.dataOut.pairsList != None:
103 103 # calculo de cross-spectra
104 104 cspc = numpy.zeros(
105 105 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
106 106 for pair in self.dataOut.pairsList:
107 107 if pair[0] not in self.dataOut.channelList:
108 108 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
109 109 str(pair), str(self.dataOut.channelList)))
110 110 if pair[1] not in self.dataOut.channelList:
111 111 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
112 112 str(pair), str(self.dataOut.channelList)))
113 113
114 114 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
115 115 numpy.conjugate(fft_volt[pair[1], :, :])
116 116 pairIndex += 1
117 117 blocksize += cspc.size
118 118
119 119 self.dataOut.data_spc = spc
120 120 self.dataOut.data_cspc = cspc
121 121 self.dataOut.data_dc = dc
122 122 self.dataOut.blockSize = blocksize
123 123 self.dataOut.flagShiftFFT = False
124 124
125 125 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
126 126
127 127 self.dataIn.runNextUnit = runNextUnit
128 128 if self.dataIn.type == "Spectra":
129 129
130 130 self.dataOut.copy(self.dataIn)
131 131 if shift_fft:
132 132 #desplaza a la derecha en el eje 2 determinadas posiciones
133 133 shift = int(self.dataOut.nFFTPoints/2)
134 134 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
135 135
136 136 if self.dataOut.data_cspc is not None:
137 137 #desplaza a la derecha en el eje 2 determinadas posiciones
138 138 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
139 139 if pairsList:
140 140 self.__selectPairs(pairsList)
141 141
142 142 elif self.dataIn.type == "Voltage":
143 143
144 144 self.dataOut.flagNoData = True
145 145
146 146 if nFFTPoints == None:
147 147 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
148 148
149 149 if nProfiles == None:
150 150 nProfiles = nFFTPoints
151 151 #print(self.dataOut.ipp)
152 152 #exit(1)
153 153 if ippFactor == None:
154 154 self.dataOut.ippFactor = 1
155 155 #if ippFactor is not None:
156 156 #self.dataOut.ippFactor = ippFactor
157 157 #print(ippFactor)
158 158 #print(self.dataOut.ippFactor)
159 159 #exit(1)
160 160
161 161 self.dataOut.nFFTPoints = nFFTPoints
162 162
163 163 if self.buffer is None:
164 164 self.buffer = numpy.zeros((self.dataIn.nChannels,
165 165 nProfiles,
166 166 self.dataIn.nHeights),
167 167 dtype='complex')
168 168
169 169 if self.dataIn.flagDataAsBlock:
170 170 nVoltProfiles = self.dataIn.data.shape[1]
171 171
172 172 if nVoltProfiles == nProfiles:
173 173 self.buffer = self.dataIn.data.copy()
174 174 self.profIndex = nVoltProfiles
175 175
176 176 elif nVoltProfiles < nProfiles:
177 177
178 178 if self.profIndex == 0:
179 179 self.id_min = 0
180 180 self.id_max = nVoltProfiles
181 181 #print(self.id_min)
182 182 #print(self.id_max)
183 183 #print(numpy.shape(self.buffer))
184 184 self.buffer[:, self.id_min:self.id_max,
185 185 :] = self.dataIn.data
186 186 self.profIndex += nVoltProfiles
187 187 self.id_min += nVoltProfiles
188 188 self.id_max += nVoltProfiles
189 189 else:
190 190 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
191 191 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
192 192 self.dataOut.flagNoData = True
193 193 else:
194 194 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
195 195 self.profIndex += 1
196 196
197 197 if self.firstdatatime == None:
198 198 self.firstdatatime = self.dataIn.utctime
199 199
200 200 if self.profIndex == nProfiles:
201 201 self.__updateSpecFromVoltage()
202 202 if pairsList == None:
203 203 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
204 204 else:
205 205 self.dataOut.pairsList = pairsList
206 206 self.__getFft()
207 207 self.dataOut.flagNoData = False
208 208 self.firstdatatime = None
209 209 self.profIndex = 0
210 210 else:
211 211 raise ValueError("The type of input object '%s' is not valid".format(
212 212 self.dataIn.type))
213 213
214 214
215 215 def __selectPairs(self, pairsList):
216 216
217 217 if not pairsList:
218 218 return
219 219
220 220 pairs = []
221 221 pairsIndex = []
222 222
223 223 for pair in pairsList:
224 224 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
225 225 continue
226 226 pairs.append(pair)
227 227 pairsIndex.append(pairs.index(pair))
228 228
229 229 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
230 230 self.dataOut.pairsList = pairs
231 231
232 232 return
233 233
234 234 def selectFFTs(self, minFFT, maxFFT ):
235 235 """
236 236 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
237 237 minFFT<= FFT <= maxFFT
238 238 """
239 239
240 240 if (minFFT > maxFFT):
241 241 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
242 242
243 243 if (minFFT < self.dataOut.getFreqRange()[0]):
244 244 minFFT = self.dataOut.getFreqRange()[0]
245 245
246 246 if (maxFFT > self.dataOut.getFreqRange()[-1]):
247 247 maxFFT = self.dataOut.getFreqRange()[-1]
248 248
249 249 minIndex = 0
250 250 maxIndex = 0
251 251 FFTs = self.dataOut.getFreqRange()
252 252
253 253 inda = numpy.where(FFTs >= minFFT)
254 254 indb = numpy.where(FFTs <= maxFFT)
255 255
256 256 try:
257 257 minIndex = inda[0][0]
258 258 except:
259 259 minIndex = 0
260 260
261 261 try:
262 262 maxIndex = indb[0][-1]
263 263 except:
264 264 maxIndex = len(FFTs)
265 265
266 266 self.selectFFTsByIndex(minIndex, maxIndex)
267 267
268 268 return 1
269 269
270 270 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
271 271 newheis = numpy.where(
272 272 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
273 273
274 274 if hei_ref != None:
275 275 newheis = numpy.where(self.dataOut.heightList > hei_ref)
276 276
277 277 minIndex = min(newheis[0])
278 278 maxIndex = max(newheis[0])
279 279 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
280 280 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
281 281
282 282 # determina indices
283 283 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
284 284 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
285 285 avg_dB = 10 * \
286 286 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
287 287 beacon_dB = numpy.sort(avg_dB)[-nheis:]
288 288 beacon_heiIndexList = []
289 289 for val in avg_dB.tolist():
290 290 if val >= beacon_dB[0]:
291 291 beacon_heiIndexList.append(avg_dB.tolist().index(val))
292 292
293 293 #data_spc = data_spc[:,:,beacon_heiIndexList]
294 294 data_cspc = None
295 295 if self.dataOut.data_cspc is not None:
296 296 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
297 297 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
298 298
299 299 data_dc = None
300 300 if self.dataOut.data_dc is not None:
301 301 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
302 302 #data_dc = data_dc[:,beacon_heiIndexList]
303 303
304 304 self.dataOut.data_spc = data_spc
305 305 self.dataOut.data_cspc = data_cspc
306 306 self.dataOut.data_dc = data_dc
307 307 self.dataOut.heightList = heightList
308 308 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
309 309
310 310 return 1
311 311
312 312 def selectFFTsByIndex(self, minIndex, maxIndex):
313 313 """
314 314
315 315 """
316 316
317 317 if (minIndex < 0) or (minIndex > maxIndex):
318 318 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
319 319
320 320 if (maxIndex >= self.dataOut.nProfiles):
321 321 maxIndex = self.dataOut.nProfiles-1
322 322
323 323 #Spectra
324 324 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
325 325
326 326 data_cspc = None
327 327 if self.dataOut.data_cspc is not None:
328 328 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
329 329
330 330 data_dc = None
331 331 if self.dataOut.data_dc is not None:
332 332 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
333 333
334 334 self.dataOut.data_spc = data_spc
335 335 self.dataOut.data_cspc = data_cspc
336 336 self.dataOut.data_dc = data_dc
337 337
338 338 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
339 339 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
340 340 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
341 341
342 342 return 1
343 343
344 344 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
345 345 # validacion de rango
346 346 print("NOISeeee")
347 347 if minHei == None:
348 348 minHei = self.dataOut.heightList[0]
349 349
350 350 if maxHei == None:
351 351 maxHei = self.dataOut.heightList[-1]
352 352
353 353 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
354 354 print('minHei: %.2f is out of the heights range' % (minHei))
355 355 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
356 356 minHei = self.dataOut.heightList[0]
357 357
358 358 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
359 359 print('maxHei: %.2f is out of the heights range' % (maxHei))
360 360 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
361 361 maxHei = self.dataOut.heightList[-1]
362 362
363 363 # validacion de velocidades
364 364 velrange = self.dataOut.getVelRange(1)
365 365
366 366 if minVel == None:
367 367 minVel = velrange[0]
368 368
369 369 if maxVel == None:
370 370 maxVel = velrange[-1]
371 371
372 372 if (minVel < velrange[0]) or (minVel > maxVel):
373 373 print('minVel: %.2f is out of the velocity range' % (minVel))
374 374 print('minVel is setting to %.2f' % (velrange[0]))
375 375 minVel = velrange[0]
376 376
377 377 if (maxVel > velrange[-1]) or (maxVel < minVel):
378 378 print('maxVel: %.2f is out of the velocity range' % (maxVel))
379 379 print('maxVel is setting to %.2f' % (velrange[-1]))
380 380 maxVel = velrange[-1]
381 381
382 382 # seleccion de indices para rango
383 383 minIndex = 0
384 384 maxIndex = 0
385 385 heights = self.dataOut.heightList
386 386
387 387 inda = numpy.where(heights >= minHei)
388 388 indb = numpy.where(heights <= maxHei)
389 389
390 390 try:
391 391 minIndex = inda[0][0]
392 392 except:
393 393 minIndex = 0
394 394
395 395 try:
396 396 maxIndex = indb[0][-1]
397 397 except:
398 398 maxIndex = len(heights)
399 399
400 400 if (minIndex < 0) or (minIndex > maxIndex):
401 401 raise ValueError("some value in (%d,%d) is not valid" % (
402 402 minIndex, maxIndex))
403 403
404 404 if (maxIndex >= self.dataOut.nHeights):
405 405 maxIndex = self.dataOut.nHeights - 1
406 406
407 407 # seleccion de indices para velocidades
408 408 indminvel = numpy.where(velrange >= minVel)
409 409 indmaxvel = numpy.where(velrange <= maxVel)
410 410 try:
411 411 minIndexVel = indminvel[0][0]
412 412 except:
413 413 minIndexVel = 0
414 414
415 415 try:
416 416 maxIndexVel = indmaxvel[0][-1]
417 417 except:
418 418 maxIndexVel = len(velrange)
419 419
420 420 # seleccion del espectro
421 421 data_spc = self.dataOut.data_spc[:,
422 422 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
423 423 # estimacion de ruido
424 424 noise = numpy.zeros(self.dataOut.nChannels)
425 425
426 426 for channel in range(self.dataOut.nChannels):
427 427 daux = data_spc[channel, :, :]
428 428 sortdata = numpy.sort(daux, axis=None)
429 429 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
430 430
431 431 self.dataOut.noise_estimation = noise.copy()
432 432
433 433 return 1
434 434
435 435 class GetSNR(Operation):
436 436 '''
437 437 Written by R. Flores
438 438 '''
439 439 """Operation to get SNR.
440 440
441 441 Parameters:
442 442 -----------
443 443
444 444 Example
445 445 --------
446 446
447 447 op = proc_unit.addOperation(name='GetSNR', optype='other')
448 448
449 449 """
450 450
451 451 def __init__(self, **kwargs):
452 452
453 453 Operation.__init__(self, **kwargs)
454 454
455 455
456 456 def run(self,dataOut):
457 457
458 #noise = dataOut.getNoise()
459 noise = dataOut.getNoise(ymin_index=-10) #Región superior donde solo debería de haber ruido
458 noise = dataOut.getNoise()
459 #noise = dataOut.getNoise(ymin_index=-10) #Región superior donde solo debería de haber ruido
460 460 #print("Noise: ", noise)
461 461 #print("Noise_dB: ", 10*numpy.log10(noise/dataOut.normFactor))
462 462 #print("Heights: ", dataOut.heightList)
463 463 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.normFactor)
464 464 ################dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints) #Before 12Jan2023
465 465 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None])/(noise[:,None])
466 466 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
467 467 dataOut.snl = numpy.log10(dataOut.data_snr)
468 468 #print("snl: ", dataOut.snl)
469 469 #exit(1)
470 470 #print(dataOut.heightList[-11])
471 471 #print(numpy.shape(dataOut.heightList))
472 472 #print(dataOut.data_snr)
473 473 #print(dataOut.data_snr[0,-11])
474 474 #exit(1)
475 475 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr)
476 476 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.1, numpy.nan, dataOut.data_snr)
477 477 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.0, numpy.nan, dataOut.data_snr)
478 478 #dataOut.data_snr = numpy.where(dataOut.data_snr<.05, numpy.nan, dataOut.data_snr)
479 479 #dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
480 480 dataOut.snl = numpy.where(dataOut.snl<-1, numpy.nan, dataOut.snl)
481 481 '''
482 482 import matplotlib.pyplot as plt
483 483 #plt.plot(10*numpy.log10(dataOut.data_snr[0]),dataOut.heightList)
484 484 plt.plot(dataOut.data_snr[0],dataOut.heightList)#,marker='*')
485 485 plt.xlim(-1,10)
486 486 plt.axvline(1,color='k')
487 487 plt.axvline(.1,color='k',linestyle='--')
488 488 plt.grid()
489 489 plt.show()
490 490 '''
491 491 #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr)
492 492 #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0)
493 493 #print(dataOut.data_snr.shape)
494 494 #exit(1)
495 495 #print("Before: ", dataOut.data_snr[0])
496 496
497 497
498 498 return dataOut
499 499
500 500 class removeDC(Operation):
501 501
502 502 def run(self, dataOut, mode=2):
503 503 self.dataOut = dataOut
504 504 jspectra = self.dataOut.data_spc
505 505 jcspectra = self.dataOut.data_cspc
506 506
507 507 num_chan = jspectra.shape[0]
508 508 num_hei = jspectra.shape[2]
509 509
510 510 if jcspectra is not None:
511 511 jcspectraExist = True
512 512 num_pairs = jcspectra.shape[0]
513 513 else:
514 514 jcspectraExist = False
515 515
516 516 freq_dc = int(jspectra.shape[1] / 2)
517 517 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
518 518 ind_vel = ind_vel.astype(int)
519 519
520 520 if ind_vel[0] < 0:
521 521 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
522 522
523 523 if mode == 1:
524 524 jspectra[:, freq_dc, :] = (
525 525 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
526 526
527 527 if jcspectraExist:
528 528 jcspectra[:, freq_dc, :] = (
529 529 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
530 530
531 531 if mode == 2:
532 532
533 533 vel = numpy.array([-2, -1, 1, 2])
534 534 xx = numpy.zeros([4, 4])
535 535
536 536 for fil in range(4):
537 537 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
538 538
539 539 xx_inv = numpy.linalg.inv(xx)
540 540 xx_aux = xx_inv[0, :]
541 541
542 542 for ich in range(num_chan):
543 543 yy = jspectra[ich, ind_vel, :]
544 544 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
545 545
546 546 junkid = jspectra[ich, freq_dc, :] <= 0
547 547 cjunkid = sum(junkid)
548 548
549 549 if cjunkid.any():
550 550 jspectra[ich, freq_dc, junkid.nonzero()] = (
551 551 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
552 552
553 553 if jcspectraExist:
554 554 for ip in range(num_pairs):
555 555 yy = jcspectra[ip, ind_vel, :]
556 556 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
557 557
558 558 self.dataOut.data_spc = jspectra
559 559 self.dataOut.data_cspc = jcspectra
560 560
561 561 return self.dataOut
562 562
563 563 class removeInterference(Operation):
564 564
565 565 def removeInterference2(self):
566 566
567 567 cspc = self.dataOut.data_cspc
568 568 spc = self.dataOut.data_spc
569 569 Heights = numpy.arange(cspc.shape[2])
570 570 realCspc = numpy.abs(cspc)
571 571
572 572 for i in range(cspc.shape[0]):
573 573 LinePower= numpy.sum(realCspc[i], axis=0)
574 574 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
575 575 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
576 576 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
577 577 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
578 578 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
579 579
580 580
581 581 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
582 582 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
583 583 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
584 584 cspc[i,InterferenceRange,:] = numpy.NaN
585 585
586 586 self.dataOut.data_cspc = cspc
587 587
588 588 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
589 589
590 590 jspectra = self.dataOut.data_spc
591 591 jcspectra = self.dataOut.data_cspc
592 592 jnoise = self.dataOut.getNoise()
593 593 num_incoh = self.dataOut.nIncohInt
594 594
595 595 num_channel = jspectra.shape[0]
596 596 num_prof = jspectra.shape[1]
597 597 num_hei = jspectra.shape[2]
598 598
599 599 # hei_interf
600 600 if hei_interf is None:
601 601 count_hei = int(num_hei / 2)
602 602 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
603 603 hei_interf = numpy.asarray(hei_interf)[0]
604 604 # nhei_interf
605 605 if (nhei_interf == None):
606 606 nhei_interf = 5
607 607 if (nhei_interf < 1):
608 608 nhei_interf = 1
609 609 if (nhei_interf > count_hei):
610 610 nhei_interf = count_hei
611 611 if (offhei_interf == None):
612 612 offhei_interf = 0
613 613
614 614 ind_hei = list(range(num_hei))
615 615 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
616 616 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
617 617 mask_prof = numpy.asarray(list(range(num_prof)))
618 618 num_mask_prof = mask_prof.size
619 619 comp_mask_prof = [0, num_prof / 2]
620 620
621 621 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
622 622 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
623 623 jnoise = numpy.nan
624 624 noise_exist = jnoise[0] < numpy.Inf
625 625
626 626 # Subrutina de Remocion de la Interferencia
627 627 for ich in range(num_channel):
628 628 # Se ordena los espectros segun su potencia (menor a mayor)
629 629 power = jspectra[ich, mask_prof, :]
630 630 power = power[:, hei_interf]
631 631 power = power.sum(axis=0)
632 632 psort = power.ravel().argsort()
633 633
634 634 # Se estima la interferencia promedio en los Espectros de Potencia empleando
635 635 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
636 636 offhei_interf, nhei_interf + offhei_interf))]]]
637 637
638 638 if noise_exist:
639 639 # tmp_noise = jnoise[ich] / num_prof
640 640 tmp_noise = jnoise[ich]
641 641 junkspc_interf = junkspc_interf - tmp_noise
642 642 #junkspc_interf[:,comp_mask_prof] = 0
643 643
644 644 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
645 645 jspc_interf = jspc_interf.transpose()
646 646 # Calculando el espectro de interferencia promedio
647 647 noiseid = numpy.where(
648 648 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
649 649 noiseid = noiseid[0]
650 650 cnoiseid = noiseid.size
651 651 interfid = numpy.where(
652 652 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
653 653 interfid = interfid[0]
654 654 cinterfid = interfid.size
655 655
656 656 if (cnoiseid > 0):
657 657 jspc_interf[noiseid] = 0
658 658
659 659 # Expandiendo los perfiles a limpiar
660 660 if (cinterfid > 0):
661 661 new_interfid = (
662 662 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
663 663 new_interfid = numpy.asarray(new_interfid)
664 664 new_interfid = {x for x in new_interfid}
665 665 new_interfid = numpy.array(list(new_interfid))
666 666 new_cinterfid = new_interfid.size
667 667 else:
668 668 new_cinterfid = 0
669 669
670 670 for ip in range(new_cinterfid):
671 671 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
672 672 jspc_interf[new_interfid[ip]
673 673 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
674 674
675 675 jspectra[ich, :, ind_hei] = jspectra[ich, :,
676 676 ind_hei] - jspc_interf # Corregir indices
677 677
678 678 # Removiendo la interferencia del punto de mayor interferencia
679 679 ListAux = jspc_interf[mask_prof].tolist()
680 680 maxid = ListAux.index(max(ListAux))
681 681
682 682 if cinterfid > 0:
683 683 for ip in range(cinterfid * (interf == 2) - 1):
684 684 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
685 685 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
686 686 cind = len(ind)
687 687
688 688 if (cind > 0):
689 689 jspectra[ich, interfid[ip], ind] = tmp_noise * \
690 690 (1 + (numpy.random.uniform(cind) - 0.5) /
691 691 numpy.sqrt(num_incoh))
692 692
693 693 ind = numpy.array([-2, -1, 1, 2])
694 694 xx = numpy.zeros([4, 4])
695 695
696 696 for id1 in range(4):
697 697 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
698 698
699 699 xx_inv = numpy.linalg.inv(xx)
700 700 xx = xx_inv[:, 0]
701 701 ind = (ind + maxid + num_mask_prof) % num_mask_prof
702 702 yy = jspectra[ich, mask_prof[ind], :]
703 703 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
704 704 yy.transpose(), xx)
705 705
706 706 indAux = (jspectra[ich, :, :] < tmp_noise *
707 707 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
708 708 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
709 709 (1 - 1 / numpy.sqrt(num_incoh))
710 710
711 711 # Remocion de Interferencia en el Cross Spectra
712 712 if jcspectra is None:
713 713 return jspectra, jcspectra
714 714 num_pairs = int(jcspectra.size / (num_prof * num_hei))
715 715 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
716 716
717 717 for ip in range(num_pairs):
718 718
719 719 #-------------------------------------------
720 720
721 721 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
722 722 cspower = cspower[:, hei_interf]
723 723 cspower = cspower.sum(axis=0)
724 724
725 725 cspsort = cspower.ravel().argsort()
726 726 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
727 727 offhei_interf, nhei_interf + offhei_interf))]]]
728 728 junkcspc_interf = junkcspc_interf.transpose()
729 729 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
730 730
731 731 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
732 732
733 733 median_real = int(numpy.median(numpy.real(
734 734 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
735 735 median_imag = int(numpy.median(numpy.imag(
736 736 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
737 737 comp_mask_prof = [int(e) for e in comp_mask_prof]
738 738 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
739 739 median_real, median_imag)
740 740
741 741 for iprof in range(num_prof):
742 742 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
743 743 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
744 744
745 745 # Removiendo la Interferencia
746 746 jcspectra[ip, :, ind_hei] = jcspectra[ip,
747 747 :, ind_hei] - jcspc_interf
748 748
749 749 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
750 750 maxid = ListAux.index(max(ListAux))
751 751
752 752 ind = numpy.array([-2, -1, 1, 2])
753 753 xx = numpy.zeros([4, 4])
754 754
755 755 for id1 in range(4):
756 756 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
757 757
758 758 xx_inv = numpy.linalg.inv(xx)
759 759 xx = xx_inv[:, 0]
760 760
761 761 ind = (ind + maxid + num_mask_prof) % num_mask_prof
762 762 yy = jcspectra[ip, mask_prof[ind], :]
763 763 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
764 764
765 765 # Guardar Resultados
766 766 self.dataOut.data_spc = jspectra
767 767 self.dataOut.data_cspc = jcspectra
768 768
769 769 return 1
770 770
771 771 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
772 772
773 773 self.dataOut = dataOut
774 774
775 775 if mode == 1:
776 776 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
777 777 elif mode == 2:
778 778 self.removeInterference2()
779 779
780 780 return self.dataOut
781 781
782 782 class removeInterferenceAtFreq(Operation):
783 783 '''
784 784 Written by R. Flores
785 785 '''
786 786 """Operation to remove interfernce at a known frequency(s).
787 787
788 788 Parameters:
789 789 -----------
790 790 None
791 791
792 792 Example
793 793 --------
794 794
795 795 op = proc_unit.addOperation(name='removeInterferenceAtFreq')
796 796
797 797 """
798 798
799 799 def __init__(self):
800 800
801 801 Operation.__init__(self)
802 802
803 803 def run(self, dataOut, freq = None, freqList = None):
804 804
805 805 VelRange = dataOut.getVelRange()
806 806 #print("VelRange: ", VelRange)
807 807
808 808 freq_ids = []
809 809
810 810 if freq is not None:
811 811 #print("freq")
812 812 #if freq < 0:
813 813 inda = numpy.where(VelRange >= freq)
814 814 minIndex = inda[0][0]
815 815 #print(numpy.shape(dataOut.dataLag_spc))
816 816 dataOut.data_spc[:,minIndex,:] = numpy.nan
817 817
818 818 #inda = numpy.where(VelRange >= ymin_noise)
819 819 #indb = numpy.where(VelRange <= ymax_noise)
820 820
821 821 #minIndex = inda[0][0]
822 822 #maxIndex = indb[0][-1]
823 823
824 824 elif freqList is not None:
825 825 #print("freqList")
826 826 for freq in freqList:
827 827 #if freq < 0:
828 828 inda = numpy.where(VelRange >= freq)
829 829 minIndex = inda[0][0]
830 830 #print(numpy.shape(dataOut.dataLag_spc))
831 831 if freq > 0:
832 832 #dataOut.data_spc[:,minIndex-1,:] = numpy.nan
833 833 freq_ids.append(minIndex-1)
834 834 else:
835 835 #dataOut.data_spc[:,minIndex,:] = numpy.nan
836 836 freq_ids.append(minIndex)
837 837 else:
838 838 raise ValueError("freq or freqList should be specified ...")
839 839
840 840 #freq_ids = numpy.array(freq_ids).flatten()
841 841
842 842 avg = numpy.mean(dataOut.data_spc[:,[t for t in range(dataOut.data_spc.shape[0]) if t not in freq_ids],:],axis=1)
843 843
844 844 for p in list(freq_ids):
845 845 dataOut.data_spc[:,p,:] = avg#numpy.nan
846 846
847 847
848 848 return dataOut
849 849
850 850 class IncohInt(Operation):
851 851
852 852 __profIndex = 0
853 853 __withOverapping = False
854 854
855 855 __byTime = False
856 856 __initime = None
857 857 __lastdatatime = None
858 858 __integrationtime = None
859 859
860 860 __buffer_spc = None
861 861 __buffer_cspc = None
862 862 __buffer_dc = None
863 863
864 864 __dataReady = False
865 865
866 866 __timeInterval = None
867 867
868 868 n = None
869 869
870 870 def __init__(self):
871 871
872 872 Operation.__init__(self)
873 873
874 874 def setup(self, n=None, timeInterval=None, overlapping=False):
875 875 """
876 876 Set the parameters of the integration class.
877 877
878 878 Inputs:
879 879
880 880 n : Number of coherent integrations
881 881 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
882 882 overlapping :
883 883
884 884 """
885 885
886 886 self.__initime = None
887 887 self.__lastdatatime = 0
888 888
889 889 self.__buffer_spc = 0
890 890 self.__buffer_cspc = 0
891 891 self.__buffer_dc = 0
892 892
893 893 self.__profIndex = 0
894 894 self.__dataReady = False
895 895 self.__byTime = False
896 896
897 897 if n is None and timeInterval is None:
898 898 raise ValueError("n or timeInterval should be specified ...")
899 899
900 900 if n is not None:
901 901 self.n = int(n)
902 902 else:
903 903
904 904 self.__integrationtime = int(timeInterval)
905 905 self.n = None
906 906 self.__byTime = True
907 907
908 908 def putData(self, data_spc, data_cspc, data_dc):
909 909 """
910 910 Add a profile to the __buffer_spc and increase in one the __profileIndex
911 911
912 912 """
913 913
914 914 self.__buffer_spc += data_spc
915 915
916 916 if data_cspc is None:
917 917 self.__buffer_cspc = None
918 918 else:
919 919 self.__buffer_cspc += data_cspc
920 920
921 921 if data_dc is None:
922 922 self.__buffer_dc = None
923 923 else:
924 924 self.__buffer_dc += data_dc
925 925
926 926 self.__profIndex += 1
927 927
928 928 return
929 929
930 930 def pushData(self):
931 931 """
932 932 Return the sum of the last profiles and the profiles used in the sum.
933 933
934 934 Affected:
935 935
936 936 self.__profileIndex
937 937
938 938 """
939 939
940 940 data_spc = self.__buffer_spc
941 941 data_cspc = self.__buffer_cspc
942 942 data_dc = self.__buffer_dc
943 943 n = self.__profIndex
944 944
945 945 self.__buffer_spc = 0
946 946 self.__buffer_cspc = 0
947 947 self.__buffer_dc = 0
948 948 self.__profIndex = 0
949 949
950 950 return data_spc, data_cspc, data_dc, n
951 951
952 952 def byProfiles(self, *args):
953 953
954 954 self.__dataReady = False
955 955 avgdata_spc = None
956 956 avgdata_cspc = None
957 957 avgdata_dc = None
958 958
959 959 self.putData(*args)
960 960
961 961 if self.__profIndex == self.n:
962 962
963 963 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
964 964 self.n = n
965 965 self.__dataReady = True
966 966
967 967 return avgdata_spc, avgdata_cspc, avgdata_dc
968 968
969 969 def byTime(self, datatime, *args):
970 970
971 971 self.__dataReady = False
972 972 avgdata_spc = None
973 973 avgdata_cspc = None
974 974 avgdata_dc = None
975 975
976 976 self.putData(*args)
977 977
978 978 if (datatime - self.__initime) >= self.__integrationtime:
979 979 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
980 980 self.n = n
981 981 self.__dataReady = True
982 982
983 983 return avgdata_spc, avgdata_cspc, avgdata_dc
984 984
985 985 def integrate(self, datatime, *args):
986 986
987 987 if self.__profIndex == 0:
988 988 self.__initime = datatime
989 989
990 990 if self.__byTime:
991 991 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
992 992 datatime, *args)
993 993 else:
994 994 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
995 995
996 996 if not self.__dataReady:
997 997 return None, None, None, None
998 998
999 999 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1000 1000
1001 1001 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1002 1002 if n == 1:
1003 1003 return dataOut
1004 1004 print("JERE")
1005 1005 dataOut.flagNoData = True
1006 1006
1007 1007 if not self.isConfig:
1008 1008 self.setup(n, timeInterval, overlapping)
1009 1009 self.isConfig = True
1010 1010
1011 1011 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1012 1012 dataOut.data_spc,
1013 1013 dataOut.data_cspc,
1014 1014 dataOut.data_dc)
1015 1015
1016 1016 if self.__dataReady:
1017 1017
1018 1018 dataOut.data_spc = avgdata_spc
1019 1019 print(numpy.sum(dataOut.data_spc))
1020 1020 exit(1)
1021 1021 dataOut.data_cspc = avgdata_cspc
1022 1022 dataOut.data_dc = avgdata_dc
1023 1023 dataOut.nIncohInt *= self.n
1024 1024 dataOut.utctime = avgdatatime
1025 1025 dataOut.flagNoData = False
1026 1026
1027 1027 return dataOut
1028 1028
1029 1029 class dopplerFlip(Operation):
1030 1030
1031 1031 def run(self, dataOut, chann = None):
1032 1032 # arreglo 1: (num_chan, num_profiles, num_heights)
1033 1033 self.dataOut = dataOut
1034 1034 # JULIA-oblicua, indice 2
1035 1035 # arreglo 2: (num_profiles, num_heights)
1036 1036 jspectra = self.dataOut.data_spc[chann]
1037 1037 jspectra_tmp = numpy.zeros(jspectra.shape)
1038 1038 num_profiles = jspectra.shape[0]
1039 1039 freq_dc = int(num_profiles / 2)
1040 1040 # Flip con for
1041 1041 for j in range(num_profiles):
1042 1042 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1043 1043 # Intercambio perfil de DC con perfil inmediato anterior
1044 1044 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1045 1045 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1046 1046 # canal modificado es re-escrito en el arreglo de canales
1047 1047 self.dataOut.data_spc[chann] = jspectra_tmp
1048 1048
1049 1049 return self.dataOut
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now