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