##// END OF EJS Templates
schain
RSS
Fix bugs in SpectraCutPlot, GenericRTIPlot and saving when using throttle
jespinoza -
r1359:7cde62c34a28
parent child
Show More
Side by Side Unified
Context file:
r1359:7cde62c34a28 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,688 +1,691
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 """Base class to create plot operations
6 6
7 7 """
8 8
9 9 import os
10 10 import sys
11 11 import zmq
12 12 import time
13 13 import numpy
14 14 import datetime
15 15 from collections import deque
16 16 from functools import wraps
17 17 from threading import Thread
18 18 import matplotlib
19 19
20 20 if 'BACKEND' in os.environ:
21 21 matplotlib.use(os.environ['BACKEND'])
22 22 elif 'linux' in sys.platform:
23 23 matplotlib.use("TkAgg")
24 24 elif 'darwin' in sys.platform:
25 25 matplotlib.use('MacOSX')
26 26 else:
27 27 from schainpy.utils import log
28 28 log.warning('Using default Backend="Agg"', 'INFO')
29 29 matplotlib.use('Agg')
30 30
31 31 import matplotlib.pyplot as plt
32 32 from matplotlib.patches import Polygon
33 33 from mpl_toolkits.axes_grid1 import make_axes_locatable
34 34 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
35 35
36 36 from schainpy.model.data.jrodata import PlotterData
37 37 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
38 38 from schainpy.utils import log
39 39
40 40 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
41 41 blu_values = matplotlib.pyplot.get_cmap(
42 42 'seismic_r', 20)(numpy.arange(20))[10:15]
43 43 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
44 44 'jro', numpy.vstack((blu_values, jet_values)))
45 45 matplotlib.pyplot.register_cmap(cmap=ncmap)
46 46
47 47 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
48 48 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
49 49
50 50 EARTH_RADIUS = 6.3710e3
51 51
52 52 def ll2xy(lat1, lon1, lat2, lon2):
53 53
54 54 p = 0.017453292519943295
55 55 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
56 56 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
57 57 r = 12742 * numpy.arcsin(numpy.sqrt(a))
58 58 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
59 59 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
60 60 theta = -theta + numpy.pi/2
61 61 return r*numpy.cos(theta), r*numpy.sin(theta)
62 62
63 63
64 64 def km2deg(km):
65 65 '''
66 66 Convert distance in km to degrees
67 67 '''
68 68
69 69 return numpy.rad2deg(km/EARTH_RADIUS)
70 70
71 71
72 72 def figpause(interval):
73 73 backend = plt.rcParams['backend']
74 74 if backend in matplotlib.rcsetup.interactive_bk:
75 75 figManager = matplotlib._pylab_helpers.Gcf.get_active()
76 76 if figManager is not None:
77 77 canvas = figManager.canvas
78 78 if canvas.figure.stale:
79 79 canvas.draw()
80 80 try:
81 81 canvas.start_event_loop(interval)
82 82 except:
83 83 pass
84 84 return
85 85
86 86 def popup(message):
87 87 '''
88 88 '''
89 89
90 90 fig = plt.figure(figsize=(12, 8), facecolor='r')
91 91 text = '\n'.join([s.strip() for s in message.split(':')])
92 92 fig.text(0.01, 0.5, text, ha='left', va='center',
93 93 size='20', weight='heavy', color='w')
94 94 fig.show()
95 95 figpause(1000)
96 96
97 97
98 98 class Throttle(object):
99 99 '''
100 100 Decorator that prevents a function from being called more than once every
101 101 time period.
102 102 To create a function that cannot be called more than once a minute, but
103 103 will sleep until it can be called:
104 104 @Throttle(minutes=1)
105 105 def foo():
106 106 pass
107 107
108 108 for i in range(10):
109 109 foo()
110 110 print "This function has run %s times." % i
111 111 '''
112 112
113 113 def __init__(self, seconds=0, minutes=0, hours=0):
114 114 self.throttle_period = datetime.timedelta(
115 115 seconds=seconds, minutes=minutes, hours=hours
116 116 )
117 117
118 118 self.time_of_last_call = datetime.datetime.min
119 119
120 120 def __call__(self, fn):
121 121 @wraps(fn)
122 122 def wrapper(*args, **kwargs):
123 123 coerce = kwargs.pop('coerce', None)
124 124 if coerce:
125 125 self.time_of_last_call = datetime.datetime.now()
126 126 return fn(*args, **kwargs)
127 127 else:
128 128 now = datetime.datetime.now()
129 129 time_since_last_call = now - self.time_of_last_call
130 130 time_left = self.throttle_period - time_since_last_call
131 131
132 132 if time_left > datetime.timedelta(seconds=0):
133 133 return
134 134
135 135 self.time_of_last_call = datetime.datetime.now()
136 136 return fn(*args, **kwargs)
137 137
138 138 return wrapper
139 139
140 140 def apply_throttle(value):
141 141
142 142 @Throttle(seconds=value)
143 143 def fnThrottled(fn):
144 144 fn()
145 145
146 146 return fnThrottled
147 147
148 148
149 149 @MPDecorator
150 150 class Plot(Operation):
151 151 """Base class for Schain plotting operations
152 152
153 153 This class should never be use directtly you must subclass a new operation,
154 154 children classes must be defined as follow:
155 155
156 156 ExamplePlot(Plot):
157 157
158 158 CODE = 'code'
159 159 colormap = 'jet'
160 160 plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
161 161
162 162 def setup(self):
163 163 pass
164 164
165 165 def plot(self):
166 166 pass
167 167
168 168 """
169 169
170 170 CODE = 'Figure'
171 171 colormap = 'jet'
172 172 bgcolor = 'white'
173 173 buffering = True
174 174 __missing = 1E30
175 175
176 176 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
177 177 'showprofile']
178 178
179 179 def __init__(self):
180 180
181 181 Operation.__init__(self)
182 182 self.isConfig = False
183 183 self.isPlotConfig = False
184 184 self.save_time = 0
185 185 self.sender_time = 0
186 186 self.data = None
187 187 self.firsttime = True
188 188 self.sender_queue = deque(maxlen=10)
189 189 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
190 190
191 191 def __fmtTime(self, x, pos):
192 192 '''
193 193 '''
194 194
195 195 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
196 196
197 197 def __setup(self, **kwargs):
198 198 '''
199 199 Initialize variables
200 200 '''
201 201
202 202 self.figures = []
203 203 self.axes = []
204 204 self.cb_axes = []
205 205 self.localtime = kwargs.pop('localtime', True)
206 206 self.show = kwargs.get('show', True)
207 207 self.save = kwargs.get('save', False)
208 208 self.save_period = kwargs.get('save_period', 0)
209 209 self.colormap = kwargs.get('colormap', self.colormap)
210 210 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
211 211 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
212 212 self.colormaps = kwargs.get('colormaps', None)
213 213 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
214 214 self.showprofile = kwargs.get('showprofile', False)
215 215 self.title = kwargs.get('wintitle', self.CODE.upper())
216 216 self.cb_label = kwargs.get('cb_label', None)
217 217 self.cb_labels = kwargs.get('cb_labels', None)
218 218 self.labels = kwargs.get('labels', None)
219 219 self.xaxis = kwargs.get('xaxis', 'frequency')
220 220 self.zmin = kwargs.get('zmin', None)
221 221 self.zmax = kwargs.get('zmax', None)
222 222 self.zlimits = kwargs.get('zlimits', None)
223 223 self.xmin = kwargs.get('xmin', None)
224 224 self.xmax = kwargs.get('xmax', None)
225 225 self.xrange = kwargs.get('xrange', 12)
226 226 self.xscale = kwargs.get('xscale', None)
227 227 self.ymin = kwargs.get('ymin', None)
228 228 self.ymax = kwargs.get('ymax', None)
229 229 self.yscale = kwargs.get('yscale', None)
230 230 self.xlabel = kwargs.get('xlabel', None)
231 231 self.attr_time = kwargs.get('attr_time', 'utctime')
232 232 self.attr_data = kwargs.get('attr_data', 'data_param')
233 233 self.decimation = kwargs.get('decimation', None)
234 self.showSNR = kwargs.get('showSNR', False)
235 234 self.oneFigure = kwargs.get('oneFigure', True)
236 235 self.width = kwargs.get('width', None)
237 236 self.height = kwargs.get('height', None)
238 237 self.colorbar = kwargs.get('colorbar', True)
239 238 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
240 239 self.channels = kwargs.get('channels', None)
241 240 self.titles = kwargs.get('titles', [])
242 241 self.polar = False
243 242 self.type = kwargs.get('type', 'iq')
244 243 self.grid = kwargs.get('grid', False)
245 244 self.pause = kwargs.get('pause', False)
246 245 self.save_code = kwargs.get('save_code', self.CODE)
247 246 self.throttle = kwargs.get('throttle', 0)
248 247 self.exp_code = kwargs.get('exp_code', None)
249 248 self.server = kwargs.get('server', False)
250 249 self.sender_period = kwargs.get('sender_period', 60)
251 250 self.tag = kwargs.get('tag', '')
252 251 self.height_index = kwargs.get('height_index', None)
253 252 self.__throttle_plot = apply_throttle(self.throttle)
254 253 code = self.attr_data if self.attr_data else self.CODE
255 254 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
256 255
257 256 if self.server:
258 257 if not self.server.startswith('tcp://'):
259 258 self.server = 'tcp://{}'.format(self.server)
260 259 log.success(
261 260 'Sending to server: {}'.format(self.server),
262 261 self.name
263 262 )
264 263
264 if isinstance(self.attr_data, str):
265 self.attr_data = [self.attr_data]
266
265 267 def __setup_plot(self):
266 268 '''
267 269 Common setup for all figures, here figures and axes are created
268 270 '''
269 271
270 272 self.setup()
271 273
272 274 self.time_label = 'LT' if self.localtime else 'UTC'
273 275
274 276 if self.width is None:
275 277 self.width = 8
276 278
277 279 self.figures = []
278 280 self.axes = []
279 281 self.cb_axes = []
280 282 self.pf_axes = []
281 283 self.cmaps = []
282 284
283 285 size = '15%' if self.ncols == 1 else '30%'
284 286 pad = '4%' if self.ncols == 1 else '8%'
285 287
286 288 if self.oneFigure:
287 289 if self.height is None:
288 290 self.height = 1.4 * self.nrows + 1
289 291 fig = plt.figure(figsize=(self.width, self.height),
290 292 edgecolor='k',
291 293 facecolor='w')
292 294 self.figures.append(fig)
293 295 for n in range(self.nplots):
294 296 ax = fig.add_subplot(self.nrows, self.ncols,
295 297 n + 1, polar=self.polar)
296 298 ax.tick_params(labelsize=8)
297 299 ax.firsttime = True
298 300 ax.index = 0
299 301 ax.press = None
300 302 self.axes.append(ax)
301 303 if self.showprofile:
302 304 cax = self.__add_axes(ax, size=size, pad=pad)
303 305 cax.tick_params(labelsize=8)
304 306 self.pf_axes.append(cax)
305 307 else:
306 308 if self.height is None:
307 309 self.height = 3
308 310 for n in range(self.nplots):
309 311 fig = plt.figure(figsize=(self.width, self.height),
310 312 edgecolor='k',
311 313 facecolor='w')
312 314 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
313 315 ax.tick_params(labelsize=8)
314 316 ax.firsttime = True
315 317 ax.index = 0
316 318 ax.press = None
317 319 self.figures.append(fig)
318 320 self.axes.append(ax)
319 321 if self.showprofile:
320 322 cax = self.__add_axes(ax, size=size, pad=pad)
321 323 cax.tick_params(labelsize=8)
322 324 self.pf_axes.append(cax)
323 325
324 326 for n in range(self.nrows):
327 print(self.nrows)
325 328 if self.colormaps is not None:
326 329 cmap = plt.get_cmap(self.colormaps[n])
327 330 else:
328 331 cmap = plt.get_cmap(self.colormap)
329 332 cmap.set_bad(self.bgcolor, 1.)
330 333 self.cmaps.append(cmap)
331 334
332 335 def __add_axes(self, ax, size='30%', pad='8%'):
333 336 '''
334 337 Add new axes to the given figure
335 338 '''
336 339 divider = make_axes_locatable(ax)
337 340 nax = divider.new_horizontal(size=size, pad=pad)
338 341 ax.figure.add_axes(nax)
339 342 return nax
340 343
341 344 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
342 345 '''
343 346 Create a masked array for missing data
344 347 '''
345 348 if x_buffer.shape[0] < 2:
346 349 return x_buffer, y_buffer, z_buffer
347 350
348 351 deltas = x_buffer[1:] - x_buffer[0:-1]
349 352 x_median = numpy.median(deltas)
350 353
351 354 index = numpy.where(deltas > 5 * x_median)
352 355
353 356 if len(index[0]) != 0:
354 357 z_buffer[::, index[0], ::] = self.__missing
355 358 z_buffer = numpy.ma.masked_inside(z_buffer,
356 359 0.99 * self.__missing,
357 360 1.01 * self.__missing)
358 361
359 362 return x_buffer, y_buffer, z_buffer
360 363
361 364 def decimate(self):
362 365
363 366 # dx = int(len(self.x)/self.__MAXNUMX) + 1
364 367 dy = int(len(self.y) / self.decimation) + 1
365 368
366 369 # x = self.x[::dx]
367 370 x = self.x
368 371 y = self.y[::dy]
369 372 z = self.z[::, ::, ::dy]
370 373
371 374 return x, y, z
372 375
373 376 def format(self):
374 377 '''
375 378 Set min and max values, labels, ticks and titles
376 379 '''
377 380
378 381 for n, ax in enumerate(self.axes):
379 382 if ax.firsttime:
380 383 if self.xaxis != 'time':
381 384 xmin = self.xmin
382 385 xmax = self.xmax
383 386 else:
384 387 xmin = self.tmin
385 388 xmax = self.tmin + self.xrange*60*60
386 389 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
387 390 ax.xaxis.set_major_locator(LinearLocator(9))
388 391 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
389 392 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
390 393 ax.set_facecolor(self.bgcolor)
391 394 if self.xscale:
392 395 ax.xaxis.set_major_formatter(FuncFormatter(
393 396 lambda x, pos: '{0:g}'.format(x*self.xscale)))
394 397 if self.yscale:
395 398 ax.yaxis.set_major_formatter(FuncFormatter(
396 399 lambda x, pos: '{0:g}'.format(x*self.yscale)))
397 400 if self.xlabel is not None:
398 401 ax.set_xlabel(self.xlabel)
399 402 if self.ylabel is not None:
400 403 ax.set_ylabel(self.ylabel)
401 404 if self.showprofile:
402 405 self.pf_axes[n].set_ylim(ymin, ymax)
403 406 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
404 407 self.pf_axes[n].set_xlabel('dB')
405 408 self.pf_axes[n].grid(b=True, axis='x')
406 409 [tick.set_visible(False)
407 410 for tick in self.pf_axes[n].get_yticklabels()]
408 411 if self.colorbar:
409 412 ax.cbar = plt.colorbar(
410 413 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
411 414 ax.cbar.ax.tick_params(labelsize=8)
412 415 ax.cbar.ax.press = None
413 416 if self.cb_label:
414 417 ax.cbar.set_label(self.cb_label, size=8)
415 418 elif self.cb_labels:
416 419 ax.cbar.set_label(self.cb_labels[n], size=8)
417 420 else:
418 421 ax.cbar = None
419 422 ax.set_xlim(xmin, xmax)
420 423 ax.set_ylim(ymin, ymax)
421 424 ax.firsttime = False
422 425 if self.grid:
423 426 ax.grid(True)
424 427 if not self.polar:
425 428 ax.set_title('{} {} {}'.format(
426 429 self.titles[n],
427 430 self.getDateTime(self.data.max_time).strftime(
428 431 '%Y-%m-%d %H:%M:%S'),
429 432 self.time_label),
430 433 size=8)
431 434 else:
432 435 ax.set_title('{}'.format(self.titles[n]), size=8)
433 436 ax.set_ylim(0, 90)
434 437 ax.set_yticks(numpy.arange(0, 90, 20))
435 438 ax.yaxis.labelpad = 40
436 439
437 440 if self.firsttime:
438 441 for n, fig in enumerate(self.figures):
439 442 fig.subplots_adjust(**self.plots_adjust)
440 443 self.firsttime = False
441 444
442 445 def clear_figures(self):
443 446 '''
444 447 Reset axes for redraw plots
445 448 '''
446 449
447 450 for ax in self.axes+self.pf_axes+self.cb_axes:
448 451 ax.clear()
449 452 ax.firsttime = True
450 453 if hasattr(ax, 'cbar') and ax.cbar:
451 454 ax.cbar.remove()
452 455
453 456 def __plot(self):
454 457 '''
455 458 Main function to plot, format and save figures
456 459 '''
457 460
458 461 self.plot()
459 462 self.format()
460 463
461 464 for n, fig in enumerate(self.figures):
462 465 if self.nrows == 0 or self.nplots == 0:
463 466 log.warning('No data', self.name)
464 467 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
465 468 fig.canvas.manager.set_window_title(self.CODE)
466 469 continue
467 470
468 471 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
469 472 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
470 473 fig.canvas.draw()
471 474 if self.show:
472 475 fig.show()
473 476 figpause(0.01)
474 477
475 478 if self.save:
476 479 self.save_figure(n)
477 480
478 481 if self.server:
479 482 self.send_to_server()
480 483
481 484 def __update(self, dataOut, timestamp):
482 485 '''
483 486 '''
484 487
485 488 metadata = {
486 489 'yrange': dataOut.heightList,
487 490 'interval': dataOut.timeInterval,
488 491 'channels': dataOut.channelList
489 492 }
490 493
491 494 data, meta = self.update(dataOut)
492 495 metadata.update(meta)
493 496 self.data.update(data, timestamp, metadata)
494 497
495 498 def save_figure(self, n):
496 499 '''
497 500 '''
498 501
499 502 if (self.data.max_time - self.save_time) <= self.save_period:
500 503 return
501 504
502 505 self.save_time = self.data.max_time
503 506
504 507 fig = self.figures[n]
505 508
509 if self.throttle == 0:
506 510 figname = os.path.join(
507 511 self.save,
508 512 self.save_code,
509 513 '{}_{}.png'.format(
510 514 self.save_code,
511 515 self.getDateTime(self.data.max_time).strftime(
512 516 '%Y%m%d_%H%M%S'
513 517 ),
514 518 )
515 519 )
516 520 log.log('Saving figure: {}'.format(figname), self.name)
517 521 if not os.path.isdir(os.path.dirname(figname)):
518 522 os.makedirs(os.path.dirname(figname))
519 523 fig.savefig(figname)
520 524
521 if self.throttle == 0:
522 525 figname = os.path.join(
523 526 self.save,
524 527 '{}_{}.png'.format(
525 528 self.save_code,
526 529 self.getDateTime(self.data.min_time).strftime(
527 530 '%Y%m%d'
528 531 ),
529 532 )
530 533 )
531 534 fig.savefig(figname)
532 535
533 536 def send_to_server(self):
534 537 '''
535 538 '''
536 539
537 540 if self.exp_code == None:
538 541 log.warning('Missing `exp_code` skipping sending to server...')
539 542
540 543 last_time = self.data.max_time
541 544 interval = last_time - self.sender_time
542 545 if interval < self.sender_period:
543 546 return
544 547
545 548 self.sender_time = last_time
546 549
547 550 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
548 551 for attr in attrs:
549 552 value = getattr(self, attr)
550 553 if value:
551 554 if isinstance(value, (numpy.float32, numpy.float64)):
552 555 value = round(float(value), 2)
553 556 self.data.meta[attr] = value
554 557 if self.colormap == 'jet':
555 558 self.data.meta['colormap'] = 'Jet'
556 559 elif 'RdBu' in self.colormap:
557 560 self.data.meta['colormap'] = 'RdBu'
558 561 else:
559 562 self.data.meta['colormap'] = 'Viridis'
560 563 self.data.meta['interval'] = int(interval)
561 564
562 565 self.sender_queue.append(last_time)
563 566
564 567 while True:
565 568 try:
566 569 tm = self.sender_queue.popleft()
567 570 except IndexError:
568 571 break
569 572 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
570 573 self.socket.send_string(msg)
571 574 socks = dict(self.poll.poll(2000))
572 575 if socks.get(self.socket) == zmq.POLLIN:
573 576 reply = self.socket.recv_string()
574 577 if reply == 'ok':
575 578 log.log("Response from server ok", self.name)
576 579 time.sleep(0.1)
577 580 continue
578 581 else:
579 582 log.warning(
580 583 "Malformed reply from server: {}".format(reply), self.name)
581 584 else:
582 585 log.warning(
583 586 "No response from server, retrying...", self.name)
584 587 self.sender_queue.appendleft(tm)
585 588 self.socket.setsockopt(zmq.LINGER, 0)
586 589 self.socket.close()
587 590 self.poll.unregister(self.socket)
588 591 self.socket = self.context.socket(zmq.REQ)
589 592 self.socket.connect(self.server)
590 593 self.poll.register(self.socket, zmq.POLLIN)
591 594 break
592 595
593 596 def setup(self):
594 597 '''
595 598 This method should be implemented in the child class, the following
596 599 attributes should be set:
597 600
598 601 self.nrows: number of rows
599 602 self.ncols: number of cols
600 603 self.nplots: number of plots (channels or pairs)
601 604 self.ylabel: label for Y axes
602 605 self.titles: list of axes title
603 606
604 607 '''
605 608 raise NotImplementedError
606 609
607 610 def plot(self):
608 611 '''
609 612 Must be defined in the child class, the actual plotting method
610 613 '''
611 614 raise NotImplementedError
612 615
613 616 def update(self, dataOut):
614 617 '''
615 618 Must be defined in the child class, update self.data with new data
616 619 '''
617 620
618 621 data = {
619 622 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
620 623 }
621 624 meta = {}
622 625
623 626 return data, meta
624 627
625 628 def run(self, dataOut, **kwargs):
626 629 '''
627 630 Main plotting routine
628 631 '''
629 632
630 633 if self.isConfig is False:
631 634 self.__setup(**kwargs)
632 635
633 636 if self.localtime:
634 637 self.getDateTime = datetime.datetime.fromtimestamp
635 638 else:
636 639 self.getDateTime = datetime.datetime.utcfromtimestamp
637 640
638 641 self.data.setup()
639 642 self.isConfig = True
640 643 if self.server:
641 644 self.context = zmq.Context()
642 645 self.socket = self.context.socket(zmq.REQ)
643 646 self.socket.connect(self.server)
644 647 self.poll = zmq.Poller()
645 648 self.poll.register(self.socket, zmq.POLLIN)
646 649
647 650 tm = getattr(dataOut, self.attr_time)
648 651
649 652 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
650 653 self.save_time = tm
651 654 self.__plot()
652 655 self.tmin += self.xrange*60*60
653 656 self.data.setup()
654 657 self.clear_figures()
655 658
656 659 self.__update(dataOut, tm)
657 660
658 661 if self.isPlotConfig is False:
659 662 self.__setup_plot()
660 663 self.isPlotConfig = True
661 664 if self.xaxis == 'time':
662 665 dt = self.getDateTime(tm)
663 666 if self.xmin is None:
664 667 self.tmin = tm
665 668 self.xmin = dt.hour
666 669 minutes = (self.xmin-int(self.xmin)) * 60
667 670 seconds = (minutes - int(minutes)) * 60
668 671 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
669 672 datetime.datetime(1970, 1, 1)).total_seconds()
670 673 if self.localtime:
671 674 self.tmin += time.timezone
672 675
673 676 if self.xmin is not None and self.xmax is not None:
674 677 self.xrange = self.xmax - self.xmin
675 678
676 679 if self.throttle == 0:
677 680 self.__plot()
678 681 else:
679 682 self.__throttle_plot(self.__plot)#, coerce=coerce)
680 683
681 684 def close(self):
682 685
683 686 if self.data and not self.data.flagNoData:
684 687 self.save_time = self.data.max_time
685 688 self.__plot()
686 689 if self.data and not self.data.flagNoData and self.pause:
687 690 figpause(10)
688 691
Show file before | Show file after | Hide comments
@@ -1,358 +1,358
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot
7 7 from schainpy.utils import log
8 8
9 9 EARTH_RADIUS = 6.3710e3
10 10
11 11
12 12 def ll2xy(lat1, lon1, lat2, lon2):
13 13
14 14 p = 0.017453292519943295
15 15 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
16 16 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
17 17 r = 12742 * numpy.arcsin(numpy.sqrt(a))
18 18 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
19 19 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
20 20 theta = -theta + numpy.pi/2
21 21 return r*numpy.cos(theta), r*numpy.sin(theta)
22 22
23 23
24 24 def km2deg(km):
25 25 '''
26 26 Convert distance in km to degrees
27 27 '''
28 28
29 29 return numpy.rad2deg(km/EARTH_RADIUS)
30 30
31 31
32 32
33 33 class SpectralMomentsPlot(SpectraPlot):
34 34 '''
35 35 Plot for Spectral Moments
36 36 '''
37 37 CODE = 'spc_moments'
38 38 colormap = 'jet'
39 39 plot_type = 'pcolor'
40 40
41 41
42 42 class SnrPlot(RTIPlot):
43 43 '''
44 44 Plot for SNR Data
45 45 '''
46 46
47 47 CODE = 'snr'
48 48 colormap = 'jet'
49 49
50 50 def update(self, dataOut):
51 51
52 52 data = {
53 53 'snr': 10*numpy.log10(dataOut.data_snr)
54 54 }
55 55
56 56 return data, {}
57 57
58 58 class DopplerPlot(RTIPlot):
59 59 '''
60 60 Plot for DOPPLER Data (1st moment)
61 61 '''
62 62
63 63 CODE = 'dop'
64 64 colormap = 'jet'
65 65
66 66 def update(self, dataOut):
67 67
68 68 data = {
69 69 'dop': 10*numpy.log10(dataOut.data_dop)
70 70 }
71 71
72 72 return data, {}
73 73
74 74 class PowerPlot(RTIPlot):
75 75 '''
76 76 Plot for Power Data (0 moment)
77 77 '''
78 78
79 79 CODE = 'pow'
80 80 colormap = 'jet'
81 81
82 82 def update(self, dataOut):
83 83
84 84 data = {
85 85 'pow': 10*numpy.log10(dataOut.data_pow)
86 86 }
87 87
88 88 return data, {}
89 89
90 90 class SpectralWidthPlot(RTIPlot):
91 91 '''
92 92 Plot for Spectral Width Data (2nd moment)
93 93 '''
94 94
95 95 CODE = 'width'
96 96 colormap = 'jet'
97 97
98 98 def update(self, dataOut):
99 99
100 100 data = {
101 101 'width': dataOut.data_width
102 102 }
103 103
104 104 return data, {}
105 105
106 106 class SkyMapPlot(Plot):
107 107 '''
108 108 Plot for meteors detection data
109 109 '''
110 110
111 111 CODE = 'param'
112 112
113 113 def setup(self):
114 114
115 115 self.ncols = 1
116 116 self.nrows = 1
117 117 self.width = 7.2
118 118 self.height = 7.2
119 119 self.nplots = 1
120 120 self.xlabel = 'Zonal Zenith Angle (deg)'
121 121 self.ylabel = 'Meridional Zenith Angle (deg)'
122 122 self.polar = True
123 123 self.ymin = -180
124 124 self.ymax = 180
125 125 self.colorbar = False
126 126
127 127 def plot(self):
128 128
129 129 arrayParameters = numpy.concatenate(self.data['param'])
130 130 error = arrayParameters[:, -1]
131 131 indValid = numpy.where(error == 0)[0]
132 132 finalMeteor = arrayParameters[indValid, :]
133 133 finalAzimuth = finalMeteor[:, 3]
134 134 finalZenith = finalMeteor[:, 4]
135 135
136 136 x = finalAzimuth * numpy.pi / 180
137 137 y = finalZenith
138 138
139 139 ax = self.axes[0]
140 140
141 141 if ax.firsttime:
142 142 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
143 143 else:
144 144 ax.plot.set_data(x, y)
145 145
146 146 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
147 147 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
148 148 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
149 149 dt2,
150 150 len(x))
151 151 self.titles[0] = title
152 152
153 153
154 154 class GenericRTIPlot(Plot):
155 155 '''
156 156 Plot for data_xxxx object
157 157 '''
158 158
159 159 CODE = 'param'
160 160 colormap = 'viridis'
161 161 plot_type = 'pcolorbuffer'
162 162
163 163 def setup(self):
164 164 self.xaxis = 'time'
165 165 self.ncols = 1
166 self.nrows = self.data.shape(self.attr_data)[0]
166 self.nrows = self.data.shape('param')[0]
167 167 self.nplots = self.nrows
168 168 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
169 169
170 170 if not self.xlabel:
171 171 self.xlabel = 'Time'
172 172
173 173 self.ylabel = 'Height [km]'
174 174 if not self.titles:
175 175 self.titles = self.data.parameters \
176 176 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
177 177
178 178 def update(self, dataOut):
179 179
180 180 data = {
181 self.attr_data : getattr(dataOut, self.attr_data)
181 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
182 182 }
183 183
184 184 meta = {}
185 185
186 186 return data, meta
187 187
188 188 def plot(self):
189 189 # self.data.normalize_heights()
190 190 self.x = self.data.times
191 191 self.y = self.data.yrange
192 self.z = self.data[self.attr_data]
192 self.z = self.data['param']
193 193
194 194 self.z = numpy.ma.masked_invalid(self.z)
195 195
196 196 if self.decimation is None:
197 197 x, y, z = self.fill_gaps(self.x, self.y, self.z)
198 198 else:
199 199 x, y, z = self.fill_gaps(*self.decimate())
200 200
201 201 for n, ax in enumerate(self.axes):
202 202
203 203 self.zmax = self.zmax if self.zmax is not None else numpy.max(
204 204 self.z[n])
205 205 self.zmin = self.zmin if self.zmin is not None else numpy.min(
206 206 self.z[n])
207 207
208 208 if ax.firsttime:
209 209 if self.zlimits is not None:
210 210 self.zmin, self.zmax = self.zlimits[n]
211 211
212 212 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
213 213 vmin=self.zmin,
214 214 vmax=self.zmax,
215 215 cmap=self.cmaps[n]
216 216 )
217 217 else:
218 218 if self.zlimits is not None:
219 219 self.zmin, self.zmax = self.zlimits[n]
220 220 ax.collections.remove(ax.collections[0])
221 221 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
222 222 vmin=self.zmin,
223 223 vmax=self.zmax,
224 224 cmap=self.cmaps[n]
225 225 )
226 226
227 227
228 228 class PolarMapPlot(Plot):
229 229 '''
230 230 Plot for weather radar
231 231 '''
232 232
233 233 CODE = 'param'
234 234 colormap = 'seismic'
235 235
236 236 def setup(self):
237 237 self.ncols = 1
238 238 self.nrows = 1
239 239 self.width = 9
240 240 self.height = 8
241 241 self.mode = self.data.meta['mode']
242 242 if self.channels is not None:
243 243 self.nplots = len(self.channels)
244 244 self.nrows = len(self.channels)
245 245 else:
246 246 self.nplots = self.data.shape(self.CODE)[0]
247 247 self.nrows = self.nplots
248 248 self.channels = list(range(self.nplots))
249 249 if self.mode == 'E':
250 250 self.xlabel = 'Longitude'
251 251 self.ylabel = 'Latitude'
252 252 else:
253 253 self.xlabel = 'Range (km)'
254 254 self.ylabel = 'Height (km)'
255 255 self.bgcolor = 'white'
256 256 self.cb_labels = self.data.meta['units']
257 257 self.lat = self.data.meta['latitude']
258 258 self.lon = self.data.meta['longitude']
259 259 self.xmin, self.xmax = float(
260 260 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
261 261 self.ymin, self.ymax = float(
262 262 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
263 263 # self.polar = True
264 264
265 265 def plot(self):
266 266
267 267 for n, ax in enumerate(self.axes):
268 268 data = self.data['param'][self.channels[n]]
269 269
270 270 zeniths = numpy.linspace(
271 271 0, self.data.meta['max_range'], data.shape[1])
272 272 if self.mode == 'E':
273 273 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
274 274 r, theta = numpy.meshgrid(zeniths, azimuths)
275 275 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
276 276 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
277 277 x = km2deg(x) + self.lon
278 278 y = km2deg(y) + self.lat
279 279 else:
280 280 azimuths = numpy.radians(self.data.yrange)
281 281 r, theta = numpy.meshgrid(zeniths, azimuths)
282 282 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
283 283 self.y = zeniths
284 284
285 285 if ax.firsttime:
286 286 if self.zlimits is not None:
287 287 self.zmin, self.zmax = self.zlimits[n]
288 288 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
289 289 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
290 290 vmin=self.zmin,
291 291 vmax=self.zmax,
292 292 cmap=self.cmaps[n])
293 293 else:
294 294 if self.zlimits is not None:
295 295 self.zmin, self.zmax = self.zlimits[n]
296 296 ax.collections.remove(ax.collections[0])
297 297 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
298 298 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
299 299 vmin=self.zmin,
300 300 vmax=self.zmax,
301 301 cmap=self.cmaps[n])
302 302
303 303 if self.mode == 'A':
304 304 continue
305 305
306 306 # plot district names
307 307 f = open('/data/workspace/schain_scripts/distrito.csv')
308 308 for line in f:
309 309 label, lon, lat = [s.strip() for s in line.split(',') if s]
310 310 lat = float(lat)
311 311 lon = float(lon)
312 312 # ax.plot(lon, lat, '.b', ms=2)
313 313 ax.text(lon, lat, label.decode('utf8'), ha='center',
314 314 va='bottom', size='8', color='black')
315 315
316 316 # plot limites
317 317 limites = []
318 318 tmp = []
319 319 for line in open('/data/workspace/schain_scripts/lima.csv'):
320 320 if '#' in line:
321 321 if tmp:
322 322 limites.append(tmp)
323 323 tmp = []
324 324 continue
325 325 values = line.strip().split(',')
326 326 tmp.append((float(values[0]), float(values[1])))
327 327 for points in limites:
328 328 ax.add_patch(
329 329 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
330 330
331 331 # plot Cuencas
332 332 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
333 333 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
334 334 values = [line.strip().split(',') for line in f]
335 335 points = [(float(s[0]), float(s[1])) for s in values]
336 336 ax.add_patch(Polygon(points, ec='b', fc='none'))
337 337
338 338 # plot grid
339 339 for r in (15, 30, 45, 60):
340 340 ax.add_artist(plt.Circle((self.lon, self.lat),
341 341 km2deg(r), color='0.6', fill=False, lw=0.2))
342 342 ax.text(
343 343 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
344 344 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
345 345 '{}km'.format(r),
346 346 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
347 347
348 348 if self.mode == 'E':
349 349 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
350 350 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
351 351 else:
352 352 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
353 353 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
354 354
355 355 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
356 356 self.titles = ['{} {}'.format(
357 357 self.data.parameters[x], title) for x in self.channels]
358 358
Show file before | Show file after | Hide comments
@@ -1,702 +1,702
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13
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 self.nplots = len(self.data.channels)
27 27 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
28 28 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
29 29 self.height = 2.6 * self.nrows
30 30 self.cb_label = 'dB'
31 31 if self.showprofile:
32 32 self.width = 4 * self.ncols
33 33 else:
34 34 self.width = 3.5 * self.ncols
35 35 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
36 36 self.ylabel = 'Range [km]'
37 37
38 38 def update(self, dataOut):
39 39
40 40 data = {}
41 41 meta = {}
42 42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 43 data['spc'] = spc
44 44 data['rti'] = dataOut.getPower()
45 45 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 46 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
47 47 if self.CODE == 'spc_moments':
48 48 data['moments'] = dataOut.moments
49 49
50 50 return data, meta
51 51
52 52 def plot(self):
53 53 if self.xaxis == "frequency":
54 54 x = self.data.xrange[0]
55 55 self.xlabel = "Frequency (kHz)"
56 56 elif self.xaxis == "time":
57 57 x = self.data.xrange[1]
58 58 self.xlabel = "Time (ms)"
59 59 else:
60 60 x = self.data.xrange[2]
61 61 self.xlabel = "Velocity (m/s)"
62 62
63 63 if self.CODE == 'spc_moments':
64 64 x = self.data.xrange[2]
65 65 self.xlabel = "Velocity (m/s)"
66 66
67 67 self.titles = []
68 68
69 69 y = self.data.yrange
70 70 self.y = y
71 71
72 72 data = self.data[-1]
73 73 z = data['spc']
74 74
75 75 for n, ax in enumerate(self.axes):
76 76 noise = data['noise'][n]
77 77 if self.CODE == 'spc_moments':
78 mean = data['moments'][n, 2]
78 mean = data['moments'][n, 1]
79 79 if ax.firsttime:
80 80 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
81 81 self.xmin = self.xmin if self.xmin else -self.xmax
82 82 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
83 83 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
84 84 ax.plt = ax.pcolormesh(x, y, z[n].T,
85 85 vmin=self.zmin,
86 86 vmax=self.zmax,
87 87 cmap=plt.get_cmap(self.colormap)
88 88 )
89 89
90 90 if self.showprofile:
91 91 ax.plt_profile = self.pf_axes[n].plot(
92 92 data['rti'][n], y)[0]
93 93 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
94 94 color="k", linestyle="dashed", lw=1)[0]
95 95 if self.CODE == 'spc_moments':
96 96 ax.plt_mean = ax.plot(mean, y, color='k')[0]
97 97 else:
98 98 ax.plt.set_array(z[n].T.ravel())
99 99 if self.showprofile:
100 100 ax.plt_profile.set_data(data['rti'][n], y)
101 101 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
102 102 if self.CODE == 'spc_moments':
103 103 ax.plt_mean.set_data(mean, y)
104 104 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
105 105
106 106
107 107 class CrossSpectraPlot(Plot):
108 108
109 109 CODE = 'cspc'
110 110 colormap = 'jet'
111 111 plot_type = 'pcolor'
112 112 zmin_coh = None
113 113 zmax_coh = None
114 114 zmin_phase = None
115 115 zmax_phase = None
116 116
117 117 def setup(self):
118 118
119 119 self.ncols = 4
120 120 self.nplots = len(self.data.pairs) * 2
121 121 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
122 122 self.width = 3.1 * self.ncols
123 123 self.height = 2.6 * self.nrows
124 124 self.ylabel = 'Range [km]'
125 125 self.showprofile = False
126 126 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
127 127
128 128 def update(self, dataOut):
129 129
130 130 data = {}
131 131 meta = {}
132 132
133 133 spc = dataOut.data_spc
134 134 cspc = dataOut.data_cspc
135 135 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
136 136 meta['pairs'] = dataOut.pairsList
137 137
138 138 tmp = []
139 139
140 140 for n, pair in enumerate(meta['pairs']):
141 141 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
142 142 coh = numpy.abs(out)
143 143 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
144 144 tmp.append(coh)
145 145 tmp.append(phase)
146 146
147 147 data['cspc'] = numpy.array(tmp)
148 148
149 149 return data, meta
150 150
151 151 def plot(self):
152 152
153 153 if self.xaxis == "frequency":
154 154 x = self.data.xrange[0]
155 155 self.xlabel = "Frequency (kHz)"
156 156 elif self.xaxis == "time":
157 157 x = self.data.xrange[1]
158 158 self.xlabel = "Time (ms)"
159 159 else:
160 160 x = self.data.xrange[2]
161 161 self.xlabel = "Velocity (m/s)"
162 162
163 163 self.titles = []
164 164
165 165 y = self.data.yrange
166 166 self.y = y
167 167
168 168 data = self.data[-1]
169 169 cspc = data['cspc']
170 170
171 171 for n in range(len(self.data.pairs)):
172 172 pair = self.data.pairs[n]
173 173 coh = cspc[n*2]
174 174 phase = cspc[n*2+1]
175 175 ax = self.axes[2 * n]
176 176 if ax.firsttime:
177 177 ax.plt = ax.pcolormesh(x, y, coh.T,
178 178 vmin=0,
179 179 vmax=1,
180 180 cmap=plt.get_cmap(self.colormap_coh)
181 181 )
182 182 else:
183 183 ax.plt.set_array(coh.T.ravel())
184 184 self.titles.append(
185 185 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
186 186
187 187 ax = self.axes[2 * n + 1]
188 188 if ax.firsttime:
189 189 ax.plt = ax.pcolormesh(x, y, phase.T,
190 190 vmin=-180,
191 191 vmax=180,
192 192 cmap=plt.get_cmap(self.colormap_phase)
193 193 )
194 194 else:
195 195 ax.plt.set_array(phase.T.ravel())
196 196 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
197 197
198 198
199 199 class RTIPlot(Plot):
200 200 '''
201 201 Plot for RTI data
202 202 '''
203 203
204 204 CODE = 'rti'
205 205 colormap = 'jet'
206 206 plot_type = 'pcolorbuffer'
207 207
208 208 def setup(self):
209 209 self.xaxis = 'time'
210 210 self.ncols = 1
211 211 self.nrows = len(self.data.channels)
212 212 self.nplots = len(self.data.channels)
213 213 self.ylabel = 'Range [km]'
214 214 self.xlabel = 'Time'
215 215 self.cb_label = 'dB'
216 216 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
217 217 self.titles = ['{} Channel {}'.format(
218 218 self.CODE.upper(), x) for x in range(self.nrows)]
219 219
220 220 def update(self, dataOut):
221 221
222 222 data = {}
223 223 meta = {}
224 224 data['rti'] = dataOut.getPower()
225 225 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
226 226
227 227 return data, meta
228 228
229 229 def plot(self):
230 230 self.x = self.data.times
231 231 self.y = self.data.yrange
232 232 self.z = self.data[self.CODE]
233 233 self.z = numpy.ma.masked_invalid(self.z)
234 234
235 235 if self.decimation is None:
236 236 x, y, z = self.fill_gaps(self.x, self.y, self.z)
237 237 else:
238 238 x, y, z = self.fill_gaps(*self.decimate())
239 239
240 240 for n, ax in enumerate(self.axes):
241 241 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
242 242 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
243 243 data = self.data[-1]
244 244 if ax.firsttime:
245 245 ax.plt = ax.pcolormesh(x, y, z[n].T,
246 246 vmin=self.zmin,
247 247 vmax=self.zmax,
248 248 cmap=plt.get_cmap(self.colormap)
249 249 )
250 250 if self.showprofile:
251 251 ax.plot_profile = self.pf_axes[n].plot(
252 252 data['rti'][n], self.y)[0]
253 253 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
254 254 color="k", linestyle="dashed", lw=1)[0]
255 255 else:
256 256 ax.collections.remove(ax.collections[0])
257 257 ax.plt = ax.pcolormesh(x, y, z[n].T,
258 258 vmin=self.zmin,
259 259 vmax=self.zmax,
260 260 cmap=plt.get_cmap(self.colormap)
261 261 )
262 262 if self.showprofile:
263 263 ax.plot_profile.set_data(data['rti'][n], self.y)
264 264 ax.plot_noise.set_data(numpy.repeat(
265 265 data['noise'][n], len(self.y)), self.y)
266 266
267 267
268 268 class CoherencePlot(RTIPlot):
269 269 '''
270 270 Plot for Coherence data
271 271 '''
272 272
273 273 CODE = 'coh'
274 274
275 275 def setup(self):
276 276 self.xaxis = 'time'
277 277 self.ncols = 1
278 278 self.nrows = len(self.data.pairs)
279 279 self.nplots = len(self.data.pairs)
280 280 self.ylabel = 'Range [km]'
281 281 self.xlabel = 'Time'
282 282 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
283 283 if self.CODE == 'coh':
284 284 self.cb_label = ''
285 285 self.titles = [
286 286 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
287 287 else:
288 288 self.cb_label = 'Degrees'
289 289 self.titles = [
290 290 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
291 291
292 292 def update(self, dataOut):
293 293
294 294 data = {}
295 295 meta = {}
296 296 data['coh'] = dataOut.getCoherence()
297 297 meta['pairs'] = dataOut.pairsList
298 298
299 299 return data, meta
300 300
301 301 class PhasePlot(CoherencePlot):
302 302 '''
303 303 Plot for Phase map data
304 304 '''
305 305
306 306 CODE = 'phase'
307 307 colormap = 'seismic'
308 308
309 309 def update(self, dataOut):
310 310
311 311 data = {}
312 312 meta = {}
313 313 data['phase'] = dataOut.getCoherence(phase=True)
314 314 meta['pairs'] = dataOut.pairsList
315 315
316 316 return data, meta
317 317
318 318 class NoisePlot(Plot):
319 319 '''
320 320 Plot for noise
321 321 '''
322 322
323 323 CODE = 'noise'
324 324 plot_type = 'scatterbuffer'
325 325
326 326 def setup(self):
327 327 self.xaxis = 'time'
328 328 self.ncols = 1
329 329 self.nrows = 1
330 330 self.nplots = 1
331 331 self.ylabel = 'Intensity [dB]'
332 332 self.xlabel = 'Time'
333 333 self.titles = ['Noise']
334 334 self.colorbar = False
335 335 self.plots_adjust.update({'right': 0.85 })
336 336
337 337 def update(self, dataOut):
338 338
339 339 data = {}
340 340 meta = {}
341 341 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
342 342 meta['yrange'] = numpy.array([])
343 343
344 344 return data, meta
345 345
346 346 def plot(self):
347 347
348 348 x = self.data.times
349 349 xmin = self.data.min_time
350 350 xmax = xmin + self.xrange * 60 * 60
351 351 Y = self.data['noise']
352 352
353 353 if self.axes[0].firsttime:
354 354 self.ymin = numpy.nanmin(Y) - 5
355 355 self.ymax = numpy.nanmax(Y) + 5
356 356 for ch in self.data.channels:
357 357 y = Y[ch]
358 358 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
359 359 plt.legend(bbox_to_anchor=(1.18, 1.0))
360 360 else:
361 361 for ch in self.data.channels:
362 362 y = Y[ch]
363 363 self.axes[0].lines[ch].set_data(x, y)
364 364
365 365
366 366 class PowerProfilePlot(Plot):
367 367
368 368 CODE = 'pow_profile'
369 369 plot_type = 'scatter'
370 370
371 371 def setup(self):
372 372
373 373 self.ncols = 1
374 374 self.nrows = 1
375 375 self.nplots = 1
376 376 self.height = 4
377 377 self.width = 3
378 378 self.ylabel = 'Range [km]'
379 379 self.xlabel = 'Intensity [dB]'
380 380 self.titles = ['Power Profile']
381 381 self.colorbar = False
382 382
383 383 def update(self, dataOut):
384 384
385 385 data = {}
386 386 meta = {}
387 387 data[self.CODE] = dataOut.getPower()
388 388
389 389 return data, meta
390 390
391 391 def plot(self):
392 392
393 393 y = self.data.yrange
394 394 self.y = y
395 395
396 396 x = self.data[-1][self.CODE]
397 397
398 398 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
399 399 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
400 400
401 401 if self.axes[0].firsttime:
402 402 for ch in self.data.channels:
403 403 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
404 404 plt.legend()
405 405 else:
406 406 for ch in self.data.channels:
407 407 self.axes[0].lines[ch].set_data(x[ch], y)
408 408
409 409
410 410 class SpectraCutPlot(Plot):
411 411
412 412 CODE = 'spc_cut'
413 413 plot_type = 'scatter'
414 414 buffering = False
415 415
416 416 def setup(self):
417 417
418 418 self.nplots = len(self.data.channels)
419 419 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
420 420 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
421 421 self.width = 3.4 * self.ncols + 1.5
422 422 self.height = 3 * self.nrows
423 423 self.ylabel = 'Power [dB]'
424 424 self.colorbar = False
425 425 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
426 426
427 427 def update(self, dataOut):
428 428
429 429 data = {}
430 430 meta = {}
431 431 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
432 432 data['spc'] = spc
433 433 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
434 434
435 435 return data, meta
436 436
437 437 def plot(self):
438 438 if self.xaxis == "frequency":
439 439 x = self.data.xrange[0][1:]
440 440 self.xlabel = "Frequency (kHz)"
441 441 elif self.xaxis == "time":
442 442 x = self.data.xrange[1]
443 443 self.xlabel = "Time (ms)"
444 444 else:
445 445 x = self.data.xrange[2]
446 446 self.xlabel = "Velocity (m/s)"
447 447
448 448 self.titles = []
449 449
450 450 y = self.data.yrange
451 451 z = self.data[-1]['spc']
452 452
453 453 if self.height_index:
454 454 index = numpy.array(self.height_index)
455 455 else:
456 456 index = numpy.arange(0, len(y), int((len(y))/9))
457 457
458 458 for n, ax in enumerate(self.axes):
459 459 if ax.firsttime:
460 460 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
461 461 self.xmin = self.xmin if self.xmin else -self.xmax
462 462 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
463 463 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
464 464 ax.plt = ax.plot(x, z[n, :, index].T)
465 465 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
466 466 self.figures[0].legend(ax.plt, labels, loc='center right')
467 467 else:
468 468 for i, line in enumerate(ax.plt):
469 line.set_data(x, z[n, :, i])
469 line.set_data(x, z[n, :, index[i]])
470 470 self.titles.append('CH {}'.format(n))
471 471
472 472
473 473 class BeaconPhase(Plot):
474 474
475 475 __isConfig = None
476 476 __nsubplots = None
477 477
478 478 PREFIX = 'beacon_phase'
479 479
480 480 def __init__(self):
481 481 Plot.__init__(self)
482 482 self.timerange = 24*60*60
483 483 self.isConfig = False
484 484 self.__nsubplots = 1
485 485 self.counter_imagwr = 0
486 486 self.WIDTH = 800
487 487 self.HEIGHT = 400
488 488 self.WIDTHPROF = 120
489 489 self.HEIGHTPROF = 0
490 490 self.xdata = None
491 491 self.ydata = None
492 492
493 493 self.PLOT_CODE = BEACON_CODE
494 494
495 495 self.FTP_WEI = None
496 496 self.EXP_CODE = None
497 497 self.SUB_EXP_CODE = None
498 498 self.PLOT_POS = None
499 499
500 500 self.filename_phase = None
501 501
502 502 self.figfile = None
503 503
504 504 self.xmin = None
505 505 self.xmax = None
506 506
507 507 def getSubplots(self):
508 508
509 509 ncol = 1
510 510 nrow = 1
511 511
512 512 return nrow, ncol
513 513
514 514 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
515 515
516 516 self.__showprofile = showprofile
517 517 self.nplots = nplots
518 518
519 519 ncolspan = 7
520 520 colspan = 6
521 521 self.__nsubplots = 2
522 522
523 523 self.createFigure(id = id,
524 524 wintitle = wintitle,
525 525 widthplot = self.WIDTH+self.WIDTHPROF,
526 526 heightplot = self.HEIGHT+self.HEIGHTPROF,
527 527 show=show)
528 528
529 529 nrow, ncol = self.getSubplots()
530 530
531 531 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
532 532
533 533 def save_phase(self, filename_phase):
534 534 f = open(filename_phase,'w+')
535 535 f.write('\n\n')
536 536 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
537 537 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
538 538 f.close()
539 539
540 540 def save_data(self, filename_phase, data, data_datetime):
541 541 f=open(filename_phase,'a')
542 542 timetuple_data = data_datetime.timetuple()
543 543 day = str(timetuple_data.tm_mday)
544 544 month = str(timetuple_data.tm_mon)
545 545 year = str(timetuple_data.tm_year)
546 546 hour = str(timetuple_data.tm_hour)
547 547 minute = str(timetuple_data.tm_min)
548 548 second = str(timetuple_data.tm_sec)
549 549 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
550 550 f.close()
551 551
552 552 def plot(self):
553 553 log.warning('TODO: Not yet implemented...')
554 554
555 555 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
556 556 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
557 557 timerange=None,
558 558 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
559 559 server=None, folder=None, username=None, password=None,
560 560 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
561 561
562 562 if dataOut.flagNoData:
563 563 return dataOut
564 564
565 565 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
566 566 return
567 567
568 568 if pairsList == None:
569 569 pairsIndexList = dataOut.pairsIndexList[:10]
570 570 else:
571 571 pairsIndexList = []
572 572 for pair in pairsList:
573 573 if pair not in dataOut.pairsList:
574 574 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
575 575 pairsIndexList.append(dataOut.pairsList.index(pair))
576 576
577 577 if pairsIndexList == []:
578 578 return
579 579
580 580 # if len(pairsIndexList) > 4:
581 581 # pairsIndexList = pairsIndexList[0:4]
582 582
583 583 hmin_index = None
584 584 hmax_index = None
585 585
586 586 if hmin != None and hmax != None:
587 587 indexes = numpy.arange(dataOut.nHeights)
588 588 hmin_list = indexes[dataOut.heightList >= hmin]
589 589 hmax_list = indexes[dataOut.heightList <= hmax]
590 590
591 591 if hmin_list.any():
592 592 hmin_index = hmin_list[0]
593 593
594 594 if hmax_list.any():
595 595 hmax_index = hmax_list[-1]+1
596 596
597 597 x = dataOut.getTimeRange()
598 598
599 599 thisDatetime = dataOut.datatime
600 600
601 601 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
602 602 xlabel = "Local Time"
603 603 ylabel = "Phase (degrees)"
604 604
605 605 update_figfile = False
606 606
607 607 nplots = len(pairsIndexList)
608 608 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
609 609 phase_beacon = numpy.zeros(len(pairsIndexList))
610 610 for i in range(nplots):
611 611 pair = dataOut.pairsList[pairsIndexList[i]]
612 612 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
613 613 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
614 614 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
615 615 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
616 616 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
617 617
618 618 if dataOut.beacon_heiIndexList:
619 619 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
620 620 else:
621 621 phase_beacon[i] = numpy.average(phase)
622 622
623 623 if not self.isConfig:
624 624
625 625 nplots = len(pairsIndexList)
626 626
627 627 self.setup(id=id,
628 628 nplots=nplots,
629 629 wintitle=wintitle,
630 630 showprofile=showprofile,
631 631 show=show)
632 632
633 633 if timerange != None:
634 634 self.timerange = timerange
635 635
636 636 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
637 637
638 638 if ymin == None: ymin = 0
639 639 if ymax == None: ymax = 360
640 640
641 641 self.FTP_WEI = ftp_wei
642 642 self.EXP_CODE = exp_code
643 643 self.SUB_EXP_CODE = sub_exp_code
644 644 self.PLOT_POS = plot_pos
645 645
646 646 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
647 647 self.isConfig = True
648 648 self.figfile = figfile
649 649 self.xdata = numpy.array([])
650 650 self.ydata = numpy.array([])
651 651
652 652 update_figfile = True
653 653
654 654 #open file beacon phase
655 655 path = '%s%03d' %(self.PREFIX, self.id)
656 656 beacon_file = os.path.join(path,'%s.txt'%self.name)
657 657 self.filename_phase = os.path.join(figpath,beacon_file)
658 658 #self.save_phase(self.filename_phase)
659 659
660 660
661 661 #store data beacon phase
662 662 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
663 663
664 664 self.setWinTitle(title)
665 665
666 666
667 667 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
668 668
669 669 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
670 670
671 671 axes = self.axesList[0]
672 672
673 673 self.xdata = numpy.hstack((self.xdata, x[0:1]))
674 674
675 675 if len(self.ydata)==0:
676 676 self.ydata = phase_beacon.reshape(-1,1)
677 677 else:
678 678 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
679 679
680 680
681 681 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
682 682 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
683 683 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
684 684 XAxisAsTime=True, grid='both'
685 685 )
686 686
687 687 self.draw()
688 688
689 689 if dataOut.ltctime >= self.xmax:
690 690 self.counter_imagwr = wr_period
691 691 self.isConfig = False
692 692 update_figfile = True
693 693
694 694 self.save(figpath=figpath,
695 695 figfile=figfile,
696 696 save=save,
697 697 ftp=ftp,
698 698 wr_period=wr_period,
699 699 thisDatetime=thisDatetime,
700 700 update_figfile=update_figfile)
701 701
702 702 return dataOut No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now