##// END OF EJS Templates
plot
joabAM -
r1749:e3f1b2c53e4b
parent child
Show More
@@ -1,1935 +1,1935
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 datetime
12 12
13 13 from schainpy.model.graphics.jroplot_base import Plot, plt, log
14 14 from itertools import combinations
15 15 from matplotlib.ticker import LinearLocator
16 16
17 17 from schainpy.model.utils.BField import BField
18 18 from scipy.interpolate import splrep
19 19 from scipy.interpolate import splev
20 20
21 21 from matplotlib import __version__ as plt_version
22 22
23 23 if plt_version >='3.3.4':
24 24 EXTRA_POINTS = 0
25 25 else:
26 26 EXTRA_POINTS = 1
27 27 class SpectraPlot(Plot):
28 28 '''
29 29 Plot for Spectra data
30 30 '''
31 31
32 32 CODE = 'spc'
33 33 colormap = 'jet'
34 34 plot_type = 'pcolor'
35 35 buffering = False
36 36 channelList = []
37 37 elevationList = []
38 38 azimuthList = []
39 39
40 40 def setup(self):
41 41
42 42 self.nplots = len(self.data.channels)
43 43 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
44 44 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
45 45 self.height = 3.4 * self.nrows
46 46 self.cb_label = 'dB'
47 47 if self.showprofile:
48 48 self.width = 5.2 * self.ncols
49 49 else:
50 50 self.width = 4.2* self.ncols
51 51 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.12})
52 52 self.ylabel = 'Range [km]'
53 53
54 54 def update_list(self,dataOut):
55 55
56 56 if len(self.channelList) == 0:
57 57 self.channelList = dataOut.channelList
58 58 if len(self.elevationList) == 0:
59 59 self.elevationList = dataOut.elevationList
60 60 if len(self.azimuthList) == 0:
61 61 self.azimuthList = dataOut.azimuthList
62 62
63 63 def update(self, dataOut):
64 64
65 65 self.update_list(dataOut)
66 66 data = {}
67 67 meta = {}
68 68 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
69 69 if dataOut.type == "Parameters":
70 70 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
71 71 spc = 10*numpy.log10(dataOut.data_spc/(dataOut.nProfiles))
72 72 else:
73 73 noise = 10*numpy.log10(dataOut.getNoise()/norm)
74 74
75 75 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
76 76 for ch in range(dataOut.nChannels):
77 77 if hasattr(dataOut.normFactor,'ndim'):
78 78 if dataOut.normFactor.ndim > 1:
79 79 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
80 80
81 81 else:
82 82 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
83 83 else:
84 84 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
85 85
86 86 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
87 87 spc = 10*numpy.log10(z)
88 88
89 89 data['spc'] = spc
90 90 data['rti'] = spc.mean(axis=1)
91 91 data['noise'] = noise
92 92 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
93 93 if self.CODE == 'spc_moments':
94 94 data['moments'] = dataOut.moments
95 95
96 96 return data, meta
97 97
98 98 def plot(self):
99 99
100 100 if self.xaxis == "frequency":
101 101 x = self.data.xrange[0]
102 102 self.xlabel = "Frequency (kHz)"
103 103 elif self.xaxis == "time":
104 104 x = self.data.xrange[1]
105 105 self.xlabel = "Time (ms)"
106 106 else:
107 107 x = self.data.xrange[2]
108 108 self.xlabel = "Velocity (m/s)"
109 109
110 110 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
111 111 x = self.data.xrange[2]
112 112 self.xlabel = "Velocity (m/s)"
113 113
114 114 self.titles = []
115 115
116 116 y = self.data.yrange
117 117 self.y = y
118 118
119 119 data = self.data[-1]
120 120 z = data['spc']
121 121
122 122 for n, ax in enumerate(self.axes):
123 123 noise = self.data['noise'][n][0]
124 124 # noise = data['noise'][n]
125 125
126 126 if self.CODE == 'spc_moments':
127 127 mean = data['moments'][n, 1]
128 128 if self.CODE == 'gaussian_fit':
129 129 gau0 = data['gaussfit'][n][2,:,0]
130 130 gau1 = data['gaussfit'][n][2,:,1]
131 131 if ax.firsttime:
132 132 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
133 133 self.xmin = self.xmin if self.xmin else -self.xmax
134 134 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
135 135 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
136 136 ax.plt = ax.pcolormesh(x, y, z[n].T,
137 137 vmin=self.zmin,
138 138 vmax=self.zmax,
139 139 cmap=plt.get_cmap(self.colormap)
140 140 )
141 141
142 142 if self.showprofile:
143 143 ax.plt_profile = self.pf_axes[n].plot(
144 144 data['rti'][n], y)[0]
145 145 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
146 146 color="k", linestyle="dashed", lw=1)[0]
147 147 if self.CODE == 'spc_moments':
148 148 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
149 149 if self.CODE == 'gaussian_fit':
150 150 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
151 151 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
152 152 else:
153 153 ax.plt.set_array(z[n].T.ravel())
154 154 if self.showprofile:
155 155 ax.plt_profile.set_data(data['rti'][n], y)
156 156 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
157 157 if self.CODE == 'spc_moments':
158 158 ax.plt_mean.set_data(mean, y)
159 159 if self.CODE == 'gaussian_fit':
160 160 ax.plt_gau0.set_data(gau0, y)
161 161 ax.plt_gau1.set_data(gau1, y)
162 162 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
163 163 self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
164 164 else:
165 165 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
166 166
167 167 class SpectraObliquePlot(Plot):
168 168 '''
169 169 Plot for Spectra data
170 170 '''
171 171
172 172 CODE = 'spc_oblique'
173 173 colormap = 'jet'
174 174 plot_type = 'pcolor'
175 175
176 176 def setup(self):
177 177 self.xaxis = "oblique"
178 178 self.nplots = len(self.data.channels)
179 179 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
180 180 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
181 181 self.height = 2.6 * self.nrows
182 182 self.cb_label = 'dB'
183 183 if self.showprofile:
184 184 self.width = 4 * self.ncols
185 185 else:
186 186 self.width = 3.5 * self.ncols
187 187 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
188 188 self.ylabel = 'Range [km]'
189 189
190 190 def update(self, dataOut):
191 191
192 192 data = {}
193 193 meta = {}
194 194
195 195 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
196 196 data['spc'] = spc
197 197 data['rti'] = dataOut.getPower()
198 198 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
199 199 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
200 200
201 201 data['shift1'] = dataOut.Dop_EEJ_T1[0]
202 202 data['shift2'] = dataOut.Dop_EEJ_T2[0]
203 203 data['max_val_2'] = dataOut.Oblique_params[0,-1,:]
204 204 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
205 205 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
206 206
207 207 return data, meta
208 208
209 209 def plot(self):
210 210
211 211 if self.xaxis == "frequency":
212 212 x = self.data.xrange[0]
213 213 self.xlabel = "Frequency (kHz)"
214 214 elif self.xaxis == "time":
215 215 x = self.data.xrange[1]
216 216 self.xlabel = "Time (ms)"
217 217 else:
218 218 x = self.data.xrange[2]
219 219 self.xlabel = "Velocity (m/s)"
220 220
221 221 self.titles = []
222 222
223 223 y = self.data.yrange
224 224 self.y = y
225 225
226 226 data = self.data[-1]
227 227 z = data['spc']
228 228
229 229 for n, ax in enumerate(self.axes):
230 230 noise = self.data['noise'][n][-1]
231 231 shift1 = data['shift1']
232 232 shift2 = data['shift2']
233 233 max_val_2 = data['max_val_2']
234 234 err1 = data['shift1_error']
235 235 err2 = data['shift2_error']
236 236 if ax.firsttime:
237 237
238 238 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
239 239 self.xmin = self.xmin if self.xmin else -self.xmax
240 240 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
241 241 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
242 242 ax.plt = ax.pcolormesh(x, y, z[n].T,
243 243 vmin=self.zmin,
244 244 vmax=self.zmax,
245 245 cmap=plt.get_cmap(self.colormap)
246 246 )
247 247
248 248 if self.showprofile:
249 249 ax.plt_profile = self.pf_axes[n].plot(
250 250 self.data['rti'][n][-1], y)[0]
251 251 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
252 252 color="k", linestyle="dashed", lw=1)[0]
253 253
254 254 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)
255 255 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)
256 256 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)
257 257
258 258 else:
259 259 self.ploterr1.remove()
260 260 self.ploterr2.remove()
261 261 self.ploterr3.remove()
262 262 ax.plt.set_array(z[n].T.ravel())
263 263 if self.showprofile:
264 264 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
265 265 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
266 266 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)
267 267 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)
268 268 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)
269 269
270 270 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
271 271
272 272
273 273 class CrossSpectraPlot(Plot):
274 274
275 275 CODE = 'cspc'
276 276 colormap = 'jet'
277 277 plot_type = 'pcolor'
278 278 zmin_coh = None
279 279 zmax_coh = None
280 280 zmin_phase = None
281 281 zmax_phase = None
282 282 realChannels = None
283 283 crossPairs = None
284 284
285 285 def setup(self):
286 286
287 287 self.ncols = 4
288 288 self.nplots = len(self.data.pairs) * 2
289 289 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
290 290 self.width = 3.1 * self.ncols
291 291 self.height = 2.6 * self.nrows
292 292 self.ylabel = 'Range [km]'
293 293 self.showprofile = False
294 294 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
295 295
296 296 def update(self, dataOut):
297 297
298 298 data = {}
299 299 meta = {}
300 300
301 301 spc = dataOut.data_spc
302 302 cspc = dataOut.data_cspc
303 303 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
304 304 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
305 305 meta['pairs'] = rawPairs
306 306 if self.crossPairs == None:
307 307 self.crossPairs = dataOut.pairsList
308 308 tmp = []
309 309
310 310 for n, pair in enumerate(meta['pairs']):
311 311 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
312 312 coh = numpy.abs(out)
313 313 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
314 314 tmp.append(coh)
315 315 tmp.append(phase)
316 316
317 317 data['cspc'] = numpy.array(tmp)
318 318
319 319 return data, meta
320 320
321 321 def plot(self):
322 322
323 323 if self.xaxis == "frequency":
324 324 x = self.data.xrange[0]
325 325 self.xlabel = "Frequency (kHz)"
326 326 elif self.xaxis == "time":
327 327 x = self.data.xrange[1]
328 328 self.xlabel = "Time (ms)"
329 329 else:
330 330 x = self.data.xrange[2]
331 331 self.xlabel = "Velocity (m/s)"
332 332
333 333 self.titles = []
334 334
335 335 y = self.data.yrange
336 336 self.y = y
337 337
338 338 data = self.data[-1]
339 339 cspc = data['cspc']
340 340
341 341 for n in range(len(self.data.pairs)):
342 342 pair = self.crossPairs[n]
343 343 coh = cspc[n*2]
344 344 phase = cspc[n*2+1]
345 345 ax = self.axes[2 * n]
346 346 if ax.firsttime:
347 347 ax.plt = ax.pcolormesh(x, y, coh.T,
348 348 vmin=self.zmin_coh,
349 349 vmax=self.zmax_coh,
350 350 cmap=plt.get_cmap(self.colormap_coh)
351 351 )
352 352 else:
353 353 ax.plt.set_array(coh.T.ravel())
354 354 self.titles.append(
355 355 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
356 356
357 357 ax = self.axes[2 * n + 1]
358 358 if ax.firsttime:
359 359 ax.plt = ax.pcolormesh(x, y, phase.T,
360 360 vmin=-180,
361 361 vmax=180,
362 362 cmap=plt.get_cmap(self.colormap_phase)
363 363 )
364 364 else:
365 365 ax.plt.set_array(phase.T.ravel())
366 366 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
367 367
368 368
369 369 class CrossSpectra4Plot(Plot):
370 370
371 371 CODE = 'cspc'
372 372 colormap = 'jet'
373 373 plot_type = 'pcolor'
374 374 zmin_coh = None
375 375 zmax_coh = None
376 376 zmin_phase = None
377 377 zmax_phase = None
378 378
379 379 def setup(self):
380 380
381 381 self.ncols = 4
382 382 self.nrows = len(self.data.pairs)
383 383 self.nplots = self.nrows * 4
384 384 self.width = 3.1 * self.ncols
385 385 self.height = 5 * self.nrows
386 386 self.ylabel = 'Range [km]'
387 387 self.showprofile = False
388 388 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
389 389
390 390 def plot(self):
391 391
392 392 if self.xaxis == "frequency":
393 393 x = self.data.xrange[0]
394 394 self.xlabel = "Frequency (kHz)"
395 395 elif self.xaxis == "time":
396 396 x = self.data.xrange[1]
397 397 self.xlabel = "Time (ms)"
398 398 else:
399 399 x = self.data.xrange[2]
400 400 self.xlabel = "Velocity (m/s)"
401 401
402 402 self.titles = []
403 403
404 404
405 405 y = self.data.heights
406 406 self.y = y
407 407 nspc = self.data['spc']
408 408 spc = self.data['cspc'][0]
409 409 cspc = self.data['cspc'][1]
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
415 415 ax = self.axes[4 * n]
416 416 if ax.firsttime:
417 417 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
418 418 self.xmin = self.xmin if self.xmin else -self.xmax
419 419 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
420 420 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
421 421 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
422 422 vmin=self.zmin,
423 423 vmax=self.zmax,
424 424 cmap=plt.get_cmap(self.colormap)
425 425 )
426 426 else:
427 427
428 428 ax.plt.set_array(nspc[pair[0]].T.ravel())
429 429 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
430 430
431 431 ax = self.axes[4 * n + 1]
432 432
433 433 if ax.firsttime:
434 434 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
435 435 vmin=self.zmin,
436 436 vmax=self.zmax,
437 437 cmap=plt.get_cmap(self.colormap)
438 438 )
439 439 else:
440 440
441 441 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
442 442 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
443 443
444 444 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
445 445 coh = numpy.abs(out)
446 446 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
447 447
448 448 ax = self.axes[4 * n + 2]
449 449 if ax.firsttime:
450 450 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
451 451 vmin=0,
452 452 vmax=1,
453 453 cmap=plt.get_cmap(self.colormap_coh)
454 454 )
455 455 else:
456 456 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
457 457 self.titles.append(
458 458 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
459 459
460 460 ax = self.axes[4 * n + 3]
461 461 if ax.firsttime:
462 462 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
463 463 vmin=-180,
464 464 vmax=180,
465 465 cmap=plt.get_cmap(self.colormap_phase)
466 466 )
467 467 else:
468 468 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
469 469 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
470 470
471 471
472 472 class CrossSpectra2Plot(Plot):
473 473
474 474 CODE = 'cspc'
475 475 colormap = 'jet'
476 476 plot_type = 'pcolor'
477 477 zmin_coh = None
478 478 zmax_coh = None
479 479 zmin_phase = None
480 480 zmax_phase = None
481 481
482 482 def setup(self):
483 483
484 484 self.ncols = 1
485 485 self.nrows = len(self.data.pairs)
486 486 self.nplots = self.nrows * 1
487 487 self.width = 3.1 * self.ncols
488 488 self.height = 5 * self.nrows
489 489 self.ylabel = 'Range [km]'
490 490 self.showprofile = False
491 491 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
492 492
493 493 def plot(self):
494 494
495 495 if self.xaxis == "frequency":
496 496 x = self.data.xrange[0]
497 497 self.xlabel = "Frequency (kHz)"
498 498 elif self.xaxis == "time":
499 499 x = self.data.xrange[1]
500 500 self.xlabel = "Time (ms)"
501 501 else:
502 502 x = self.data.xrange[2]
503 503 self.xlabel = "Velocity (m/s)"
504 504
505 505 self.titles = []
506 506
507 507
508 508 y = self.data.heights
509 509 self.y = y
510 510 cspc = self.data['cspc'][1]
511 511
512 512 for n in range(self.nrows):
513 513 noise = self.data['noise'][:,-1]
514 514 pair = self.data.pairs[n]
515 515 out = cspc[n]
516 516 cross = numpy.abs(out)
517 517 z = cross/self.data.nFactor
518 518 cross = 10*numpy.log10(z)
519 519
520 520 ax = self.axes[1 * n]
521 521 if ax.firsttime:
522 522 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
523 523 self.xmin = self.xmin if self.xmin else -self.xmax
524 524 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
525 525 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
526 526 ax.plt = ax.pcolormesh(x, y, cross.T,
527 527 vmin=self.zmin,
528 528 vmax=self.zmax,
529 529 cmap=plt.get_cmap(self.colormap)
530 530 )
531 531 else:
532 532 ax.plt.set_array(cross.T.ravel())
533 533 self.titles.append(
534 534 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
535 535
536 536
537 537 class CrossSpectra3Plot(Plot):
538 538
539 539 CODE = 'cspc'
540 540 colormap = 'jet'
541 541 plot_type = 'pcolor'
542 542 zmin_coh = None
543 543 zmax_coh = None
544 544 zmin_phase = None
545 545 zmax_phase = None
546 546
547 547 def setup(self):
548 548
549 549 self.ncols = 3
550 550 self.nrows = len(self.data.pairs)
551 551 self.nplots = self.nrows * 3
552 552 self.width = 3.1 * self.ncols
553 553 self.height = 5 * self.nrows
554 554 self.ylabel = 'Range [km]'
555 555 self.showprofile = False
556 556 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
557 557
558 558 def plot(self):
559 559
560 560 if self.xaxis == "frequency":
561 561 x = self.data.xrange[0]
562 562 self.xlabel = "Frequency (kHz)"
563 563 elif self.xaxis == "time":
564 564 x = self.data.xrange[1]
565 565 self.xlabel = "Time (ms)"
566 566 else:
567 567 x = self.data.xrange[2]
568 568 self.xlabel = "Velocity (m/s)"
569 569
570 570 self.titles = []
571 571
572 572
573 573 y = self.data.heights
574 574 self.y = y
575 575
576 576 cspc = self.data['cspc'][1]
577 577
578 578 for n in range(self.nrows):
579 579 noise = self.data['noise'][:,-1]
580 580 pair = self.data.pairs[n]
581 581 out = cspc[n]
582 582
583 583 cross = numpy.abs(out)
584 584 z = cross/self.data.nFactor
585 585 cross = 10*numpy.log10(z)
586 586
587 587 out_r= out.real/self.data.nFactor
588 588
589 589 out_i= out.imag/self.data.nFactor
590 590
591 591 ax = self.axes[3 * n]
592 592 if ax.firsttime:
593 593 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
594 594 self.xmin = self.xmin if self.xmin else -self.xmax
595 595 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
596 596 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
597 597 ax.plt = ax.pcolormesh(x, y, cross.T,
598 598 vmin=self.zmin,
599 599 vmax=self.zmax,
600 600 cmap=plt.get_cmap(self.colormap)
601 601 )
602 602 else:
603 603 ax.plt.set_array(cross.T.ravel())
604 604 self.titles.append(
605 605 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
606 606
607 607 ax = self.axes[3 * n + 1]
608 608 if ax.firsttime:
609 609 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
610 610 self.xmin = self.xmin if self.xmin else -self.xmax
611 611 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
612 612 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
613 613 ax.plt = ax.pcolormesh(x, y, out_r.T,
614 614 vmin=-1.e6,
615 615 vmax=0,
616 616 cmap=plt.get_cmap(self.colormap)
617 617 )
618 618 else:
619 619 ax.plt.set_array(out_r.T.ravel())
620 620 self.titles.append(
621 621 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
622 622
623 623 ax = self.axes[3 * n + 2]
624 624
625 625
626 626 if ax.firsttime:
627 627 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
628 628 self.xmin = self.xmin if self.xmin else -self.xmax
629 629 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
630 630 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
631 631 ax.plt = ax.pcolormesh(x, y, out_i.T,
632 632 vmin=-1.e6,
633 633 vmax=1.e6,
634 634 cmap=plt.get_cmap(self.colormap)
635 635 )
636 636 else:
637 637 ax.plt.set_array(out_i.T.ravel())
638 638 self.titles.append(
639 639 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
640 640
641 641 class RTIPlot(Plot):
642 642 '''
643 643 Plot for RTI data
644 644 '''
645 645
646 646 CODE = 'rti'
647 647 colormap = 'jet'
648 648 plot_type = 'pcolorbuffer'
649 649 titles = None
650 650 channelList = []
651 651 elevationList = []
652 652 azimuthList = []
653 653
654 654 def setup(self):
655 655 self.xaxis = 'time'
656 656 self.ncols = 1
657 657 self.nrows = len(self.data.channels)
658 658 self.nplots = len(self.data.channels)
659 659 self.ylabel = 'Range [km]'
660 660 #self.xlabel = 'Time'
661 661 self.cb_label = 'dB'
662 662 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
663 663 self.titles = ['{} Channel {}'.format(
664 664 self.CODE.upper(), x) for x in range(self.nplots)]
665 665
666 666 def update_list(self,dataOut):
667 667
668 668 if len(self.channelList) == 0:
669 669 self.channelList = dataOut.channelList
670 670 if len(self.elevationList) == 0:
671 671 self.elevationList = dataOut.elevationList
672 672 if len(self.azimuthList) == 0:
673 673 self.azimuthList = dataOut.azimuthList
674 674
675 675
676 676 def update(self, dataOut):
677 677
678 678 if len(self.channelList) == 0:
679 679 self.update_list(dataOut)
680 680 data = {}
681 681 meta = {}
682 682 data['rti'] = dataOut.getPower()
683 683 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
684 684 noise = 10*numpy.log10(dataOut.getNoise()/norm)
685 685 data['noise'] = noise
686 686
687 687 return data, meta
688 688
689 689 def plot(self):
690 690
691 691 self.x = self.data.times
692 692 self.y = self.data.yrange
693 693 self.z = self.data[self.CODE]
694 694 self.z = numpy.array(self.z, dtype=float)
695 695 self.z = numpy.ma.masked_invalid(self.z)
696 696
697 697 try:
698 698 if self.channelList != None:
699 699 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
700 700 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
701 701 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
702 702 else:
703 703 self.titles = ['{} Channel {}'.format(
704 704 self.CODE.upper(), x) for x in self.channelList]
705 705 except:
706 706 if self.channelList.any() != None:
707 707 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
708 708 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
709 709 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
710 710 else:
711 711 self.titles = ['{} Channel {}'.format(
712 712 self.CODE.upper(), x) for x in self.channelList]
713 713
714 714 if self.decimation is None:
715 715 x, y, z = self.fill_gaps(self.x, self.y, self.z)
716 716 else:
717 717 x, y, z = self.fill_gaps(*self.decimate())
718 718
719 719 for n, ax in enumerate(self.axes):
720 720
721 721 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
722 722 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
723 723 data = self.data[-1]
724 724 if ax.firsttime:
725 725 ax.plt = ax.pcolormesh(x, y, z[n].T,
726 726 vmin=self.zmin,
727 727 vmax=self.zmax,
728 728 cmap=plt.get_cmap(self.colormap)
729 729 )
730 730 if self.showprofile:
731 731 ax.plot_profile = self.pf_axes[n].plot(
732 732 data[self.CODE][n], self.y)[0]
733 733 if "noise" in self.data:
734 734 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
735 735 color="k", linestyle="dashed", lw=1)[0]
736 736 else:
737 # ax.collections.remove(ax.collections[0]) # error while running
737 ax.collections.remove(ax.collections[0]) # error while running in 3.12
738 738 ax.plt = ax.pcolormesh(x, y, z[n].T,
739 739 vmin=self.zmin,
740 740 vmax=self.zmax,
741 741 cmap=plt.get_cmap(self.colormap)
742 742 )
743 743 if self.showprofile:
744 744 ax.plot_profile.set_data(data[self.CODE][n], self.y)
745 745 if "noise" in self.data:
746 746 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
747 747 color="k", linestyle="dashed", lw=1)[0]
748 748
749 749 class SpectrogramPlot(Plot):
750 750 '''
751 751 Plot for Spectrogram data
752 752 '''
753 753
754 754 CODE = 'Spectrogram_Profile'
755 755 colormap = 'binary'
756 756 plot_type = 'pcolorbuffer'
757 757
758 758 def setup(self):
759 759 self.xaxis = 'time'
760 760 self.ncols = 1
761 761 self.nrows = len(self.data.channels)
762 762 self.nplots = len(self.data.channels)
763 763 self.xlabel = 'Time'
764 764 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
765 765 self.titles = []
766 766
767 767 self.titles = ['{} Channel {}'.format(
768 768 self.CODE.upper(), x) for x in range(self.nrows)]
769 769
770 770
771 771 def update(self, dataOut):
772 772 data = {}
773 773 meta = {}
774 774
775 775 maxHei = 1620#+12000
776 776 indb = numpy.where(dataOut.heightList <= maxHei)
777 777 hei = indb[0][-1]
778 778
779 779 factor = dataOut.nIncohInt
780 780 z = dataOut.data_spc[:,:,hei] / factor
781 781 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
782 782
783 783 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
784 784 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
785 785
786 786 data['hei'] = hei
787 787 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
788 788 data['nProfiles'] = dataOut.nProfiles
789 789
790 790 return data, meta
791 791
792 792 def plot(self):
793 793
794 794 self.x = self.data.times
795 795 self.z = self.data[self.CODE]
796 796 self.y = self.data.xrange[0]
797 797
798 798 hei = self.data['hei'][-1]
799 799 DH = self.data['DH'][-1]
800 800 nProfiles = self.data['nProfiles'][-1]
801 801
802 802 self.ylabel = "Frequency (kHz)"
803 803
804 804 self.z = numpy.ma.masked_invalid(self.z)
805 805
806 806 if self.decimation is None:
807 807 x, y, z = self.fill_gaps(self.x, self.y, self.z)
808 808 else:
809 809 x, y, z = self.fill_gaps(*self.decimate())
810 810
811 811 for n, ax in enumerate(self.axes):
812 812 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
813 813 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
814 814 data = self.data[-1]
815 815 if ax.firsttime:
816 816 ax.plt = ax.pcolormesh(x, y, z[n].T,
817 817 vmin=self.zmin,
818 818 vmax=self.zmax,
819 819 cmap=plt.get_cmap(self.colormap)
820 820 )
821 821 else:
822 822 # ax.collections.remove(ax.collections[0]) # error while running
823 823 ax.plt = ax.pcolormesh(x, y, z[n].T,
824 824 vmin=self.zmin,
825 825 vmax=self.zmax,
826 826 cmap=plt.get_cmap(self.colormap)
827 827 )
828 828
829 829
830 830
831 831 class CoherencePlot(RTIPlot):
832 832 '''
833 833 Plot for Coherence data
834 834 '''
835 835
836 836 CODE = 'coh'
837 837 titles = None
838 838
839 839 def setup(self):
840 840 self.xaxis = 'time'
841 841 self.ncols = 1
842 842 self.nrows = len(self.data.pairs)
843 843 self.nplots = len(self.data.pairs)
844 844 self.ylabel = 'Range [km]'
845 845 self.xlabel = 'Time'
846 846 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
847 847 if self.CODE == 'coh':
848 848 self.cb_label = ''
849 849 self.titles = [
850 850 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
851 851 else:
852 852 self.cb_label = 'Degrees'
853 853 self.titles = [
854 854 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
855 855
856 856 def update(self, dataOut):
857 857
858 858 data = {}
859 859 meta = {}
860 860 data['coh'] = dataOut.getCoherence()
861 861 meta['pairs'] = dataOut.pairsList
862 862
863 863 return data, meta
864 864
865 865 class PhasePlot(CoherencePlot):
866 866 '''
867 867 Plot for Phase map data
868 868 '''
869 869
870 870 CODE = 'phase'
871 871 colormap = 'seismic'
872 872
873 873 def update(self, dataOut):
874 874
875 875 data = {}
876 876 meta = {}
877 877 data['phase'] = dataOut.getCoherence(phase=True)
878 878 meta['pairs'] = dataOut.pairsList
879 879
880 880 return data, meta
881 881
882 882 class NoisePlot(Plot):
883 883 '''
884 884 Plot for noise
885 885 '''
886 886
887 887 CODE = 'noise'
888 888 plot_type = 'scatterbuffer'
889 889
890 890 def setup(self):
891 891 self.xaxis = 'time'
892 892 self.ncols = 1
893 893 self.nrows = 1
894 894 self.nplots = 1
895 895 self.ylabel = 'Intensity [dB]'
896 896 self.xlabel = 'Time'
897 897 self.titles = ['Noise']
898 898 self.colorbar = False
899 899 self.plots_adjust.update({'right': 0.85 })
900 900 self.titles = ['Noise Plot']
901 901
902 902 def update(self, dataOut):
903 903
904 904 data = {}
905 905 meta = {}
906 906 noise = 10*numpy.log10(dataOut.getNoise())
907 907 noise = noise.reshape(dataOut.nChannels, 1)
908 908 data['noise'] = noise
909 909 meta['yrange'] = numpy.array([])
910 910
911 911 return data, meta
912 912
913 913 def plot(self):
914 914
915 915 x = self.data.times
916 916 xmin = self.data.min_time
917 917 xmax = xmin + self.xrange * 60 * 60
918 918 Y = self.data['noise']
919 919
920 920 if self.axes[0].firsttime:
921 921 self.ymin = numpy.nanmin(Y) - 5
922 922 self.ymax = numpy.nanmax(Y) + 5
923 923 for ch in self.data.channels:
924 924 y = Y[ch]
925 925 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
926 926 plt.legend(bbox_to_anchor=(1.18, 1.0))
927 927 else:
928 928 for ch in self.data.channels:
929 929 y = Y[ch]
930 930 self.axes[0].lines[ch].set_data(x, y)
931 931
932 932 class PowerProfilePlot(Plot):
933 933
934 934 CODE = 'pow_profile'
935 935 plot_type = 'scatter'
936 936
937 937 def setup(self):
938 938
939 939 self.ncols = 1
940 940 self.nrows = 1
941 941 self.nplots = 1
942 942 self.height = 4
943 943 self.width = 3
944 944 self.ylabel = 'Range [km]'
945 945 self.xlabel = 'Intensity [dB]'
946 946 self.titles = ['Power Profile']
947 947 self.colorbar = False
948 948
949 949 def update(self, dataOut):
950 950
951 951 data = {}
952 952 meta = {}
953 953 data[self.CODE] = dataOut.getPower()
954 954
955 955 return data, meta
956 956
957 957 def plot(self):
958 958
959 959 y = self.data.yrange
960 960 self.y = y
961 961
962 962 x = self.data[-1][self.CODE]
963 963
964 964 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
965 965 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
966 966
967 967 if self.axes[0].firsttime:
968 968 for ch in self.data.channels:
969 969 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
970 970 plt.legend()
971 971 else:
972 972 for ch in self.data.channels:
973 973 self.axes[0].lines[ch].set_data(x[ch], y)
974 974
975 975
976 976 class SpectraCutPlot(Plot):
977 977
978 978 CODE = 'spc_cut'
979 979 plot_type = 'scatter'
980 980 buffering = False
981 981 heights = []
982 982 channelList = []
983 983 maintitle = "Spectra Cuts"
984 984 flag_setIndex = False
985 985
986 986 def setup(self):
987 987
988 988 self.nplots = len(self.data.channels)
989 989 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
990 990 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
991 991 self.width = 4.5 * self.ncols + 2.5
992 992 self.height = 4.8 * self.nrows
993 993 self.ylabel = 'Power [dB]'
994 994 self.colorbar = False
995 995 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
996 996
997 997 if len(self.selectedHeightsList) > 0:
998 998 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
999 999
1000 1000
1001 1001
1002 1002 def update(self, dataOut):
1003 1003 if len(self.channelList) == 0:
1004 1004 self.channelList = dataOut.channelList
1005 1005
1006 1006 self.heights = dataOut.heightList
1007 1007 #print("sels: ",self.selectedHeightsList)
1008 1008 if len(self.selectedHeightsList)>0 and not self.flag_setIndex:
1009 1009
1010 1010 for sel_height in self.selectedHeightsList:
1011 1011 index_list = numpy.where(self.heights >= sel_height)
1012 1012 index_list = index_list[0]
1013 1013 self.height_index.append(index_list[0])
1014 1014 #print("sels i:"", self.height_index)
1015 1015 self.flag_setIndex = True
1016 1016 #print(self.height_index)
1017 1017 data = {}
1018 1018 meta = {}
1019 1019
1020 1020 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
1021 1021 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1022 1022 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1023 1023
1024 1024
1025 1025 z = []
1026 1026 for ch in range(dataOut.nChannels):
1027 1027 if hasattr(dataOut.normFactor,'shape'):
1028 1028 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1029 1029 else:
1030 1030 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1031 1031
1032 1032 z = numpy.asarray(z)
1033 1033 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1034 1034 spc = 10*numpy.log10(z)
1035 1035
1036 1036
1037 1037 data['spc'] = spc - noise
1038 1038 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1039 1039
1040 1040 return data, meta
1041 1041
1042 1042 def plot(self):
1043 1043 if self.xaxis == "frequency":
1044 1044 x = self.data.xrange[0][0:]
1045 1045 self.xlabel = "Frequency (kHz)"
1046 1046 elif self.xaxis == "time":
1047 1047 x = self.data.xrange[1]
1048 1048 self.xlabel = "Time (ms)"
1049 1049 else:
1050 1050 x = self.data.xrange[2]
1051 1051 self.xlabel = "Velocity (m/s)"
1052 1052
1053 1053 self.titles = []
1054 1054
1055 1055 y = self.data.yrange
1056 1056 z = self.data[-1]['spc']
1057 1057 #print(z.shape)
1058 1058 if len(self.height_index) > 0:
1059 1059 index = self.height_index
1060 1060 else:
1061 1061 index = numpy.arange(0, len(y), int((len(y))/9))
1062 1062 #print("inde x ", index, self.axes)
1063 1063
1064 1064 for n, ax in enumerate(self.axes):
1065 1065
1066 1066 if ax.firsttime:
1067 1067
1068 1068
1069 1069 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1070 1070 self.xmin = self.xmin if self.xmin else -self.xmax
1071 1071 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
1072 1072 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
1073 1073
1074 1074
1075 1075 ax.plt = ax.plot(x, z[n, :, index].T)
1076 1076 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1077 1077 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
1078 1078 ax.minorticks_on()
1079 1079 ax.grid(which='major', axis='both')
1080 1080 ax.grid(which='minor', axis='x')
1081 1081 else:
1082 1082 for i, line in enumerate(ax.plt):
1083 1083 line.set_data(x, z[n, :, index[i]])
1084 1084
1085 1085
1086 1086 self.titles.append('CH {}'.format(self.channelList[n]))
1087 1087 plt.suptitle(self.maintitle, fontsize=10)
1088 1088
1089 1089
1090 1090 class BeaconPhase(Plot):
1091 1091
1092 1092 __isConfig = None
1093 1093 __nsubplots = None
1094 1094
1095 1095 PREFIX = 'beacon_phase'
1096 1096
1097 1097 def __init__(self):
1098 1098 Plot.__init__(self)
1099 1099 self.timerange = 24*60*60
1100 1100 self.isConfig = False
1101 1101 self.__nsubplots = 1
1102 1102 self.counter_imagwr = 0
1103 1103 self.WIDTH = 800
1104 1104 self.HEIGHT = 400
1105 1105 self.WIDTHPROF = 120
1106 1106 self.HEIGHTPROF = 0
1107 1107 self.xdata = None
1108 1108 self.ydata = None
1109 1109
1110 1110 self.PLOT_CODE = BEACON_CODE
1111 1111
1112 1112 self.FTP_WEI = None
1113 1113 self.EXP_CODE = None
1114 1114 self.SUB_EXP_CODE = None
1115 1115 self.PLOT_POS = None
1116 1116
1117 1117 self.filename_phase = None
1118 1118
1119 1119 self.figfile = None
1120 1120
1121 1121 self.xmin = None
1122 1122 self.xmax = None
1123 1123
1124 1124 def getSubplots(self):
1125 1125
1126 1126 ncol = 1
1127 1127 nrow = 1
1128 1128
1129 1129 return nrow, ncol
1130 1130
1131 1131 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1132 1132
1133 1133 self.__showprofile = showprofile
1134 1134 self.nplots = nplots
1135 1135
1136 1136 ncolspan = 7
1137 1137 colspan = 6
1138 1138 self.__nsubplots = 2
1139 1139
1140 1140 self.createFigure(id = id,
1141 1141 wintitle = wintitle,
1142 1142 widthplot = self.WIDTH+self.WIDTHPROF,
1143 1143 heightplot = self.HEIGHT+self.HEIGHTPROF,
1144 1144 show=show)
1145 1145
1146 1146 nrow, ncol = self.getSubplots()
1147 1147
1148 1148 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1149 1149
1150 1150 def save_phase(self, filename_phase):
1151 1151 f = open(filename_phase,'w+')
1152 1152 f.write('\n\n')
1153 1153 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1154 1154 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1155 1155 f.close()
1156 1156
1157 1157 def save_data(self, filename_phase, data, data_datetime):
1158 1158 f=open(filename_phase,'a')
1159 1159 timetuple_data = data_datetime.timetuple()
1160 1160 day = str(timetuple_data.tm_mday)
1161 1161 month = str(timetuple_data.tm_mon)
1162 1162 year = str(timetuple_data.tm_year)
1163 1163 hour = str(timetuple_data.tm_hour)
1164 1164 minute = str(timetuple_data.tm_min)
1165 1165 second = str(timetuple_data.tm_sec)
1166 1166 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1167 1167 f.close()
1168 1168
1169 1169 def plot(self):
1170 1170 log.warning('TODO: Not yet implemented...')
1171 1171
1172 1172 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1173 1173 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1174 1174 timerange=None,
1175 1175 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1176 1176 server=None, folder=None, username=None, password=None,
1177 1177 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1178 1178
1179 1179 if dataOut.flagNoData:
1180 1180 return dataOut
1181 1181
1182 1182 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1183 1183 return
1184 1184
1185 1185 if pairsList == None:
1186 1186 pairsIndexList = dataOut.pairsIndexList[:10]
1187 1187 else:
1188 1188 pairsIndexList = []
1189 1189 for pair in pairsList:
1190 1190 if pair not in dataOut.pairsList:
1191 1191 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1192 1192 pairsIndexList.append(dataOut.pairsList.index(pair))
1193 1193
1194 1194 if pairsIndexList == []:
1195 1195 return
1196 1196
1197 1197 # if len(pairsIndexList) > 4:
1198 1198 # pairsIndexList = pairsIndexList[0:4]
1199 1199
1200 1200 hmin_index = None
1201 1201 hmax_index = None
1202 1202
1203 1203 if hmin != None and hmax != None:
1204 1204 indexes = numpy.arange(dataOut.nHeights)
1205 1205 hmin_list = indexes[dataOut.heightList >= hmin]
1206 1206 hmax_list = indexes[dataOut.heightList <= hmax]
1207 1207
1208 1208 if hmin_list.any():
1209 1209 hmin_index = hmin_list[0]
1210 1210
1211 1211 if hmax_list.any():
1212 1212 hmax_index = hmax_list[-1]+1
1213 1213
1214 1214 x = dataOut.getTimeRange()
1215 1215
1216 1216 thisDatetime = dataOut.datatime
1217 1217
1218 1218 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1219 1219 xlabel = "Local Time"
1220 1220 ylabel = "Phase (degrees)"
1221 1221
1222 1222 update_figfile = False
1223 1223
1224 1224 nplots = len(pairsIndexList)
1225 1225 phase_beacon = numpy.zeros(len(pairsIndexList))
1226 1226 for i in range(nplots):
1227 1227 pair = dataOut.pairsList[pairsIndexList[i]]
1228 1228 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1229 1229 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1230 1230 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1231 1231 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1232 1232 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1233 1233
1234 1234 if dataOut.beacon_heiIndexList:
1235 1235 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1236 1236 else:
1237 1237 phase_beacon[i] = numpy.average(phase)
1238 1238
1239 1239 if not self.isConfig:
1240 1240
1241 1241 nplots = len(pairsIndexList)
1242 1242
1243 1243 self.setup(id=id,
1244 1244 nplots=nplots,
1245 1245 wintitle=wintitle,
1246 1246 showprofile=showprofile,
1247 1247 show=show)
1248 1248
1249 1249 if timerange != None:
1250 1250 self.timerange = timerange
1251 1251
1252 1252 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1253 1253
1254 1254 if ymin == None: ymin = 0
1255 1255 if ymax == None: ymax = 360
1256 1256
1257 1257 self.FTP_WEI = ftp_wei
1258 1258 self.EXP_CODE = exp_code
1259 1259 self.SUB_EXP_CODE = sub_exp_code
1260 1260 self.PLOT_POS = plot_pos
1261 1261
1262 1262 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1263 1263 self.isConfig = True
1264 1264 self.figfile = figfile
1265 1265 self.xdata = numpy.array([])
1266 1266 self.ydata = numpy.array([])
1267 1267
1268 1268 update_figfile = True
1269 1269
1270 1270 #open file beacon phase
1271 1271 path = '%s%03d' %(self.PREFIX, self.id)
1272 1272 beacon_file = os.path.join(path,'%s.txt'%self.name)
1273 1273 self.filename_phase = os.path.join(figpath,beacon_file)
1274 1274
1275 1275 self.setWinTitle(title)
1276 1276
1277 1277
1278 1278 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1279 1279
1280 1280 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1281 1281
1282 1282 axes = self.axesList[0]
1283 1283
1284 1284 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1285 1285
1286 1286 if len(self.ydata)==0:
1287 1287 self.ydata = phase_beacon.reshape(-1,1)
1288 1288 else:
1289 1289 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1290 1290
1291 1291
1292 1292 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1293 1293 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1294 1294 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1295 1295 XAxisAsTime=True, grid='both'
1296 1296 )
1297 1297
1298 1298 self.draw()
1299 1299
1300 1300 if dataOut.ltctime >= self.xmax:
1301 1301 self.counter_imagwr = wr_period
1302 1302 self.isConfig = False
1303 1303 update_figfile = True
1304 1304
1305 1305 self.save(figpath=figpath,
1306 1306 figfile=figfile,
1307 1307 save=save,
1308 1308 ftp=ftp,
1309 1309 wr_period=wr_period,
1310 1310 thisDatetime=thisDatetime,
1311 1311 update_figfile=update_figfile)
1312 1312
1313 1313 return dataOut
1314 1314
1315 1315 #####################################
1316 1316 class NoiselessSpectraPlot(Plot):
1317 1317 '''
1318 1318 Plot for Spectra data, subtracting
1319 1319 the noise in all channels, using for
1320 1320 amisr-14 data
1321 1321 '''
1322 1322
1323 1323 CODE = 'noiseless_spc'
1324 1324 colormap = 'jet'
1325 1325 plot_type = 'pcolor'
1326 1326 buffering = False
1327 1327 channelList = []
1328 1328 last_noise = None
1329 1329
1330 1330 def setup(self):
1331 1331
1332 1332 self.nplots = len(self.data.channels)
1333 1333 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
1334 1334 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
1335 1335 self.height = 3.5 * self.nrows
1336 1336
1337 1337 self.cb_label = 'dB'
1338 1338 if self.showprofile:
1339 1339 self.width = 5.8 * self.ncols
1340 1340 else:
1341 1341 self.width = 4.8* self.ncols
1342 1342 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12})
1343 1343
1344 1344 self.ylabel = 'Range [km]'
1345 1345
1346 1346
1347 1347 def update_list(self,dataOut):
1348 1348 if len(self.channelList) == 0:
1349 1349 self.channelList = dataOut.channelList
1350 1350
1351 1351 def update(self, dataOut):
1352 1352
1353 1353 self.update_list(dataOut)
1354 1354 data = {}
1355 1355 meta = {}
1356 1356
1357 1357 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1358 1358 n0 = (dataOut.getNoise()/norm)
1359 1359 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1360 1360 noise = 10*numpy.log10(noise)
1361 1361
1362 1362 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
1363 1363 for ch in range(dataOut.nChannels):
1364 1364 if hasattr(dataOut.normFactor,'ndim'):
1365 1365 if dataOut.normFactor.ndim > 1:
1366 1366 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1367 1367 else:
1368 1368 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1369 1369 else:
1370 1370 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1371 1371
1372 1372 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1373 1373 spc = 10*numpy.log10(z)
1374 1374
1375 1375
1376 1376 data['spc'] = spc - noise
1377 1377 #print(spc.shape)
1378 1378 data['rti'] = spc.mean(axis=1)
1379 1379 data['noise'] = noise
1380 1380
1381 1381
1382 1382
1383 1383 # data['noise'] = noise
1384 1384 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1385 1385
1386 1386 return data, meta
1387 1387
1388 1388 def plot(self):
1389 1389 if self.xaxis == "frequency":
1390 1390 x = self.data.xrange[0]
1391 1391 self.xlabel = "Frequency (kHz)"
1392 1392 elif self.xaxis == "time":
1393 1393 x = self.data.xrange[1]
1394 1394 self.xlabel = "Time (ms)"
1395 1395 else:
1396 1396 x = self.data.xrange[2]
1397 1397 self.xlabel = "Velocity (m/s)"
1398 1398
1399 1399 self.titles = []
1400 1400 y = self.data.yrange
1401 1401 self.y = y
1402 1402
1403 1403 data = self.data[-1]
1404 1404 z = data['spc']
1405 1405
1406 1406 for n, ax in enumerate(self.axes):
1407 1407 #noise = data['noise'][n]
1408 1408
1409 1409 if ax.firsttime:
1410 1410 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1411 1411 self.xmin = self.xmin if self.xmin else -self.xmax
1412 1412 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1413 1413 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1414 1414 ax.plt = ax.pcolormesh(x, y, z[n].T,
1415 1415 vmin=self.zmin,
1416 1416 vmax=self.zmax,
1417 1417 cmap=plt.get_cmap(self.colormap)
1418 1418 )
1419 1419
1420 1420 if self.showprofile:
1421 1421 ax.plt_profile = self.pf_axes[n].plot(
1422 1422 data['rti'][n], y)[0]
1423 1423
1424 1424
1425 1425 else:
1426 1426 ax.plt.set_array(z[n].T.ravel())
1427 1427 if self.showprofile:
1428 1428 ax.plt_profile.set_data(data['rti'][n], y)
1429 1429
1430 1430
1431 1431 self.titles.append('CH {}'.format(self.channelList[n]))
1432 1432
1433 1433
1434 1434 class NoiselessRTIPlot(RTIPlot):
1435 1435 '''
1436 1436 Plot for RTI data
1437 1437 '''
1438 1438
1439 1439 CODE = 'noiseless_rti'
1440 1440 colormap = 'jet'
1441 1441 plot_type = 'pcolorbuffer'
1442 1442 titles = None
1443 1443 channelList = []
1444 1444 elevationList = []
1445 1445 azimuthList = []
1446 1446 last_noise = None
1447 1447
1448 1448 def setup(self):
1449 1449 self.xaxis = 'time'
1450 1450 self.ncols = 1
1451 1451 #print("dataChannels ",self.data.channels)
1452 1452 self.nrows = len(self.data.channels)
1453 1453 self.nplots = len(self.data.channels)
1454 1454 self.ylabel = 'Range [km]'
1455 1455 #self.xlabel = 'Time'
1456 1456 self.cb_label = 'dB'
1457 1457 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1458 1458 self.titles = ['{} Channel {}'.format(
1459 1459 self.CODE.upper(), x) for x in range(self.nplots)]
1460 1460
1461 1461 def update_list(self,dataOut):
1462 1462 if len(self.channelList) == 0:
1463 1463 self.channelList = dataOut.channelList
1464 1464 if len(self.elevationList) == 0:
1465 1465 self.elevationList = dataOut.elevationList
1466 1466 if len(self.azimuthList) == 0:
1467 1467 self.azimuthList = dataOut.azimuthList
1468 1468
1469 1469 def update(self, dataOut):
1470 1470 if len(self.channelList) == 0:
1471 1471 self.update_list(dataOut)
1472 1472
1473 1473 data = {}
1474 1474 meta = {}
1475 1475 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1476 1476 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt
1477 1477 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1478 1478 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1479 1479 data['noise'] = n0
1480 1480 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1481 1481 noiseless_data = dataOut.getPower() - noise
1482 1482
1483 1483 #print("power, noise:", dataOut.getPower(), n0)
1484 1484 #print(noise)
1485 1485 #print(noiseless_data)
1486 1486
1487 1487 data['noiseless_rti'] = noiseless_data
1488 1488
1489 1489 return data, meta
1490 1490
1491 1491 def plot(self):
1492 1492 from matplotlib import pyplot as plt
1493 1493 self.x = self.data.times
1494 1494 self.y = self.data.yrange
1495 1495 self.z = self.data['noiseless_rti']
1496 1496 self.z = numpy.array(self.z, dtype=float)
1497 1497 self.z = numpy.ma.masked_invalid(self.z)
1498 1498
1499 1499
1500 1500 try:
1501 1501 if self.channelList != None:
1502 1502 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1503 1503 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1504 1504 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1505 1505 else:
1506 1506 self.titles = ['{} Channel {}'.format(
1507 1507 self.CODE.upper(), x) for x in self.channelList]
1508 1508 except:
1509 1509 if self.channelList.any() != None:
1510 1510 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1511 1511 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1512 1512 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1513 1513 else:
1514 1514 self.titles = ['{} Channel {}'.format(
1515 1515 self.CODE.upper(), x) for x in self.channelList]
1516 1516
1517 1517
1518 1518 if self.decimation is None:
1519 1519 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1520 1520 else:
1521 1521 x, y, z = self.fill_gaps(*self.decimate())
1522 1522
1523 1523 dummy_var = self.axes #ExtraΓ±amente esto actualiza el valor axes
1524 1524 #print("plot shapes ", z.shape, x.shape, y.shape)
1525 1525 #print(self.axes)
1526 1526 for n, ax in enumerate(self.axes):
1527 1527
1528 1528
1529 1529 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
1530 1530 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1531 1531 data = self.data[-1]
1532 1532 if ax.firsttime:
1533 1533 if (n+1) == len(self.channelList):
1534 1534 ax.set_xlabel('Time')
1535 1535 ax.plt = ax.pcolormesh(x, y, z[n].T,
1536 1536 vmin=self.zmin,
1537 1537 vmax=self.zmax,
1538 1538 cmap=plt.get_cmap(self.colormap)
1539 1539 )
1540 1540 if self.showprofile:
1541 1541 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1542 1542
1543 1543 else:
1544 1544 # ax.collections.remove(ax.collections[0]) # error while running
1545 1545 ax.plt = ax.pcolormesh(x, y, z[n].T,
1546 1546 vmin=self.zmin,
1547 1547 vmax=self.zmax,
1548 1548 cmap=plt.get_cmap(self.colormap)
1549 1549 )
1550 1550 if self.showprofile:
1551 1551 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1552 1552 # if "noise" in self.data:
1553 1553 # #ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
1554 1554 # ax.plot_noise.set_data(data['noise'][n], self.y)
1555 1555
1556 1556
1557 1557 class OutliersRTIPlot(Plot):
1558 1558 '''
1559 1559 Plot for data_xxxx object
1560 1560 '''
1561 1561
1562 1562 CODE = 'outlier_rtc' # Range Time Counts
1563 1563 colormap = 'cool'
1564 1564 plot_type = 'pcolorbuffer'
1565 1565
1566 1566 def setup(self):
1567 1567 self.xaxis = 'time'
1568 1568 self.ncols = 1
1569 1569 self.nrows = self.data.shape('outlier_rtc')[0]
1570 1570 self.nplots = self.nrows
1571 1571 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1572 1572
1573 1573
1574 1574 if not self.xlabel:
1575 1575 self.xlabel = 'Time'
1576 1576
1577 1577 self.ylabel = 'Height [km]'
1578 1578 if not self.titles:
1579 1579 self.titles = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)]
1580 1580
1581 1581 def update(self, dataOut):
1582 1582
1583 1583 data = {}
1584 1584 data['outlier_rtc'] = dataOut.data_outlier
1585 1585
1586 1586 meta = {}
1587 1587
1588 1588 return data, meta
1589 1589
1590 1590 def plot(self):
1591 1591 # self.data.normalize_heights()
1592 1592 self.x = self.data.times
1593 1593 self.y = self.data.yrange
1594 1594 self.z = self.data['outlier_rtc']
1595 1595
1596 1596 #self.z = numpy.ma.masked_invalid(self.z)
1597 1597
1598 1598 if self.decimation is None:
1599 1599 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1600 1600 else:
1601 1601 x, y, z = self.fill_gaps(*self.decimate())
1602 1602
1603 1603 for n, ax in enumerate(self.axes):
1604 1604
1605 1605 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1606 1606 self.z[n])
1607 1607 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1608 1608 self.z[n])
1609 1609 data = self.data[-1]
1610 1610 if ax.firsttime:
1611 1611 if self.zlimits is not None:
1612 1612 self.zmin, self.zmax = self.zlimits[n]
1613 1613
1614 1614 ax.plt = ax.pcolormesh(x, y, z[n].T,
1615 1615 vmin=self.zmin,
1616 1616 vmax=self.zmax,
1617 1617 cmap=self.cmaps[n]
1618 1618 )
1619 1619 if self.showprofile:
1620 1620 ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0]
1621 1621 self.pf_axes[n].set_xlabel('')
1622 1622 else:
1623 1623 if self.zlimits is not None:
1624 1624 self.zmin, self.zmax = self.zlimits[n]
1625 1625 # ax.collections.remove(ax.collections[0]) # error while running
1626 1626 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1627 1627 vmin=self.zmin,
1628 1628 vmax=self.zmax,
1629 1629 cmap=self.cmaps[n]
1630 1630 )
1631 1631 if self.showprofile:
1632 1632 ax.plot_profile.set_data(data['outlier_rtc'][n], self.y)
1633 1633 self.pf_axes[n].set_xlabel('')
1634 1634
1635 1635 class NIncohIntRTIPlot(Plot):
1636 1636 '''
1637 1637 Plot for data_xxxx object
1638 1638 '''
1639 1639
1640 1640 CODE = 'integrations_rtc' # Range Time Counts
1641 1641 colormap = 'BuGn'
1642 1642 plot_type = 'pcolorbuffer'
1643 1643
1644 1644 def setup(self):
1645 1645 self.xaxis = 'time'
1646 1646 self.ncols = 1
1647 1647 self.nrows = self.data.shape('integrations_rtc')[0]
1648 1648 self.nplots = self.nrows
1649 1649 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1650 1650
1651 1651
1652 1652 if not self.xlabel:
1653 1653 self.xlabel = 'Time'
1654 1654
1655 1655 self.ylabel = 'Height [km]'
1656 1656 if not self.titles:
1657 1657 self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)]
1658 1658
1659 1659 def update(self, dataOut):
1660 1660
1661 1661 data = {}
1662 1662 data['integrations_rtc'] = dataOut.nIncohInt
1663 1663
1664 1664 meta = {}
1665 1665
1666 1666 return data, meta
1667 1667
1668 1668 def plot(self):
1669 1669 # self.data.normalize_heights()
1670 1670 self.x = self.data.times
1671 1671 self.y = self.data.yrange
1672 1672 self.z = self.data['integrations_rtc']
1673 1673
1674 1674 #self.z = numpy.ma.masked_invalid(self.z)
1675 1675
1676 1676 if self.decimation is None:
1677 1677 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1678 1678 else:
1679 1679 x, y, z = self.fill_gaps(*self.decimate())
1680 1680
1681 1681 for n, ax in enumerate(self.axes):
1682 1682
1683 1683 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1684 1684 self.z[n])
1685 1685 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1686 1686 self.z[n])
1687 1687 data = self.data[-1]
1688 1688 if ax.firsttime:
1689 1689 if self.zlimits is not None:
1690 1690 self.zmin, self.zmax = self.zlimits[n]
1691 1691
1692 1692 ax.plt = ax.pcolormesh(x, y, z[n].T,
1693 1693 vmin=self.zmin,
1694 1694 vmax=self.zmax,
1695 1695 cmap=self.cmaps[n]
1696 1696 )
1697 1697 if self.showprofile:
1698 1698 ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0]
1699 1699 self.pf_axes[n].set_xlabel('')
1700 1700 else:
1701 1701 if self.zlimits is not None:
1702 1702 self.zmin, self.zmax = self.zlimits[n]
1703 1703 # ax.collections.remove(ax.collections[0]) # error while running
1704 1704 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1705 1705 vmin=self.zmin,
1706 1706 vmax=self.zmax,
1707 1707 cmap=self.cmaps[n]
1708 1708 )
1709 1709 if self.showprofile:
1710 1710 ax.plot_profile.set_data(data['integrations_rtc'][n], self.y)
1711 1711 self.pf_axes[n].set_xlabel('')
1712 1712
1713 1713
1714 1714
1715 1715 class RTIMapPlot(Plot):
1716 1716 '''
1717 1717 Plot for RTI data
1718 1718
1719 1719 Example:
1720 1720
1721 1721 controllerObj = Project()
1722 1722 controllerObj.setup(id = '11', name='eej_proc', description=desc)
1723 1723 ##.......................................................................................
1724 1724 ##.......................................................................................
1725 1725 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', path=inPath, startDate='2023/05/24',endDate='2023/05/24',
1726 1726 startTime='12:00:00',endTime='12:45:59',walk=1,timezone='lt',margin_days=1,code = code,nCode = nCode,
1727 1727 nBaud = nBaud,nOsamp = nosamp,nChannels=nChannels,nFFT=NFFT,
1728 1728 syncronization=False,shiftChannels=0)
1729 1729
1730 1730 volts_proc = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
1731 1731
1732 1732 opObj01 = volts_proc.addOperation(name='Decoder', optype='other')
1733 1733 opObj01.addParameter(name='code', value=code, format='floatlist')
1734 1734 opObj01.addParameter(name='nCode', value=1, format='int')
1735 1735 opObj01.addParameter(name='nBaud', value=nBaud, format='int')
1736 1736 opObj01.addParameter(name='osamp', value=nosamp, format='int')
1737 1737
1738 1738 opObj12 = volts_proc.addOperation(name='selectHeights', optype='self')
1739 1739 opObj12.addParameter(name='minHei', value='90', format='float')
1740 1740 opObj12.addParameter(name='maxHei', value='150', format='float')
1741 1741
1742 1742 proc_spc = controllerObj.addProcUnit(datatype='SpectraProc', inputId=volts_proc.getId())
1743 1743 proc_spc.addParameter(name='nFFTPoints', value='8', format='int')
1744 1744
1745 1745 opObj11 = proc_spc.addOperation(name='IncohInt', optype='other')
1746 1746 opObj11.addParameter(name='n', value='1', format='int')
1747 1747
1748 1748 beamMapFile = "/home/japaza/Documents/AMISR_sky_mapper/UMET_beamcodes.csv"
1749 1749
1750 1750 opObj12 = proc_spc.addOperation(name='RTIMapPlot', optype='external')
1751 1751 opObj12.addParameter(name='selectedHeightsList', value='95, 100, 105, 110 ', format='int')
1752 1752 opObj12.addParameter(name='bField', value='100', format='int')
1753 1753 opObj12.addParameter(name='filename', value=beamMapFile, format='str')
1754 1754
1755 1755 '''
1756 1756
1757 1757 CODE = 'rti_skymap'
1758 1758
1759 1759 plot_type = 'scatter'
1760 1760 titles = None
1761 1761 colormap = 'jet'
1762 1762 channelList = []
1763 1763 elevationList = []
1764 1764 azimuthList = []
1765 1765 last_noise = None
1766 1766 flag_setIndex = False
1767 1767 heights = []
1768 1768 dcosx = []
1769 1769 dcosy = []
1770 1770 fullDcosy = None
1771 1771 fullDcosy = None
1772 1772 hindex = []
1773 1773 mapFile = False
1774 1774 ##### BField ####
1775 1775 flagBField = False
1776 1776 dcosxB = []
1777 1777 dcosyB = []
1778 1778 Bmarker = ['+','*','D','x','s','>','o','^']
1779 1779
1780 1780
1781 1781 def setup(self):
1782 1782
1783 1783 self.xaxis = 'Range (Km)'
1784 1784 if len(self.selectedHeightsList) > 0:
1785 1785 self.nplots = len(self.selectedHeightsList)
1786 1786 else:
1787 1787 self.nplots = 4
1788 1788 self.ncols = int(numpy.ceil(self.nplots/2))
1789 1789 self.nrows = int(numpy.ceil(self.nplots/self.ncols))
1790 1790 self.ylabel = 'dcosy'
1791 1791 self.xlabel = 'dcosx'
1792 1792 self.colorbar = True
1793 1793 self.width = 6 + 4.1*self.nrows
1794 1794 self.height = 3 + 3.5*self.ncols
1795 1795
1796 1796
1797 1797 if self.extFile!=None:
1798 1798 try:
1799 1799 pointings = numpy.genfromtxt(self.extFile, delimiter=',')
1800 1800 full_azi = pointings[:,1]
1801 1801 full_elev = pointings[:,2]
1802 1802 self.fullDcosx = numpy.cos(numpy.radians(full_elev))*numpy.sin(numpy.radians(full_azi))
1803 1803 self.fullDcosy = numpy.cos(numpy.radians(full_elev))*numpy.cos(numpy.radians(full_azi))
1804 1804 mapFile = True
1805 1805 except Exception as e:
1806 1806 self.extFile = None
1807 1807 print(e)
1808 1808
1809 1809
1810 1810 def update_list(self,dataOut):
1811 1811 if len(self.channelList) == 0:
1812 1812 self.channelList = dataOut.channelList
1813 1813 if len(self.elevationList) == 0:
1814 1814 self.elevationList = dataOut.elevationList
1815 1815 if len(self.azimuthList) == 0:
1816 1816 self.azimuthList = dataOut.azimuthList
1817 1817 a = numpy.radians(numpy.asarray(self.azimuthList))
1818 1818 e = numpy.radians(numpy.asarray(self.elevationList))
1819 1819 self.heights = dataOut.heightList
1820 1820 self.dcosx = numpy.cos(e)*numpy.sin(a)
1821 1821 self.dcosy = numpy.cos(e)*numpy.cos(a)
1822 1822
1823 1823 if len(self.bFieldList)>0:
1824 1824 datetObj = datetime.datetime.fromtimestamp(dataOut.utctime)
1825 1825 doy = datetObj.timetuple().tm_yday
1826 1826 year = datetObj.year
1827 1827 # self.dcosxB, self.dcosyB
1828 1828 ObjB = BField(year=year,doy=doy,site=2,heights=self.bFieldList)
1829 1829 [dcos, alpha, nlon, nlat] = ObjB.getBField()
1830 1830
1831 1831 alpha_location = numpy.zeros((nlon,2,len(self.bFieldList)))
1832 1832 for ih in range(len(self.bFieldList)):
1833 1833 alpha_location[:,0,ih] = dcos[:,0,ih,0]
1834 1834 for ilon in numpy.arange(nlon):
1835 1835 myx = (alpha[ilon,:,ih])[::-1]
1836 1836 myy = (dcos[ilon,:,ih,0])[::-1]
1837 1837 tck = splrep(myx,myy,s=0)
1838 1838 mydcosx = splev(ObjB.alpha_i,tck,der=0)
1839 1839
1840 1840 myx = (alpha[ilon,:,ih])[::-1]
1841 1841 myy = (dcos[ilon,:,ih,1])[::-1]
1842 1842 tck = splrep(myx,myy,s=0)
1843 1843 mydcosy = splev(ObjB.alpha_i,tck,der=0)
1844 1844 alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy])
1845 1845 self.dcosxB.append(alpha_location[:,0,ih])
1846 1846 self.dcosyB.append(alpha_location[:,1,ih])
1847 1847 self.flagBField = True
1848 1848
1849 1849 if len(self.celestialList)>0:
1850 1850 #getBField(self.bFieldList, date)
1851 1851 #pass = kwargs.get('celestial', [])
1852 1852 pass
1853 1853
1854 1854
1855 1855 def update(self, dataOut):
1856 1856
1857 1857 if len(self.channelList) == 0:
1858 1858 self.update_list(dataOut)
1859 1859
1860 1860 if not self.flag_setIndex:
1861 1861 if len(self.selectedHeightsList)>0:
1862 1862 for sel_height in self.selectedHeightsList:
1863 1863 index_list = numpy.where(self.heights >= sel_height)
1864 1864 index_list = index_list[0]
1865 1865 self.hindex.append(index_list[0])
1866 1866 self.flag_setIndex = True
1867 1867
1868 1868 data = {}
1869 1869 meta = {}
1870 1870
1871 1871 data['rti_skymap'] = dataOut.getPower()
1872 1872 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1873 1873 noise = 10*numpy.log10(dataOut.getNoise()/norm)
1874 1874 data['noise'] = noise
1875 1875
1876 1876 return data, meta
1877 1877
1878 1878 def plot(self):
1879 1879
1880 1880 self.x = self.dcosx
1881 1881 self.y = self.dcosy
1882 1882 self.z = self.data[-1]['rti_skymap']
1883 1883 self.z = numpy.array(self.z, dtype=float)
1884 1884
1885 1885 if len(self.hindex) > 0:
1886 1886 index = self.hindex
1887 1887 else:
1888 1888 index = numpy.arange(0, len(self.heights), int((len(self.heights))/4.2))
1889 1889
1890 1890 self.titles = ['Height {:.2f} km '.format(self.heights[i])+" " for i in index]
1891 1891 for n, ax in enumerate(self.axes):
1892 1892
1893 1893 if ax.firsttime:
1894 1894
1895 1895 self.xmax = self.xmax if self.xmax else numpy.nanmax(self.x)
1896 1896 self.xmin = self.xmin if self.xmin else numpy.nanmin(self.x)
1897 1897 self.ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
1898 1898 self.ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
1899 1899 self.zmax = self.zmax if self.zmax else numpy.nanmax(self.z)
1900 1900 self.zmin = self.zmin if self.zmin else numpy.nanmin(self.z)
1901 1901
1902 1902 if self.extFile!=None:
1903 1903 ax.scatter(self.fullDcosx, self.fullDcosy, marker="+", s=20)
1904 1904
1905 1905 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1906 1906 s=60, marker="s", vmax = self.zmax)
1907 1907
1908 1908
1909 1909 ax.minorticks_on()
1910 1910 ax.grid(which='major', axis='both')
1911 1911 ax.grid(which='minor', axis='x')
1912 1912
1913 1913 if self.flagBField :
1914 1914
1915 1915 for ih in range(len(self.bFieldList)):
1916 1916 label = str(self.bFieldList[ih]) + ' km'
1917 1917 ax.plot(self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1918 1918 label=label, linestyle='--', ms=4.0,lw=0.5)
1919 1919 handles, labels = ax.get_legend_handles_labels()
1920 1920 a = -0.05
1921 1921 b = 1.15 - 1.19*(self.nrows)
1922 1922 self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field βŠ₯')
1923 1923
1924 1924 else:
1925 1925
1926 1926 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1927 1927 s=80, marker="s", vmax = self.zmax)
1928 1928
1929 1929 if self.flagBField :
1930 1930 for ih in range(len(self.bFieldList)):
1931 1931 ax.plot (self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1932 1932 linestyle='--', ms=4.0,lw=0.5)
1933 1933
1934 1934
1935 1935
@@ -1,779 +1,778
1 1 ''''
2 2 Created on Set 9, 2015
3 3
4 4 @author: roj-idl71 Karim Kuyeng
5 5
6 6 @upgrade: 2021, Joab Apaza
7 7
8 8 '''
9 9
10 10 import os
11 11 import sys
12 12 import glob
13 13 import fnmatch
14 14 import datetime
15 15 import time
16 16 import re
17 17 import h5py
18 18 import numpy
19 19
20 20 try:
21 21 from gevent import sleep
22 22 except:
23 23 from time import sleep
24 24
25 25 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader,ProcessingHeader
26 26 from schainpy.model.data.jrodata import Voltage
27 27 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
28 28 from numpy import imag
29 29 from schainpy.utils import log
30 30
31 31
32 32 class AMISRReader(ProcessingUnit):
33 33 '''
34 34 classdocs
35 35 '''
36 36
37 37 def __init__(self):
38 38 '''
39 39 Constructor
40 40 '''
41 41
42 42 ProcessingUnit.__init__(self)
43 43
44 44 self.set = None
45 45 self.subset = None
46 46 self.extension_file = '.h5'
47 47 self.dtc_str = 'dtc'
48 48 self.dtc_id = 0
49 49 self.status = True
50 50 self.isConfig = False
51 51 self.dirnameList = []
52 52 self.filenameList = []
53 53 self.fileIndex = None
54 54 self.flagNoMoreFiles = False
55 55 self.flagIsNewFile = 0
56 56 self.filename = ''
57 57 self.amisrFilePointer = None
58 58
59 59 self.beamCodeMap = None
60 60 self.azimuthList = []
61 61 self.elevationList = []
62 62 self.dataShape = None
63 63 self.flag_old_beams = False
64 64
65 65 self.flagAsync = False #Use when the experiment has no syncronization
66 66 self.shiftChannels = 0
67 67 self.profileIndex = 0
68 68
69 69
70 70 self.beamCodeByFrame = None
71 71 self.radacTimeByFrame = None
72 72
73 73 self.dataset = None
74 74
75 75 self.__firstFile = True
76 76
77 77 self.buffer = None
78 78
79 79 self.timezone = 'ut'
80 80
81 81 self.__waitForNewFile = 20
82 82 self.__filename_online = None
83 83 #Is really necessary create the output object in the initializer
84 84 self.dataOut = Voltage()
85 85 self.dataOut.error=False
86 86 self.margin_days = 1
87 87 self.flag_ignoreFiles = False #to activate the ignoring Files flag
88 88 self.flag_standby = False # just keep waiting, use when ignoring files
89 89 self.ignStartDateTime=None
90 90 self.ignEndDateTime=None
91 91
92 92 def setup(self,path=None,
93 93 startDate=None,
94 94 endDate=None,
95 95 startTime=None,
96 96 endTime=None,
97 97 walk=True,
98 98 timezone='ut',
99 99 all=0,
100 100 code = 1,
101 101 nCode = 1,
102 102 nBaud = 0,
103 103 nOsamp = 0,
104 104 online=False,
105 105 old_beams=False,
106 106 margin_days=1,
107 107 nFFT = None,
108 108 nChannels = None,
109 109 ignStartDate=None,
110 110 ignEndDate=None,
111 111 ignStartTime=None,
112 112 ignEndTime=None,
113 113 syncronization=True,
114 114 shiftChannels=0
115 115 ):
116 116
117 117
118 118
119 119 self.timezone = timezone
120 120 self.all = all
121 121 self.online = online
122 122 self.flag_old_beams = old_beams
123 123 self.code = code
124 124 self.nCode = int(nCode)
125 125 self.nBaud = int(nBaud)
126 126 self.nOsamp = int(nOsamp)
127 127 self.margin_days = margin_days
128 128 self.__sampleRate = None
129 129 self.flagAsync = not syncronization
130 130 self.shiftChannels = shiftChannels
131 131 self.nFFT = nFFT
132 132 self.nChannels = nChannels
133 133 if ignStartTime!=None and ignEndTime!=None:
134 134 if ignStartDate!=None and ignEndDate!=None:
135 135 self.ignStartDateTime=datetime.datetime.combine(ignStartDate,ignStartTime)
136 136 self.ignEndDateTime=datetime.datetime.combine(ignEndDate,ignEndTime)
137 137 else:
138 138 self.ignStartDateTime=datetime.datetime.combine(startDate,ignStartTime)
139 139 self.ignEndDateTime=datetime.datetime.combine(endDate,ignEndTime)
140 140 self.flag_ignoreFiles = True
141 141
142 142 #self.findFiles()
143 143 if not(online):
144 144 #Busqueda de archivos offline
145 145 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk,)
146 146 else:
147 147 self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
148 148
149 149 if not(self.filenameList):
150 150 raise schainpy.admin.SchainWarning("There is no files into the folder: %s"%(path))
151 151 #sys.exit(0)
152 152 self.dataOut.error = True
153 153
154 154 self.fileIndex = 0
155 155
156 156 self.readNextFile(online)
157 157
158 158 '''
159 159 Add code
160 160 '''
161 161 self.isConfig = True
162 162 # print("Setup Done")
163 163 pass
164 164
165 165
166 166 def readAMISRHeader(self,fp):
167 167
168 168 if self.isConfig and (not self.flagNoMoreFiles):
169 169 newShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
170 170 if self.dataShape != newShape and newShape != None and not self.flag_standby:
171 171 raise schainpy.admin.SchainError("NEW FILE HAS A DIFFERENT SHAPE: ")
172 172 print(self.dataShape,newShape,"\n")
173 173 return 0
174 174 else:
175 175 self.dataShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
176 176
177 177
178 178 header = 'Raw11/Data/RadacHeader'
179 179 if self.nChannels == None:
180 180 expFile = fp['Setup/Experimentfile'][()].decode()
181 181 linesExp = expFile.split("\n")
182 182 a = [line for line in linesExp if "nbeamcodes" in line]
183 183 self.nChannels = int(a[0][11:])
184 184
185 185 if not self.flagAsync: #for experiments with no syncronization
186 186 self.shiftChannels = 0
187 187
188 188
189 189
190 190 self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
191 191
192 192
193 193 if (self.startDate > datetime.date(2021, 7, 15)) or self.flag_old_beams: #Se cambiΓ³ la forma de extracciΓ³n de Apuntes el 17 o forzar con flag de reorganizaciΓ³n
194 194 self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode()
195 195 self.trueBeams = self.beamcodeFile.split("\n")
196 196 self.trueBeams.pop()#remove last
197 197 if self.nFFT == None:
198 198 log.error("FFT or number of repetitions per channels is needed",self.name)
199 199 beams_idx = [k*self.nFFT for k in range(self.nChannels)]
200 200 beams = [self.trueBeams[b] for b in beams_idx]
201 201 self.beamCode = [int(x, 16) for x in beams]
202 202
203 203 if(self.flagAsync and self.shiftChannels == 0):
204 204 initBeam = self.beamCodeByPulse[0, 0]
205 205 self.shiftChannels = numpy.argwhere(self.beamCode ==initBeam)[0,0]
206 206
207 207 else:
208 208 _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes
209 209 self.beamCode = _beamCode[0,:]
210 210
211 211
212 212
213 213
214 214 if self.beamCodeMap == None:
215 215 self.beamCodeMap = fp['Setup/BeamcodeMap']
216 216 for beam in self.beamCode:
217 217 beamAziElev = numpy.where(self.beamCodeMap[:,0]==beam)
218 218 beamAziElev = beamAziElev[0].squeeze()
219 219 self.azimuthList.append(self.beamCodeMap[beamAziElev,1])
220 220 self.elevationList.append(self.beamCodeMap[beamAziElev,2])
221 221 #print("Beamssss: ",self.beamCodeMap[beamAziElev,1],self.beamCodeMap[beamAziElev,2])
222 222 #print(self.beamCode)
223 223 #self.code = fp.get(header+'/Code') # NOT USE FOR THIS
224 224 self.frameCount = fp.get(header+'/FrameCount')# NOT USE FOR THIS
225 225 self.modeGroup = fp.get(header+'/ModeGroup')# NOT USE FOR THIS
226 226 self.nsamplesPulse = fp.get(header+'/NSamplesPulse')# TO GET NSA OR USING DATA FOR THAT
227 227 self.pulseCount = fp.get(header+'/PulseCount')# NOT USE FOR THIS
228 228 self.radacTime = fp.get(header+'/RadacTime')# 1st TIME ON FILE ANDE CALCULATE THE REST WITH IPP*nindexprofile
229 229 self.timeCount = fp.get(header+'/TimeCount')# NOT USE FOR THIS
230 230 self.timeStatus = fp.get(header+'/TimeStatus')# NOT USE FOR THIS
231 231 self.rangeFromFile = fp.get('Raw11/Data/Samples/Range')
232 232 self.frequency = fp.get('Rx/Frequency')
233 233 txAus = fp.get('Raw11/Data/Pulsewidth') #seconds
234 234 self.baud = fp.get('Raw11/Data/TxBaud')
235 235 sampleRate = fp.get('Rx/SampleRate')
236 236 self.__sampleRate = sampleRate[()]
237 237 self.nblocks = self.pulseCount.shape[0] #nblocks
238 238 self.profPerBlockRAW = self.pulseCount.shape[1] #profiles per block in raw data
239 239 self.nprofiles = self.pulseCount.shape[1] #nprofile
240 240 #self.nsa = self.nsamplesPulse[0,0] #ngates
241 241 self.nsa = len(self.rangeFromFile[0])
242 242 self.nchannels = len(self.beamCode)
243 243 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
244 244 #print("IPPS secs: ",self.ippSeconds)
245 245 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
246 246 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
247 247
248 248 #filling radar controller header parameters
249 249 self.__ippKm = self.ippSeconds *.15*1e6 # in km
250 250 #self.__txA = txAus[()]*.15 #(ipp[us]*.15km/1us) in km
251 251 self.__txA = txAus[()] #seconds
252 252 self.__txAKm = self.__txA*1e6*.15
253 253 self.__txB = 0
254 254 nWindows=1
255 255 self.__nSamples = self.nsa
256 256 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
257 257 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
258 258 #print("amisr-ipp:",self.ippSeconds, self.__ippKm)
259 259 #for now until understand why the code saved is different (code included even though code not in tuf file)
260 260 #self.__codeType = 0
261 261 # self.__nCode = None
262 262 # self.__nBaud = None
263 263 self.__code = self.code
264 264 self.__codeType = 0
265 265 if self.code != None:
266 266 self.__codeType = 1
267 267 self.__nCode = self.nCode
268 268 self.__nBaud = self.nBaud
269 269 self.__baudTX = self.__txA/(self.nBaud)
270 270 #self.__code = 0
271 271
272 272 #filling system header parameters
273 273 self.__nSamples = self.nsa
274 274 self.newProfiles = self.nprofiles/self.nchannels
275 275 self.__channelList = [n for n in range(self.nchannels)]
276 276
277 277 self.__frequency = self.frequency[0][0]
278 278
279 279
280 280 return 1
281 281
282 282
283 283 def createBuffers(self):
284 284
285 285 pass
286 286
287 287 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
288 288 self.path = path
289 289 self.startDate = startDate
290 290 self.endDate = endDate
291 291 self.startTime = startTime
292 292 self.endTime = endTime
293 293 self.walk = walk
294 294
295 295
296 296 def __checkPath(self):
297 297 if os.path.exists(self.path):
298 298 self.status = 1
299 299 else:
300 300 self.status = 0
301 301 print('Path:%s does not exists'%self.path)
302 302
303 303 return
304 304
305 305
306 306 def __selDates(self, amisr_dirname_format):
307 307 try:
308 308 year = int(amisr_dirname_format[0:4])
309 309 month = int(amisr_dirname_format[4:6])
310 310 dom = int(amisr_dirname_format[6:8])
311 311 thisDate = datetime.date(year,month,dom)
312 312 #margen de un dΓ­a extra, igual luego se filtra for fecha y hora
313 313 if (thisDate>=(self.startDate - datetime.timedelta(days=self.margin_days)) and thisDate <= (self.endDate)+ datetime.timedelta(days=1)):
314 314 return amisr_dirname_format
315 315 except:
316 316 return None
317 317
318 318
319 319 def __findDataForDates(self,online=False):
320 320
321 321 if not(self.status):
322 322 return None
323 323
324 324 pat = '\d+.\d+'
325 325 dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
326 326 dirnameList = [x for x in dirnameList if x!=None]
327 327 dirnameList = [x.string for x in dirnameList]
328 328 if not(online):
329 329 dirnameList = [self.__selDates(x) for x in dirnameList]
330 330 dirnameList = [x for x in dirnameList if x!=None]
331 331 if len(dirnameList)>0:
332 332 self.status = 1
333 333 self.dirnameList = dirnameList
334 334 self.dirnameList.sort()
335 335 else:
336 336 self.status = 0
337 337 return None
338 338
339 339 def __getTimeFromData(self):
340 340 startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
341 341 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
342 342
343 343 print('Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader))
344 344 print('........................................')
345 345 filter_filenameList = []
346 346 self.filenameList.sort()
347 347 total_files = len(self.filenameList)
348 348 #for i in range(len(self.filenameList)-1):
349 349 for i in range(total_files):
350 350 filename = self.filenameList[i]
351 351 #print("file-> ",filename)
352 352 try:
353 353 fp = h5py.File(filename,'r')
354 354 time_str = fp.get('Time/RadacTimeString')
355 355
356 356 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
357 357 #startDateTimeStr_File = "2019-12-16 09:21:11"
358 358 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
359 359 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
360 360
361 361 #endDateTimeStr_File = "2019-12-16 11:10:11"
362 362 endDateTimeStr_File = time_str[-1][-1].decode('UTF-8').split('.')[0]
363 363 junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
364 364 endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
365 365
366 366 fp.close()
367 367
368 368 #print("check time", startDateTime_File)
369 369 if self.timezone == 'lt':
370 370 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
371 371 endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
372 372 if (startDateTime_File >=startDateTime_Reader and endDateTime_File<=endDateTime_Reader):
373 373 filter_filenameList.append(filename)
374 374
375 375 if (startDateTime_File>endDateTime_Reader):
376 376 break
377 377 except Exception as e:
378 378 log.warning("Error opening file {} -> {}".format(os.path.split(filename)[1],e))
379 379
380 380 filter_filenameList.sort()
381 381 self.filenameList = filter_filenameList
382 382
383 383 return 1
384 384
385 385 def __filterByGlob1(self, dirName):
386 386 filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
387 387 filter_files.sort()
388 388 filterDict = {}
389 389 filterDict.setdefault(dirName)
390 390 filterDict[dirName] = filter_files
391 391 return filterDict
392 392
393 393 def __getFilenameList(self, fileListInKeys, dirList):
394 394 for value in fileListInKeys:
395 395 dirName = list(value.keys())[0]
396 396 for file in value[dirName]:
397 397 filename = os.path.join(dirName, file)
398 398 self.filenameList.append(filename)
399 399
400 400
401 401 def __selectDataForTimes(self, online=False):
402 402 #aun no esta implementado el filtro for tiempo-> implementado en readNextFile
403 403 if not(self.status):
404 404 return None
405 405
406 406 dirList = [os.path.join(self.path,x) for x in self.dirnameList]
407 407 fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
408 408 self.__getFilenameList(fileListInKeys, dirList)
409 409 if not(online):
410 410 #filtro por tiempo
411 411 if not(self.all):
412 412 self.__getTimeFromData()
413 413
414 414 if len(self.filenameList)>0:
415 415 self.status = 1
416 416 self.filenameList.sort()
417 417 else:
418 418 self.status = 0
419 419 return None
420 420
421 421 else:
422 422 #get the last file - 1
423 423 self.filenameList = [self.filenameList[-2]]
424 424 new_dirnameList = []
425 425 for dirname in self.dirnameList:
426 426 junk = numpy.array([dirname in x for x in self.filenameList])
427 427 junk_sum = junk.sum()
428 428 if junk_sum > 0:
429 429 new_dirnameList.append(dirname)
430 430 self.dirnameList = new_dirnameList
431 431 return 1
432 432
433 433 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
434 434 endTime=datetime.time(23,59,59),walk=True):
435 435
436 436 if endDate ==None:
437 437 startDate = datetime.datetime.utcnow().date()
438 438 endDate = datetime.datetime.utcnow().date()
439 439
440 440 self.__setParameters(path=path, startDate=startDate, endDate=endDate,startTime = startTime,endTime=endTime, walk=walk)
441 441
442 442 self.__checkPath()
443 443
444 444 self.__findDataForDates(online=True)
445 445
446 446 self.dirnameList = [self.dirnameList[-1]]
447 447
448 448 self.__selectDataForTimes(online=True)
449 449
450 450 return
451 451
452 452
453 453 def searchFilesOffLine(self,
454 454 path,
455 455 startDate,
456 456 endDate,
457 457 startTime=datetime.time(0,0,0),
458 458 endTime=datetime.time(23,59,59),
459 459 walk=True):
460 460
461 461 self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
462 462
463 463 self.__checkPath()
464 464
465 465 self.__findDataForDates()
466 466
467 467 self.__selectDataForTimes()
468 468
469 469 for i in range(len(self.filenameList)):
470 470 print("%s" %(self.filenameList[i]))
471 471
472 472 return
473 473
474 474 def __setNextFileOffline(self):
475 475
476 476 try:
477 477 self.filename = self.filenameList[self.fileIndex]
478 478 self.amisrFilePointer = h5py.File(self.filename,'r')
479 479 self.fileIndex += 1
480 480 except:
481 481 self.flagNoMoreFiles = 1
482 482 raise schainpy.admin.SchainError('No more files to read')
483 483 return 0
484 484
485 485 self.flagIsNewFile = 1
486 486 print("Setting the file: %s"%self.filename)
487 487
488 488 return 1
489 489
490 490
491 491 def __setNextFileOnline(self):
492 492 filename = self.filenameList[0]
493 493 if self.__filename_online != None:
494 494 self.__selectDataForTimes(online=True)
495 495 filename = self.filenameList[0]
496 496 wait = 0
497 497 self.__waitForNewFile=300 ## DEBUG:
498 498 while self.__filename_online == filename:
499 499 print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
500 500 if wait == 5:
501 501 self.flagNoMoreFiles = 1
502 502 return 0
503 503 sleep(self.__waitForNewFile)
504 504 self.__selectDataForTimes(online=True)
505 505 filename = self.filenameList[0]
506 506 wait += 1
507 507
508 508 self.__filename_online = filename
509 509
510 510 self.amisrFilePointer = h5py.File(filename,'r')
511 511 self.flagIsNewFile = 1
512 512 self.filename = filename
513 513 print("Setting the file: %s"%self.filename)
514 514 return 1
515 515
516 516
517 517 def readData(self):
518 518 buffer = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
519 519 re = buffer[:,:,:,0]
520 520 im = buffer[:,:,:,1]
521 521 dataset = re + im*1j
522 522
523 523 self.radacTime = self.amisrFilePointer.get('Raw11/Data/RadacHeader/RadacTime')
524 524 timeset = self.radacTime[:,0]
525 525
526 526 return dataset,timeset
527 527
528 528 def reshapeData(self):
529 529 #print(self.beamCodeByPulse, self.beamCode, self.nblocks, self.nprofiles, self.nsa)
530 530 channels = self.beamCodeByPulse[0,:]
531 531 nchan = self.nchannels
532 532 #self.newProfiles = self.nprofiles/nchan #must be defined on filljroheader
533 533 nblocks = self.nblocks
534 534 nsamples = self.nsa
535 535 #print("Channels: ",self.nChannels)
536 536 #Dimensions : nChannels, nProfiles, nSamples
537 537 new_block = numpy.empty((nblocks, nchan, numpy.int_(self.newProfiles), nsamples), dtype="complex64")
538 538 ############################################
539 539 profPerCH = int(self.profPerBlockRAW / (self.nFFT* self.nChannels))
540 540 #profPerCH = int(self.profPerBlockRAW / self.nChannels)
541 541 for thisChannel in range(nchan):
542 542
543 543 ich = thisChannel
544 544
545 545 idx_ch = [self.nFFT*(ich + nchan*k) for k in range(profPerCH)]
546 546 #print(idx_ch)
547 547 if self.nFFT > 1:
548 548 aux = [numpy.arange(i, i+self.nFFT) for i in idx_ch]
549 549 idx_ch = None
550 550 idx_ch =aux
551 551 idx_ch = numpy.array(idx_ch, dtype=int).flatten()
552 552 else:
553 553 idx_ch = numpy.array(idx_ch, dtype=int)
554 554
555 555 #print(ich,profPerCH,idx_ch)
556 556 #print(numpy.where(channels==self.beamCode[ich])[0])
557 557 #new_block[:,ich,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[ich])[0],:]
558 558 new_block[:,ich,:,:] = self.dataset[:,idx_ch,:]
559 559
560 560 new_block = numpy.transpose(new_block, (1,0,2,3))
561 561 new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
562 562 if self.flagAsync:
563 563 new_block = numpy.roll(new_block, self.shiftChannels, axis=0)
564 564 return new_block
565 565
566 566 def updateIndexes(self):
567 567
568 568 pass
569 569
570 570 def fillJROHeader(self):
571 571
572 572 #fill radar controller header
573 573
574 574 #fill system header
575 575 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
576 576 nProfiles=self.newProfiles,
577 577 nChannels=len(self.__channelList),
578 578 adcResolution=14,
579 579 pciDioBusWidth=32)
580 580
581 581 self.dataOut.type = "Voltage"
582 582 self.dataOut.data = None
583 583 self.dataOut.dtype = numpy.dtype([('real','<i8'),('imag','<i8')])
584 584 # self.dataOut.nChannels = 0
585 585
586 586 # self.dataOut.nHeights = 0
587 587
588 588 self.dataOut.nProfiles = self.newProfiles*self.nblocks
589 589 #self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
590 590 ranges = numpy.reshape(self.rangeFromFile[()],(-1))
591 591 self.dataOut.heightList = ranges/1000.0 #km
592 592 self.dataOut.channelList = self.__channelList
593 593
594 594 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
595 595
596 596 # self.dataOut.channelIndexList = None
597 597
598 598
599 599 # #self.dataOut.azimuthList = numpy.roll( numpy.array(self.azimuthList) ,self.shiftChannels)
600 600 # #self.dataOut.elevationList = numpy.roll(numpy.array(self.elevationList) ,self.shiftChannels)
601 601 # #self.dataOut.codeList = numpy.roll(numpy.array(self.beamCode), self.shiftChannels)
602 602
603 603 self.dataOut.azimuthList = self.azimuthList
604 604 self.dataOut.elevationList = self.elevationList
605 605 self.dataOut.codeList = self.beamCode
606 606
607 607
608 608
609 609 #print(self.dataOut.elevationList)
610 610 self.dataOut.flagNoData = True
611 611
612 612 #Set to TRUE if the data is discontinuous
613 613 self.dataOut.flagDiscontinuousBlock = False
614 614
615 615 self.dataOut.utctime = None
616 616
617 617 #self.dataOut.timeZone = -5 #self.__timezone/60 #timezone like jroheader, difference in minutes between UTC and localtime
618 618 if self.timezone == 'lt':
619 619 self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
620 620 else:
621 621 self.dataOut.timeZone = 0 #by default time is UTC
622 622
623 623 self.dataOut.dstFlag = 0
624 624 self.dataOut.errorCount = 0
625 625 self.dataOut.nCohInt = 1
626 626 self.dataOut.flagDecodeData = False #asumo que la data esta decodificada
627 627 self.dataOut.flagDeflipData = False #asumo que la data esta sin flip
628 628 self.dataOut.flagShiftFFT = False
629 629 self.dataOut.ippSeconds = self.ippSeconds
630 630 self.dataOut.ipp = self.__ippKm
631 631 self.dataOut.nCode = self.__nCode
632 632 self.dataOut.code = self.__code
633 633 self.dataOut.nBaud = self.__nBaud
634 634
635 635
636 636 self.dataOut.frequency = self.__frequency
637 637 self.dataOut.realtime = self.online
638 638
639 639 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(ipp=self.__ippKm,
640 640 txA=self.__txAKm,
641 641 txB=0,
642 642 nWindows=1,
643 643 nHeights=self.__nSamples,
644 644 firstHeight=self.__firstHeight,
645 645 codeType=self.__codeType,
646 646 nCode=self.__nCode, nBaud=self.__nBaud,
647 647 code = self.__code,
648 648 nOsamp=self.nOsamp,
649 649 frequency = self.__frequency,
650 650 sampleRate= self.__sampleRate,
651 651 fClock=self.__sampleRate)
652 652
653 653
654 654 self.dataOut.radarControllerHeaderObj.heightList = ranges/1000.0 #km
655 655 self.dataOut.radarControllerHeaderObj.heightResolution = self.__deltaHeight
656 656 self.dataOut.radarControllerHeaderObj.rangeIpp = self.__ippKm #km
657 657 self.dataOut.radarControllerHeaderObj.rangeTxA = self.__txA*1e6*.15 #km
658 658 self.dataOut.radarControllerHeaderObj.nChannels = self.nchannels
659 659 self.dataOut.radarControllerHeaderObj.channelList = self.__channelList
660 660 self.dataOut.radarControllerHeaderObj.azimuthList = self.azimuthList
661 661 self.dataOut.radarControllerHeaderObj.elevationList = self.elevationList
662 662 self.dataOut.radarControllerHeaderObj.dtype = "Voltage"
663 663 self.dataOut.ippSeconds = self.ippSeconds
664 664 self.dataOut.ippFactor = self.nFFT
665 665 pass
666 666
667 667 def readNextFile(self,online=False):
668 668
669 669 if not(online):
670 670 newFile = self.__setNextFileOffline()
671 671 else:
672 672 newFile = self.__setNextFileOnline()
673 673
674 674 if not(newFile):
675 675 self.dataOut.error = True
676 676 return 0
677 677
678 678 if not self.readAMISRHeader(self.amisrFilePointer):
679 679 self.dataOut.error = True
680 680 return 0
681 681
682 682 #self.createBuffers()
683 683 self.fillJROHeader()
684 684
685 685 #self.__firstFile = False
686 686
687 687 self.dataset,self.timeset = self.readData()
688
689 688 if self.endDate!=None:
690 689 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
691 690 time_str = self.amisrFilePointer.get('Time/RadacTimeString')
692 691 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
693 692 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
694 693 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
695 694 if self.timezone == 'lt':
696 695 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
697 696 if (startDateTime_File>endDateTime_Reader):
698 697 self.flag_standby = False
699 698 return 0
700 699 if self.flag_ignoreFiles and (startDateTime_File >= self.ignStartDateTime and startDateTime_File <= self.ignEndDateTime):
701 700 print("Ignoring...")
702 701 self.flag_standby = True
702 self.profileIndex = 99999999999999999
703 703 return 1
704 704 self.flag_standby = False
705 705
706 706 self.jrodataset = self.reshapeData()
707 707 #----self.updateIndexes()
708 708 self.profileIndex = 0
709 709
710 710 return 1
711 711
712 712
713 713 def __hasNotDataInBuffer(self):
714 714 if self.profileIndex >= (self.newProfiles*self.nblocks):
715 715 return 1
716 716 return 0
717 717
718 718
719 719 def getData(self):
720 720
721 721 if self.flagNoMoreFiles:
722 722 self.dataOut.flagNoData = True
723 723 return 0
724
725 724 if self.profileIndex >= (self.newProfiles*self.nblocks): #
726 725 #if self.__hasNotDataInBuffer():
727 726 if not (self.readNextFile(self.online)):
728 727 print("Profile Index break...")
729 728 return 0
730 729
731 730 if self.flag_standby: #Standby mode, if files are being ignoring, just return with no error flag
732 731 return 0
733 732
734 733 if self.dataset is None: # setear esta condicion cuando no hayan datos por leer
735 734 self.dataOut.flagNoData = True
736 735 print("No more data break...")
737 736 return 0
738 737
739 738 #self.dataOut.data = numpy.reshape(self.jrodataset[self.profileIndex,:],(1,-1))
740 739
741 740 self.dataOut.data = self.jrodataset[:,self.profileIndex,:]
742 741
743 742 #print("R_t",self.timeset)
744 743
745 744 #self.dataOut.utctime = self.jrotimeset[self.profileIndex]
746 745 #verificar basic header de jro data y ver si es compatible con este valor
747 746 #self.dataOut.utctime = self.timeset + (self.profileIndex * self.ippSeconds * self.nchannels)
748 747 indexprof = numpy.mod(self.profileIndex, self.newProfiles)
749 748 indexblock = self.profileIndex/self.newProfiles
750 749 #print (indexblock, indexprof)
751 750 diffUTC = 0
752 751 t_comp = (indexprof * self.ippSeconds * self.nchannels) + diffUTC #
753 752
754 753 #print("utc :",indexblock," __ ",t_comp)
755 754 #print(numpy.shape(self.timeset))
756 755 self.dataOut.utctime = self.timeset[numpy.int_(indexblock)] + t_comp
757 756 #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
758 757
759 758 self.dataOut.profileIndex = self.profileIndex
760 759 #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
761 760 self.dataOut.flagNoData = False
762 761 # if indexprof == 0:
763 762 # print("kamisr: ",self.dataOut.utctime)
764 763
765 764 self.profileIndex += 1
766
765
767 766 return self.dataOut.data #retorno necesario??
768 767
769 768
770 769 def run(self, **kwargs):
771 770 '''
772 771 This method will be called many times so here you should put all your code
773 772 '''
774 773 #print("running kamisr")
775 774 if not self.isConfig:
776 775 self.setup(**kwargs)
777 776 self.isConfig = True
778 777
779 778 self.getData()
@@ -1,1735 +1,1737
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.model.data import _noise
20 20 from schainpy.utils import log
21 21 import matplotlib.pyplot as plt
22 22 from schainpy.model.io.utilsIO import getHei_index
23 23 import datetime
24 24
25 25 class SpectraProc(ProcessingUnit):
26 26
27 27 def __init__(self):
28 28
29 29 ProcessingUnit.__init__(self)
30 30
31 31 self.buffer = None
32 32 self.firstdatatime = None
33 33 self.profIndex = 0
34 34 self.dataOut = Spectra()
35 self.dataOut.error=False
35 36 self.id_min = None
36 37 self.id_max = None
37 38 self.setupReq = False #Agregar a todas las unidades de proc
38 39 self.nsamplesFFT = 0
39 40
40 41 def __updateSpecFromVoltage(self):
41 42
42 43 self.dataOut.timeZone = self.dataIn.timeZone
43 44 self.dataOut.dstFlag = self.dataIn.dstFlag
44 45 self.dataOut.errorCount = self.dataIn.errorCount
45 46 self.dataOut.useLocalTime = self.dataIn.useLocalTime
46 47 try:
47 48 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
48 49 except:
49 50 pass
50 51 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
51 52 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
52 53 self.dataOut.ippSeconds = self.dataIn.ippSeconds
53 54 self.dataOut.ipp = self.dataIn.ipp
54 55 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
55 56 self.dataOut.channelList = self.dataIn.channelList
56 57 self.dataOut.heightList = self.dataIn.heightList
57 58 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
58 59 self.dataOut.nProfiles = self.dataOut.nFFTPoints
59 60 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
60 61 self.dataOut.utctime = self.firstdatatime
61 62 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
62 63 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
63 64 self.dataOut.flagShiftFFT = False
64 65 self.dataOut.nCohInt = self.dataIn.nCohInt
65 66 self.dataOut.nIncohInt = 1
66 67 self.dataOut.deltaHeight = self.dataIn.deltaHeight
67 68 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
68 69 self.dataOut.frequency = self.dataIn.frequency
69 70 self.dataOut.realtime = self.dataIn.realtime
70 71 self.dataOut.azimuth = self.dataIn.azimuth
71 72 self.dataOut.zenith = self.dataIn.zenith
72 73 self.dataOut.codeList = self.dataIn.codeList
73 74 self.dataOut.azimuthList = self.dataIn.azimuthList
74 75 self.dataOut.elevationList = self.dataIn.elevationList
75 76 self.dataOut.code = self.dataIn.code
76 77 self.dataOut.nCode = self.dataIn.nCode
77 78 self.dataOut.flagProfilesByRange = self.dataIn.flagProfilesByRange
78 79 self.dataOut.nProfilesByRange = self.dataIn.nProfilesByRange
79 80 self.dataOut.runNextUnit = self.dataIn.runNextUnit
80 81 try:
81 82 self.dataOut.step = self.dataIn.step
82 83 except:
83 84 pass
84 85
85 86 def __getFft(self):
86 87 """
87 88 Convierte valores de Voltaje a Spectra
88 89
89 90 Affected:
90 91 self.dataOut.data_spc
91 92 self.dataOut.data_cspc
92 93 self.dataOut.data_dc
93 94 self.dataOut.heightList
94 95 self.profIndex
95 96 self.buffer
96 97 self.dataOut.flagNoData
97 98 """
98 99 fft_volt = numpy.fft.fft(
99 100 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
100 101 fft_volt = fft_volt.astype(numpy.dtype('complex'))
101 102 dc = fft_volt[:, 0, :]
102 103
103 104 # calculo de self-spectra
104 105 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
105 106 spc = fft_volt * numpy.conjugate(fft_volt)
106 107 spc = spc.real
107 108
108 109 blocksize = 0
109 110 blocksize += dc.size
110 111 blocksize += spc.size
111 112
112 113 cspc = None
113 114 pairIndex = 0
114 115 if self.dataOut.pairsList != None:
115 116 # calculo de cross-spectra
116 117 cspc = numpy.zeros(
117 118 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
118 119 for pair in self.dataOut.pairsList:
119 120 if pair[0] not in self.dataOut.channelList:
120 121 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
121 122 str(pair), str(self.dataOut.channelList)))
122 123 if pair[1] not in self.dataOut.channelList:
123 124 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
124 125 str(pair), str(self.dataOut.channelList)))
125 126
126 127 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
127 128 numpy.conjugate(fft_volt[pair[1], :, :])
128 129 pairIndex += 1
129 130 blocksize += cspc.size
130 131
131 132 self.dataOut.data_spc = spc
132 133 self.dataOut.data_cspc = cspc
133 134 self.dataOut.data_dc = dc
134 135 self.dataOut.blockSize = blocksize
135 136 self.dataOut.flagShiftFFT = False
136 137
137 138 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False,
138 zeroPad=False, zeroPoints=0, runNextUnit = 0):
139
139 zeroPad=False, zeroPoints=0, runNextUnit=0):
140 140 self.dataIn.runNextUnit = runNextUnit
141 141 try:
142 142 type = self.dataIn.type.decode("utf-8")
143 143 self.dataIn.type = type
144 except:
144 except Exception as e:
145 # print("spc -> ",e)
145 146 pass
146 147
147 148 if self.dataIn.type == "Spectra":
149 #print("AQUI")
148 150 try:
149 151 self.dataOut.copy(self.dataIn)
150 152 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
151 153 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
152 154 self.dataOut.nProfiles = self.dataOut.nFFTPoints
153 155 #self.dataOut.nHeights = len(self.dataOut.heightList)
154 156 except Exception as e:
155 157 print("Error dataIn ",e)
156 158
157 159 if shift_fft:
158 160 #desplaza a la derecha en el eje 2 determinadas posiciones
159 161 shift = int(self.dataOut.nFFTPoints/2)
160 162 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
161 163
162 164 if self.dataOut.data_cspc is not None:
163 165 #desplaza a la derecha en el eje 2 determinadas posiciones
164 166 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
165 167 if pairsList:
166 168 self.__selectPairs(pairsList)
167 169
168 170 elif self.dataIn.type == "Voltage":
169 171
170 172 self.dataOut.flagNoData = True
171 173 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
172 174 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
173 175
174 176 if nFFTPoints == None:
175 177 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
176 178
177 179 if nProfiles == None:
178 180 nProfiles = nFFTPoints
179 181
180 182 if ippFactor == None:
181 183 self.dataOut.ippFactor = self.dataIn.ippFactor
182 184 else:
183 185 self.dataOut.ippFactor = ippFactor
184 186
185 187 if self.buffer is None:
186 188 if not zeroPad:
187 189 self.buffer = numpy.zeros((self.dataIn.nChannels,
188 190 nProfiles,
189 191 self.dataIn.nHeights),
190 192 dtype='complex')
191 193 zeroPoints = 0
192 194 else:
193 195 self.buffer = numpy.zeros((self.dataIn.nChannels,
194 196 nFFTPoints+int(zeroPoints),
195 197 self.dataIn.nHeights),
196 198 dtype='complex')
197 199
198 200 self.dataOut.nFFTPoints = nFFTPoints
199 201
200 202 if self.buffer is None:
201 203 self.buffer = numpy.zeros((self.dataIn.nChannels,
202 204 nProfiles,
203 205 self.dataIn.nHeights),
204 206 dtype='complex')
205 207
206 208 if self.dataIn.flagDataAsBlock:
207 209 nVoltProfiles = self.dataIn.data.shape[1]
208 210 zeroPoints = 0
209 211 if nVoltProfiles == nProfiles or zeroPad:
210 212 self.buffer = self.dataIn.data.copy()
211 213 self.profIndex = nVoltProfiles
212 214
213 215 elif nVoltProfiles < nProfiles:
214 216
215 217 if self.profIndex == 0:
216 218 self.id_min = 0
217 219 self.id_max = nVoltProfiles
218 220
219 221 self.buffer[:, self.id_min:self.id_max,
220 222 :] = self.dataIn.data
221 223 self.profIndex += nVoltProfiles
222 224 self.id_min += nVoltProfiles
223 225 self.id_max += nVoltProfiles
224 226 elif nVoltProfiles > nProfiles:
225 227 self.reader.bypass = True
226 228 if self.profIndex == 0:
227 229 self.id_min = 0
228 230 self.id_max = nProfiles
229 231
230 232 self.buffer = self.dataIn.data[:, self.id_min:self.id_max,:]
231 233 self.profIndex += nProfiles
232 234 self.id_min += nProfiles
233 235 self.id_max += nProfiles
234 236 if self.id_max == nVoltProfiles:
235 237 self.reader.bypass = False
236 238
237 239 else:
238 240 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
239 241 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
240 242 self.dataOut.flagNoData = True
241 243 else:
242 244 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
243 245 self.profIndex += 1
244 246
245 247 if self.firstdatatime == None:
246 248 self.firstdatatime = self.dataIn.utctime
247 249
248 250 if self.profIndex == nProfiles or (zeroPad and zeroPoints==0):
249 251
250 252 self.__updateSpecFromVoltage()
251 253 if pairsList == None:
252 254 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
253 255 else:
254 256 self.dataOut.pairsList = pairsList
255 257 self.__getFft()
256 258 self.dataOut.flagNoData = False
257 259 self.firstdatatime = None
258 260 self.nsamplesFFT = self.profIndex
259 261 #if not self.reader.bypass:
260 262 self.profIndex = 0
261 263 #update Processing Header:
262 264 self.dataOut.processingHeaderObj.dtype = "Spectra"
263 265 self.dataOut.processingHeaderObj.nFFTPoints = self.dataOut.nFFTPoints
264 266 self.dataOut.processingHeaderObj.nSamplesFFT = self.nsamplesFFT
265 267 self.dataOut.processingHeaderObj.nIncohInt = 1
266 268
267 269 elif self.dataIn.type == "Parameters": #when get data from h5 spc file
268 270
269 271 self.dataOut.data_spc = self.dataIn.data_spc
270 272 self.dataOut.data_cspc = self.dataIn.data_cspc
271 273 self.dataOut.data_outlier = self.dataIn.data_outlier
272 274 self.dataOut.nProfiles = self.dataIn.nProfiles
273 275 self.dataOut.nIncohInt = self.dataIn.nIncohInt
274 276 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
275 277 self.dataOut.ippFactor = self.dataIn.ippFactor
276 278 self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
277 279 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
278 280 self.dataOut.ProcessingHeader = self.dataIn.ProcessingHeader.copy()
279 281 self.dataOut.ippSeconds = self.dataIn.ippSeconds
280 282 self.dataOut.ipp = self.dataIn.ipp
281 283 #self.dataOut.abscissaList = self.dataIn.getVelRange(1)
282 284 #self.dataOut.spc_noise = self.dataIn.getNoise()
283 285 #self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
284 286 # self.dataOut.normFactor = self.dataIn.normFactor
285 287 if hasattr(self.dataIn, 'channelList'):
286 288 self.dataOut.channelList = self.dataIn.channelList
287 289 if hasattr(self.dataIn, 'pairsList'):
288 290 self.dataOut.pairsList = self.dataIn.pairsList
289 291 self.dataOut.groupList = self.dataIn.pairsList
290 292
291 293 self.dataOut.flagNoData = False
292 294
293 295 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
294 296 self.dataOut.ChanDist = self.dataIn.ChanDist
295 297 else: self.dataOut.ChanDist = None
296 298
297 299 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
298 300 # self.dataOut.VelRange = self.dataIn.VelRange
299 301 #else: self.dataOut.VelRange = None
300 302
301 303 else:
302 304 raise ValueError("The type of input object '%s' is not valid".format(
303 305 self.dataIn.type))
304
306 # print("SPC done")
305 307
306 308 def __selectPairs(self, pairsList):
307 309
308 310 if not pairsList:
309 311 return
310 312
311 313 pairs = []
312 314 pairsIndex = []
313 315
314 316 for pair in pairsList:
315 317 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
316 318 continue
317 319 pairs.append(pair)
318 320 pairsIndex.append(pairs.index(pair))
319 321
320 322 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
321 323 self.dataOut.pairsList = pairs
322 324
323 325 return
324 326
325 327 def selectFFTs(self, minFFT, maxFFT ):
326 328 """
327 329 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
328 330 minFFT<= FFT <= maxFFT
329 331 """
330 332
331 333 if (minFFT > maxFFT):
332 334 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
333 335
334 336 if (minFFT < self.dataOut.getFreqRange()[0]):
335 337 minFFT = self.dataOut.getFreqRange()[0]
336 338
337 339 if (maxFFT > self.dataOut.getFreqRange()[-1]):
338 340 maxFFT = self.dataOut.getFreqRange()[-1]
339 341
340 342 minIndex = 0
341 343 maxIndex = 0
342 344 FFTs = self.dataOut.getFreqRange()
343 345
344 346 inda = numpy.where(FFTs >= minFFT)
345 347 indb = numpy.where(FFTs <= maxFFT)
346 348
347 349 try:
348 350 minIndex = inda[0][0]
349 351 except:
350 352 minIndex = 0
351 353
352 354 try:
353 355 maxIndex = indb[0][-1]
354 356 except:
355 357 maxIndex = len(FFTs)
356 358
357 359 self.selectFFTsByIndex(minIndex, maxIndex)
358 360
359 361 return 1
360 362
361 363 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
362 364 newheis = numpy.where(
363 365 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
364 366
365 367 if hei_ref != None:
366 368 newheis = numpy.where(self.dataOut.heightList > hei_ref)
367 369
368 370 minIndex = min(newheis[0])
369 371 maxIndex = max(newheis[0])
370 372 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
371 373 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
372 374
373 375 # determina indices
374 376 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
375 377 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
376 378 avg_dB = 10 * \
377 379 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
378 380 beacon_dB = numpy.sort(avg_dB)[-nheis:]
379 381 beacon_heiIndexList = []
380 382 for val in avg_dB.tolist():
381 383 if val >= beacon_dB[0]:
382 384 beacon_heiIndexList.append(avg_dB.tolist().index(val))
383 385
384 386 data_cspc = None
385 387 if self.dataOut.data_cspc is not None:
386 388 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
387 389
388 390 data_dc = None
389 391 if self.dataOut.data_dc is not None:
390 392 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
391 393
392 394 self.dataOut.data_spc = data_spc
393 395 self.dataOut.data_cspc = data_cspc
394 396 self.dataOut.data_dc = data_dc
395 397 self.dataOut.heightList = heightList
396 398 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
397 399
398 400 return 1
399 401
400 402 def selectFFTsByIndex(self, minIndex, maxIndex):
401 403 """
402 404
403 405 """
404 406
405 407 if (minIndex < 0) or (minIndex > maxIndex):
406 408 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
407 409
408 410 if (maxIndex >= self.dataOut.nProfiles):
409 411 maxIndex = self.dataOut.nProfiles-1
410 412
411 413 #Spectra
412 414 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
413 415
414 416 data_cspc = None
415 417 if self.dataOut.data_cspc is not None:
416 418 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
417 419
418 420 data_dc = None
419 421 if self.dataOut.data_dc is not None:
420 422 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
421 423
422 424 self.dataOut.data_spc = data_spc
423 425 self.dataOut.data_cspc = data_cspc
424 426 self.dataOut.data_dc = data_dc
425 427
426 428 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
427 429 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
428 430 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
429 431
430 432 return 1
431 433
432 434 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
433 435 # validacion de rango
434 436 if minHei == None:
435 437 minHei = self.dataOut.heightList[0]
436 438
437 439 if maxHei == None:
438 440 maxHei = self.dataOut.heightList[-1]
439 441
440 442 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
441 443 print('minHei: %.2f is out of the heights range' % (minHei))
442 444 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
443 445 minHei = self.dataOut.heightList[0]
444 446
445 447 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
446 448 print('maxHei: %.2f is out of the heights range' % (maxHei))
447 449 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
448 450 maxHei = self.dataOut.heightList[-1]
449 451
450 452 # validacion de velocidades
451 453 velrange = self.dataOut.getVelRange(1)
452 454
453 455 if minVel == None:
454 456 minVel = velrange[0]
455 457
456 458 if maxVel == None:
457 459 maxVel = velrange[-1]
458 460
459 461 if (minVel < velrange[0]) or (minVel > maxVel):
460 462 print('minVel: %.2f is out of the velocity range' % (minVel))
461 463 print('minVel is setting to %.2f' % (velrange[0]))
462 464 minVel = velrange[0]
463 465
464 466 if (maxVel > velrange[-1]) or (maxVel < minVel):
465 467 print('maxVel: %.2f is out of the velocity range' % (maxVel))
466 468 print('maxVel is setting to %.2f' % (velrange[-1]))
467 469 maxVel = velrange[-1]
468 470
469 471 # seleccion de indices para rango
470 472 minIndex = 0
471 473 maxIndex = 0
472 474 heights = self.dataOut.heightList
473 475
474 476 inda = numpy.where(heights >= minHei)
475 477 indb = numpy.where(heights <= maxHei)
476 478
477 479 try:
478 480 minIndex = inda[0][0]
479 481 except:
480 482 minIndex = 0
481 483
482 484 try:
483 485 maxIndex = indb[0][-1]
484 486 except:
485 487 maxIndex = len(heights)
486 488
487 489 if (minIndex < 0) or (minIndex > maxIndex):
488 490 raise ValueError("some value in (%d,%d) is not valid" % (
489 491 minIndex, maxIndex))
490 492
491 493 if (maxIndex >= self.dataOut.nHeights):
492 494 maxIndex = self.dataOut.nHeights - 1
493 495
494 496 # seleccion de indices para velocidades
495 497 indminvel = numpy.where(velrange >= minVel)
496 498 indmaxvel = numpy.where(velrange <= maxVel)
497 499 try:
498 500 minIndexVel = indminvel[0][0]
499 501 except:
500 502 minIndexVel = 0
501 503
502 504 try:
503 505 maxIndexVel = indmaxvel[0][-1]
504 506 except:
505 507 maxIndexVel = len(velrange)
506 508
507 509 # seleccion del espectro
508 510 data_spc = self.dataOut.data_spc[:,
509 511 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
510 512 # estimacion de ruido
511 513 noise = numpy.zeros(self.dataOut.nChannels)
512 514
513 515 for channel in range(self.dataOut.nChannels):
514 516 daux = data_spc[channel, :, :]
515 517 sortdata = numpy.sort(daux, axis=None)
516 518 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
517 519
518 520 self.dataOut.noise_estimation = noise.copy()
519 521
520 522 return 1
521 523
522 524 class GetSNR(Operation):
523 525 '''
524 526 Written by R. Flores
525 527 '''
526 528 """Operation to get SNR.
527 529
528 530 Parameters:
529 531 -----------
530 532
531 533 Example
532 534 --------
533 535
534 536 op = proc_unit.addOperation(name='GetSNR', optype='other')
535 537
536 538 """
537 539
538 540 def __init__(self, **kwargs):
539 541
540 542 Operation.__init__(self, **kwargs)
541 543
542 544 def run(self,dataOut):
543 545
544 546 noise = dataOut.getNoise(ymin_index=-10) #RegiΓ³n superior donde solo deberΓ­a de haber ruido
545 547 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
546 548 dataOut.snl = numpy.log10(dataOut.data_snr)
547 549 dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
548 550
549 551 return dataOut
550 552
551 553 class removeDC(Operation):
552 554
553 555 def run(self, dataOut, mode=2):
554 556 self.dataOut = dataOut
555 557 jspectra = self.dataOut.data_spc
556 558 jcspectra = self.dataOut.data_cspc
557 559
558 560 num_chan = jspectra.shape[0]
559 561 num_hei = jspectra.shape[2]
560 562
561 563 if jcspectra is not None:
562 564 jcspectraExist = True
563 565 num_pairs = jcspectra.shape[0]
564 566 else:
565 567 jcspectraExist = False
566 568
567 569 freq_dc = int(jspectra.shape[1] / 2)
568 570 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
569 571 ind_vel = ind_vel.astype(int)
570 572
571 573 if ind_vel[0] < 0:
572 574 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
573 575
574 576 if mode == 1:
575 577 jspectra[:, freq_dc, :] = (
576 578 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
577 579
578 580 if jcspectraExist:
579 581 jcspectra[:, freq_dc, :] = (
580 582 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
581 583
582 584 if mode == 2:
583 585
584 586 vel = numpy.array([-2, -1, 1, 2])
585 587 xx = numpy.zeros([4, 4])
586 588
587 589 for fil in range(4):
588 590 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
589 591
590 592 xx_inv = numpy.linalg.inv(xx)
591 593 xx_aux = xx_inv[0, :]
592 594
593 595 for ich in range(num_chan):
594 596 yy = jspectra[ich, ind_vel, :]
595 597 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
596 598
597 599 junkid = jspectra[ich, freq_dc, :] <= 0
598 600 cjunkid = sum(junkid)
599 601
600 602 if cjunkid.any():
601 603 jspectra[ich, freq_dc, junkid.nonzero()] = (
602 604 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
603 605
604 606 if jcspectraExist:
605 607 for ip in range(num_pairs):
606 608 yy = jcspectra[ip, ind_vel, :]
607 609 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
608 610
609 611 self.dataOut.data_spc = jspectra
610 612 self.dataOut.data_cspc = jcspectra
611 613
612 614 return self.dataOut
613 615 class getNoiseB(Operation):
614 616 """
615 617 Get noise from custom heights and frequency ranges,
616 618 offset for additional manual correction
617 619 J. Apaza -> developed to amisr isr spectra
618 620
619 621 """
620 622 __slots__ =('offset','warnings', 'isConfig', 'minIndex','maxIndex','minIndexFFT','maxIndexFFT')
621 623 def __init__(self):
622 624
623 625 Operation.__init__(self)
624 626 self.isConfig = False
625 627
626 628 def setup(self, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
627 629
628 630 self.warnings = warnings
629 631 if minHei == None:
630 632 minHei = self.dataOut.heightList[0]
631 633
632 634 if maxHei == None:
633 635 maxHei = self.dataOut.heightList[-1]
634 636
635 637 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
636 638 if self.warnings:
637 639 print('minHei: %.2f is out of the heights range' % (minHei))
638 640 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
639 641 minHei = self.dataOut.heightList[0]
640 642
641 643 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
642 644 if self.warnings:
643 645 print('maxHei: %.2f is out of the heights range' % (maxHei))
644 646 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
645 647 maxHei = self.dataOut.heightList[-1]
646 648
647 649
648 650 #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
649 651 minIndexFFT = 0
650 652 maxIndexFFT = 0
651 653 # validacion de velocidades
652 654 indminPoint = None
653 655 indmaxPoint = None
654 656 if self.dataOut.type == 'Spectra':
655 657 if minVel == None and maxVel == None :
656 658
657 659 freqrange = self.dataOut.getFreqRange(1)
658 660
659 661 if minFreq == None:
660 662 minFreq = freqrange[0]
661 663
662 664 if maxFreq == None:
663 665 maxFreq = freqrange[-1]
664 666
665 667 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
666 668 if self.warnings:
667 669 print('minFreq: %.2f is out of the frequency range' % (minFreq))
668 670 print('minFreq is setting to %.2f' % (freqrange[0]))
669 671 minFreq = freqrange[0]
670 672
671 673 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
672 674 if self.warnings:
673 675 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
674 676 print('maxFreq is setting to %.2f' % (freqrange[-1]))
675 677 maxFreq = freqrange[-1]
676 678
677 679 indminPoint = numpy.where(freqrange >= minFreq)
678 680 indmaxPoint = numpy.where(freqrange <= maxFreq)
679 681
680 682 else:
681 683
682 684 velrange = self.dataOut.getVelRange(1)
683 685
684 686 if minVel == None:
685 687 minVel = velrange[0]
686 688
687 689 if maxVel == None:
688 690 maxVel = velrange[-1]
689 691
690 692 if (minVel < velrange[0]) or (minVel > maxVel):
691 693 if self.warnings:
692 694 print('minVel: %.2f is out of the velocity range' % (minVel))
693 695 print('minVel is setting to %.2f' % (velrange[0]))
694 696 minVel = velrange[0]
695 697
696 698 if (maxVel > velrange[-1]) or (maxVel < minVel):
697 699 if self.warnings:
698 700 print('maxVel: %.2f is out of the velocity range' % (maxVel))
699 701 print('maxVel is setting to %.2f' % (velrange[-1]))
700 702 maxVel = velrange[-1]
701 703
702 704 indminPoint = numpy.where(velrange >= minVel)
703 705 indmaxPoint = numpy.where(velrange <= maxVel)
704 706
705 707
706 708 # seleccion de indices para rango REEMPLAZAR FOR FUNCION EXTERNA LUEGO
707 709 # minIndex = 0
708 710 # maxIndex = 0
709 711 # heights = self.dataOut.heightList
710 712 # inda = numpy.where(heights >= minHei)
711 713 # indb = numpy.where(heights <= maxHei)
712 714 # try:
713 715 # minIndex = inda[0][0]
714 716 # except:
715 717 # minIndex = 0
716 718 # try:
717 719 # maxIndex = indb[0][-1]
718 720 # except:
719 721 # maxIndex = len(heights)
720 722 # if (minIndex < 0) or (minIndex > maxIndex):
721 723 # raise ValueError("some value in (%d,%d) is not valid" % (
722 724 # minIndex, maxIndex))
723 725 # if (maxIndex >= self.dataOut.nHeights):
724 726 # maxIndex = self.dataOut.nHeights - 1
725 727
726 728 minIndex, maxIndex = getHei_index(minHei,maxHei,self.dataOut.heightList)
727 729
728 730
729 731 #############################################################3
730 732 # seleccion de indices para velocidades
731 733 if self.dataOut.type == 'Spectra':
732 734 try:
733 735 minIndexFFT = indminPoint[0][0]
734 736 except:
735 737 minIndexFFT = 0
736 738
737 739 try:
738 740 maxIndexFFT = indmaxPoint[0][-1]
739 741 except:
740 742 maxIndexFFT = len( self.dataOut.getFreqRange(1))
741 743
742 744 self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT
743 745 self.isConfig = True
744 746 self.offset = 1
745 747 if offset!=None:
746 748 self.offset = 10**(offset/10)
747 749
748 750
749 751 def run(self, dataOut, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
750 752 self.dataOut = dataOut
751 753
752 754 if not self.isConfig:
753 755 self.setup(offset, minHei, maxHei,minVel, maxVel, minFreq, maxFreq, warnings)
754 756
755 757 self.dataOut.noise_estimation = None
756 758 noise = None
757 759 if self.dataOut.type == 'Voltage':
758 760 noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
759 761 elif self.dataOut.type == 'Spectra':
760 762 noise = numpy.zeros( self.dataOut.nChannels)
761 763 norm = 1
762 764
763 765 for channel in range( self.dataOut.nChannels):
764 766 if not hasattr(self.dataOut.nIncohInt,'__len__'):
765 767 norm = 1
766 768 else:
767 769 norm = self.dataOut.max_nIncohInt[channel]/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
768 770
769 771 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
770 772 daux = numpy.multiply(daux, norm)
771 773 sortdata = numpy.sort(daux, axis=None)
772 774 noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt[channel])/self.offset
773 775
774 776 else:
775 777 noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
776 778
777 779 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
778 780
779 781 return self.dataOut
780 782
781 783 def getNoiseByMean(self,data):
782 784 #data debe estar ordenado
783 785 data = numpy.mean(data,axis=1)
784 786 sortdata = numpy.sort(data, axis=None)
785 787 pnoise = None
786 788 j = 0
787 789
788 790 mean = numpy.mean(sortdata)
789 791 min = numpy.min(sortdata)
790 792 delta = mean - min
791 793 indexes = numpy.where(sortdata > (mean+delta))[0] #only array of indexes
792 794 #print(len(indexes))
793 795 if len(indexes)==0:
794 796 pnoise = numpy.mean(sortdata)
795 797 else:
796 798 j = indexes[0]
797 799 pnoise = numpy.mean(sortdata[0:j])
798 800
799 801 return pnoise
800 802
801 803 def getNoiseByHS(self,data, navg):
802 804 #data debe estar ordenado
803 805 #data = numpy.mean(data,axis=1)
804 806 sortdata = numpy.sort(data, axis=None)
805 807
806 808 lenOfData = len(sortdata)
807 809 nums_min = lenOfData*0.2
808 810
809 811 if nums_min <= 5:
810 812
811 813 nums_min = 5
812 814
813 815 sump = 0.
814 816 sumq = 0.
815 817
816 818 j = 0
817 819 cont = 1
818 820
819 821 while((cont == 1)and(j < lenOfData)):
820 822
821 823 sump += sortdata[j]
822 824 sumq += sortdata[j]**2
823 825 #sumq -= sump**2
824 826 if j > nums_min:
825 827 rtest = float(j)/(j-1) + 1.0/navg
826 828 #if ((sumq*j) > (sump**2)):
827 829 if ((sumq*j) > (rtest*sump**2)):
828 830 j = j - 1
829 831 sump = sump - sortdata[j]
830 832 sumq = sumq - sortdata[j]**2
831 833 cont = 0
832 834
833 835 j += 1
834 836
835 837 lnoise = sump / j
836 838
837 839 return lnoise
838 840
839 841 class removeInterference(Operation):
840 842
841 843 def removeInterference2(self):
842 844
843 845 cspc = self.dataOut.data_cspc
844 846 spc = self.dataOut.data_spc
845 847 Heights = numpy.arange(cspc.shape[2])
846 848 realCspc = numpy.abs(cspc)
847 849
848 850 for i in range(cspc.shape[0]):
849 851 LinePower= numpy.sum(realCspc[i], axis=0)
850 852 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
851 853 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
852 854 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
853 855 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
854 856 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
855 857
856 858
857 859 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
858 860 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
859 861 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
860 862 cspc[i,InterferenceRange,:] = numpy.NaN
861 863
862 864 self.dataOut.data_cspc = cspc
863 865
864 866 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
865 867
866 868 jspectra = self.dataOut.data_spc
867 869 jcspectra = self.dataOut.data_cspc
868 870 jnoise = self.dataOut.getNoise()
869 871 num_incoh = self.dataOut.nIncohInt
870 872
871 873 num_channel = jspectra.shape[0]
872 874 num_prof = jspectra.shape[1]
873 875 num_hei = jspectra.shape[2]
874 876
875 877 # hei_interf
876 878 if hei_interf is None:
877 879 count_hei = int(num_hei / 2)
878 880 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
879 881 hei_interf = numpy.asarray(hei_interf)[0]
880 882 # nhei_interf
881 883 if (nhei_interf == None):
882 884 nhei_interf = 5
883 885 if (nhei_interf < 1):
884 886 nhei_interf = 1
885 887 if (nhei_interf > count_hei):
886 888 nhei_interf = count_hei
887 889 if (offhei_interf == None):
888 890 offhei_interf = 0
889 891
890 892 ind_hei = list(range(num_hei))
891 893 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
892 894 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
893 895 mask_prof = numpy.asarray(list(range(num_prof)))
894 896 num_mask_prof = mask_prof.size
895 897 comp_mask_prof = [0, num_prof / 2]
896 898
897 899 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
898 900 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
899 901 jnoise = numpy.nan
900 902 noise_exist = jnoise[0] < numpy.Inf
901 903
902 904 # Subrutina de Remocion de la Interferencia
903 905 for ich in range(num_channel):
904 906 # Se ordena los espectros segun su potencia (menor a mayor)
905 907 power = jspectra[ich, mask_prof, :]
906 908 power = power[:, hei_interf]
907 909 power = power.sum(axis=0)
908 910 psort = power.ravel().argsort()
909 911
910 912 # Se estima la interferencia promedio en los Espectros de Potencia empleando
911 913 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
912 914 offhei_interf, nhei_interf + offhei_interf))]]]
913 915
914 916 if noise_exist:
915 917 # tmp_noise = jnoise[ich] / num_prof
916 918 tmp_noise = jnoise[ich]
917 919 junkspc_interf = junkspc_interf - tmp_noise
918 920 #junkspc_interf[:,comp_mask_prof] = 0
919 921
920 922 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
921 923 jspc_interf = jspc_interf.transpose()
922 924 # Calculando el espectro de interferencia promedio
923 925 noiseid = numpy.where(
924 926 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
925 927 noiseid = noiseid[0]
926 928 cnoiseid = noiseid.size
927 929 interfid = numpy.where(
928 930 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
929 931 interfid = interfid[0]
930 932 cinterfid = interfid.size
931 933
932 934 if (cnoiseid > 0):
933 935 jspc_interf[noiseid] = 0
934 936
935 937 # Expandiendo los perfiles a limpiar
936 938 if (cinterfid > 0):
937 939 new_interfid = (
938 940 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
939 941 new_interfid = numpy.asarray(new_interfid)
940 942 new_interfid = {x for x in new_interfid}
941 943 new_interfid = numpy.array(list(new_interfid))
942 944 new_cinterfid = new_interfid.size
943 945 else:
944 946 new_cinterfid = 0
945 947
946 948 for ip in range(new_cinterfid):
947 949 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
948 950 jspc_interf[new_interfid[ip]
949 951 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
950 952
951 953 jspectra[ich, :, ind_hei] = jspectra[ich, :,
952 954 ind_hei] - jspc_interf # Corregir indices
953 955
954 956 # Removiendo la interferencia del punto de mayor interferencia
955 957 ListAux = jspc_interf[mask_prof].tolist()
956 958 maxid = ListAux.index(max(ListAux))
957 959
958 960 if cinterfid > 0:
959 961 for ip in range(cinterfid * (interf == 2) - 1):
960 962 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
961 963 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
962 964 cind = len(ind)
963 965
964 966 if (cind > 0):
965 967 jspectra[ich, interfid[ip], ind] = tmp_noise * \
966 968 (1 + (numpy.random.uniform(cind) - 0.5) /
967 969 numpy.sqrt(num_incoh))
968 970
969 971 ind = numpy.array([-2, -1, 1, 2])
970 972 xx = numpy.zeros([4, 4])
971 973
972 974 for id1 in range(4):
973 975 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
974 976
975 977 xx_inv = numpy.linalg.inv(xx)
976 978 xx = xx_inv[:, 0]
977 979 ind = (ind + maxid + num_mask_prof) % num_mask_prof
978 980 yy = jspectra[ich, mask_prof[ind], :]
979 981 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
980 982 yy.transpose(), xx)
981 983
982 984 indAux = (jspectra[ich, :, :] < tmp_noise *
983 985 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
984 986 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
985 987 (1 - 1 / numpy.sqrt(num_incoh))
986 988
987 989 # Remocion de Interferencia en el Cross Spectra
988 990 if jcspectra is None:
989 991 return jspectra, jcspectra
990 992 num_pairs = int(jcspectra.size / (num_prof * num_hei))
991 993 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
992 994
993 995 for ip in range(num_pairs):
994 996
995 997 #-------------------------------------------
996 998
997 999 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
998 1000 cspower = cspower[:, hei_interf]
999 1001 cspower = cspower.sum(axis=0)
1000 1002
1001 1003 cspsort = cspower.ravel().argsort()
1002 1004 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1003 1005 offhei_interf, nhei_interf + offhei_interf))]]]
1004 1006 junkcspc_interf = junkcspc_interf.transpose()
1005 1007 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1006 1008
1007 1009 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1008 1010
1009 1011 median_real = int(numpy.median(numpy.real(
1010 1012 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1011 1013 median_imag = int(numpy.median(numpy.imag(
1012 1014 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1013 1015 comp_mask_prof = [int(e) for e in comp_mask_prof]
1014 1016 junkcspc_interf[comp_mask_prof, :] = numpy.complex_(
1015 1017 median_real, median_imag)
1016 1018
1017 1019 for iprof in range(num_prof):
1018 1020 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1019 1021 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1020 1022
1021 1023 # Removiendo la Interferencia
1022 1024 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1023 1025 :, ind_hei] - jcspc_interf
1024 1026
1025 1027 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1026 1028 maxid = ListAux.index(max(ListAux))
1027 1029
1028 1030 ind = numpy.array([-2, -1, 1, 2])
1029 1031 xx = numpy.zeros([4, 4])
1030 1032
1031 1033 for id1 in range(4):
1032 1034 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1033 1035
1034 1036 xx_inv = numpy.linalg.inv(xx)
1035 1037 xx = xx_inv[:, 0]
1036 1038
1037 1039 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1038 1040 yy = jcspectra[ip, mask_prof[ind], :]
1039 1041 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1040 1042
1041 1043 # Guardar Resultados
1042 1044 self.dataOut.data_spc = jspectra
1043 1045 self.dataOut.data_cspc = jcspectra
1044 1046
1045 1047 return 1
1046 1048
1047 1049
1048 1050 def run(self, dataOut, interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None, mode=1):
1049 1051
1050 1052 self.dataOut = dataOut
1051 1053
1052 1054 if mode == 1:
1053 1055 self.removeInterference(interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None)
1054 1056 elif mode == 2:
1055 1057 self.removeInterference2()
1056 1058
1057 1059 return self.dataOut
1058 1060
1059 1061
1060 1062 class deflip(Operation):
1061 1063
1062 1064 def run(self, dataOut):
1063 1065 # arreglo 1: (num_chan, num_profiles, num_heights)
1064 1066 self.dataOut = dataOut
1065 1067
1066 1068 # JULIA-oblicua, indice 2
1067 1069 # arreglo 2: (num_profiles, num_heights)
1068 1070 jspectra = self.dataOut.data_spc[2]
1069 1071 jspectra_tmp=numpy.zeros(jspectra.shape)
1070 1072 num_profiles=jspectra.shape[0]
1071 1073 freq_dc = int(num_profiles / 2)
1072 1074 # Flip con for
1073 1075 for j in range(num_profiles):
1074 1076 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1075 1077 # Intercambio perfil de DC con perfil inmediato anterior
1076 1078 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1077 1079 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1078 1080 # canal modificado es re-escrito en el arreglo de canales
1079 1081 self.dataOut.data_spc[2] = jspectra_tmp
1080 1082
1081 1083 return self.dataOut
1082 1084
1083 1085
1084 1086 class IncohInt(Operation):
1085 1087
1086 1088 __profIndex = 0
1087 1089 __withOverapping = False
1088 1090
1089 1091 __byTime = False
1090 1092 __initime = None
1091 1093 __lastdatatime = None
1092 1094 __integrationtime = None
1093 1095
1094 1096 __buffer_spc = None
1095 1097 __buffer_cspc = None
1096 1098 __buffer_dc = None
1097 1099
1098 1100 __dataReady = False
1099 1101
1100 1102 __timeInterval = None
1101 1103 incohInt = 0
1102 1104 nOutliers = 0
1103 1105 n = None
1104 1106
1105 1107 _flagProfilesByRange = False
1106 1108 _nProfilesByRange = 0
1107 1109 def __init__(self):
1108 1110
1109 1111 Operation.__init__(self)
1110 1112
1111 1113 def setup(self, n=None, timeInterval=None, overlapping=False):
1112 1114 """
1113 1115 Set the parameters of the integration class.
1114 1116
1115 1117 Inputs:
1116 1118
1117 1119 n : Number of coherent integrations
1118 1120 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1119 1121 overlapping :
1120 1122
1121 1123 """
1122 1124
1123 1125 self.__initime = None
1124 1126 self.__lastdatatime = 0
1125 1127
1126 1128 self.__buffer_spc = 0
1127 1129 self.__buffer_cspc = 0
1128 1130 self.__buffer_dc = 0
1129 1131
1130 1132 self.__profIndex = 0
1131 1133 self.__dataReady = False
1132 1134 self.__byTime = False
1133 1135 self.incohInt = 0
1134 1136 self.nOutliers = 0
1135 1137 if n is None and timeInterval is None:
1136 1138 raise ValueError("n or timeInterval should be specified ...")
1137 1139
1138 1140 if n is not None:
1139 1141 self.n = int(n)
1140 1142 else:
1141 1143
1142 1144 self.__integrationtime = int(timeInterval)
1143 1145 self.n = None
1144 1146 self.__byTime = True
1145 1147
1146 1148
1147 1149
1148 1150 def putData(self, data_spc, data_cspc, data_dc):
1149 1151 """
1150 1152 Add a profile to the __buffer_spc and increase in one the __profileIndex
1151 1153
1152 1154 """
1153 1155 if data_spc.all() == numpy.nan :
1154 1156 print("nan ")
1155 1157 return
1156 1158 self.__buffer_spc += data_spc
1157 1159
1158 1160 if data_cspc is None:
1159 1161 self.__buffer_cspc = None
1160 1162 else:
1161 1163 self.__buffer_cspc += data_cspc
1162 1164
1163 1165 if data_dc is None:
1164 1166 self.__buffer_dc = None
1165 1167 else:
1166 1168 self.__buffer_dc += data_dc
1167 1169
1168 1170 self.__profIndex += 1
1169 1171
1170 1172 return
1171 1173
1172 1174 def pushData(self):
1173 1175 """
1174 1176 Return the sum of the last profiles and the profiles used in the sum.
1175 1177
1176 1178 Affected:
1177 1179
1178 1180 self.__profileIndex
1179 1181
1180 1182 """
1181 1183
1182 1184 data_spc = self.__buffer_spc
1183 1185 data_cspc = self.__buffer_cspc
1184 1186 data_dc = self.__buffer_dc
1185 1187 n = self.__profIndex
1186 1188
1187 1189 self.__buffer_spc = 0
1188 1190 self.__buffer_cspc = 0
1189 1191 self.__buffer_dc = 0
1190 1192
1191 1193
1192 1194 return data_spc, data_cspc, data_dc, n
1193 1195
1194 1196 def byProfiles(self, *args):
1195 1197
1196 1198 self.__dataReady = False
1197 1199 avgdata_spc = None
1198 1200 avgdata_cspc = None
1199 1201 avgdata_dc = None
1200 1202
1201 1203 self.putData(*args)
1202 1204
1203 1205 if self.__profIndex == self.n:
1204 1206
1205 1207 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1206 1208 self.n = n
1207 1209 self.__dataReady = True
1208 1210
1209 1211 return avgdata_spc, avgdata_cspc, avgdata_dc
1210 1212
1211 1213 def byTime(self, datatime, *args):
1212 1214
1213 1215 self.__dataReady = False
1214 1216 avgdata_spc = None
1215 1217 avgdata_cspc = None
1216 1218 avgdata_dc = None
1217 1219
1218 1220 self.putData(*args)
1219 1221
1220 1222 if (datatime - self.__initime) >= self.__integrationtime:
1221 1223 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1222 1224 self.n = n
1223 1225 self.__dataReady = True
1224 1226
1225 1227 return avgdata_spc, avgdata_cspc, avgdata_dc
1226 1228
1227 1229 def integrate(self, datatime, *args):
1228 1230
1229 1231 if self.__profIndex == 0:
1230 1232 self.__initime = datatime
1231 1233
1232 1234 if self.__byTime:
1233 1235 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1234 1236 datatime, *args)
1235 1237 else:
1236 1238 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1237 1239
1238 1240 if not self.__dataReady:
1239 1241 return None, None, None, None
1240 1242
1241 1243 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1242 1244
1243 1245 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1244 1246 if n == 1:
1245 1247 return dataOut
1246 1248
1247 1249 if dataOut.flagNoData == True:
1248 1250 return dataOut
1249 1251
1250 1252 if dataOut.flagProfilesByRange == True:
1251 1253 self._flagProfilesByRange = True
1252 1254
1253 1255 dataOut.flagNoData = True
1254 1256 dataOut.processingHeaderObj.timeIncohInt = timeInterval
1255 1257 if not self.isConfig:
1256 1258 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
1257 1259 self.setup(n, timeInterval, overlapping)
1258 1260 self.isConfig = True
1259 1261
1260 1262
1261 1263 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1262 1264 dataOut.data_spc,
1263 1265 dataOut.data_cspc,
1264 1266 dataOut.data_dc)
1265 1267
1266 1268 self.incohInt += dataOut.nIncohInt
1267 1269
1268 1270
1269 1271 if isinstance(dataOut.data_outlier,numpy.ndarray) or isinstance(dataOut.data_outlier,int) or isinstance(dataOut.data_outlier, float):
1270 1272 self.nOutliers += dataOut.data_outlier
1271 1273
1272 1274 if self._flagProfilesByRange:
1273 1275 dataOut.flagProfilesByRange = True
1274 1276 self._nProfilesByRange += dataOut.nProfilesByRange
1275 1277
1276 1278 if self.__dataReady:
1277 1279 #print("prof: ",dataOut.max_nIncohInt,self.__profIndex)
1278 1280 dataOut.data_spc = avgdata_spc
1279 1281 dataOut.data_cspc = avgdata_cspc
1280 1282 dataOut.data_dc = avgdata_dc
1281 1283 dataOut.nIncohInt = self.incohInt
1282 1284 dataOut.data_outlier = self.nOutliers
1283 1285 dataOut.utctime = avgdatatime
1284 1286 dataOut.flagNoData = False
1285 1287 self.incohInt = 0
1286 1288 self.nOutliers = 0
1287 1289 self.__profIndex = 0
1288 1290 dataOut.nProfilesByRange = self._nProfilesByRange
1289 1291 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
1290 1292 self._flagProfilesByRange = False
1291 1293 # print("IncohInt Done")
1292 1294 return dataOut
1293 1295
1294 1296
1295 1297 class IntegrationFaradaySpectra(Operation):
1296 1298
1297 1299 __profIndex = 0
1298 1300 __withOverapping = False
1299 1301
1300 1302 __byTime = False
1301 1303 __initime = None
1302 1304 __lastdatatime = None
1303 1305 __integrationtime = None
1304 1306
1305 1307 __buffer_spc = None
1306 1308 __buffer_cspc = None
1307 1309 __buffer_dc = None
1308 1310
1309 1311 __dataReady = False
1310 1312
1311 1313 __timeInterval = None
1312 1314 n_ints = None #matriz de numero de integracions (CH,HEI)
1313 1315 n = None
1314 1316 minHei_ind = None
1315 1317 maxHei_ind = None
1316 1318 navg = 1.0
1317 1319 factor = 0.0
1318 1320 dataoutliers = None # (CHANNELS, HEIGHTS)
1319 1321
1320 1322 _flagProfilesByRange = False
1321 1323 _nProfilesByRange = 0
1322 1324
1323 1325 def __init__(self):
1324 1326
1325 1327 Operation.__init__(self)
1326 1328
1327 1329 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, avg=1,factor=0.75):
1328 1330 """
1329 1331 Set the parameters of the integration class.
1330 1332
1331 1333 Inputs:
1332 1334
1333 1335 n : Number of coherent integrations
1334 1336 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1335 1337 overlapping :
1336 1338
1337 1339 """
1338 1340
1339 1341 self.__initime = None
1340 1342 self.__lastdatatime = 0
1341 1343
1342 1344 self.__buffer_spc = []
1343 1345 self.__buffer_cspc = []
1344 1346 self.__buffer_dc = 0
1345 1347
1346 1348 self.__profIndex = 0
1347 1349 self.__dataReady = False
1348 1350 self.__byTime = False
1349 1351
1350 1352 self.factor = factor
1351 1353 self.navg = avg
1352 1354 #self.ByLags = dataOut.ByLags ###REDEFINIR
1353 1355 self.ByLags = False
1354 1356 self.maxProfilesInt = 0
1355 1357 self.__nChannels = dataOut.nChannels
1356 1358 if DPL != None:
1357 1359 self.DPL=DPL
1358 1360 else:
1359 1361 #self.DPL=dataOut.DPL ###REDEFINIR
1360 1362 self.DPL=0
1361 1363
1362 1364 if n is None and timeInterval is None:
1363 1365 raise ValueError("n or timeInterval should be specified ...")
1364 1366
1365 1367 if n is not None:
1366 1368 self.n = int(n)
1367 1369 else:
1368 1370 self.__integrationtime = int(timeInterval)
1369 1371 self.n = None
1370 1372 self.__byTime = True
1371 1373
1372 1374
1373 1375 if minHei == None:
1374 1376 minHei = self.dataOut.heightList[0]
1375 1377
1376 1378 if maxHei == None:
1377 1379 maxHei = self.dataOut.heightList[-1]
1378 1380
1379 1381 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1380 1382 print('minHei: %.2f is out of the heights range' % (minHei))
1381 1383 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
1382 1384 minHei = self.dataOut.heightList[0]
1383 1385
1384 1386 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
1385 1387 print('maxHei: %.2f is out of the heights range' % (maxHei))
1386 1388 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
1387 1389 maxHei = self.dataOut.heightList[-1]
1388 1390
1389 1391 ind_list1 = numpy.where(self.dataOut.heightList >= minHei)
1390 1392 ind_list2 = numpy.where(self.dataOut.heightList <= maxHei)
1391 1393 self.minHei_ind = ind_list1[0][0]
1392 1394 self.maxHei_ind = ind_list2[0][-1]
1393 1395
1394 1396 def putData(self, data_spc, data_cspc, data_dc):
1395 1397 """
1396 1398 Add a profile to the __buffer_spc and increase in one the __profileIndex
1397 1399
1398 1400 """
1399 1401
1400 1402 self.__buffer_spc.append(data_spc)
1401 1403
1402 1404 if self.__nChannels < 2:
1403 1405 self.__buffer_cspc = None
1404 1406 else:
1405 1407 self.__buffer_cspc.append(data_cspc)
1406 1408
1407 1409 if data_dc is None:
1408 1410 self.__buffer_dc = None
1409 1411 else:
1410 1412 self.__buffer_dc += data_dc
1411 1413
1412 1414 self.__profIndex += 1
1413 1415
1414 1416 return
1415 1417
1416 1418 def hildebrand_sekhon_Integration(self,sortdata,navg, factor):
1417 1419 #data debe estar ordenado
1418 1420 #sortdata = numpy.sort(data, axis=None)
1419 1421 #sortID=data.argsort()
1420 1422 lenOfData = len(sortdata)
1421 1423 nums_min = lenOfData*factor
1422 1424 if nums_min <= 5:
1423 1425 nums_min = 5
1424 1426 sump = 0.
1425 1427 sumq = 0.
1426 1428 j = 0
1427 1429 cont = 1
1428 1430 while((cont == 1)and(j < lenOfData)):
1429 1431 sump += sortdata[j]
1430 1432 sumq += sortdata[j]**2
1431 1433 if j > nums_min:
1432 1434 rtest = float(j)/(j-1) + 1.0/navg
1433 1435 if ((sumq*j) > (rtest*sump**2)):
1434 1436 j = j - 1
1435 1437 sump = sump - sortdata[j]
1436 1438 sumq = sumq - sortdata[j]**2
1437 1439 cont = 0
1438 1440 j += 1
1439 1441 #lnoise = sump / j
1440 1442 #print("H S done")
1441 1443 #return j,sortID
1442 1444 return j
1443 1445
1444 1446
1445 1447 def pushData(self):
1446 1448 """
1447 1449 Return the sum of the last profiles and the profiles used in the sum.
1448 1450
1449 1451 Affected:
1450 1452
1451 1453 self.__profileIndex
1452 1454
1453 1455 """
1454 1456 bufferH=None
1455 1457 buffer=None
1456 1458 buffer1=None
1457 1459 buffer_cspc=None
1458 1460 #print("aes: ", self.__buffer_cspc)
1459 1461 self.__buffer_spc=numpy.array(self.__buffer_spc)
1460 1462 if self.__nChannels > 1 :
1461 1463 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1462 1464
1463 1465 #print("FREQ_DC",self.__buffer_spc.shape,self.__buffer_cspc.shape)
1464 1466
1465 1467 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1466 1468 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1467 1469
1468 1470 self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
1469 1471
1470 1472 for k in range(self.minHei_ind,self.maxHei_ind):
1471 1473 if self.__nChannels > 1:
1472 1474 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1473 1475
1474 1476 outliers_IDs_cspc=[]
1475 1477 cspc_outliers_exist=False
1476 1478 for i in range(self.nChannels):#dataOut.nChannels):
1477 1479
1478 1480 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1479 1481 indexes=[]
1480 1482 #sortIDs=[]
1481 1483 outliers_IDs=[]
1482 1484
1483 1485 for j in range(self.nProfiles): #frecuencias en el tiempo
1484 1486 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1485 1487 # continue
1486 1488 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1487 1489 # continue
1488 1490 buffer=buffer1[:,j]
1489 1491 sortdata = numpy.sort(buffer, axis=None)
1490 1492
1491 1493 sortID=buffer.argsort()
1492 1494 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1493 1495
1494 1496 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1495 1497
1496 1498 # fig,ax = plt.subplots()
1497 1499 # ax.set_title(str(k)+" "+str(j))
1498 1500 # x=range(len(sortdata))
1499 1501 # ax.scatter(x,sortdata)
1500 1502 # ax.axvline(index)
1501 1503 # plt.show()
1502 1504
1503 1505 indexes.append(index)
1504 1506 #sortIDs.append(sortID)
1505 1507 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1506 1508
1507 1509 #print("Outliers: ",outliers_IDs)
1508 1510 outliers_IDs=numpy.array(outliers_IDs)
1509 1511 outliers_IDs=outliers_IDs.ravel()
1510 1512 outliers_IDs=numpy.unique(outliers_IDs)
1511 1513 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1512 1514 indexes=numpy.array(indexes)
1513 1515 indexmin=numpy.min(indexes)
1514 1516
1515 1517
1516 1518 #print(indexmin,buffer1.shape[0], k)
1517 1519
1518 1520 # fig,ax = plt.subplots()
1519 1521 # ax.plot(sortdata)
1520 1522 # ax2 = ax.twinx()
1521 1523 # x=range(len(indexes))
1522 1524 # #plt.scatter(x,indexes)
1523 1525 # ax2.scatter(x,indexes)
1524 1526 # plt.show()
1525 1527
1526 1528 if indexmin != buffer1.shape[0]:
1527 1529 if self.__nChannels > 1:
1528 1530 cspc_outliers_exist= True
1529 1531
1530 1532 lt=outliers_IDs
1531 1533 #avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1532 1534
1533 1535 for p in list(outliers_IDs):
1534 1536 #buffer1[p,:]=avg
1535 1537 buffer1[p,:] = numpy.NaN
1536 1538
1537 1539 self.dataOutliers[i,k] = len(outliers_IDs)
1538 1540
1539 1541
1540 1542 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1541 1543
1542 1544
1543 1545 if self.__nChannels > 1:
1544 1546 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1545 1547
1546 1548
1547 1549 if self.__nChannels > 1:
1548 1550 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1549 1551 if cspc_outliers_exist:
1550 1552
1551 1553 lt=outliers_IDs_cspc
1552 1554
1553 1555 #avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1554 1556 for p in list(outliers_IDs_cspc):
1555 1557 #buffer_cspc[p,:]=avg
1556 1558 buffer_cspc[p,:] = numpy.NaN
1557 1559
1558 1560 if self.__nChannels > 1:
1559 1561 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1560 1562
1561 1563
1562 1564
1563 1565
1564 1566 nOutliers = len(outliers_IDs)
1565 1567 #print("Outliers n: ",self.dataOutliers,nOutliers)
1566 1568 buffer=None
1567 1569 bufferH=None
1568 1570 buffer1=None
1569 1571 buffer_cspc=None
1570 1572
1571 1573
1572 1574 buffer=None
1573 1575
1574 1576 #data_spc = numpy.sum(self.__buffer_spc,axis=0)
1575 1577 data_spc = numpy.nansum(self.__buffer_spc,axis=0)
1576 1578 if self.__nChannels > 1:
1577 1579 #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1578 1580 data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
1579 1581 else:
1580 1582 data_cspc = None
1581 1583 data_dc = self.__buffer_dc
1582 1584 #(CH, HEIGH)
1583 1585 self.maxProfilesInt = self.__profIndex - 1
1584 1586 n = self.__profIndex - self.dataOutliers # n becomes a matrix
1585 1587
1586 1588 self.__buffer_spc = []
1587 1589 self.__buffer_cspc = []
1588 1590 self.__buffer_dc = 0
1589 1591 self.__profIndex = 0
1590 1592 #print("cleaned ",data_cspc)
1591 1593 return data_spc, data_cspc, data_dc, n
1592 1594
1593 1595 def byProfiles(self, *args):
1594 1596
1595 1597 self.__dataReady = False
1596 1598 avgdata_spc = None
1597 1599 avgdata_cspc = None
1598 1600 avgdata_dc = None
1599 1601
1600 1602 self.putData(*args)
1601 1603
1602 1604 if self.__profIndex >= self.n:
1603 1605
1604 1606 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1605 1607 self.n_ints = n
1606 1608 self.__dataReady = True
1607 1609
1608 1610 return avgdata_spc, avgdata_cspc, avgdata_dc
1609 1611
1610 1612 def byTime(self, datatime, *args):
1611 1613
1612 1614 self.__dataReady = False
1613 1615 avgdata_spc = None
1614 1616 avgdata_cspc = None
1615 1617 avgdata_dc = None
1616 1618
1617 1619 self.putData(*args)
1618 1620
1619 1621 if (datatime - self.__initime) >= self.__integrationtime:
1620 1622 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1621 1623 self.n_ints = n
1622 1624 self.__dataReady = True
1623 1625
1624 1626 return avgdata_spc, avgdata_cspc, avgdata_dc
1625 1627
1626 1628 def integrate(self, datatime, *args):
1627 1629
1628 1630 if self.__profIndex == 0:
1629 1631 self.__initime = datatime
1630 1632
1631 1633 if self.__byTime:
1632 1634 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1633 1635 datatime, *args)
1634 1636 else:
1635 1637 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1636 1638
1637 1639 if not self.__dataReady:
1638 1640 return None, None, None, None
1639 1641
1640 1642 #print("integrate", avgdata_cspc)
1641 1643 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1642 1644
1643 1645 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
1644 1646 self.dataOut = dataOut
1645 1647 if n == 1:
1646 1648 return self.dataOut
1647 1649 self.dataOut.processingHeaderObj.timeIncohInt = timeInterval
1648 1650
1649 1651 if dataOut.flagProfilesByRange:
1650 1652 self._flagProfilesByRange = True
1651 1653
1652 1654 if self.dataOut.nChannels == 1:
1653 1655 self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
1654 1656 #print("IN spc:", self.dataOut.data_spc.shape, self.dataOut.data_cspc)
1655 1657 if not self.isConfig:
1656 1658 self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor)
1657 1659 self.isConfig = True
1658 1660
1659 1661 if not self.ByLags:
1660 1662 self.nProfiles=self.dataOut.nProfiles
1661 1663 self.nChannels=self.dataOut.nChannels
1662 1664 self.nHeights=self.dataOut.nHeights
1663 1665 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1664 1666 self.dataOut.data_spc,
1665 1667 self.dataOut.data_cspc,
1666 1668 self.dataOut.data_dc)
1667 1669 else:
1668 1670 self.nProfiles=self.dataOut.nProfiles
1669 1671 self.nChannels=self.dataOut.nChannels
1670 1672 self.nHeights=self.dataOut.nHeights
1671 1673 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1672 1674 self.dataOut.dataLag_spc,
1673 1675 self.dataOut.dataLag_cspc,
1674 1676 self.dataOut.dataLag_dc)
1675 1677 self.dataOut.flagNoData = True
1676 1678
1677 1679 if self._flagProfilesByRange:
1678 1680 dataOut.flagProfilesByRange = True
1679 1681 self._nProfilesByRange += dataOut.nProfilesByRange
1680 1682
1681 1683 if self.__dataReady:
1682 1684
1683 1685 if not self.ByLags:
1684 1686 if self.nChannels == 1:
1685 1687 #print("f int", avgdata_spc.shape)
1686 1688 self.dataOut.data_spc = avgdata_spc
1687 1689 self.dataOut.data_cspc = None
1688 1690 else:
1689 1691 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
1690 1692 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1691 1693 self.dataOut.data_dc = avgdata_dc
1692 1694 self.dataOut.data_outlier = self.dataOutliers
1693 1695
1694 1696
1695 1697 else:
1696 1698 self.dataOut.dataLag_spc = avgdata_spc
1697 1699 self.dataOut.dataLag_cspc = avgdata_cspc
1698 1700 self.dataOut.dataLag_dc = avgdata_dc
1699 1701
1700 1702 self.dataOut.data_spc=self.dataOut.dataLag_spc[:,:,:,self.dataOut.LagPlot]
1701 1703 self.dataOut.data_cspc=self.dataOut.dataLag_cspc[:,:,:,self.dataOut.LagPlot]
1702 1704 self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot]
1703 1705
1704 1706 self.dataOut.nIncohInt *= self.n_ints
1705 1707
1706 1708 self.dataOut.utctime = avgdatatime
1707 1709 self.dataOut.flagNoData = False
1708 1710
1709 1711 dataOut.nProfilesByRange = self._nProfilesByRange
1710 1712 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
1711 1713 self._flagProfilesByRange = False
1712 1714
1713 1715 return self.dataOut
1714 1716
1715 1717 class dopplerFlip(Operation):
1716 1718
1717 1719 def run(self, dataOut, chann = None):
1718 1720 # arreglo 1: (num_chan, num_profiles, num_heights)
1719 1721 self.dataOut = dataOut
1720 1722 # JULIA-oblicua, indice 2
1721 1723 # arreglo 2: (num_profiles, num_heights)
1722 1724 jspectra = self.dataOut.data_spc[chann]
1723 1725 jspectra_tmp = numpy.zeros(jspectra.shape)
1724 1726 num_profiles = jspectra.shape[0]
1725 1727 freq_dc = int(num_profiles / 2)
1726 1728 # Flip con for
1727 1729 for j in range(num_profiles):
1728 1730 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1729 1731 # Intercambio perfil de DC con perfil inmediato anterior
1730 1732 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1731 1733 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1732 1734 # canal modificado es re-escrito en el arreglo de canales
1733 1735 self.dataOut.data_spc[chann] = jspectra_tmp
1734 1736
1735 1737 return self.dataOut No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now