##// END OF EJS Templates
schain
RSS
New Updates
rflores -
r1561:3446c343bac9
parent child
Show More
Side by Side Unified
Context file:
r1561:3446c343bac9 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,702 +1,704
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.xlimits = kwargs.get('xlimits', None)
224 224 self.xstep_given = kwargs.get('xstep_given', None)
225 225 self.ystep_given = kwargs.get('ystep_given', None)
226 226 self.autoxticks = kwargs.get('autoxticks', True)
227 227 self.xmin = kwargs.get('xmin', None)
228 228 self.xmax = kwargs.get('xmax', None)
229 229 self.xrange = kwargs.get('xrange', 12)
230 230 self.xscale = kwargs.get('xscale', None)
231 231 self.ymin = kwargs.get('ymin', None)
232 232 self.ymax = kwargs.get('ymax', None)
233 233 self.yscale = kwargs.get('yscale', None)
234 234 self.xlabel = kwargs.get('xlabel', None)
235 235 self.attr_time = kwargs.get('attr_time', 'utctime')
236 236 self.attr_data = kwargs.get('attr_data', 'data_param')
237 237 self.decimation = kwargs.get('decimation', None)
238 238 self.oneFigure = kwargs.get('oneFigure', True)
239 239 self.width = kwargs.get('width', None)
240 240 self.height = kwargs.get('height', None)
241 241 self.colorbar = kwargs.get('colorbar', True)
242 242 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
243 243 self.channels = kwargs.get('channels', None)
244 244 self.titles = kwargs.get('titles', [])
245 245 self.polar = False
246 246 self.type = kwargs.get('type', 'iq')
247 247 self.grid = kwargs.get('grid', False)
248 248 self.pause = kwargs.get('pause', False)
249 249 self.save_code = kwargs.get('save_code', self.CODE)
250 250 self.throttle = kwargs.get('throttle', 0)
251 251 self.exp_code = kwargs.get('exp_code', None)
252 252 self.server = kwargs.get('server', False)
253 253 self.sender_period = kwargs.get('sender_period', 60)
254 254 self.tag = kwargs.get('tag', '')
255 255 self.height_index = kwargs.get('height_index', None)
256 256 self.__throttle_plot = apply_throttle(self.throttle)
257 257 code = self.attr_data if self.attr_data else self.CODE
258 258 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
259 259 #self.EEJtype = kwargs.get('EEJtype', 2)
260 260
261 261 if self.server:
262 262 if not self.server.startswith('tcp://'):
263 263 self.server = 'tcp://{}'.format(self.server)
264 264 log.success(
265 265 'Sending to server: {}'.format(self.server),
266 266 self.name
267 267 )
268 268
269 269 if isinstance(self.attr_data, str):
270 270 self.attr_data = [self.attr_data]
271 271
272 272 def __setup_plot(self):
273 273 '''
274 274 Common setup for all figures, here figures and axes are created
275 275 '''
276 276
277 277 self.setup()
278 278
279 279 self.time_label = 'LT' if self.localtime else 'UTC'
280 280
281 281 if self.width is None:
282 282 self.width = 8
283 283
284 284 self.figures = []
285 285 self.axes = []
286 286 self.cb_axes = []
287 287 self.pf_axes = []
288 288 self.cmaps = []
289 289
290 290 size = '15%' if self.ncols == 1 else '30%'
291 291 pad = '4%' if self.ncols == 1 else '8%'
292 292
293 293 if self.oneFigure:
294 294 if self.height is None:
295 295 self.height = 1.4 * self.nrows + 1
296 296 fig = plt.figure(figsize=(self.width, self.height),
297 297 edgecolor='k',
298 298 facecolor='w')
299 299 self.figures.append(fig)
300 300 for n in range(self.nplots):
301 301 ax = fig.add_subplot(self.nrows, self.ncols,
302 302 n + 1, polar=self.polar)
303 303 ax.tick_params(labelsize=8)
304 304 ax.firsttime = True
305 305 ax.index = 0
306 306 ax.press = None
307 307 self.axes.append(ax)
308 308 if self.showprofile:
309 309 cax = self.__add_axes(ax, size=size, pad=pad)
310 310 cax.tick_params(labelsize=8)
311 311 self.pf_axes.append(cax)
312 312 else:
313 313 if self.height is None:
314 314 self.height = 3
315 315 for n in range(self.nplots):
316 316 fig = plt.figure(figsize=(self.width, self.height),
317 317 edgecolor='k',
318 318 facecolor='w')
319 319 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
320 320 ax.tick_params(labelsize=8)
321 321 ax.firsttime = True
322 322 ax.index = 0
323 323 ax.press = None
324 324 self.figures.append(fig)
325 325 self.axes.append(ax)
326 326 if self.showprofile:
327 327 cax = self.__add_axes(ax, size=size, pad=pad)
328 328 cax.tick_params(labelsize=8)
329 329 self.pf_axes.append(cax)
330 330
331 331 for n in range(self.nrows):
332 332 if self.colormaps is not None:
333 333 cmap = plt.get_cmap(self.colormaps[n])
334 334 else:
335 335 cmap = plt.get_cmap(self.colormap)
336 336 cmap.set_bad(self.bgcolor, 1.)
337 337 self.cmaps.append(cmap)
338 338
339 339 def __add_axes(self, ax, size='30%', pad='8%'):
340 340 '''
341 341 Add new axes to the given figure
342 342 '''
343 343 divider = make_axes_locatable(ax)
344 344 nax = divider.new_horizontal(size=size, pad=pad)
345 345 ax.figure.add_axes(nax)
346 346 return nax
347 347
348 348 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
349 349 '''
350 350 Create a masked array for missing data
351 351 '''
352 352 if x_buffer.shape[0] < 2:
353 353 return x_buffer, y_buffer, z_buffer
354 354
355 355 deltas = x_buffer[1:] - x_buffer[0:-1]
356 356 x_median = numpy.median(deltas)
357 357
358 358 index = numpy.where(deltas > 5 * x_median)
359 359
360 360 if len(index[0]) != 0:
361 361 z_buffer[::, index[0], ::] = self.__missing
362 362 z_buffer = numpy.ma.masked_inside(z_buffer,
363 363 0.99 * self.__missing,
364 364 1.01 * self.__missing)
365 365
366 366 return x_buffer, y_buffer, z_buffer
367 367
368 368 def decimate(self):
369 369
370 370 # dx = int(len(self.x)/self.__MAXNUMX) + 1
371 371 dy = int(len(self.y) / self.decimation) + 1
372 372
373 373 # x = self.x[::dx]
374 374 x = self.x
375 375 y = self.y[::dy]
376 376 z = self.z[::, ::, ::dy]
377 377
378 378 return x, y, z
379 379
380 380 def format(self):
381 381 '''
382 382 Set min and max values, labels, ticks and titles
383 383 '''
384 384
385 385 for n, ax in enumerate(self.axes):
386 386 if ax.firsttime:
387 387 if self.xaxis != 'time':
388 388 xmin = self.xmin
389 389 xmax = self.xmax
390 390 else:
391 391 xmin = self.tmin
392 392 xmax = self.tmin + self.xrange*60*60
393 393 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
394 394 ax.xaxis.set_major_locator(LinearLocator(9))
395 395 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
396 396 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
397 397 ax.set_facecolor(self.bgcolor)
398 398 if self.xscale:
399 399 ax.xaxis.set_major_formatter(FuncFormatter(
400 400 lambda x, pos: '{0:g}'.format(x*self.xscale)))
401 401 if self.yscale:
402 402 ax.yaxis.set_major_formatter(FuncFormatter(
403 403 lambda x, pos: '{0:g}'.format(x*self.yscale)))
404 404 if self.xlabel is not None:
405 405 ax.set_xlabel(self.xlabel)
406 406 if self.ylabel is not None:
407 407 ax.set_ylabel(self.ylabel)
408 408 if self.showprofile:
409 409 if self.zlimits is not None:
410 410 self.zmin, self.zmax = self.zlimits[n]
411 411 self.pf_axes[n].set_ylim(ymin, ymax)
412 412 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
413 413 self.pf_axes[n].set_xlabel('dB')
414 414 self.pf_axes[n].grid(b=True, axis='x')
415 415 [tick.set_visible(False)
416 416 for tick in self.pf_axes[n].get_yticklabels()]
417 417 if self.colorbar:
418 418 ax.cbar = plt.colorbar(
419 419 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
420 420 ax.cbar.ax.tick_params(labelsize=8)
421 421 ax.cbar.ax.press = None
422 422 if self.cb_label:
423 423 ax.cbar.set_label(self.cb_label, size=8)
424 424 elif self.cb_labels:
425 425 ax.cbar.set_label(self.cb_labels[n], size=8)
426 426 else:
427 427 ax.cbar = None
428 428 ax.set_xlim(xmin, xmax)
429 429 ax.set_ylim(ymin, ymax)
430 430 ax.firsttime = False
431 431 if self.grid:
432 432 ax.grid(True)
433 433 if not self.polar:
434 434 ax.set_title('{} {} {}'.format(
435 435 self.titles[n],
436 436 self.getDateTime(self.data.max_time).strftime(
437 437 '%Y-%m-%d %H:%M:%S'),
438 438 self.time_label),
439 439 size=8)
440 440 else:
441 441 ax.set_title('{}'.format(self.titles[n]), size=8)
442 442 ax.set_ylim(0, 90)
443 443 ax.set_yticks(numpy.arange(0, 90, 20))
444 444 ax.yaxis.labelpad = 40
445 445
446 446 if self.firsttime:
447 447 for n, fig in enumerate(self.figures):
448 448 fig.subplots_adjust(**self.plots_adjust)
449 449 self.firsttime = False
450 450
451 451 def clear_figures(self):
452 452 '''
453 453 Reset axes for redraw plots
454 454 '''
455 455
456 456 for ax in self.axes+self.pf_axes+self.cb_axes:
457 457 ax.clear()
458 458 ax.firsttime = True
459 459 if hasattr(ax, 'cbar') and ax.cbar:
460 460 ax.cbar.remove()
461 461
462 462 def __plot(self):
463 463 '''
464 464 Main function to plot, format and save figures
465 465 '''
466 466
467 467 self.plot()
468 468 self.format()
469 469
470 470 for n, fig in enumerate(self.figures):
471 471 if self.nrows == 0 or self.nplots == 0:
472 472 log.warning('No data', self.name)
473 473 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
474 474 fig.canvas.manager.set_window_title(self.CODE)
475 475 continue
476 476
477 477 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
478 478 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
479 479 fig.canvas.draw()
480 480 if self.show:
481 481 fig.show()
482 482 figpause(0.01)
483 483
484 484 if self.save:
485 485 self.save_figure(n)
486 486
487 487 if self.server:
488 488 self.send_to_server()
489 489
490 490 def __update(self, dataOut, timestamp):
491 491 '''
492 492 '''
493 493
494 494 metadata = {
495 495 'yrange': dataOut.heightList,
496 496 'interval': dataOut.timeInterval,
497 497 'channels': dataOut.channelList
498 498 }
499 499
500 500 data, meta = self.update(dataOut)
501 501 metadata.update(meta)
502 502 self.data.update(data, timestamp, metadata)
503 503
504 504 def save_figure(self, n):
505 505 '''
506 506 '''
507 507
508 508 if (self.data.max_time - self.save_time) <= self.save_period:
509 509 return
510 510
511 511 self.save_time = self.data.max_time
512 512
513 513 fig = self.figures[n]
514 514
515 515 if self.throttle == 0:
516 516 figname = os.path.join(
517 517 self.save,
518 518 self.save_code,
519 519 '{}_{}.png'.format(
520 520 self.save_code,
521 521 self.getDateTime(self.data.max_time).strftime(
522 522 '%Y%m%d_%H%M%S'
523 523 ),
524 524 )
525 525 )
526 526 log.log('Saving figure: {}'.format(figname), self.name)
527 527 if not os.path.isdir(os.path.dirname(figname)):
528 528 os.makedirs(os.path.dirname(figname))
529 529 fig.savefig(figname)
530 530
531 531 figname = os.path.join(
532 532 self.save,
533 533 #self.save_code,
534 534 '{}_{}.png'.format(
535 535 self.save_code,
536 536 self.getDateTime(self.data.min_time).strftime(
537 537 '%Y%m%d'
538 538 ),
539 539 )
540 540 )
541 541 log.log('Saving figure: {}'.format(figname), self.name)
542 542 if not os.path.isdir(os.path.dirname(figname)):
543 543 os.makedirs(os.path.dirname(figname))
544 544 fig.savefig(figname)
545 545
546 546 def send_to_server(self):
547 547 '''
548 548 '''
549 549
550 550 if self.exp_code == None:
551 551 log.warning('Missing `exp_code` skipping sending to server...')
552 552
553 553 last_time = self.data.max_time
554 554 interval = last_time - self.sender_time
555 555 if interval < self.sender_period:
556 556 return
557 557
558 558 self.sender_time = last_time
559 559
560 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
560 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax', 'zlimits']
561 561 for attr in attrs:
562 562 value = getattr(self, attr)
563 563 if value:
564 564 if isinstance(value, (numpy.float32, numpy.float64)):
565 565 value = round(float(value), 2)
566 566 self.data.meta[attr] = value
567 567 if self.colormap == 'jet':
568 568 self.data.meta['colormap'] = 'Jet'
569 569 elif 'RdBu' in self.colormap:
570 570 self.data.meta['colormap'] = 'RdBu'
571 571 else:
572 572 self.data.meta['colormap'] = 'Viridis'
573 573 self.data.meta['interval'] = int(interval)
574 574 #print(last_time)
575 575 #print(time.time())
576 576 #exit(1)
577 577 self.sender_queue.append(last_time)
578 578
579 579 while True:
580 580 try:
581 581 tm = self.sender_queue.popleft()
582 582 except IndexError:
583 583 break
584 584 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
585 585 self.socket.send_string(msg)
586 586 socks = dict(self.poll.poll(2000))
587 587 if socks.get(self.socket) == zmq.POLLIN:
588 588 reply = self.socket.recv_string()
589 589 if reply == 'ok':
590 590 log.log("Response from server ok", self.name)
591 591 time.sleep(0.1)
592 592 continue
593 593 else:
594 594 log.warning(
595 595 "Malformed reply from server: {}".format(reply), self.name)
596 596 else:
597 597 log.warning(
598 598 "No response from server, retrying...", self.name)
599 599 self.sender_queue.appendleft(tm)
600 600 self.socket.setsockopt(zmq.LINGER, 0)
601 601 self.socket.close()
602 602 self.poll.unregister(self.socket)
603 603 self.socket = self.context.socket(zmq.REQ)
604 604 self.socket.connect(self.server)
605 605 self.poll.register(self.socket, zmq.POLLIN)
606 606 break
607 607
608 608 def setup(self):
609 609 '''
610 610 This method should be implemented in the child class, the following
611 611 attributes should be set:
612 612
613 613 self.nrows: number of rows
614 614 self.ncols: number of cols
615 615 self.nplots: number of plots (channels or pairs)
616 616 self.ylabel: label for Y axes
617 617 self.titles: list of axes title
618 618
619 619 '''
620 620 raise NotImplementedError
621 621
622 622 def plot(self):
623 623 '''
624 624 Must be defined in the child class, the actual plotting method
625 625 '''
626 626 raise NotImplementedError
627 627
628 628 def update(self, dataOut):
629 629 '''
630 630 Must be defined in the child class, update self.data with new data
631 631 '''
632 632
633 633 data = {
634 634 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
635 635 }
636 636 meta = {}
637 637
638 638 return data, meta
639 639
640 640 def run(self, dataOut, **kwargs):
641 641 '''
642 642 Main plotting routine
643 643 '''
644 644
645 645 if self.isConfig is False:
646 646 self.__setup(**kwargs)
647 647
648 648 if self.localtime:
649 649 self.getDateTime = datetime.datetime.fromtimestamp
650 650 else:
651 651 self.getDateTime = datetime.datetime.utcfromtimestamp
652 652
653 653 self.data.setup()
654 654 self.isConfig = True
655 655 if self.server:
656 656 self.context = zmq.Context()
657 657 self.socket = self.context.socket(zmq.REQ)
658 658 self.socket.connect(self.server)
659 659 self.poll = zmq.Poller()
660 660 self.poll.register(self.socket, zmq.POLLIN)
661 661
662 662 tm = getattr(dataOut, self.attr_time)
663
664 663 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
665 664 self.save_time = tm
666 665 self.__plot()
667 self.tmin += self.xrange*60*60
666 #self.tmin += self.xrange*60*60 #Modified by R. Flores
667 self.tmin += 24*60*60 #Modified by R. Flores
668
668 669 self.data.setup()
669 670 self.clear_figures()
670 671
671 672 self.__update(dataOut, tm)
672 673
673 674 if self.isPlotConfig is False:
674 675 self.__setup_plot()
675 676 self.isPlotConfig = True
676 677 if self.xaxis == 'time':
677 678 dt = self.getDateTime(tm)
679
678 680 if self.xmin is None:
679 681 self.tmin = tm
680 682 self.xmin = dt.hour
681 683 minutes = (self.xmin-int(self.xmin)) * 60
682 684 seconds = (minutes - int(minutes)) * 60
683 685 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
684 686 datetime.datetime(1970, 1, 1)).total_seconds()
685 687 if self.localtime:
686 688 self.tmin += time.timezone
687 689
688 690 if self.xmin is not None and self.xmax is not None:
689 691 self.xrange = self.xmax - self.xmin
690 692
691 693 if self.throttle == 0:
692 694 self.__plot()
693 695 else:
694 696 self.__throttle_plot(self.__plot)#, coerce=coerce)
695 697
696 698 def close(self):
697 699
698 700 if self.data and not self.data.flagNoData:
699 701 self.save_time = 0
700 702 self.__plot()
701 703 if self.data and not self.data.flagNoData and self.pause:
702 704 figpause(10)
Show file before | Show file after | Hide comments
@@ -1,1337 +1,1343
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 import collections.abc
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("Spc: ",spc[0])
46 46 #exit(1)
47 47 data['spc'] = spc
48 48 data['rti'] = dataOut.getPower()
49 49 #print(data['rti'][0])
50 50 #exit(1)
51 51 #print("NormFactor: ",dataOut.normFactor)
52 52 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
53 53 if hasattr(dataOut, 'LagPlot'): #Double Pulse
54 54 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
55 55 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=46,ymax_index=max_hei_id)/dataOut.normFactor)
56 56 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=40,ymax_index=max_hei_id)/dataOut.normFactor)
57 57 data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)/dataOut.normFactor)
58 58 data['noise'][0] = 10*numpy.log10(dataOut.getNoise(ymin_index=53)[0]/dataOut.normFactor)
59 59 #data['noise'][1] = 22.035507
60 60 else:
61 61 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
62 62 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=26,ymax_index=44)/dataOut.normFactor)
63 63 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
64 64
65 65 if self.CODE == 'spc_moments':
66 66 data['moments'] = dataOut.moments
67 67 if self.CODE == 'gaussian_fit':
68 68 data['gaussfit'] = dataOut.DGauFitParams
69 69
70 70 return data, meta
71 71
72 72 def plot(self):
73 73
74 74 if self.xaxis == "frequency":
75 75 x = self.data.xrange[0]
76 76 self.xlabel = "Frequency (kHz)"
77 77 elif self.xaxis == "time":
78 78 x = self.data.xrange[1]
79 79 self.xlabel = "Time (ms)"
80 80 else:
81 81 x = self.data.xrange[2]
82 82 self.xlabel = "Velocity (m/s)"
83 83
84 84 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
85 85 x = self.data.xrange[2]
86 86 self.xlabel = "Velocity (m/s)"
87 87
88 88 self.titles = []
89 89
90 90 y = self.data.yrange
91 91 self.y = y
92 92
93 93 data = self.data[-1]
94 94 z = data['spc']
95 95
96 96 self.CODE2 = 'spc_oblique'
97 97
98 98 for n, ax in enumerate(self.axes):
99 99 noise = data['noise'][n]
100 100 if self.CODE == 'spc_moments':
101 101 mean = data['moments'][n, 1]
102 102 if self.CODE == 'gaussian_fit':
103 103 gau0 = data['gaussfit'][n][2,:,0]
104 104 gau1 = data['gaussfit'][n][2,:,1]
105 105 if ax.firsttime:
106 106 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
107 107 self.xmin = self.xmin if self.xmin else numpy.nanmin(x)#-self.xmax
108 108 #self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
109 109 #self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
110 110 if self.zlimits is not None:
111 111 self.zmin, self.zmax = self.zlimits[n]
112 112
113 113 ax.plt = ax.pcolormesh(x, y, z[n].T,
114 114 vmin=self.zmin,
115 115 vmax=self.zmax,
116 116 cmap=plt.get_cmap(self.colormap),
117 117 )
118 118
119 119 if self.showprofile:
120 120 ax.plt_profile = self.pf_axes[n].plot(
121 121 data['rti'][n], y)[0]
122 122 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
123 123 color="k", linestyle="dashed", lw=1)[0]
124 124 if self.CODE == 'spc_moments':
125 125 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
126 126 if self.CODE == 'gaussian_fit':
127 127 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
128 128 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
129 129 else:
130 130 if self.zlimits is not None:
131 131 self.zmin, self.zmax = self.zlimits[n]
132 132 ax.plt.set_array(z[n].T.ravel())
133 133 if self.showprofile:
134 134 ax.plt_profile.set_data(data['rti'][n], y)
135 135 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
136 136 if self.CODE == 'spc_moments':
137 137 ax.plt_mean.set_data(mean, y)
138 138 if self.CODE == 'gaussian_fit':
139 139 ax.plt_gau0.set_data(gau0, y)
140 140 ax.plt_gau1.set_data(gau1, y)
141 141 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
142 142
143 143 class SpectraObliquePlot(Plot):
144 144 '''
145 145 Plot for Spectra data
146 146 '''
147 147
148 148 CODE = 'spc_oblique'
149 149 colormap = 'jet'
150 150 plot_type = 'pcolor'
151 151
152 152 def setup(self):
153 153 self.xaxis = "oblique"
154 154 self.nplots = len(self.data.channels)
155 155 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
156 156 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
157 157 self.height = 2.6 * self.nrows
158 158 self.cb_label = 'dB'
159 159 if self.showprofile:
160 160 self.width = 4 * self.ncols
161 161 else:
162 162 self.width = 3.5 * self.ncols
163 163 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
164 164 self.ylabel = 'Range [km]'
165 165
166 166 def update(self, dataOut):
167 167
168 168 data = {}
169 169 meta = {}
170 170 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
171 171 data['spc'] = spc
172 172 data['rti'] = dataOut.getPower()
173 173 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
174 174 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
175 175 '''
176 176 data['shift1'] = dataOut.Oblique_params[0,-2,:]
177 177 data['shift2'] = dataOut.Oblique_params[0,-1,:]
178 178 data['shift1_error'] = dataOut.Oblique_param_errors[0,-2,:]
179 179 data['shift2_error'] = dataOut.Oblique_param_errors[0,-1,:]
180 180 '''
181 181 '''
182 182 data['shift1'] = dataOut.Oblique_params[0,1,:]
183 183 data['shift2'] = dataOut.Oblique_params[0,4,:]
184 184 data['shift1_error'] = dataOut.Oblique_param_errors[0,1,:]
185 185 data['shift2_error'] = dataOut.Oblique_param_errors[0,4,:]
186 186 '''
187 187 data['shift1'] = dataOut.Dop_EEJ_T1[0]
188 188 data['shift2'] = dataOut.Dop_EEJ_T2[0]
189 data['max_val_2'] = dataOut.Oblique_params[0,-1,:]
189 190 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
190 191 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
191 192
192 193 return data, meta
193 194
194 195 def plot(self):
195 196
196 197 if self.xaxis == "frequency":
197 198 x = self.data.xrange[0]
198 199 self.xlabel = "Frequency (kHz)"
199 200 elif self.xaxis == "time":
200 201 x = self.data.xrange[1]
201 202 self.xlabel = "Time (ms)"
202 203 else:
203 204 x = self.data.xrange[2]
204 205 self.xlabel = "Velocity (m/s)"
205 206
206 207 self.titles = []
207 208
208 209 y = self.data.yrange
209 210 self.y = y
210 211
211 212 data = self.data[-1]
212 213 z = data['spc']
213 214
214 215 for n, ax in enumerate(self.axes):
215 216 noise = self.data['noise'][n][-1]
216 217 shift1 = data['shift1']
217 218 #print(shift1)
218 219 shift2 = data['shift2']
220 max_val_2 = data['max_val_2']
219 221 err1 = data['shift1_error']
220 222 err2 = data['shift2_error']
221 223 if ax.firsttime:
222 224
223 225 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
224 226 self.xmin = self.xmin if self.xmin else -self.xmax
225 227 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
226 228 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
227 229 ax.plt = ax.pcolormesh(x, y, z[n].T,
228 230 vmin=self.zmin,
229 231 vmax=self.zmax,
230 232 cmap=plt.get_cmap(self.colormap)
231 233 )
232 234
233 235 if self.showprofile:
234 236 ax.plt_profile = self.pf_axes[n].plot(
235 237 self.data['rti'][n][-1], y)[0]
236 238 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
237 239 color="k", linestyle="dashed", lw=1)[0]
238 240
239 241 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)
240 242 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)
243 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)
244
241 245 #print("plotter1: ", self.ploterr1,shift1)
242 246
243 247 else:
244 248 #print("else plotter1: ", self.ploterr1,shift1)
245 249 self.ploterr1.remove()
246 250 self.ploterr2.remove()
251 self.ploterr3.remove()
247 252 ax.plt.set_array(z[n].T.ravel())
248 253 if self.showprofile:
249 254 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
250 255 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
251 256 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)
252 257 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)
258 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)
253 259
254 260 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
255 261
256 262
257 263 class CrossSpectraPlot(Plot):
258 264
259 265 CODE = 'cspc'
260 266 colormap = 'jet'
261 267 plot_type = 'pcolor'
262 268 zmin_coh = None
263 269 zmax_coh = None
264 270 zmin_phase = None
265 271 zmax_phase = None
266 272
267 273 def setup(self):
268 274
269 275 self.ncols = 4
270 276 self.nplots = len(self.data.pairs) * 2
271 277 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
272 278 self.width = 3.1 * self.ncols
273 279 self.height = 5 * self.nrows
274 280 self.ylabel = 'Range [km]'
275 281 self.showprofile = False
276 282 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
277 283
278 284 def update(self, dataOut):
279 285
280 286 data = {}
281 287 meta = {}
282 288
283 289 spc = dataOut.data_spc
284 290 cspc = dataOut.data_cspc
285 291 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
286 292 meta['pairs'] = dataOut.pairsList
287 293
288 294 tmp = []
289 295
290 296 for n, pair in enumerate(meta['pairs']):
291 297 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
292 298 coh = numpy.abs(out)
293 299 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
294 300 tmp.append(coh)
295 301 tmp.append(phase)
296 302
297 303 data['cspc'] = numpy.array(tmp)
298 304
299 305 return data, meta
300 306
301 307 def plot(self):
302 308
303 309 if self.xaxis == "frequency":
304 310 x = self.data.xrange[0]
305 311 self.xlabel = "Frequency (kHz)"
306 312 elif self.xaxis == "time":
307 313 x = self.data.xrange[1]
308 314 self.xlabel = "Time (ms)"
309 315 else:
310 316 x = self.data.xrange[2]
311 317 self.xlabel = "Velocity (m/s)"
312 318
313 319 self.titles = []
314 320
315 321 y = self.data.yrange
316 322 self.y = y
317 323
318 324 data = self.data[-1]
319 325 cspc = data['cspc']
320 326
321 327 for n in range(len(self.data.pairs)):
322 328 pair = self.data.pairs[n]
323 329 coh = cspc[n*2]
324 330 phase = cspc[n*2+1]
325 331 ax = self.axes[2 * n]
326 332 if ax.firsttime:
327 333 ax.plt = ax.pcolormesh(x, y, coh.T,
328 334 vmin=0,
329 335 vmax=1,
330 336 cmap=plt.get_cmap(self.colormap_coh)
331 337 )
332 338 else:
333 339 ax.plt.set_array(coh.T.ravel())
334 340 self.titles.append(
335 341 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
336 342
337 343 ax = self.axes[2 * n + 1]
338 344 if ax.firsttime:
339 345 ax.plt = ax.pcolormesh(x, y, phase.T,
340 346 vmin=-180,
341 347 vmax=180,
342 348 cmap=plt.get_cmap(self.colormap_phase)
343 349 )
344 350 else:
345 351 ax.plt.set_array(phase.T.ravel())
346 352 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
347 353
348 354
349 355 class CrossSpectra4Plot(Plot):
350 356
351 357 CODE = 'cspc'
352 358 colormap = 'jet'
353 359 plot_type = 'pcolor'
354 360 zmin_coh = None
355 361 zmax_coh = None
356 362 zmin_phase = None
357 363 zmax_phase = None
358 364
359 365 def setup(self):
360 366
361 367 self.ncols = 4
362 368 self.nrows = len(self.data.pairs)
363 369 self.nplots = self.nrows * 4
364 370 self.width = 3.1 * self.ncols
365 371 self.height = 5 * self.nrows
366 372 self.ylabel = 'Range [km]'
367 373 self.showprofile = False
368 374 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
369 375
370 376 def plot(self):
371 377
372 378 if self.xaxis == "frequency":
373 379 x = self.data.xrange[0]
374 380 self.xlabel = "Frequency (kHz)"
375 381 elif self.xaxis == "time":
376 382 x = self.data.xrange[1]
377 383 self.xlabel = "Time (ms)"
378 384 else:
379 385 x = self.data.xrange[2]
380 386 self.xlabel = "Velocity (m/s)"
381 387
382 388 self.titles = []
383 389
384 390
385 391 y = self.data.heights
386 392 self.y = y
387 393 nspc = self.data['spc']
388 394 #print(numpy.shape(self.data['spc']))
389 395 spc = self.data['cspc'][0]
390 396 #print(numpy.shape(nspc))
391 397 #exit()
392 398 #nspc[1,:,:] = numpy.flip(nspc[1,:,:],axis=0)
393 399 #print(numpy.shape(spc))
394 400 #exit()
395 401 cspc = self.data['cspc'][1]
396 402
397 403 #xflip=numpy.flip(x)
398 404 #print(numpy.shape(cspc))
399 405 #exit()
400 406
401 407 for n in range(self.nrows):
402 408 noise = self.data['noise'][:,-1]
403 409 pair = self.data.pairs[n]
404 410 #print(pair)
405 411 #exit()
406 412 ax = self.axes[4 * n]
407 413 if ax.firsttime:
408 414 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
409 415 self.xmin = self.xmin if self.xmin else -self.xmax
410 416 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
411 417 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
412 418 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
413 419 vmin=self.zmin,
414 420 vmax=self.zmax,
415 421 cmap=plt.get_cmap(self.colormap)
416 422 )
417 423 else:
418 424 #print(numpy.shape(nspc[pair[0]].T))
419 425 #exit()
420 426 ax.plt.set_array(nspc[pair[0]].T.ravel())
421 427 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
422 428
423 429 ax = self.axes[4 * n + 1]
424 430
425 431 if ax.firsttime:
426 432 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
427 433 vmin=self.zmin,
428 434 vmax=self.zmax,
429 435 cmap=plt.get_cmap(self.colormap)
430 436 )
431 437 else:
432 438
433 439 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
434 440 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
435 441
436 442 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
437 443 coh = numpy.abs(out)
438 444 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
439 445
440 446 ax = self.axes[4 * n + 2]
441 447 if ax.firsttime:
442 448 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
443 449 vmin=0,
444 450 vmax=1,
445 451 cmap=plt.get_cmap(self.colormap_coh)
446 452 )
447 453 else:
448 454 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
449 455 self.titles.append(
450 456 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
451 457
452 458 ax = self.axes[4 * n + 3]
453 459 if ax.firsttime:
454 460 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
455 461 vmin=-180,
456 462 vmax=180,
457 463 cmap=plt.get_cmap(self.colormap_phase)
458 464 )
459 465 else:
460 466 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
461 467 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
462 468
463 469
464 470 class CrossSpectra2Plot(Plot):
465 471
466 472 CODE = 'cspc'
467 473 colormap = 'jet'
468 474 plot_type = 'pcolor'
469 475 zmin_coh = None
470 476 zmax_coh = None
471 477 zmin_phase = None
472 478 zmax_phase = None
473 479
474 480 def setup(self):
475 481
476 482 self.ncols = 1
477 483 self.nrows = len(self.data.pairs)
478 484 self.nplots = self.nrows * 1
479 485 self.width = 3.1 * self.ncols
480 486 self.height = 5 * self.nrows
481 487 self.ylabel = 'Range [km]'
482 488 self.showprofile = False
483 489 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
484 490
485 491 def plot(self):
486 492
487 493 if self.xaxis == "frequency":
488 494 x = self.data.xrange[0]
489 495 self.xlabel = "Frequency (kHz)"
490 496 elif self.xaxis == "time":
491 497 x = self.data.xrange[1]
492 498 self.xlabel = "Time (ms)"
493 499 else:
494 500 x = self.data.xrange[2]
495 501 self.xlabel = "Velocity (m/s)"
496 502
497 503 self.titles = []
498 504
499 505
500 506 y = self.data.heights
501 507 self.y = y
502 508 #nspc = self.data['spc']
503 509 #print(numpy.shape(self.data['spc']))
504 510 #spc = self.data['cspc'][0]
505 511 #print(numpy.shape(spc))
506 512 #exit()
507 513 cspc = self.data['cspc'][1]
508 514 #print(numpy.shape(cspc))
509 515 #exit()
510 516
511 517 for n in range(self.nrows):
512 518 noise = self.data['noise'][:,-1]
513 519 pair = self.data.pairs[n]
514 520 #print(pair) #exit()
515 521
516 522
517 523
518 524 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
519 525
520 526 #print(out[:,53])
521 527 #exit()
522 528 cross = numpy.abs(out)
523 529 z = cross/self.data.nFactor
524 530 #print("here")
525 531 #print(dataOut.data_spc[0,0,0])
526 532 #exit()
527 533
528 534 cross = 10*numpy.log10(z)
529 535 #print(numpy.shape(cross))
530 536 #print(cross[0,:])
531 537 #print(self.data.nFactor)
532 538 #exit()
533 539 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
534 540
535 541 ax = self.axes[1 * n]
536 542 if ax.firsttime:
537 543 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
538 544 self.xmin = self.xmin if self.xmin else -self.xmax
539 545 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
540 546 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
541 547 ax.plt = ax.pcolormesh(x, y, cross.T,
542 548 vmin=self.zmin,
543 549 vmax=self.zmax,
544 550 cmap=plt.get_cmap(self.colormap)
545 551 )
546 552 else:
547 553 ax.plt.set_array(cross.T.ravel())
548 554 self.titles.append(
549 555 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
550 556
551 557
552 558 class CrossSpectra3Plot(Plot):
553 559
554 560 CODE = 'cspc'
555 561 colormap = 'jet'
556 562 plot_type = 'pcolor'
557 563 zmin_coh = None
558 564 zmax_coh = None
559 565 zmin_phase = None
560 566 zmax_phase = None
561 567
562 568 def setup(self):
563 569
564 570 self.ncols = 3
565 571 self.nrows = len(self.data.pairs)
566 572 self.nplots = self.nrows * 3
567 573 self.width = 3.1 * self.ncols
568 574 self.height = 5 * self.nrows
569 575 self.ylabel = 'Range [km]'
570 576 self.showprofile = False
571 577 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
572 578
573 579 def plot(self):
574 580
575 581 if self.xaxis == "frequency":
576 582 x = self.data.xrange[0]
577 583 self.xlabel = "Frequency (kHz)"
578 584 elif self.xaxis == "time":
579 585 x = self.data.xrange[1]
580 586 self.xlabel = "Time (ms)"
581 587 else:
582 588 x = self.data.xrange[2]
583 589 self.xlabel = "Velocity (m/s)"
584 590
585 591 self.titles = []
586 592
587 593
588 594 y = self.data.heights
589 595 self.y = y
590 596 #nspc = self.data['spc']
591 597 #print(numpy.shape(self.data['spc']))
592 598 #spc = self.data['cspc'][0]
593 599 #print(numpy.shape(spc))
594 600 #exit()
595 601 cspc = self.data['cspc'][1]
596 602 #print(numpy.shape(cspc))
597 603 #exit()
598 604
599 605 for n in range(self.nrows):
600 606 noise = self.data['noise'][:,-1]
601 607 pair = self.data.pairs[n]
602 608 #print(pair) #exit()
603 609
604 610
605 611
606 612 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
607 613
608 614 #print(out[:,53])
609 615 #exit()
610 616 cross = numpy.abs(out)
611 617 z = cross/self.data.nFactor
612 618 cross = 10*numpy.log10(z)
613 619
614 620 out_r= out.real/self.data.nFactor
615 621 #out_r = 10*numpy.log10(out_r)
616 622
617 623 out_i= out.imag/self.data.nFactor
618 624 #out_i = 10*numpy.log10(out_i)
619 625 #print(numpy.shape(cross))
620 626 #print(cross[0,:])
621 627 #print(self.data.nFactor)
622 628 #exit()
623 629 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
624 630
625 631 ax = self.axes[3 * n]
626 632 if ax.firsttime:
627 633 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
628 634 self.xmin = self.xmin if self.xmin else -self.xmax
629 635 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
630 636 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
631 637 ax.plt = ax.pcolormesh(x, y, cross.T,
632 638 vmin=self.zmin,
633 639 vmax=self.zmax,
634 640 cmap=plt.get_cmap(self.colormap)
635 641 )
636 642 else:
637 643 ax.plt.set_array(cross.T.ravel())
638 644 self.titles.append(
639 645 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
640 646
641 647 ax = self.axes[3 * n + 1]
642 648 if ax.firsttime:
643 649 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
644 650 self.xmin = self.xmin if self.xmin else -self.xmax
645 651 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
646 652 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
647 653 ax.plt = ax.pcolormesh(x, y, out_r.T,
648 654 vmin=-1.e6,
649 655 vmax=0,
650 656 cmap=plt.get_cmap(self.colormap)
651 657 )
652 658 else:
653 659 ax.plt.set_array(out_r.T.ravel())
654 660 self.titles.append(
655 661 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
656 662
657 663 ax = self.axes[3 * n + 2]
658 664
659 665
660 666 if ax.firsttime:
661 667 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
662 668 self.xmin = self.xmin if self.xmin else -self.xmax
663 669 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
664 670 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
665 671 ax.plt = ax.pcolormesh(x, y, out_i.T,
666 672 vmin=-1.e6,
667 673 vmax=1.e6,
668 674 cmap=plt.get_cmap(self.colormap)
669 675 )
670 676 else:
671 677 ax.plt.set_array(out_i.T.ravel())
672 678 self.titles.append(
673 679 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
674 680
675 681 class RTIPlot(Plot):
676 682 '''
677 683 Plot for RTI data
678 684 '''
679 685
680 686 CODE = 'rti'
681 687 colormap = 'jet'
682 688 plot_type = 'pcolorbuffer'
683 689
684 690 def setup(self):
685 691 self.xaxis = 'time'
686 692 self.ncols = 1
687 693 self.nrows = len(self.data.channels)
688 694 self.nplots = len(self.data.channels)
689 695 self.ylabel = 'Range [km]'
690 696 self.xlabel = 'Time'
691 697 self.cb_label = 'dB'
692 698 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
693 699 self.titles = ['{} Channel {}'.format(
694 700 self.CODE.upper(), x) for x in range(self.nrows)]
695 701
696 702 def update(self, dataOut):
697 703
698 704 data = {}
699 705 meta = {}
700 706 data['rti'] = dataOut.getPower()
701 707 #print(numpy.shape(data['rti']))
702 708
703 709 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
704 710
705 711 return data, meta
706 712
707 713 def plot(self):
708 714
709 715 self.x = self.data.times
710 716 self.y = self.data.yrange
711 717 self.z = self.data[self.CODE]
712 718
713 719 self.z = numpy.ma.masked_invalid(self.z)
714 720
715 721 if self.decimation is None:
716 722 x, y, z = self.fill_gaps(self.x, self.y, self.z)
717 723 else:
718 724 x, y, z = self.fill_gaps(*self.decimate())
719 725
720 726 '''
721 727 if not isinstance(self.zmin, collections.abc.Sequence):
722 728 if not self.zmin:
723 729 self.zmin = [numpy.min(self.z)]*len(self.axes)
724 730 else:
725 731 self.zmin = [self.zmin]*len(self.axes)
726 732
727 733 if not isinstance(self.zmax, collections.abc.Sequence):
728 734 if not self.zmax:
729 735 self.zmax = [numpy.max(self.z)]*len(self.axes)
730 736 else:
731 737 self.zmax = [self.zmax]*len(self.axes)
732 738 '''
733 739 for n, ax in enumerate(self.axes):
734 740
735 741 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
736 742 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
737 743
738 744 if ax.firsttime:
739 745 if self.zlimits is not None:
740 746 self.zmin, self.zmax = self.zlimits[n]
741 747 ax.plt = ax.pcolormesh(x, y, z[n].T,
742 748 vmin=self.zmin,
743 749 vmax=self.zmax,
744 750 cmap=plt.get_cmap(self.colormap)
745 751 )
746 752 if self.showprofile:
747 753 ax.plot_profile = self.pf_axes[n].plot(
748 754 self.data['rti'][n][-1], self.y)[0]
749 755 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
750 756 color="k", linestyle="dashed", lw=1)[0]
751 757 else:
752 758 if self.zlimits is not None:
753 759 self.zmin, self.zmax = self.zlimits[n]
754 760 ax.collections.remove(ax.collections[0])
755 761 ax.plt = ax.pcolormesh(x, y, z[n].T,
756 762 vmin=self.zmin,
757 763 vmax=self.zmax,
758 764 cmap=plt.get_cmap(self.colormap)
759 765 )
760 766 if self.showprofile:
761 767 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
762 768 ax.plot_noise.set_data(numpy.repeat(
763 769 self.data['noise'][n][-1], len(self.y)), self.y)
764 770
765 771
766 772 class SpectrogramPlot(Plot):
767 773 '''
768 774 Plot for Spectrogram data
769 775 '''
770 776
771 777 CODE = 'Spectrogram_Profile'
772 778 colormap = 'binary'
773 779 plot_type = 'pcolorbuffer'
774 780
775 781 def setup(self):
776 782 self.xaxis = 'time'
777 783 self.ncols = 1
778 784 self.nrows = len(self.data.channels)
779 785 self.nplots = len(self.data.channels)
780 786 self.xlabel = 'Time'
781 787 #self.cb_label = 'dB'
782 788 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
783 789 self.titles = []
784 790
785 791 #self.titles = ['{} Channel {} \n H = {} km ({} - {})'.format(
786 792 #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)]
787 793
788 794 self.titles = ['{} Channel {}'.format(
789 795 self.CODE.upper(), x) for x in range(self.nrows)]
790 796
791 797
792 798 def update(self, dataOut):
793 799 data = {}
794 800 meta = {}
795 801
796 802 maxHei = 1620#+12000
797 803 indb = numpy.where(dataOut.heightList <= maxHei)
798 804 hei = indb[0][-1]
799 805 #print(dataOut.heightList)
800 806
801 807 factor = dataOut.nIncohInt
802 808 z = dataOut.data_spc[:,:,hei] / factor
803 809 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
804 810 #buffer = 10 * numpy.log10(z)
805 811
806 812 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
807 813
808 814
809 815 #self.hei = hei
810 816 #self.heightList = dataOut.heightList
811 817 #self.DH = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
812 818 #self.nProfiles = dataOut.nProfiles
813 819
814 820 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
815 821
816 822 data['hei'] = hei
817 823 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
818 824 data['nProfiles'] = dataOut.nProfiles
819 825 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
820 826 '''
821 827 import matplotlib.pyplot as plt
822 828 plt.plot(10 * numpy.log10(z[0,:]))
823 829 plt.show()
824 830
825 831 from time import sleep
826 832 sleep(10)
827 833 '''
828 834 return data, meta
829 835
830 836 def plot(self):
831 837
832 838 self.x = self.data.times
833 839 self.z = self.data[self.CODE]
834 840 self.y = self.data.xrange[0]
835 841
836 842 hei = self.data['hei'][-1]
837 843 DH = self.data['DH'][-1]
838 844 nProfiles = self.data['nProfiles'][-1]
839 845
840 846 self.ylabel = "Frequency (kHz)"
841 847
842 848 self.z = numpy.ma.masked_invalid(self.z)
843 849
844 850 if self.decimation is None:
845 851 x, y, z = self.fill_gaps(self.x, self.y, self.z)
846 852 else:
847 853 x, y, z = self.fill_gaps(*self.decimate())
848 854
849 855 for n, ax in enumerate(self.axes):
850 856 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
851 857 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
852 858 data = self.data[-1]
853 859 if ax.firsttime:
854 860 ax.plt = ax.pcolormesh(x, y, z[n].T,
855 861 vmin=self.zmin,
856 862 vmax=self.zmax,
857 863 cmap=plt.get_cmap(self.colormap)
858 864 )
859 865 else:
860 866 ax.collections.remove(ax.collections[0])
861 867 ax.plt = ax.pcolormesh(x, y, z[n].T,
862 868 vmin=self.zmin,
863 869 vmax=self.zmax,
864 870 cmap=plt.get_cmap(self.colormap)
865 871 )
866 872
867 873 #self.titles.append('Spectrogram')
868 874
869 875 #self.titles.append('{} Channel {} \n H = {} km ({} - {})'.format(
870 876 #self.CODE.upper(), x, y[hei], y[hei],y[hei]+(DH*nProfiles)))
871 877
872 878
873 879
874 880
875 881 class CoherencePlot(RTIPlot):
876 882 '''
877 883 Plot for Coherence data
878 884 '''
879 885
880 886 CODE = 'coh'
881 887
882 888 def setup(self):
883 889 self.xaxis = 'time'
884 890 self.ncols = 1
885 891 self.nrows = len(self.data.pairs)
886 892 self.nplots = len(self.data.pairs)
887 893 self.ylabel = 'Range [km]'
888 894 self.xlabel = 'Time'
889 895 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
890 896 if self.CODE == 'coh':
891 897 self.cb_label = ''
892 898 self.titles = [
893 899 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
894 900 else:
895 901 self.cb_label = 'Degrees'
896 902 self.titles = [
897 903 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
898 904
899 905 def update(self, dataOut):
900 906
901 907 data = {}
902 908 meta = {}
903 909 data['coh'] = dataOut.getCoherence()
904 910 meta['pairs'] = dataOut.pairsList
905 911
906 912 return data, meta
907 913
908 914 class PhasePlot(CoherencePlot):
909 915 '''
910 916 Plot for Phase map data
911 917 '''
912 918
913 919 CODE = 'phase'
914 920 colormap = 'seismic'
915 921
916 922 def update(self, dataOut):
917 923
918 924 data = {}
919 925 meta = {}
920 926 data['phase'] = dataOut.getCoherence(phase=True)
921 927 meta['pairs'] = dataOut.pairsList
922 928
923 929 return data, meta
924 930
925 931 class NoisePlot(Plot):
926 932 '''
927 933 Plot for noise
928 934 '''
929 935
930 936 CODE = 'noise'
931 937 plot_type = 'scatterbuffer'
932 938
933 939 def setup(self):
934 940 self.xaxis = 'time'
935 941 self.ncols = 1
936 942 self.nrows = 1
937 943 self.nplots = 1
938 944 self.ylabel = 'Intensity [dB]'
939 945 self.xlabel = 'Time'
940 946 self.titles = ['Noise']
941 947 self.colorbar = False
942 948 self.plots_adjust.update({'right': 0.85 })
943 949
944 950 def update(self, dataOut):
945 951
946 952 data = {}
947 953 meta = {}
948 954 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
949 955 meta['yrange'] = numpy.array([])
950 956
951 957 return data, meta
952 958
953 959 def plot(self):
954 960
955 961 x = self.data.times
956 962 xmin = self.data.min_time
957 963 xmax = xmin + self.xrange * 60 * 60
958 964 Y = self.data['noise']
959 965
960 966 if self.axes[0].firsttime:
961 967 self.ymin = numpy.nanmin(Y) - 5
962 968 self.ymax = numpy.nanmax(Y) + 5
963 969 for ch in self.data.channels:
964 970 y = Y[ch]
965 971 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
966 972 plt.legend(bbox_to_anchor=(1.18, 1.0))
967 973 else:
968 974 for ch in self.data.channels:
969 975 y = Y[ch]
970 976 self.axes[0].lines[ch].set_data(x, y)
971 977
972 978 self.ymin = numpy.nanmin(Y) - 5
973 979 self.ymax = numpy.nanmax(Y) + 10
974 980
975 981
976 982 class PowerProfilePlot(Plot):
977 983
978 984 CODE = 'pow_profile'
979 985 plot_type = 'scatter'
980 986
981 987 def setup(self):
982 988
983 989 self.ncols = 1
984 990 self.nrows = 1
985 991 self.nplots = 1
986 992 self.height = 4
987 993 self.width = 3
988 994 self.ylabel = 'Range [km]'
989 995 self.xlabel = 'Intensity [dB]'
990 996 self.titles = ['Power Profile']
991 997 self.colorbar = False
992 998
993 999 def update(self, dataOut):
994 1000
995 1001 data = {}
996 1002 meta = {}
997 1003 data[self.CODE] = dataOut.getPower()
998 1004
999 1005 return data, meta
1000 1006
1001 1007 def plot(self):
1002 1008
1003 1009 y = self.data.yrange
1004 1010 self.y = y
1005 1011
1006 1012 x = self.data[-1][self.CODE]
1007 1013
1008 1014 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
1009 1015 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
1010 1016
1011 1017 if self.axes[0].firsttime:
1012 1018 for ch in self.data.channels:
1013 1019 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
1014 1020 plt.legend()
1015 1021 else:
1016 1022 for ch in self.data.channels:
1017 1023 self.axes[0].lines[ch].set_data(x[ch], y)
1018 1024
1019 1025
1020 1026 class SpectraCutPlot(Plot):
1021 1027
1022 1028 CODE = 'spc_cut'
1023 1029 plot_type = 'scatter'
1024 1030 buffering = False
1025 1031
1026 1032 def setup(self):
1027 1033
1028 1034 self.nplots = len(self.data.channels)
1029 1035 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
1030 1036 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
1031 1037 self.width = 3.4 * self.ncols + 1.5
1032 1038 self.height = 3 * self.nrows
1033 1039 self.ylabel = 'Power [dB]'
1034 1040 self.colorbar = False
1035 1041 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
1036 1042
1037 1043 def update(self, dataOut):
1038 1044
1039 1045 data = {}
1040 1046 meta = {}
1041 1047 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
1042 1048 data['spc'] = spc
1043 1049 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
1044 1050 if self.CODE == 'cut_gaussian_fit':
1045 1051 data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
1046 1052 data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
1047 1053 return data, meta
1048 1054
1049 1055 def plot(self):
1050 1056 if self.xaxis == "frequency":
1051 1057 x = self.data.xrange[0][1:]
1052 1058 self.xlabel = "Frequency (kHz)"
1053 1059 elif self.xaxis == "time":
1054 1060 x = self.data.xrange[1]
1055 1061 self.xlabel = "Time (ms)"
1056 1062 else:
1057 1063 x = self.data.xrange[2][:-1]
1058 1064 self.xlabel = "Velocity (m/s)"
1059 1065
1060 1066 if self.CODE == 'cut_gaussian_fit':
1061 1067 x = self.data.xrange[2][:-1]
1062 1068 self.xlabel = "Velocity (m/s)"
1063 1069
1064 1070 self.titles = []
1065 1071
1066 1072 y = self.data.yrange
1067 1073 data = self.data[-1]
1068 1074 z = data['spc']
1069 1075
1070 1076 if self.height_index:
1071 1077 index = numpy.array(self.height_index)
1072 1078 else:
1073 1079 index = numpy.arange(0, len(y), int((len(y))/9))
1074 1080
1075 1081 for n, ax in enumerate(self.axes):
1076 1082 if self.CODE == 'cut_gaussian_fit':
1077 1083 gau0 = data['gauss_fit0']
1078 1084 gau1 = data['gauss_fit1']
1079 1085 if ax.firsttime:
1080 1086 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1081 1087 self.xmin = self.xmin if self.xmin else -self.xmax
1082 1088 self.ymin = self.ymin if self.ymin else numpy.nanmin(z[:,:,index])
1083 1089 self.ymax = self.ymax if self.ymax else numpy.nanmax(z[:,:,index])
1084 1090 #print(self.ymax)
1085 1091 #print(z[n, :, index])
1086 1092 ax.plt = ax.plot(x, z[n, :, index].T, lw=0.25)
1087 1093 if self.CODE == 'cut_gaussian_fit':
1088 1094 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
1089 1095 for i, line in enumerate(ax.plt_gau0):
1090 1096 line.set_color(ax.plt[i].get_color())
1091 1097 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
1092 1098 for i, line in enumerate(ax.plt_gau1):
1093 1099 line.set_color(ax.plt[i].get_color())
1094 1100 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1095 1101 self.figures[0].legend(ax.plt, labels, loc='center right')
1096 1102 else:
1097 1103 for i, line in enumerate(ax.plt):
1098 1104 line.set_data(x, z[n, :, index[i]].T)
1099 1105 for i, line in enumerate(ax.plt_gau0):
1100 1106 line.set_data(x, gau0[n, :, index[i]].T)
1101 1107 line.set_color(ax.plt[i].get_color())
1102 1108 for i, line in enumerate(ax.plt_gau1):
1103 1109 line.set_data(x, gau1[n, :, index[i]].T)
1104 1110 line.set_color(ax.plt[i].get_color())
1105 1111 self.titles.append('CH {}'.format(n))
1106 1112
1107 1113
1108 1114 class BeaconPhase(Plot):
1109 1115
1110 1116 __isConfig = None
1111 1117 __nsubplots = None
1112 1118
1113 1119 PREFIX = 'beacon_phase'
1114 1120
1115 1121 def __init__(self):
1116 1122 Plot.__init__(self)
1117 1123 self.timerange = 24*60*60
1118 1124 self.isConfig = False
1119 1125 self.__nsubplots = 1
1120 1126 self.counter_imagwr = 0
1121 1127 self.WIDTH = 800
1122 1128 self.HEIGHT = 400
1123 1129 self.WIDTHPROF = 120
1124 1130 self.HEIGHTPROF = 0
1125 1131 self.xdata = None
1126 1132 self.ydata = None
1127 1133
1128 1134 self.PLOT_CODE = BEACON_CODE
1129 1135
1130 1136 self.FTP_WEI = None
1131 1137 self.EXP_CODE = None
1132 1138 self.SUB_EXP_CODE = None
1133 1139 self.PLOT_POS = None
1134 1140
1135 1141 self.filename_phase = None
1136 1142
1137 1143 self.figfile = None
1138 1144
1139 1145 self.xmin = None
1140 1146 self.xmax = None
1141 1147
1142 1148 def getSubplots(self):
1143 1149
1144 1150 ncol = 1
1145 1151 nrow = 1
1146 1152
1147 1153 return nrow, ncol
1148 1154
1149 1155 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1150 1156
1151 1157 self.__showprofile = showprofile
1152 1158 self.nplots = nplots
1153 1159
1154 1160 ncolspan = 7
1155 1161 colspan = 6
1156 1162 self.__nsubplots = 2
1157 1163
1158 1164 self.createFigure(id = id,
1159 1165 wintitle = wintitle,
1160 1166 widthplot = self.WIDTH+self.WIDTHPROF,
1161 1167 heightplot = self.HEIGHT+self.HEIGHTPROF,
1162 1168 show=show)
1163 1169
1164 1170 nrow, ncol = self.getSubplots()
1165 1171
1166 1172 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1167 1173
1168 1174 def save_phase(self, filename_phase):
1169 1175 f = open(filename_phase,'w+')
1170 1176 f.write('\n\n')
1171 1177 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1172 1178 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1173 1179 f.close()
1174 1180
1175 1181 def save_data(self, filename_phase, data, data_datetime):
1176 1182 f=open(filename_phase,'a')
1177 1183 timetuple_data = data_datetime.timetuple()
1178 1184 day = str(timetuple_data.tm_mday)
1179 1185 month = str(timetuple_data.tm_mon)
1180 1186 year = str(timetuple_data.tm_year)
1181 1187 hour = str(timetuple_data.tm_hour)
1182 1188 minute = str(timetuple_data.tm_min)
1183 1189 second = str(timetuple_data.tm_sec)
1184 1190 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1185 1191 f.close()
1186 1192
1187 1193 def plot(self):
1188 1194 log.warning('TODO: Not yet implemented...')
1189 1195
1190 1196 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1191 1197 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1192 1198 timerange=None,
1193 1199 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1194 1200 server=None, folder=None, username=None, password=None,
1195 1201 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1196 1202
1197 1203 if dataOut.flagNoData:
1198 1204 return dataOut
1199 1205
1200 1206 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1201 1207 return
1202 1208
1203 1209 if pairsList == None:
1204 1210 pairsIndexList = dataOut.pairsIndexList[:10]
1205 1211 else:
1206 1212 pairsIndexList = []
1207 1213 for pair in pairsList:
1208 1214 if pair not in dataOut.pairsList:
1209 1215 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1210 1216 pairsIndexList.append(dataOut.pairsList.index(pair))
1211 1217
1212 1218 if pairsIndexList == []:
1213 1219 return
1214 1220
1215 1221 # if len(pairsIndexList) > 4:
1216 1222 # pairsIndexList = pairsIndexList[0:4]
1217 1223
1218 1224 hmin_index = None
1219 1225 hmax_index = None
1220 1226
1221 1227 if hmin != None and hmax != None:
1222 1228 indexes = numpy.arange(dataOut.nHeights)
1223 1229 hmin_list = indexes[dataOut.heightList >= hmin]
1224 1230 hmax_list = indexes[dataOut.heightList <= hmax]
1225 1231
1226 1232 if hmin_list.any():
1227 1233 hmin_index = hmin_list[0]
1228 1234
1229 1235 if hmax_list.any():
1230 1236 hmax_index = hmax_list[-1]+1
1231 1237
1232 1238 x = dataOut.getTimeRange()
1233 1239
1234 1240 thisDatetime = dataOut.datatime
1235 1241
1236 1242 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1237 1243 xlabel = "Local Time"
1238 1244 ylabel = "Phase (degrees)"
1239 1245
1240 1246 update_figfile = False
1241 1247
1242 1248 nplots = len(pairsIndexList)
1243 1249 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1244 1250 phase_beacon = numpy.zeros(len(pairsIndexList))
1245 1251 for i in range(nplots):
1246 1252 pair = dataOut.pairsList[pairsIndexList[i]]
1247 1253 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1248 1254 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1249 1255 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1250 1256 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1251 1257 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1252 1258
1253 1259 if dataOut.beacon_heiIndexList:
1254 1260 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1255 1261 else:
1256 1262 phase_beacon[i] = numpy.average(phase)
1257 1263
1258 1264 if not self.isConfig:
1259 1265
1260 1266 nplots = len(pairsIndexList)
1261 1267
1262 1268 self.setup(id=id,
1263 1269 nplots=nplots,
1264 1270 wintitle=wintitle,
1265 1271 showprofile=showprofile,
1266 1272 show=show)
1267 1273
1268 1274 if timerange != None:
1269 1275 self.timerange = timerange
1270 1276
1271 1277 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1272 1278
1273 1279 if ymin == None: ymin = 0
1274 1280 if ymax == None: ymax = 360
1275 1281
1276 1282 self.FTP_WEI = ftp_wei
1277 1283 self.EXP_CODE = exp_code
1278 1284 self.SUB_EXP_CODE = sub_exp_code
1279 1285 self.PLOT_POS = plot_pos
1280 1286
1281 1287 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1282 1288 self.isConfig = True
1283 1289 self.figfile = figfile
1284 1290 self.xdata = numpy.array([])
1285 1291 self.ydata = numpy.array([])
1286 1292
1287 1293 update_figfile = True
1288 1294
1289 1295 #open file beacon phase
1290 1296 path = '%s%03d' %(self.PREFIX, self.id)
1291 1297 beacon_file = os.path.join(path,'%s.txt'%self.name)
1292 1298 self.filename_phase = os.path.join(figpath,beacon_file)
1293 1299 #self.save_phase(self.filename_phase)
1294 1300
1295 1301
1296 1302 #store data beacon phase
1297 1303 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1298 1304
1299 1305 self.setWinTitle(title)
1300 1306
1301 1307
1302 1308 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1303 1309
1304 1310 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1305 1311
1306 1312 axes = self.axesList[0]
1307 1313
1308 1314 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1309 1315
1310 1316 if len(self.ydata)==0:
1311 1317 self.ydata = phase_beacon.reshape(-1,1)
1312 1318 else:
1313 1319 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1314 1320
1315 1321
1316 1322 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1317 1323 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1318 1324 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1319 1325 XAxisAsTime=True, grid='both'
1320 1326 )
1321 1327
1322 1328 self.draw()
1323 1329
1324 1330 if dataOut.ltctime >= self.xmax:
1325 1331 self.counter_imagwr = wr_period
1326 1332 self.isConfig = False
1327 1333 update_figfile = True
1328 1334
1329 1335 self.save(figpath=figpath,
1330 1336 figfile=figfile,
1331 1337 save=save,
1332 1338 ftp=ftp,
1333 1339 wr_period=wr_period,
1334 1340 thisDatetime=thisDatetime,
1335 1341 update_figfile=update_figfile)
1336 1342
1337 1343 return dataOut
Show file before | Show file after | Hide comments
@@ -1,1287 +1,1287
1 1
2 2 import os
3 3 import time
4 4 import math
5 5 import datetime
6 6 import numpy
7 import collections.abc
7
8 8 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator #YONG
9 9
10 10 from .jroplot_spectra import RTIPlot, NoisePlot
11 11
12 12 from schainpy.utils import log
13 13 from .plotting_codes import *
14 14
15 15 from schainpy.model.graphics.jroplot_base import Plot, plt
16 16
17 17 import matplotlib.pyplot as plt
18 18 import matplotlib.colors as colors
19 19 from matplotlib.ticker import MultipleLocator
20 20
21 21
22 22 class RTIDPPlot(RTIPlot):
23 23
24 24 '''Plot for RTI Double Pulse Experiment
25 25 '''
26 26
27 27 CODE = 'RTIDP'
28 28 colormap = 'jet'
29 29 plot_name = 'RTI'
30 30 plot_type = 'pcolorbuffer'
31 31
32 32 def setup(self):
33 33 self.xaxis = 'time'
34 34 self.ncols = 1
35 35 self.nrows = 3
36 36 self.nplots = self.nrows
37 37
38 38 self.ylabel = 'Range [km]'
39 39 self.xlabel = 'Time (LT)'
40 40
41 41 self.cb_label = 'Intensity (dB)'
42 42
43 43 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
44 44
45 45 self.titles = ['{} Channel {}'.format(
46 46 self.plot_name.upper(), '0x1'),'{} Channel {}'.format(
47 47 self.plot_name.upper(), '0'),'{} Channel {}'.format(
48 48 self.plot_name.upper(), '1')]
49 49
50 50 def update(self, dataOut):
51 51
52 52 data = {}
53 53 meta = {}
54 54 data['rti'] = dataOut.data_for_RTI_DP
55 55 data['NDP'] = dataOut.NDP
56 56
57 57 return data, meta
58 58
59 59 def plot(self):
60 60
61 61 NDP = self.data['NDP'][-1]
62 62 self.x = self.data.times
63 63 self.y = self.data.yrange[0:NDP]
64 64 self.z = self.data['rti']
65 65 self.z = numpy.ma.masked_invalid(self.z)
66 66
67 67 if self.decimation is None:
68 68 x, y, z = self.fill_gaps(self.x, self.y, self.z)
69 69 else:
70 70 x, y, z = self.fill_gaps(*self.decimate())
71 71
72 72 for n, ax in enumerate(self.axes):
73 73
74 74 self.zmax = self.zmax if self.zmax is not None else numpy.max(
75 75 self.z[1][0,12:40])
76 76 self.zmin = self.zmin if self.zmin is not None else numpy.min(
77 77 self.z[1][0,12:40])
78 78
79 79 if ax.firsttime:
80 80
81 81 if self.zlimits is not None:
82 82 self.zmin, self.zmax = self.zlimits[n]
83 83
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 else:
90 90 #if self.zlimits is not None:
91 91 #self.zmin, self.zmax = self.zlimits[n]
92 92 ax.collections.remove(ax.collections[0])
93 93 ax.plt = ax.pcolormesh(x, y, z[n].T,
94 94 vmin=self.zmin,
95 95 vmax=self.zmax,
96 96 cmap=plt.get_cmap(self.colormap)
97 97 )
98 98
99 99
100 100 class RTILPPlot(RTIPlot):
101 101
102 102 '''
103 103 Plot for RTI Long Pulse
104 104 '''
105 105
106 106 CODE = 'RTILP'
107 107 colormap = 'jet'
108 108 plot_name = 'RTI LP'
109 109 plot_type = 'pcolorbuffer'
110 110
111 111 def setup(self):
112 112 self.xaxis = 'time'
113 113 self.ncols = 1
114 114 self.nrows = 2
115 115 self.nplots = self.nrows
116 116
117 117 self.ylabel = 'Range [km]'
118 118 self.xlabel = 'Time (LT)'
119 119
120 120 self.cb_label = 'Intensity (dB)'
121 121
122 122 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
123 123
124 124
125 125 self.titles = ['{} Channel {}'.format(
126 126 self.plot_name.upper(), '0'),'{} Channel {}'.format(
127 127 self.plot_name.upper(), '1'),'{} Channel {}'.format(
128 128 self.plot_name.upper(), '2'),'{} Channel {}'.format(
129 129 self.plot_name.upper(), '3')]
130 130
131 131
132 132 def update(self, dataOut):
133 133
134 134 data = {}
135 135 meta = {}
136 136 data['rti'] = dataOut.data_for_RTI_LP
137 137 data['NRANGE'] = dataOut.NRANGE
138 138
139 139 return data, meta
140 140
141 141 def plot(self):
142 142
143 143 NRANGE = self.data['NRANGE'][-1]
144 144 self.x = self.data.times
145 145 self.y = self.data.yrange[0:NRANGE]
146 146
147 147 self.z = self.data['rti']
148 148
149 149 self.z = numpy.ma.masked_invalid(self.z)
150 150
151 151 if self.decimation is None:
152 152 x, y, z = self.fill_gaps(self.x, self.y, self.z)
153 153 else:
154 154 x, y, z = self.fill_gaps(*self.decimate())
155 155
156 156 for n, ax in enumerate(self.axes):
157 157
158 158 self.zmax = self.zmax if self.zmax is not None else numpy.max(
159 159 self.z[1][0,12:40])
160 160 self.zmin = self.zmin if self.zmin is not None else numpy.min(
161 161 self.z[1][0,12:40])
162 162
163 163 if ax.firsttime:
164 164
165 165 if self.zlimits is not None:
166 166 self.zmin, self.zmax = self.zlimits[n]
167 167
168 168
169 169 ax.plt = ax.pcolormesh(x, y, z[n].T,
170 170 vmin=self.zmin,
171 171 vmax=self.zmax,
172 172 cmap=plt.get_cmap(self.colormap)
173 173 )
174 174
175 175 else:
176 176 if self.zlimits is not None:
177 177 self.zmin, self.zmax = self.zlimits[n]
178 178 ax.collections.remove(ax.collections[0])
179 179 ax.plt = ax.pcolormesh(x, y, z[n].T,
180 180 vmin=self.zmin,
181 181 vmax=self.zmax,
182 182 cmap=plt.get_cmap(self.colormap)
183 183 )
184 184
185 185
186 186 class DenRTIPlot(RTIPlot):
187 187
188 188 '''
189 189 Plot for Den
190 190 '''
191 191
192 192 CODE = 'denrti'
193 193 colormap = 'jet'
194 194
195 195 def setup(self):
196 196 self.xaxis = 'time'
197 197 self.ncols = 1
198 198 self.nrows = self.data.shape(self.CODE)[0]
199 199 self.nplots = self.nrows
200 200
201 201 self.ylabel = 'Range [km]'
202 202 self.xlabel = 'Time (LT)'
203 203
204 204 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
205 205
206 206 if self.CODE == 'denrti':
207 207 self.cb_label = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
208 208
209 209
210 210 self.titles = ['Electron Density RTI']
211 211
212 212 def update(self, dataOut):
213 213
214 214 data = {}
215 215 meta = {}
216 216
217 217 data['denrti'] = dataOut.DensityFinal*1.e-6 #To Plot in cm^-3
218 218
219 219 return data, meta
220 220
221 221 def plot(self):
222 222
223 223 self.x = self.data.times
224 224 self.y = self.data.yrange
225 225
226 226 self.z = self.data[self.CODE]
227 227
228 228 self.z = numpy.ma.masked_invalid(self.z)
229 229
230 230 if self.decimation is None:
231 231 x, y, z = self.fill_gaps(self.x, self.y, self.z)
232 232 else:
233 233 x, y, z = self.fill_gaps(*self.decimate())
234 234
235 235 for n, ax in enumerate(self.axes):
236 236
237 237 self.zmax = self.zmax if self.zmax is not None else numpy.max(
238 238 self.z[n])
239 239 self.zmin = self.zmin if self.zmin is not None else numpy.min(
240 240 self.z[n])
241 241
242 242 if ax.firsttime:
243 243
244 244 if self.zlimits is not None:
245 245 self.zmin, self.zmax = self.zlimits[n]
246 246 if numpy.log10(self.zmin)<0:
247 247 self.zmin=1
248 248 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
249 249 vmin=self.zmin,
250 250 vmax=self.zmax,
251 251 cmap=self.cmaps[n],
252 252 norm=colors.LogNorm()
253 253 )
254 254
255 255 else:
256 256 if self.zlimits is not None:
257 257 self.zmin, self.zmax = self.zlimits[n]
258 258 ax.collections.remove(ax.collections[0])
259 259 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
260 260 vmin=self.zmin,
261 261 vmax=self.zmax,
262 262 cmap=self.cmaps[n],
263 263 norm=colors.LogNorm()
264 264 )
265 265
266 266
267 267 class ETempRTIPlot(RTIPlot):
268 268
269 269 '''
270 270 Plot for Electron Temperature
271 271 '''
272 272
273 273 CODE = 'ETemp'
274 274 colormap = 'jet'
275 275
276 276 def setup(self):
277 277 self.xaxis = 'time'
278 278 self.ncols = 1
279 279 self.nrows = self.data.shape(self.CODE)[0]
280 280 self.nplots = self.nrows
281 281
282 282 self.ylabel = 'Range [km]'
283 283 self.xlabel = 'Time (LT)'
284 284 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
285 285 if self.CODE == 'ETemp':
286 286 self.cb_label = 'Electron Temperature (K)'
287 287 self.titles = ['Electron Temperature RTI']
288 288 if self.CODE == 'ITemp':
289 289 self.cb_label = 'Ion Temperature (K)'
290 290 self.titles = ['Ion Temperature RTI']
291 291 if self.CODE == 'HeFracLP':
292 292 self.cb_label='He+ Fraction'
293 293 self.titles = ['He+ Fraction RTI']
294 294 self.zmax=0.16
295 295 if self.CODE== 'HFracLP':
296 296 self.cb_label='H+ Fraction'
297 297 self.titles = ['H+ Fraction RTI']
298 298
299 299 def update(self, dataOut):
300 300
301 301 data = {}
302 302 meta = {}
303 303
304 304 data['ETemp'] = dataOut.ElecTempFinal
305 305
306 306 return data, meta
307 307
308 308 def plot(self):
309 309
310 310 self.x = self.data.times
311 311 self.y = self.data.yrange
312 312
313 313
314 314 self.z = self.data[self.CODE]
315 315
316 316 self.z = numpy.ma.masked_invalid(self.z)
317 317
318 318 if self.decimation is None:
319 319 x, y, z = self.fill_gaps(self.x, self.y, self.z)
320 320 else:
321 321 x, y, z = self.fill_gaps(*self.decimate())
322 322
323 323 for n, ax in enumerate(self.axes):
324 324
325 325 self.zmax = self.zmax if self.zmax is not None else numpy.max(
326 326 self.z[n])
327 327 self.zmin = self.zmin if self.zmin is not None else numpy.min(
328 328 self.z[n])
329 329
330 330 if ax.firsttime:
331 331
332 332 if self.zlimits is not None:
333 333 self.zmin, self.zmax = self.zlimits[n]
334 334
335 335 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
336 336 vmin=self.zmin,
337 337 vmax=self.zmax,
338 338 cmap=self.cmaps[n]
339 339 )
340 340 #plt.tight_layout()
341 341
342 342 else:
343 343 if self.zlimits is not None:
344 344 self.zmin, self.zmax = self.zlimits[n]
345 345 ax.collections.remove(ax.collections[0])
346 346 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
347 347 vmin=self.zmin,
348 348 vmax=self.zmax,
349 349 cmap=self.cmaps[n]
350 350 )
351 351
352 352
353 353 class ITempRTIPlot(ETempRTIPlot):
354 354
355 355 '''
356 356 Plot for Ion Temperature
357 357 '''
358 358
359 359 CODE = 'ITemp'
360 360 colormap = 'jet'
361 361 plot_name = 'Ion Temperature'
362 362
363 363 def update(self, dataOut):
364 364
365 365 data = {}
366 366 meta = {}
367 367
368 368 data['ITemp'] = dataOut.IonTempFinal
369 369
370 370 return data, meta
371 371
372 372
373 373 class HFracRTIPlot(ETempRTIPlot):
374 374
375 375 '''
376 376 Plot for H+ LP
377 377 '''
378 378
379 379 CODE = 'HFracLP'
380 380 colormap = 'jet'
381 381 plot_name = 'H+ Frac'
382 382
383 383 def update(self, dataOut):
384 384
385 385 data = {}
386 386 meta = {}
387 387 data['HFracLP'] = dataOut.PhyFinal
388 388
389 389 return data, meta
390 390
391 391
392 392 class HeFracRTIPlot(ETempRTIPlot):
393 393
394 394 '''
395 395 Plot for He+ LP
396 396 '''
397 397
398 398 CODE = 'HeFracLP'
399 399 colormap = 'jet'
400 400 plot_name = 'He+ Frac'
401 401
402 402 def update(self, dataOut):
403 403
404 404 data = {}
405 405 meta = {}
406 406 data['HeFracLP'] = dataOut.PheFinal
407 407
408 408 return data, meta
409 409
410 410
411 411 class TempsDPPlot(Plot):
412 412 '''
413 413 Plot for Electron - Ion Temperatures
414 414 '''
415 415
416 416 CODE = 'tempsDP'
417 417 #plot_name = 'Temperatures'
418 418 plot_type = 'scatterbuffer'
419 419
420 420 def setup(self):
421 421
422 422 self.ncols = 1
423 423 self.nrows = 1
424 424 self.nplots = 1
425 425 self.ylabel = 'Range [km]'
426 426 self.xlabel = 'Temperature (K)'
427 427 self.titles = ['Electron/Ion Temperatures']
428 428 self.width = 3.5
429 429 self.height = 5.5
430 430 self.colorbar = False
431 431 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
432 432
433 433 def update(self, dataOut):
434 434 data = {}
435 435 meta = {}
436 436
437 437 data['Te'] = dataOut.te2
438 438 data['Ti'] = dataOut.ti2
439 439 data['Te_error'] = dataOut.ete2
440 440 data['Ti_error'] = dataOut.eti2
441 441
442 442 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
443 443
444 444 return data, meta
445 445
446 446 def plot(self):
447 447
448 448 y = self.data.yrange
449 449
450 450 self.xmin = -100
451 451 self.xmax = 5000
452 452
453 453 ax = self.axes[0]
454 454
455 455 data = self.data[-1]
456 456
457 457 Te = data['Te']
458 458 Ti = data['Ti']
459 459 errTe = data['Te_error']
460 460 errTi = data['Ti_error']
461 461
462 462 if ax.firsttime:
463 463 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
464 464 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
465 465 plt.legend(loc='lower right')
466 466 self.ystep_given = 50
467 467 ax.yaxis.set_minor_locator(MultipleLocator(15))
468 468 ax.grid(which='minor')
469 469
470 470 else:
471 471 self.clear_figures()
472 472 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
473 473 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
474 474 plt.legend(loc='lower right')
475 475 ax.yaxis.set_minor_locator(MultipleLocator(15))
476 476
477 477
478 478 class TempsHPPlot(Plot):
479 479 '''
480 480 Plot for Temperatures Hybrid Experiment
481 481 '''
482 482
483 483 CODE = 'temps_LP'
484 484 #plot_name = 'Temperatures'
485 485 plot_type = 'scatterbuffer'
486 486
487 487
488 488 def setup(self):
489 489
490 490 self.ncols = 1
491 491 self.nrows = 1
492 492 self.nplots = 1
493 493 self.ylabel = 'Range [km]'
494 494 self.xlabel = 'Temperature (K)'
495 495 self.titles = ['Electron/Ion Temperatures']
496 496 self.width = 3.5
497 497 self.height = 6.5
498 498 self.colorbar = False
499 499 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
500 500
501 501 def update(self, dataOut):
502 502 data = {}
503 503 meta = {}
504 504
505 505
506 506 data['Te'] = numpy.concatenate((dataOut.te2[:dataOut.cut],dataOut.te[dataOut.cut:]))
507 507 data['Ti'] = numpy.concatenate((dataOut.ti2[:dataOut.cut],dataOut.ti[dataOut.cut:]))
508 508 data['Te_error'] = numpy.concatenate((dataOut.ete2[:dataOut.cut],dataOut.ete[dataOut.cut:]))
509 509 data['Ti_error'] = numpy.concatenate((dataOut.eti2[:dataOut.cut],dataOut.eti[dataOut.cut:]))
510 510
511 511 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
512 512
513 513 return data, meta
514 514
515 515 def plot(self):
516 516
517 517
518 518 self.y = self.data.yrange
519 519 self.xmin = -100
520 520 self.xmax = 4500
521 521 ax = self.axes[0]
522 522
523 523 data = self.data[-1]
524 524
525 525 Te = data['Te']
526 526 Ti = data['Ti']
527 527 errTe = data['Te_error']
528 528 errTi = data['Ti_error']
529 529
530 530 if ax.firsttime:
531 531
532 532 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
533 533 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
534 534 plt.legend(loc='lower right')
535 535 self.ystep_given = 200
536 536 ax.yaxis.set_minor_locator(MultipleLocator(15))
537 537 ax.grid(which='minor')
538 538
539 539 else:
540 540 self.clear_figures()
541 541 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
542 542 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
543 543 plt.legend(loc='lower right')
544 544 ax.yaxis.set_minor_locator(MultipleLocator(15))
545 545 ax.grid(which='minor')
546 546
547 547
548 548 class FracsHPPlot(Plot):
549 549 '''
550 550 Plot for Composition LP
551 551 '''
552 552
553 553 CODE = 'fracs_LP'
554 554 plot_type = 'scatterbuffer'
555 555
556 556
557 557 def setup(self):
558 558
559 559 self.ncols = 1
560 560 self.nrows = 1
561 561 self.nplots = 1
562 562 self.ylabel = 'Range [km]'
563 563 self.xlabel = 'Frac'
564 564 self.titles = ['Composition']
565 565 self.width = 3.5
566 566 self.height = 6.5
567 567 self.colorbar = False
568 568 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
569 569
570 570 def update(self, dataOut):
571 571 data = {}
572 572 meta = {}
573 573
574 574 #aux_nan=numpy.zeros(dataOut.cut,'float32')
575 575 #aux_nan[:]=numpy.nan
576 576 #data['ph'] = numpy.concatenate((aux_nan,dataOut.ph[dataOut.cut:]))
577 577 #data['eph'] = numpy.concatenate((aux_nan,dataOut.eph[dataOut.cut:]))
578 578
579 579 data['ph'] = dataOut.ph[dataOut.cut:]
580 580 data['eph'] = dataOut.eph[dataOut.cut:]
581 581 data['phe'] = dataOut.phe[dataOut.cut:]
582 582 data['ephe'] = dataOut.ephe[dataOut.cut:]
583 583
584 584 data['cut'] = dataOut.cut
585 585
586 586 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
587 587
588 588
589 589 return data, meta
590 590
591 591 def plot(self):
592 592
593 593 data = self.data[-1]
594 594
595 595 ph = data['ph']
596 596 eph = data['eph']
597 597 phe = data['phe']
598 598 ephe = data['ephe']
599 599 cut = data['cut']
600 600 self.y = self.data.yrange
601 601
602 602 self.xmin = 0
603 603 self.xmax = 1
604 604 ax = self.axes[0]
605 605
606 606 if ax.firsttime:
607 607
608 608 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
609 609 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
610 610 plt.legend(loc='lower right')
611 611 self.xstep_given = 0.2
612 612 self.ystep_given = 200
613 613 ax.yaxis.set_minor_locator(MultipleLocator(15))
614 614 ax.grid(which='minor')
615 615
616 616 else:
617 617 self.clear_figures()
618 618 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
619 619 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
620 620 plt.legend(loc='lower right')
621 621 ax.yaxis.set_minor_locator(MultipleLocator(15))
622 622 ax.grid(which='minor')
623 623
624 624 class EDensityPlot(Plot):
625 625 '''
626 626 Plot for electron density
627 627 '''
628 628
629 629 CODE = 'den'
630 630 #plot_name = 'Electron Density'
631 631 plot_type = 'scatterbuffer'
632 632
633 633 def setup(self):
634 634
635 635 self.ncols = 1
636 636 self.nrows = 1
637 637 self.nplots = 1
638 638 self.ylabel = 'Range [km]'
639 639 self.xlabel = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
640 640 self.titles = ['Electron Density']
641 641 self.width = 3.5
642 642 self.height = 5.5
643 643 self.colorbar = False
644 644 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
645 645
646 646 def update(self, dataOut):
647 647 data = {}
648 648 meta = {}
649 649
650 650 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
651 651 data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS]
652 652 data['den_error'] = dataOut.sdp2[:dataOut.NSHTS]
653 653 #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS]
654 654
655 655 data['NSHTS'] = dataOut.NSHTS
656 656
657 657 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
658 658
659 659 return data, meta
660 660
661 661 def plot(self):
662 662
663 663 y = self.data.yrange
664 664
665 665 self.xmin = 1e3
666 666 self.xmax = 1e7
667 667
668 668 ax = self.axes[0]
669 669
670 670 data = self.data[-1]
671 671
672 672 DenPow = data['den_power']
673 673 DenFar = data['den_Faraday']
674 674 errDenPow = data['den_error']
675 675 #errFaraday = data['err_Faraday']
676 676
677 677 NSHTS = data['NSHTS']
678 678
679 679 if self.CODE == 'denLP':
680 680 DenPowLP = data['den_LP']
681 681 errDenPowLP = data['den_LP_error']
682 682 cut = data['cut']
683 683
684 684 if ax.firsttime:
685 685 self.autoxticks=False
686 686 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
687 687 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2)
688 688 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
689 689 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2)
690 690
691 691 if self.CODE=='denLP':
692 692 ax.errorbar(DenPowLP[cut:], y[cut:], xerr=errDenPowLP[cut:], fmt='r^-',elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
693 693
694 694 plt.legend(loc='upper left',fontsize=8.5)
695 695 #plt.legend(loc='lower left',fontsize=8.5)
696 696 ax.set_xscale("log", nonposx='clip')
697 697 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
698 698 self.ystep_given=100
699 699 if self.CODE=='denLP':
700 700 self.ystep_given=200
701 701 ax.set_yticks(grid_y_ticks,minor=True)
702 702 ax.grid(which='minor')
703 703
704 704 else:
705 705 dataBefore = self.data[-2]
706 706 DenPowBefore = dataBefore['den_power']
707 707 self.clear_figures()
708 708 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
709 709 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2)
710 710 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
711 711 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2)
712 712 ax.errorbar(DenPowBefore, y[:NSHTS], elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed")
713 713
714 714 if self.CODE=='denLP':
715 715 ax.errorbar(DenPowLP[cut:], y[cut:], fmt='r^-', xerr=errDenPowLP[cut:],elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
716 716
717 717 ax.set_xscale("log", nonposx='clip')
718 718 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
719 719 ax.set_yticks(grid_y_ticks,minor=True)
720 720 ax.grid(which='minor')
721 721 plt.legend(loc='upper left',fontsize=8.5)
722 722 #plt.legend(loc='lower left',fontsize=8.5)
723 723
724 724 class FaradayAnglePlot(Plot):
725 725 '''
726 726 Plot for electron density
727 727 '''
728 728
729 729 CODE = 'angle'
730 730 plot_name = 'Faraday Angle'
731 731 plot_type = 'scatterbuffer'
732 732
733 733 def setup(self):
734 734
735 735 self.ncols = 1
736 736 self.nrows = 1
737 737 self.nplots = 1
738 738 self.ylabel = 'Range [km]'
739 739 self.xlabel = 'Faraday Angle (º)'
740 740 self.titles = ['Electron Density']
741 741 self.width = 3.5
742 742 self.height = 5.5
743 743 self.colorbar = False
744 744 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
745 745
746 746 def update(self, dataOut):
747 747 data = {}
748 748 meta = {}
749 749
750 750 data['angle'] = numpy.degrees(dataOut.phi)
751 751 #'''
752 print(dataOut.phi_uwrp)
753 print(data['angle'])
754 exit(1)
752 #print(dataOut.phi_uwrp)
753 #print(data['angle'])
754 #exit(1)
755 755 #'''
756 756 data['dphi'] = dataOut.dphi_uc*10
757 757 #print(dataOut.dphi)
758 758
759 759 #data['NSHTS'] = dataOut.NSHTS
760 760
761 761 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
762 762
763 763 return data, meta
764 764
765 765 def plot(self):
766 766
767 767 data = self.data[-1]
768 768 self.x = data[self.CODE]
769 769 dphi = data['dphi']
770 770 self.y = self.data.yrange
771 771 self.xmin = -360#-180
772 772 self.xmax = 360#180
773 773 ax = self.axes[0]
774 774
775 775 if ax.firsttime:
776 776 self.autoxticks=False
777 777 #if self.CODE=='den':
778 778 ax.plot(self.x, self.y,marker='o',color='g',linewidth=1.0,markersize=2)
779 779 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
780 780
781 781 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
782 782 self.ystep_given=100
783 783 if self.CODE=='denLP':
784 784 self.ystep_given=200
785 785 ax.set_yticks(grid_y_ticks,minor=True)
786 786 ax.grid(which='minor')
787 787 #plt.tight_layout()
788 788 else:
789 789
790 790 self.clear_figures()
791 791 #if self.CODE=='den':
792 792 #print(numpy.shape(self.x))
793 793 ax.plot(self.x, self.y, marker='o',color='g',linewidth=1.0, markersize=2)
794 794 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
795 795
796 796 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
797 797 ax.set_yticks(grid_y_ticks,minor=True)
798 798 ax.grid(which='minor')
799 799
800 800 class EDensityHPPlot(EDensityPlot):
801 801
802 802 '''
803 803 Plot for Electron Density Hybrid Experiment
804 804 '''
805 805
806 806 CODE = 'denLP'
807 807 plot_name = 'Electron Density'
808 808 plot_type = 'scatterbuffer'
809 809
810 810 def update(self, dataOut):
811 811 data = {}
812 812 meta = {}
813 813
814 814 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
815 815 data['den_Faraday']=dataOut.dphi[:dataOut.NSHTS]
816 816 data['den_error']=dataOut.sdp2[:dataOut.NSHTS]
817 817 data['den_LP']=dataOut.ne[:dataOut.NACF]
818 818 data['den_LP_error']=dataOut.ene[:dataOut.NACF]*dataOut.ne[:dataOut.NACF]*0.434
819 819 #self.ene=10**dataOut.ene[:dataOut.NACF]
820 820 data['NSHTS']=dataOut.NSHTS
821 821 data['cut']=dataOut.cut
822 822
823 823 return data, meta
824 824
825 825
826 826 class ACFsPlot(Plot):
827 827 '''
828 828 Plot for ACFs Double Pulse Experiment
829 829 '''
830 830
831 831 CODE = 'acfs'
832 832 #plot_name = 'ACF'
833 833 plot_type = 'scatterbuffer'
834 834
835 835
836 836 def setup(self):
837 837 self.ncols = 1
838 838 self.nrows = 1
839 839 self.nplots = 1
840 840 self.ylabel = 'Range [km]'
841 841 self.xlabel = 'Lag (ms)'
842 842 self.titles = ['ACFs']
843 843 self.width = 3.5
844 844 self.height = 5.5
845 845 self.colorbar = False
846 846 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
847 847
848 848 def update(self, dataOut):
849 849 data = {}
850 850 meta = {}
851 851
852 852 data['ACFs'] = dataOut.acfs_to_plot
853 853 data['ACFs_error'] = dataOut.acfs_error_to_plot
854 854 data['lags'] = dataOut.lags_to_plot
855 855 data['Lag_contaminated_1'] = dataOut.x_igcej_to_plot
856 856 data['Lag_contaminated_2'] = dataOut.x_ibad_to_plot
857 857 data['Height_contaminated_1'] = dataOut.y_igcej_to_plot
858 858 data['Height_contaminated_2'] = dataOut.y_ibad_to_plot
859 859
860 860 meta['yrange'] = numpy.array([])
861 861 #meta['NSHTS'] = dataOut.NSHTS
862 862 #meta['DPL'] = dataOut.DPL
863 863 data['NSHTS'] = dataOut.NSHTS #This is metadata
864 864 data['DPL'] = dataOut.DPL #This is metadata
865 865
866 866 return data, meta
867 867
868 868 def plot(self):
869 869
870 870 data = self.data[-1]
871 871 #NSHTS = self.meta['NSHTS']
872 872 #DPL = self.meta['DPL']
873 873 NSHTS = data['NSHTS'] #This is metadata
874 874 DPL = data['DPL'] #This is metadata
875 875
876 876 lags = data['lags']
877 877 ACFs = data['ACFs']
878 878 errACFs = data['ACFs_error']
879 879 BadLag1 = data['Lag_contaminated_1']
880 880 BadLag2 = data['Lag_contaminated_2']
881 881 BadHei1 = data['Height_contaminated_1']
882 882 BadHei2 = data['Height_contaminated_2']
883 883
884 884 self.xmin = 0.0
885 885 self.xmax = 2.0
886 886 self.y = ACFs
887 887
888 888 ax = self.axes[0]
889 889
890 890 if ax.firsttime:
891 891
892 892 for i in range(NSHTS):
893 893 x_aux = numpy.isfinite(lags[i,:])
894 894 y_aux = numpy.isfinite(ACFs[i,:])
895 895 yerr_aux = numpy.isfinite(errACFs[i,:])
896 896 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
897 897 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
898 898 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
899 899 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
900 900 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
901 901 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',marker='o',linewidth=1.0,markersize=2)
902 902 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
903 903 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
904 904
905 905 self.xstep_given = (self.xmax-self.xmin)/(DPL-1)
906 906 self.ystep_given = 50
907 907 ax.yaxis.set_minor_locator(MultipleLocator(15))
908 908 ax.grid(which='minor')
909 909
910 910 else:
911 911 self.clear_figures()
912 912 for i in range(NSHTS):
913 913 x_aux = numpy.isfinite(lags[i,:])
914 914 y_aux = numpy.isfinite(ACFs[i,:])
915 915 yerr_aux = numpy.isfinite(errACFs[i,:])
916 916 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
917 917 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
918 918 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
919 919 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
920 920 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
921 921 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],linewidth=1.0,markersize=2,color='b',marker='o')
922 922 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
923 923 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
924 924 ax.yaxis.set_minor_locator(MultipleLocator(15))
925 925
926 926 class ACFsLPPlot(Plot):
927 927 '''
928 928 Plot for ACFs Double Pulse Experiment
929 929 '''
930 930
931 931 CODE = 'acfs_LP'
932 932 #plot_name = 'ACF'
933 933 plot_type = 'scatterbuffer'
934 934
935 935
936 936 def setup(self):
937 937 self.ncols = 1
938 938 self.nrows = 1
939 939 self.nplots = 1
940 940 self.ylabel = 'Range [km]'
941 941 self.xlabel = 'Lag (ms)'
942 942 self.titles = ['ACFs']
943 943 self.width = 3.5
944 944 self.height = 5.5
945 945 self.colorbar = False
946 946 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
947 947
948 948 def update(self, dataOut):
949 949 data = {}
950 950 meta = {}
951 951
952 952 aux=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
953 953 errors=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
954 954 lags_LP_to_plot=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
955 955
956 956 for i in range(dataOut.NACF):
957 957 for j in range(dataOut.IBITS):
958 958 if numpy.abs(dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0])<1.0:
959 959 aux[i,j]=dataOut.output_LP_integrated.real[j,i,0]/dataOut.output_LP_integrated.real[0,i,0]
960 960 aux[i,j]=max(min(aux[i,j],1.0),-1.0)*dataOut.DH+dataOut.heightList[i]
961 961 lags_LP_to_plot[i,j]=dataOut.lags_LP[j]
962 962 errors[i,j]=dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0]*dataOut.DH
963 963 else:
964 964 aux[i,j]=numpy.nan
965 965 lags_LP_to_plot[i,j]=numpy.nan
966 966 errors[i,j]=numpy.nan
967 967
968 968 data['ACFs'] = aux
969 969 data['ACFs_error'] = errors
970 970 data['lags'] = lags_LP_to_plot
971 971
972 972 meta['yrange'] = numpy.array([])
973 973 #meta['NACF'] = dataOut.NACF
974 974 #meta['NLAG'] = dataOut.NLAG
975 975 data['NACF'] = dataOut.NACF #This is metadata
976 976 data['NLAG'] = dataOut.NLAG #This is metadata
977 977
978 978 return data, meta
979 979
980 980 def plot(self):
981 981
982 982 data = self.data[-1]
983 983 #NACF = self.meta['NACF']
984 984 #NLAG = self.meta['NLAG']
985 985 NACF = data['NACF'] #This is metadata
986 986 NLAG = data['NLAG'] #This is metadata
987 987
988 988 lags = data['lags']
989 989 ACFs = data['ACFs']
990 990 errACFs = data['ACFs_error']
991 991
992 992 self.xmin = 0.0
993 993 self.xmax = 1.5
994 994
995 995 self.y = ACFs
996 996
997 997 ax = self.axes[0]
998 998
999 999 if ax.firsttime:
1000 1000
1001 1001 for i in range(NACF):
1002 1002 x_aux = numpy.isfinite(lags[i,:])
1003 1003 y_aux = numpy.isfinite(ACFs[i,:])
1004 1004 yerr_aux = numpy.isfinite(errACFs[i,:])
1005 1005
1006 1006 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1007 1007 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1008 1008
1009 1009 #self.xstep_given = (self.xmax-self.xmin)/(self.data.NLAG-1)
1010 1010 self.xstep_given=0.3
1011 1011 self.ystep_given = 200
1012 1012 ax.yaxis.set_minor_locator(MultipleLocator(15))
1013 1013 ax.grid(which='minor')
1014 1014
1015 1015 else:
1016 1016 self.clear_figures()
1017 1017
1018 1018 for i in range(NACF):
1019 1019 x_aux = numpy.isfinite(lags[i,:])
1020 1020 y_aux = numpy.isfinite(ACFs[i,:])
1021 1021 yerr_aux = numpy.isfinite(errACFs[i,:])
1022 1022
1023 1023 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1024 1024 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1025 1025
1026 1026 ax.yaxis.set_minor_locator(MultipleLocator(15))
1027 1027
1028 1028
1029 1029 class CrossProductsPlot(Plot):
1030 1030 '''
1031 1031 Plot for cross products
1032 1032 '''
1033 1033
1034 1034 CODE = 'crossprod'
1035 1035 plot_name = 'Cross Products'
1036 1036 plot_type = 'scatterbuffer'
1037 1037
1038 1038 def setup(self):
1039 1039
1040 1040 self.ncols = 3
1041 1041 self.nrows = 1
1042 1042 self.nplots = 3
1043 1043 self.ylabel = 'Range [km]'
1044 1044 self.titles = []
1045 1045 self.width = 3.5*self.nplots
1046 1046 self.height = 5.5
1047 1047 self.colorbar = False
1048 1048 self.plots_adjust.update({'wspace':.3, 'left': 0.12, 'right': 0.92, 'bottom': 0.1})
1049 1049
1050 1050
1051 1051 def update(self, dataOut):
1052 1052
1053 1053 data = {}
1054 1054 meta = {}
1055 1055
1056 1056 data['crossprod'] = dataOut.crossprods
1057 1057 data['NDP'] = dataOut.NDP
1058 1058
1059 1059 return data, meta
1060 1060
1061 1061 def plot(self):
1062 1062
1063 1063 NDP = self.data['NDP'][-1]
1064 1064 x = self.data['crossprod'][:,-1,:,:,:,:]
1065 1065 y = self.data.yrange[0:NDP]
1066 1066
1067 1067 for n, ax in enumerate(self.axes):
1068 1068
1069 1069 self.xmin=numpy.min(numpy.concatenate((x[n][0,20:30,0,0],x[n][1,20:30,0,0],x[n][2,20:30,0,0],x[n][3,20:30,0,0])))
1070 1070 self.xmax=numpy.max(numpy.concatenate((x[n][0,20:30,0,0],x[n][1,20:30,0,0],x[n][2,20:30,0,0],x[n][3,20:30,0,0])))
1071 1071
1072 1072 if ax.firsttime:
1073 1073
1074 1074 self.autoxticks=False
1075 1075 if n==0:
1076 1076 label1='kax'
1077 1077 label2='kay'
1078 1078 label3='kbx'
1079 1079 label4='kby'
1080 1080 self.xlimits=[(self.xmin,self.xmax)]
1081 1081 elif n==1:
1082 1082 label1='kax2'
1083 1083 label2='kay2'
1084 1084 label3='kbx2'
1085 1085 label4='kby2'
1086 1086 self.xlimits.append((self.xmin,self.xmax))
1087 1087 elif n==2:
1088 1088 label1='kaxay'
1089 1089 label2='kbxby'
1090 1090 label3='kaxbx'
1091 1091 label4='kaxby'
1092 1092 self.xlimits.append((self.xmin,self.xmax))
1093 1093
1094 1094 ax.plotline1 = ax.plot(x[n][0,:,0,0], y, color='r',linewidth=2.0, label=label1)
1095 1095 ax.plotline2 = ax.plot(x[n][1,:,0,0], y, color='k',linewidth=2.0, label=label2)
1096 1096 ax.plotline3 = ax.plot(x[n][2,:,0,0], y, color='b',linewidth=2.0, label=label3)
1097 1097 ax.plotline4 = ax.plot(x[n][3,:,0,0], y, color='m',linewidth=2.0, label=label4)
1098 1098 ax.legend(loc='upper right')
1099 1099 ax.set_xlim(self.xmin, self.xmax)
1100 1100 self.titles.append('{}'.format(self.plot_name.upper()))
1101 1101
1102 1102 else:
1103 1103
1104 1104 if n==0:
1105 1105 self.xlimits=[(self.xmin,self.xmax)]
1106 1106 else:
1107 1107 self.xlimits.append((self.xmin,self.xmax))
1108 1108
1109 1109 ax.set_xlim(self.xmin, self.xmax)
1110 1110
1111 1111 ax.plotline1[0].set_data(x[n][0,:,0,0],y)
1112 1112 ax.plotline2[0].set_data(x[n][1,:,0,0],y)
1113 1113 ax.plotline3[0].set_data(x[n][2,:,0,0],y)
1114 1114 ax.plotline4[0].set_data(x[n][3,:,0,0],y)
1115 1115 self.titles.append('{}'.format(self.plot_name.upper()))
1116 1116
1117 1117
1118 1118 class CrossProductsLPPlot(Plot):
1119 1119 '''
1120 1120 Plot for cross products LP
1121 1121 '''
1122 1122
1123 1123 CODE = 'crossprodslp'
1124 1124 plot_name = 'Cross Products LP'
1125 1125 plot_type = 'scatterbuffer'
1126 1126
1127 1127
1128 1128 def setup(self):
1129 1129
1130 1130 self.ncols = 2
1131 1131 self.nrows = 1
1132 1132 self.nplots = 2
1133 1133 self.ylabel = 'Range [km]'
1134 1134 self.xlabel = 'dB'
1135 1135 self.width = 3.5*self.nplots
1136 1136 self.height = 5.5
1137 1137 self.colorbar = False
1138 1138 self.titles = []
1139 1139 self.plots_adjust.update({'wspace': .8 ,'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1140 1140
1141 1141 def update(self, dataOut):
1142 1142 data = {}
1143 1143 meta = {}
1144 1144
1145 1145 data['crossprodslp'] = 10*numpy.log10(numpy.abs(dataOut.output_LP))
1146 1146
1147 1147 data['NRANGE'] = dataOut.NRANGE #This is metadata
1148 1148 data['NLAG'] = dataOut.NLAG #This is metadata
1149 1149
1150 1150 return data, meta
1151 1151
1152 1152 def plot(self):
1153 1153
1154 1154 NRANGE = self.data['NRANGE'][-1]
1155 1155 NLAG = self.data['NLAG'][-1]
1156 1156
1157 1157 x = self.data[self.CODE][:,-1,:,:]
1158 1158 self.y = self.data.yrange[0:NRANGE]
1159 1159
1160 1160 label_array=numpy.array(['lag '+ str(x) for x in range(NLAG)])
1161 1161 color_array=['r','k','g','b','c','m','y','orange','steelblue','purple','peru','darksalmon','grey','limegreen','olive','midnightblue']
1162 1162
1163 1163
1164 1164 for n, ax in enumerate(self.axes):
1165 1165
1166 1166 self.xmin=28#30
1167 1167 self.xmax=70#70
1168 1168 #self.xmin=numpy.min(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1169 1169 #self.xmax=numpy.max(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1170 1170
1171 1171 if ax.firsttime:
1172 1172
1173 1173 self.autoxticks=False
1174 1174 if n == 0:
1175 1175 self.plotline_array=numpy.zeros((2,NLAG),dtype=object)
1176 1176
1177 1177 for i in range(NLAG):
1178 1178 self.plotline_array[n,i], = ax.plot(x[i,:,n], self.y, color=color_array[i],linewidth=1.0, label=label_array[i])
1179 1179
1180 1180 ax.legend(loc='upper right')
1181 1181 ax.set_xlim(self.xmin, self.xmax)
1182 1182 if n==0:
1183 1183 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1184 1184 if n==1:
1185 1185 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1186 1186 else:
1187 1187 for i in range(NLAG):
1188 1188 self.plotline_array[n,i].set_data(x[i,:,n],self.y)
1189 1189
1190 1190 if n==0:
1191 1191 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1192 1192 if n==1:
1193 1193 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1194 1194
1195 1195
1196 1196 class NoiseDPPlot(NoisePlot):
1197 1197 '''
1198 1198 Plot for noise Double Pulse
1199 1199 '''
1200 1200
1201 1201 CODE = 'noise'
1202 1202 #plot_name = 'Noise'
1203 1203 #plot_type = 'scatterbuffer'
1204 1204
1205 1205 def update(self, dataOut):
1206 1206
1207 1207 data = {}
1208 1208 meta = {}
1209 1209 data['noise'] = 10*numpy.log10(dataOut.noise_final)
1210 1210
1211 1211 return data, meta
1212 1212
1213 1213
1214 1214 class XmitWaveformPlot(Plot):
1215 1215 '''
1216 1216 Plot for xmit waveform
1217 1217 '''
1218 1218
1219 1219 CODE = 'xmit'
1220 1220 plot_name = 'Xmit Waveform'
1221 1221 plot_type = 'scatterbuffer'
1222 1222
1223 1223
1224 1224 def setup(self):
1225 1225
1226 1226 self.ncols = 1
1227 1227 self.nrows = 1
1228 1228 self.nplots = 1
1229 1229 self.ylabel = ''
1230 1230 self.xlabel = 'Number of Lag'
1231 1231 self.width = 5.5
1232 1232 self.height = 3.5
1233 1233 self.colorbar = False
1234 1234 self.plots_adjust.update({'right': 0.85 })
1235 1235 self.titles = [self.plot_name]
1236 1236 #self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1237 1237
1238 1238 #if not self.titles:
1239 1239 #self.titles = self.data.parameters \
1240 1240 #if self.data.parameters else ['{}'.format(self.plot_name.upper())]
1241 1241
1242 1242 def update(self, dataOut):
1243 1243
1244 1244 data = {}
1245 1245 meta = {}
1246 1246
1247 1247 y_1=numpy.arctan2(dataOut.output_LP[:,0,2].imag,dataOut.output_LP[:,0,2].real)* 180 / (numpy.pi*10)
1248 1248 y_2=numpy.abs(dataOut.output_LP[:,0,2])
1249 1249 norm=numpy.max(y_2)
1250 1250 norm=max(norm,0.1)
1251 1251 y_2=y_2/norm
1252 1252
1253 1253 meta['yrange'] = numpy.array([])
1254 1254
1255 1255 data['xmit'] = numpy.vstack((y_1,y_2))
1256 1256 data['NLAG'] = dataOut.NLAG
1257 1257
1258 1258 return data, meta
1259 1259
1260 1260 def plot(self):
1261 1261
1262 1262 data = self.data[-1]
1263 1263 NLAG = data['NLAG']
1264 1264 x = numpy.arange(0,NLAG,1,'float32')
1265 1265 y = data['xmit']
1266 1266
1267 1267 self.xmin = 0
1268 1268 self.xmax = NLAG-1
1269 1269 self.ymin = -1.0
1270 1270 self.ymax = 1.0
1271 1271 ax = self.axes[0]
1272 1272
1273 1273 if ax.firsttime:
1274 1274 ax.plotline0=ax.plot(x,y[0,:],color='blue')
1275 1275 ax.plotline1=ax.plot(x,y[1,:],color='red')
1276 1276 secax=ax.secondary_xaxis(location=0.5)
1277 1277 secax.xaxis.tick_bottom()
1278 1278 secax.tick_params( labelleft=False, labeltop=False,
1279 1279 labelright=False, labelbottom=False)
1280 1280
1281 1281 self.xstep_given = 3
1282 1282 self.ystep_given = .25
1283 1283 secax.set_xticks(numpy.linspace(self.xmin, self.xmax, 6)) #only works on matplotlib.version>3.2
1284 1284
1285 1285 else:
1286 1286 ax.plotline0[0].set_data(x,y[0,:])
1287 1287 ax.plotline1[0].set_data(x,y[1,:])
Show file before | Show file after | Hide comments
@@ -1,1615 +1,1614
1 1 """
2 2 Created on Jul 2, 2014
3 3
4 4 @author: roj-idl71
5 5 """
6 6 import os
7 7 import sys
8 8 import glob
9 9 import time
10 10 import numpy
11 11 import fnmatch
12 12 import inspect
13 13 import time
14 14 import datetime
15 15 import zmq
16 16
17 17 from schainpy.model.proc.jroproc_base import Operation, MPDecorator
18 18 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
19 19 from schainpy.model.data.jroheaderIO import get_dtype_index, get_numpy_dtype, get_procflag_dtype, get_dtype_width
20 20 from schainpy.utils import log
21 21 import schainpy.admin
22 22
23 23 LOCALTIME = True
24 24 DT_DIRECTIVES = {
25 25 '%Y': 4,
26 26 '%y': 2,
27 27 '%m': 2,
28 28 '%d': 2,
29 29 '%j': 3,
30 30 '%H': 2,
31 31 '%M': 2,
32 32 '%S': 2,
33 33 '%f': 6
34 34 }
35 35
36 36
37 37 def isNumber(cad):
38 38 """
39 39 Chequea si el conjunto de caracteres que componen un string puede ser convertidos a un numero.
40 40
41 41 Excepciones:
42 42 Si un determinado string no puede ser convertido a numero
43 43 Input:
44 44 str, string al cual se le analiza para determinar si convertible a un numero o no
45 45
46 46 Return:
47 47 True : si el string es uno numerico
48 48 False : no es un string numerico
49 49 """
50 50 try:
51 51 float(cad)
52 52 return True
53 53 except:
54 54 return False
55 55
56 56
57 57 def isFileInEpoch(filename, startUTSeconds, endUTSeconds):
58 58 """
59 59 Esta funcion determina si un archivo de datos se encuentra o no dentro del rango de fecha especificado.
60 60
61 61 Inputs:
62 62 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
63 63
64 64 startUTSeconds : fecha inicial del rango seleccionado. La fecha esta dada en
65 65 segundos contados desde 01/01/1970.
66 66 endUTSeconds : fecha final del rango seleccionado. La fecha esta dada en
67 67 segundos contados desde 01/01/1970.
68 68
69 69 Return:
70 70 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
71 71 fecha especificado, de lo contrario retorna False.
72 72
73 73 Excepciones:
74 74 Si el archivo no existe o no puede ser abierto
75 75 Si la cabecera no puede ser leida.
76 76
77 77 """
78 78 basicHeaderObj = BasicHeader(LOCALTIME)
79 79
80 80 try:
81 81
82 82 fp = open(filename, 'rb')
83 83 except IOError:
84 84 print("The file %s can't be opened" % (filename))
85 85 return 0
86 86
87 87 sts = basicHeaderObj.read(fp)
88 88 fp.close()
89 89
90 90 if not(sts):
91 91 print("Skipping the file %s because it has not a valid header" % (filename))
92 92 return 0
93 93
94 94 if not ((startUTSeconds <= basicHeaderObj.utc) and (endUTSeconds > basicHeaderObj.utc)):
95 95 return 0
96 96
97 97 return 1
98 98
99 99
100 100 def isTimeInRange(thisTime, startTime, endTime):
101 101 if endTime >= startTime:
102 102 if (thisTime < startTime) or (thisTime > endTime):
103 103 return 0
104 104 return 1
105 105 else:
106 106 if (thisTime < startTime) and (thisTime > endTime):
107 107 return 0
108 108 return 1
109 109
110 110
111 111 def isFileInTimeRange(filename, startDate, endDate, startTime, endTime):
112 112 """
113 113 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
114 114
115 115 Inputs:
116 116 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
117 117
118 118 startDate : fecha inicial del rango seleccionado en formato datetime.date
119 119
120 120 endDate : fecha final del rango seleccionado en formato datetime.date
121 121
122 122 startTime : tiempo inicial del rango seleccionado en formato datetime.time
123 123
124 124 endTime : tiempo final del rango seleccionado en formato datetime.time
125 125
126 126 Return:
127 127 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
128 128 fecha especificado, de lo contrario retorna False.
129 129
130 130 Excepciones:
131 131 Si el archivo no existe o no puede ser abierto
132 132 Si la cabecera no puede ser leida.
133 133
134 134 """
135 135
136 136 try:
137 137 fp = open(filename, 'rb')
138 138 except IOError:
139 139 print("The file %s can't be opened" % (filename))
140 140 return None
141 141
142 142 firstBasicHeaderObj = BasicHeader(LOCALTIME)
143 143 systemHeaderObj = SystemHeader()
144 144
145 145 radarControllerHeaderObj = RadarControllerHeader()
146 146 processingHeaderObj = ProcessingHeader()
147 147
148 148 lastBasicHeaderObj = BasicHeader(LOCALTIME)
149 149
150 150 sts = firstBasicHeaderObj.read(fp)
151 151
152 152 if not(sts):
153 153 print("[Reading] Skipping the file %s because it has not a valid header" % (filename))
154 154 return None
155 155
156 156 if not systemHeaderObj.read(fp):
157 157 return None
158 158
159 159 if not radarControllerHeaderObj.read(fp):
160 160 return None
161 161
162 162 if not processingHeaderObj.read(fp):
163 163 return None
164 164
165 165 filesize = os.path.getsize(filename)
166 166
167 167 offset = processingHeaderObj.blockSize + 24 # header size
168 168
169 169 if filesize <= offset:
170 170 print("[Reading] %s: This file has not enough data" % filename)
171 171 return None
172 172
173 173 fp.seek(-offset, 2)
174 174
175 175 sts = lastBasicHeaderObj.read(fp)
176 176
177 177 fp.close()
178 178
179 179 thisDatetime = lastBasicHeaderObj.datatime
180 180 thisTime_last_block = thisDatetime.time()
181 181
182 182 thisDatetime = firstBasicHeaderObj.datatime
183 183 thisDate = thisDatetime.date()
184 184 thisTime_first_block = thisDatetime.time()
185 185
186 186 # General case
187 187 # o>>>>>>>>>>>>>><<<<<<<<<<<<<<o
188 188 #-----------o----------------------------o-----------
189 189 # startTime endTime
190 190
191 191 if endTime >= startTime:
192 192 if (thisTime_last_block < startTime) or (thisTime_first_block > endTime):
193 193 return None
194 194
195 195 return thisDatetime
196 196
197 197 # If endTime < startTime then endTime belongs to the next day
198 198
199 199 #<<<<<<<<<<<o o>>>>>>>>>>>
200 200 #-----------o----------------------------o-----------
201 201 # endTime startTime
202 202
203 203 if (thisDate == startDate) and (thisTime_last_block < startTime):
204 204 return None
205 205
206 206 if (thisDate == endDate) and (thisTime_first_block > endTime):
207 207 return None
208 208
209 209 if (thisTime_last_block < startTime) and (thisTime_first_block > endTime):
210 210 return None
211 211
212 212 return thisDatetime
213 213
214 214
215 215 def isFolderInDateRange(folder, startDate=None, endDate=None):
216 216 """
217 217 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
218 218
219 219 Inputs:
220 220 folder : nombre completo del directorio.
221 221 Su formato deberia ser "/path_root/?YYYYDDD"
222 222
223 223 siendo:
224 224 YYYY : Anio (ejemplo 2015)
225 225 DDD : Dia del anio (ejemplo 305)
226 226
227 227 startDate : fecha inicial del rango seleccionado en formato datetime.date
228 228
229 229 endDate : fecha final del rango seleccionado en formato datetime.date
230 230
231 231 Return:
232 232 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
233 233 fecha especificado, de lo contrario retorna False.
234 234 Excepciones:
235 235 Si el directorio no tiene el formato adecuado
236 236 """
237 237
238 238 basename = os.path.basename(folder)
239 239
240 240 if not isRadarFolder(basename):
241 241 print("The folder %s has not the rigth format" % folder)
242 242 return 0
243 243
244 244 if startDate and endDate:
245 245 thisDate = getDateFromRadarFolder(basename)
246 246
247 247 if thisDate < startDate:
248 248 return 0
249 249
250 250 if thisDate > endDate:
251 251 return 0
252 252
253 253 return 1
254 254
255 255
256 256 def isFileInDateRange(filename, startDate=None, endDate=None):
257 257 """
258 258 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
259 259
260 260 Inputs:
261 261 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
262 262
263 263 Su formato deberia ser "?YYYYDDDsss"
264 264
265 265 siendo:
266 266 YYYY : Anio (ejemplo 2015)
267 267 DDD : Dia del anio (ejemplo 305)
268 268 sss : set
269 269
270 270 startDate : fecha inicial del rango seleccionado en formato datetime.date
271 271
272 272 endDate : fecha final del rango seleccionado en formato datetime.date
273 273
274 274 Return:
275 275 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
276 276 fecha especificado, de lo contrario retorna False.
277 277 Excepciones:
278 278 Si el archivo no tiene el formato adecuado
279 279 """
280 280
281 281 basename = os.path.basename(filename)
282 282
283 283 if not isRadarFile(basename):
284 284 print("The filename %s has not the rigth format" % filename)
285 285 return 0
286 286
287 287 if startDate and endDate:
288 288 thisDate = getDateFromRadarFile(basename)
289 289
290 290 if thisDate < startDate:
291 291 return 0
292 292
293 293 if thisDate > endDate:
294 294 return 0
295 295
296 296 return 1
297 297
298 298
299 299 def getFileFromSet(path, ext, set):
300 300 validFilelist = []
301 301 fileList = os.listdir(path)
302 302
303 303 # 0 1234 567 89A BCDE
304 304 # H YYYY DDD SSS .ext
305 305
306 306 for thisFile in fileList:
307 307 try:
308 308 year = int(thisFile[1:5])
309 309 doy = int(thisFile[5:8])
310 310 except:
311 311 continue
312 312
313 313 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
314 314 continue
315 315
316 316 validFilelist.append(thisFile)
317 317
318 318 myfile = fnmatch.filter(
319 319 validFilelist, '*%4.4d%3.3d%3.3d*' % (year, doy, set))
320 320
321 321 if len(myfile) != 0:
322 322 return myfile[0]
323 323 else:
324 324 filename = '*%4.4d%3.3d%3.3d%s' % (year, doy, set, ext.lower())
325 325 print('the filename %s does not exist' % filename)
326 326 print('...going to the last file: ')
327 327
328 328 if validFilelist:
329 329 validFilelist = sorted(validFilelist, key=str.lower)
330 330 return validFilelist[-1]
331 331
332 332 return None
333 333
334 334
335 335 def getlastFileFromPath(path, ext):
336 336 """
337 337 Depura el fileList dejando solo los que cumplan el formato de "PYYYYDDDSSS.ext"
338 338 al final de la depuracion devuelve el ultimo file de la lista que quedo.
339 339
340 340 Input:
341 341 fileList : lista conteniendo todos los files (sin path) que componen una determinada carpeta
342 342 ext : extension de los files contenidos en una carpeta
343 343
344 344 Return:
345 345 El ultimo file de una determinada carpeta, no se considera el path.
346 346 """
347 347 validFilelist = []
348 348 fileList = os.listdir(path)
349 349
350 350 # 0 1234 567 89A BCDE
351 351 # H YYYY DDD SSS .ext
352 352
353 353 for thisFile in fileList:
354 354
355 355 year = thisFile[1:5]
356 356 if not isNumber(year):
357 357 continue
358 358
359 359 doy = thisFile[5:8]
360 360 if not isNumber(doy):
361 361 continue
362 362
363 363 year = int(year)
364 364 doy = int(doy)
365 365
366 366 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
367 367 continue
368 368
369 369 validFilelist.append(thisFile)
370 370
371 371 if validFilelist:
372 372 validFilelist = sorted(validFilelist, key=str.lower)
373 373 return validFilelist[-1]
374 374
375 375 return None
376 376
377 377
378 378 def isRadarFolder(folder):
379 379 try:
380 380 year = int(folder[1:5])
381 381 doy = int(folder[5:8])
382 382 except:
383 383 return 0
384 384
385 385 return 1
386 386
387 387
388 388 def isRadarFile(file):
389 389 try:
390 390 year = int(file[1:5])
391 391 doy = int(file[5:8])
392 392 set = int(file[8:11])
393 393 except:
394 394 return 0
395 395
396 396 return 1
397 397
398 398
399 399 def getDateFromRadarFile(file):
400 400 try:
401 401 year = int(file[1:5])
402 402 doy = int(file[5:8])
403 403 set = int(file[8:11])
404 404 except:
405 405 return None
406 406
407 407 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
408 408 return thisDate
409 409
410 410
411 411 def getDateFromRadarFolder(folder):
412 412 try:
413 413 year = int(folder[1:5])
414 414 doy = int(folder[5:8])
415 415 except:
416 416 return None
417 417
418 418 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
419 419 return thisDate
420 420
421 421 def parse_format(s, fmt):
422 422
423 423 for i in range(fmt.count('%')):
424 424 x = fmt.index('%')
425 425 d = DT_DIRECTIVES[fmt[x:x+2]]
426 426 fmt = fmt.replace(fmt[x:x+2], s[x:x+d])
427 427 return fmt
428 428
429 429 class Reader(object):
430 430
431 431 c = 3E8
432 432 isConfig = False
433 433 dtype = None
434 434 pathList = []
435 435 filenameList = []
436 436 datetimeList = []
437 437 filename = None
438 438 ext = None
439 439 flagIsNewFile = 1
440 440 flagDiscontinuousBlock = 0
441 441 flagIsNewBlock = 0
442 442 flagNoMoreFiles = 0
443 443 fp = None
444 444 firstHeaderSize = 0
445 445 basicHeaderSize = 24
446 446 versionFile = 1103
447 447 fileSize = None
448 448 fileSizeByHeader = None
449 449 fileIndex = -1
450 450 profileIndex = None
451 451 blockIndex = 0
452 452 nTotalBlocks = 0
453 453 maxTimeStep = 30
454 454 lastUTTime = None
455 455 datablock = None
456 456 dataOut = None
457 457 getByBlock = False
458 458 path = None
459 459 startDate = None
460 460 endDate = None
461 461 startTime = datetime.time(0, 0, 0)
462 462 endTime = datetime.time(23, 59, 59)
463 463 set = None
464 464 expLabel = ""
465 465 online = False
466 466 delay = 60
467 467 nTries = 3 # quantity tries
468 468 nFiles = 3 # number of files for searching
469 469 walk = True
470 470 getblock = False
471 471 nTxs = 1
472 472 realtime = False
473 473 blocksize = 0
474 474 blocktime = None
475 475 warnings = True
476 476 verbose = True
477 477 server = None
478 478 format = None
479 479 oneDDict = None
480 480 twoDDict = None
481 481 independentParam = None
482 482 filefmt = None
483 483 folderfmt = None
484 484 open_file = open
485 485 open_mode = 'rb'
486 486
487 487 def run(self):
488 488
489 489 raise NotImplementedError
490 490
491 491 def getAllowedArgs(self):
492 492 if hasattr(self, '__attrs__'):
493 493 return self.__attrs__
494 494 else:
495 495 return inspect.getargspec(self.run).args
496 496
497 497 def set_kwargs(self, **kwargs):
498 498
499 499 for key, value in kwargs.items():
500 500 setattr(self, key, value)
501 501
502 502 def find_folders(self, path, startDate, endDate, folderfmt, last=False):
503 503
504 504 folders = [x for f in path.split(',')
505 505 for x in os.listdir(f) if os.path.isdir(os.path.join(f, x))]
506 506 folders.sort()
507 507
508 508 if last:
509 509 folders = [folders[-1]]
510 510
511 511 for folder in folders:
512 512 try:
513 513 dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date()
514 514 if dt >= startDate and dt <= endDate:
515 515 yield os.path.join(path, folder)
516 516 else:
517 517 log.log('Skiping folder {}'.format(folder), self.name)
518 518 except Exception as e:
519 519 log.log('Skiping folder {}'.format(folder), self.name)
520 520 continue
521 521 return
522 522
523 523 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
524 524 expLabel='', last=False):
525
526 525 for path in folders:
527 files = glob.glob1(path, '*{}'.format(ext))
526 files = glob.glob1(path+'/'+expLabel, '*{}'.format(ext))
528 527 files.sort()
529 528 if last:
530 529 if files:
531 530 fo = files[-1]
532 531 try:
533 532 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
534 533 yield os.path.join(path, expLabel, fo)
535 534 except Exception as e:
536 535 pass
537 536 return
538 537 else:
539 538 return
540 539
541 540 for fo in files:
542 541 try:
543 542 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
544 543 #print(dt)
545 544 #print(startDate)
546 545 #print(endDate)
547 546 if dt >= startDate and dt <= endDate:
548 547
549 548 yield os.path.join(path, expLabel, fo)
550 549
551 550 else:
552 551
553 552 log.log('Skiping file {}'.format(fo), self.name)
554 553 except Exception as e:
555 554 log.log('Skiping file {}'.format(fo), self.name)
556 555 continue
557 556
558 557 def searchFilesOffLine(self, path, startDate, endDate,
559 558 expLabel, ext, walk,
560 559 filefmt, folderfmt):
561 560 """Search files in offline mode for the given arguments
562 561
563 562 Return:
564 563 Generator of files
565 564 """
566 565
567 566 if walk:
568 567 folders = self.find_folders(
569 568 path, startDate, endDate, folderfmt)
569 #print("folders: ", folders)
570 570 else:
571 571 folders = path.split(',')
572 572
573 573 return self.find_files(
574 574 folders, ext, filefmt, startDate, endDate, expLabel)
575 575
576 576 def searchFilesOnLine(self, path, startDate, endDate,
577 577 expLabel, ext, walk,
578 578 filefmt, folderfmt):
579 579 """Search for the last file of the last folder
580 580
581 581 Arguments:
582 582 path : carpeta donde estan contenidos los files que contiene data
583 583 expLabel : Nombre del subexperimento (subfolder)
584 584 ext : extension de los files
585 585 walk : Si es habilitado no realiza busquedas dentro de los ubdirectorios (doypath)
586 586
587 587 Return:
588 588 generator with the full path of last filename
589 589 """
590 590
591 591 if walk:
592 592 folders = self.find_folders(
593 593 path, startDate, endDate, folderfmt, last=True)
594 594 else:
595 595 folders = path.split(',')
596 596
597 597 return self.find_files(
598 598 folders, ext, filefmt, startDate, endDate, expLabel, last=True)
599 599
600 600 def setNextFile(self):
601 601 """Set the next file to be readed open it and parse de file header"""
602 602
603 603 #print("fp: ",self.fp)
604 604 while True:
605 605
606 606 #print(self.fp)
607 607 if self.fp != None:
608 608 self.fp.close()
609 609
610 610 #print("setNextFile")
611 611 #print("BEFORE OPENING",self.filename)
612 612 if self.online:
613 613 newFile = self.setNextFileOnline()
614 614
615 615 else:
616 616
617 617 newFile = self.setNextFileOffline()
618 618
619 619 #print("newFile: ",newFile)
620 620 if not(newFile):
621 621
622 622 if self.online:
623 623 raise schainpy.admin.SchainError('Time to wait for new files reach')
624 624 else:
625 625 if self.fileIndex == -1:
626 626 #print("OKK")
627 627 raise schainpy.admin.SchainWarning('No files found in the given path')
628 628 else:
629 629
630 630 raise schainpy.admin.SchainWarning('No more files to read')
631 631
632 632 if self.verifyFile(self.filename):
633 633
634 634 break
635 635
636 636 ##print("BEFORE OPENING",self.filename)
637 637
638 638 log.log('Opening file: %s' % self.filename, self.name)
639 639
640 640 self.readFirstHeader()
641 641 self.nReadBlocks = 0
642 642
643 643 def setNextFileOnline(self):
644 644 """Check for the next file to be readed in online mode.
645 645
646 646 Set:
647 647 self.filename
648 648 self.fp
649 649 self.filesize
650 650
651 651 Return:
652 652 boolean
653 653
654 654 """
655 655
656 656 nextFile = True
657 657 nextDay = False
658 658
659 659 for nFiles in range(self.nFiles+1):
660 660 for nTries in range(self.nTries):
661 661 fullfilename, filename = self.checkForRealPath(nextFile, nextDay)
662 662 if fullfilename is not None:
663 663 break
664 664 log.warning(
665 665 "Waiting %0.2f sec for the next file: \"%s\" , try %02d ..." % (self.delay, filename, nTries + 1),
666 666 self.name)
667 667 time.sleep(self.delay)
668 668 nextFile = False
669 669 continue
670 670
671 671 if fullfilename is not None:
672 672 break
673 673
674 674 #self.nTries = 1
675 675 nextFile = True
676 676
677 677 if nFiles == (self.nFiles - 1):
678 678 log.log('Trying with next day...', self.name)
679 679 nextDay = True
680 680 self.nTries = 3
681 681
682 682 if fullfilename:
683 683 self.fileSize = os.path.getsize(fullfilename)
684 684 self.filename = fullfilename
685 685 self.flagIsNewFile = 1
686 686 if self.fp != None:
687 687 self.fp.close()
688 688 #print(fullfilename)
689 689 self.fp = self.open_file(fullfilename, self.open_mode)
690 690
691 691 self.flagNoMoreFiles = 0
692 692 self.fileIndex += 1
693 693 return 1
694 694 else:
695 695 return 0
696 696
697 697 def setNextFileOffline(self):
698 698 """Open the next file to be readed in offline mode"""
699 699
700 700 try:
701 701 filename = next(self.filenameList)
702 702 self.fileIndex +=1
703 703 except StopIteration:
704 704 self.flagNoMoreFiles = 1
705 705 return 0
706 706 #print(self.fileIndex)
707 707 #print(filename)
708 708 self.filename = filename
709 709 self.fileSize = os.path.getsize(filename)
710 710 self.fp = self.open_file(filename, self.open_mode)
711 711 self.flagIsNewFile = 1
712 712
713 713 return 1
714 714
715 715 @staticmethod
716 716 def isDateTimeInRange(dt, startDate, endDate, startTime, endTime):
717 717 """Check if the given datetime is in range"""
718 718
719 719 if startDate <= dt.date() <= endDate:
720 720 if startTime <= dt.time() <= endTime:
721 721 return True
722 722 return False
723 723
724 724 def verifyFile(self, filename):
725 725 """Check for a valid file
726 726
727 727 Arguments:
728 728 filename -- full path filename
729 729
730 730 Return:
731 731 boolean
732 732 """
733 733
734 734 return True
735 735
736 736 def checkForRealPath(self, nextFile, nextDay):
737 737 """Check if the next file to be readed exists"""
738 738
739 739 raise NotImplementedError
740 740
741 741 def readFirstHeader(self):
742 742 """Parse the file header"""
743 743
744 744
745 745 pass
746 746
747 747 def waitDataBlock(self, pointer_location, blocksize=None):
748 748 """
749 749 """
750 750
751 751 currentPointer = pointer_location
752 752 if blocksize is None:
753 753 neededSize = self.processingHeaderObj.blockSize # + self.basicHeaderSize
754 754 else:
755 755 neededSize = blocksize
756 756
757 757 for nTries in range(self.nTries):
758 758 self.fp.close()
759 759 self.fp = open(self.filename, 'rb')
760 760 self.fp.seek(currentPointer)
761 761
762 762 self.fileSize = os.path.getsize(self.filename)
763 763 currentSize = self.fileSize - currentPointer
764 764
765 765 if (currentSize >= neededSize):
766 766 return 1
767 767
768 768 log.warning(
769 769 "Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1),
770 770 self.name
771 771 )
772 772 time.sleep(self.delay)
773 773
774 774 return 0
775 775
776 776 class JRODataReader(Reader):
777 777
778 778 utc = 0
779 779 nReadBlocks = 0
780 780 foldercounter = 0
781 781 firstHeaderSize = 0
782 782 basicHeaderSize = 24
783 783 __isFirstTimeOnline = 1
784 784 filefmt = "*%Y%j***"
785 785 folderfmt = "*%Y%j"
786 786 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'online', 'delay', 'walk']
787 787
788 788 def getDtypeWidth(self):
789 789
790 790 dtype_index = get_dtype_index(self.dtype)
791 791 dtype_width = get_dtype_width(dtype_index)
792 792
793 793 return dtype_width
794 794
795 795 def checkForRealPath(self, nextFile, nextDay):
796 796 """Check if the next file to be readed exists.
797 797
798 798 Example :
799 799 nombre correcto del file es .../.../D2009307/P2009307367.ext
800 800
801 801 Entonces la funcion prueba con las siguientes combinaciones
802 802 .../.../y2009307367.ext
803 803 .../.../Y2009307367.ext
804 804 .../.../x2009307/y2009307367.ext
805 805 .../.../x2009307/Y2009307367.ext
806 806 .../.../X2009307/y2009307367.ext
807 807 .../.../X2009307/Y2009307367.ext
808 808 siendo para este caso, la ultima combinacion de letras, identica al file buscado
809 809
810 810 Return:
811 811 str -- fullpath of the file
812 812 """
813 813
814 814
815 815 if nextFile:
816 816 self.set += 1
817 817 if nextDay:
818 818 self.set = 0
819 819 self.doy += 1
820 820 foldercounter = 0
821 821 prefixDirList = [None, 'd', 'D']
822 822 if self.ext.lower() == ".r": # voltage
823 823 prefixFileList = ['d', 'D']
824 824 elif self.ext.lower() == ".pdata": # spectra
825 825 prefixFileList = ['p', 'P']
826 826
827 827 ##############DP##############
828 828
829 829 elif self.ext.lower() == ".dat": # dat
830 830 prefixFileList = ['z', 'Z']
831 831
832 832
833 833
834 834 ##############DP##############
835 835 # barrido por las combinaciones posibles
836 836 for prefixDir in prefixDirList:
837 837 thispath = self.path
838 838 if prefixDir != None:
839 839 # formo el nombre del directorio xYYYYDDD (x=d o x=D)
840 840 if foldercounter == 0:
841 841 thispath = os.path.join(self.path, "%s%04d%03d" %
842 842 (prefixDir, self.year, self.doy))
843 843 else:
844 844 thispath = os.path.join(self.path, "%s%04d%03d_%02d" % (
845 845 prefixDir, self.year, self.doy, foldercounter))
846 846 for prefixFile in prefixFileList: # barrido por las dos combinaciones posibles de "D"
847 847 # formo el nombre del file xYYYYDDDSSS.ext
848 848 filename = "%s%04d%03d%03d%s" % (prefixFile, self.year, self.doy, self.set, self.ext)
849 849 fullfilename = os.path.join(
850 850 thispath, filename)
851 851
852 852 if os.path.exists(fullfilename):
853 853 return fullfilename, filename
854 854
855 855 return None, filename
856 856
857 857 def __waitNewBlock(self):
858 858 """
859 859 Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
860 860
861 861 Si el modo de lectura es OffLine siempre retorn 0
862 862 """
863 863 if not self.online:
864 864 return 0
865 865
866 866 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
867 867 return 0
868 868
869 869 currentPointer = self.fp.tell()
870 870
871 871 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
872 872
873 873 for nTries in range(self.nTries):
874 874
875 875 self.fp.close()
876 876 self.fp = open(self.filename, 'rb')
877 877 self.fp.seek(currentPointer)
878 878
879 879 self.fileSize = os.path.getsize(self.filename)
880 880 currentSize = self.fileSize - currentPointer
881 881
882 882 if (currentSize >= neededSize):
883 883 self.basicHeaderObj.read(self.fp)
884 884 return 1
885 885
886 886 if self.fileSize == self.fileSizeByHeader:
887 887 # self.flagEoF = True
888 888 return 0
889 889
890 890 print("[Reading] Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1))
891 891 #print(self.filename)
892 892 time.sleep(self.delay)
893 893
894 894 return 0
895 895
896 896 def __setNewBlock(self):
897 897
898 898 if self.fp == None:
899 899 return 0
900 900
901 901 if self.flagIsNewFile:
902 902 self.lastUTTime = self.basicHeaderObj.utc
903 903 return 1
904 904
905 905 if self.realtime:
906 906 self.flagDiscontinuousBlock = 1
907 907 if not(self.setNextFile()):
908 908 return 0
909 909 else:
910 910 return 1
911 911
912 912 currentSize = self.fileSize - self.fp.tell()
913 913 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
914 914
915 915 if (currentSize >= neededSize):
916 916 self.basicHeaderObj.read(self.fp)
917 917 self.lastUTTime = self.basicHeaderObj.utc
918 918 return 1
919 919
920 920 if self.__waitNewBlock():
921 921 self.lastUTTime = self.basicHeaderObj.utc
922 922 return 1
923 923
924 924 if not(self.setNextFile()):
925 925 return 0
926 926
927 927 deltaTime = self.basicHeaderObj.utc - self.lastUTTime
928 928 self.lastUTTime = self.basicHeaderObj.utc
929 929
930 930 self.flagDiscontinuousBlock = 0
931
932 931 if deltaTime > self.maxTimeStep:
933 932 self.flagDiscontinuousBlock = 1
934 933
935 934 return 1
936 935
937 936 def readNextBlock(self):
938 937
939 938 while True:
940 939 if not(self.__setNewBlock()):
941 940 continue
942 941
943 942 if not(self.readBlock()):
944 943 return 0
945 944
946 945 self.getBasicHeader()
947 946
948 947 if not self.isDateTimeInRange(self.dataOut.datatime, self.startDate, self.endDate, self.startTime, self.endTime):
949 948 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.nReadBlocks,
950 949 self.processingHeaderObj.dataBlocksPerFile,
951 950 self.dataOut.datatime.ctime()))
952 951 continue
953 952
954 953 break
955 954
956 955 if self.verbose:
957 956 print("[Reading] Block No. %d/%d -> %s" % (self.nReadBlocks,
958 957 self.processingHeaderObj.dataBlocksPerFile,
959 958 self.dataOut.datatime.ctime()))
960 959 #################DP#################
961 960 self.dataOut.TimeBlockDate=self.dataOut.datatime.ctime()
962 961 self.dataOut.TimeBlockSeconds=time.mktime(time.strptime(self.dataOut.datatime.ctime()))
963 962 #################DP#################
964 963 return 1
965 964
966 965 def readFirstHeader(self):
967 966
968 967 self.basicHeaderObj.read(self.fp)
969 968 self.systemHeaderObj.read(self.fp)
970 969 self.radarControllerHeaderObj.read(self.fp)
971 970 self.processingHeaderObj.read(self.fp)
972 971 self.firstHeaderSize = self.basicHeaderObj.size
973 972
974 973 datatype = int(numpy.log2((self.processingHeaderObj.processFlags &
975 974 PROCFLAG.DATATYPE_MASK)) - numpy.log2(PROCFLAG.DATATYPE_CHAR))
976 975 if datatype == 0:
977 976 datatype_str = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
978 977 elif datatype == 1:
979 978 datatype_str = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
980 979 elif datatype == 2:
981 980 datatype_str = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
982 981 elif datatype == 3:
983 982 datatype_str = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
984 983 elif datatype == 4:
985 984 datatype_str = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
986 985 elif datatype == 5:
987 986 datatype_str = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
988 987 else:
989 988 raise ValueError('Data type was not defined')
990 989
991 990 self.dtype = datatype_str
992 991 #self.ippSeconds = 2 * 1000 * self.radarControllerHeaderObj.ipp / self.c
993 992 self.fileSizeByHeader = self.processingHeaderObj.dataBlocksPerFile * self.processingHeaderObj.blockSize + \
994 993 self.firstHeaderSize + self.basicHeaderSize * \
995 994 (self.processingHeaderObj.dataBlocksPerFile - 1)
996 995 # self.dataOut.channelList = numpy.arange(self.systemHeaderObj.numChannels)
997 996 # self.dataOut.channelIndexList = numpy.arange(self.systemHeaderObj.numChannels)
998 997 self.getBlockDimension()
999 998
1000 999 def verifyFile(self, filename):
1001 1000
1002 1001 flag = True
1003 1002
1004 1003 try:
1005 1004 fp = open(filename, 'rb')
1006 1005 except IOError:
1007 1006 log.error("File {} can't be opened".format(filename), self.name)
1008 1007 return False
1009 1008
1010 1009 if self.online and self.waitDataBlock(0):
1011 1010 pass
1012 1011
1013 1012 basicHeaderObj = BasicHeader(LOCALTIME)
1014 1013 systemHeaderObj = SystemHeader()
1015 1014 radarControllerHeaderObj = RadarControllerHeader()
1016 1015 processingHeaderObj = ProcessingHeader()
1017 1016
1018 1017 if not(basicHeaderObj.read(fp)):
1019 1018 flag = False
1020 1019 if not(systemHeaderObj.read(fp)):
1021 1020 flag = False
1022 1021 if not(radarControllerHeaderObj.read(fp)):
1023 1022 flag = False
1024 1023 if not(processingHeaderObj.read(fp)):
1025 1024 flag = False
1026 1025 if not self.online:
1027 1026 dt1 = basicHeaderObj.datatime
1028 1027 pos = self.fileSize-processingHeaderObj.blockSize-24
1029 1028 if pos<0:
1030 1029 flag = False
1031 1030 log.error('Invalid size for file: {}'.format(self.filename), self.name)
1032 1031 else:
1033 1032 fp.seek(pos)
1034 1033 if not(basicHeaderObj.read(fp)):
1035 1034 flag = False
1036 1035 dt2 = basicHeaderObj.datatime
1037 1036 if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \
1038 1037 self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime):
1039 1038 flag = False
1040 1039
1041 1040 fp.close()
1042 1041 return flag
1043 1042
1044 1043 def findDatafiles(self, path, startDate=None, endDate=None, expLabel='', ext='.r', walk=True, include_path=False):
1045 1044
1046 1045 path_empty = True
1047 1046
1048 1047 dateList = []
1049 1048 pathList = []
1050 1049
1051 1050 multi_path = path.split(',')
1052 1051
1053 1052 if not walk:
1054 1053
1055 1054 for single_path in multi_path:
1056 1055
1057 1056 if not os.path.isdir(single_path):
1058 1057 continue
1059 1058
1060 1059 fileList = glob.glob1(single_path, "*" + ext)
1061 1060
1062 1061 if not fileList:
1063 1062 continue
1064 1063
1065 1064 path_empty = False
1066 1065
1067 1066 fileList.sort()
1068 1067
1069 1068 for thisFile in fileList:
1070 1069
1071 1070 if not os.path.isfile(os.path.join(single_path, thisFile)):
1072 1071 continue
1073 1072
1074 1073 if not isRadarFile(thisFile):
1075 1074 continue
1076 1075
1077 1076 if not isFileInDateRange(thisFile, startDate, endDate):
1078 1077 continue
1079 1078
1080 1079 thisDate = getDateFromRadarFile(thisFile)
1081 1080
1082 1081 if thisDate in dateList or single_path in pathList:
1083 1082 continue
1084 1083
1085 1084 dateList.append(thisDate)
1086 1085 pathList.append(single_path)
1087 1086
1088 1087 else:
1089 1088 for single_path in multi_path:
1090 1089
1091 1090 if not os.path.isdir(single_path):
1092 1091 continue
1093 1092
1094 1093 dirList = []
1095 1094
1096 1095 for thisPath in os.listdir(single_path):
1097 1096
1098 1097 if not os.path.isdir(os.path.join(single_path, thisPath)):
1099 1098 continue
1100 1099
1101 1100 if not isRadarFolder(thisPath):
1102 1101 continue
1103 1102
1104 1103 if not isFolderInDateRange(thisPath, startDate, endDate):
1105 1104 continue
1106 1105
1107 1106 dirList.append(thisPath)
1108 1107
1109 1108 if not dirList:
1110 1109 continue
1111 1110
1112 1111 dirList.sort()
1113 1112
1114 1113 for thisDir in dirList:
1115 1114
1116 1115 datapath = os.path.join(single_path, thisDir, expLabel)
1117 1116 fileList = glob.glob1(datapath, "*" + ext)
1118 1117
1119 1118 if not fileList:
1120 1119 continue
1121 1120
1122 1121 path_empty = False
1123 1122
1124 1123 thisDate = getDateFromRadarFolder(thisDir)
1125 1124
1126 1125 pathList.append(datapath)
1127 1126 dateList.append(thisDate)
1128 1127
1129 1128 dateList.sort()
1130 1129
1131 1130 if walk:
1132 1131 pattern_path = os.path.join(multi_path[0], "[dYYYYDDD]", expLabel)
1133 1132 else:
1134 1133 pattern_path = multi_path[0]
1135 1134
1136 1135 if path_empty:
1137 1136 raise schainpy.admin.SchainError("[Reading] No *%s files in %s for %s to %s" % (ext, pattern_path, startDate, endDate))
1138 1137 else:
1139 1138 if not dateList:
1140 1139 raise schainpy.admin.SchainError("[Reading] Date range selected invalid [%s - %s]: No *%s files in %s)" % (startDate, endDate, ext, path))
1141 1140
1142 1141 if include_path:
1143 1142 return dateList, pathList
1144 1143
1145 1144 return dateList
1146 1145
1147 1146 def setup(self, **kwargs):
1148 1147
1149 1148 self.set_kwargs(**kwargs)
1150 1149 if not self.ext.startswith('.'):
1151 1150 self.ext = '.{}'.format(self.ext)
1152 1151
1153 1152 if self.server is not None:
1154 1153 if 'tcp://' in self.server:
1155 1154 address = server
1156 1155 else:
1157 1156 address = 'ipc:///tmp/%s' % self.server
1158 1157 self.server = address
1159 1158 self.context = zmq.Context()
1160 1159 self.receiver = self.context.socket(zmq.PULL)
1161 1160 self.receiver.connect(self.server)
1162 1161 time.sleep(0.5)
1163 1162 print('[Starting] ReceiverData from {}'.format(self.server))
1164 1163 else:
1165 1164 self.server = None
1166 1165 if self.path == None:
1167 1166 raise ValueError("[Reading] The path is not valid")
1168 1167
1169 1168 if self.online:
1170 1169 log.log("[Reading] Searching files in online mode...", self.name)
1171 1170
1172 1171 for nTries in range(self.nTries):
1173 1172 fullpath = self.searchFilesOnLine(self.path, self.startDate,
1174 1173 self.endDate, self.expLabel, self.ext, self.walk,
1175 1174 self.filefmt, self.folderfmt)
1176 1175
1177 1176 try:
1178 1177 fullpath = next(fullpath)
1179 1178 except:
1180 1179 fullpath = None
1181 1180
1182 1181 if fullpath:
1183 1182 break
1184 1183
1185 1184 log.warning(
1186 1185 'Waiting {} sec for a valid file in {}: try {} ...'.format(
1187 1186 self.delay, self.path, nTries + 1),
1188 1187 self.name)
1189 1188 time.sleep(self.delay)
1190 1189
1191 1190 if not(fullpath):
1192 1191 raise schainpy.admin.SchainError(
1193 1192 'There isn\'t any valid file in {}'.format(self.path))
1194 1193
1195 1194 pathname, filename = os.path.split(fullpath)
1196 1195 self.year = int(filename[1:5])
1197 1196 self.doy = int(filename[5:8])
1198 1197 self.set = int(filename[8:11]) - 1
1199 1198 else:
1200 1199 log.log("Searching files in {}".format(self.path), self.name)
1201 1200 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
1202 1201 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
1203 1202
1204 1203 self.setNextFile()
1205 1204
1206 1205 return
1207 1206
1208 1207 def getBasicHeader(self):
1209 1208
1210 1209 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
1211 1210 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
1212 1211
1213 1212 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
1214 1213
1215 1214 self.dataOut.timeZone = self.basicHeaderObj.timeZone
1216 1215
1217 1216 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
1218 1217
1219 1218 self.dataOut.errorCount = self.basicHeaderObj.errorCount
1220 1219
1221 1220 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
1222 1221
1223 1222 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
1224 1223
1225 1224 def getFirstHeader(self):
1226 1225
1227 1226 raise NotImplementedError
1228 1227
1229 1228 def getData(self):
1230 1229
1231 1230 raise NotImplementedError
1232 1231
1233 1232 def hasNotDataInBuffer(self):
1234 1233
1235 1234 raise NotImplementedError
1236 1235
1237 1236 def readBlock(self):
1238 1237
1239 1238 raise NotImplementedError
1240 1239
1241 1240 def isEndProcess(self):
1242 1241
1243 1242 return self.flagNoMoreFiles
1244 1243
1245 1244 def printReadBlocks(self):
1246 1245
1247 1246 print("[Reading] Number of read blocks per file %04d" % self.nReadBlocks)
1248 1247
1249 1248 def printTotalBlocks(self):
1250 1249
1251 1250 print("[Reading] Number of read blocks %04d" % self.nTotalBlocks)
1252 1251
1253 1252 def run(self, **kwargs):
1254 1253 """
1255 1254
1256 1255 Arguments:
1257 1256 path :
1258 1257 startDate :
1259 1258 endDate :
1260 1259 startTime :
1261 1260 endTime :
1262 1261 set :
1263 1262 expLabel :
1264 1263 ext :
1265 1264 online :
1266 1265 delay :
1267 1266 walk :
1268 1267 getblock :
1269 1268 nTxs :
1270 1269 realtime :
1271 1270 blocksize :
1272 1271 blocktime :
1273 1272 skip :
1274 1273 cursor :
1275 1274 warnings :
1276 1275 server :
1277 1276 verbose :
1278 1277 format :
1279 1278 oneDDict :
1280 1279 twoDDict :
1281 1280 independentParam :
1282 1281 """
1283 1282
1284 1283 if not(self.isConfig):
1285 1284 self.setup(**kwargs)
1286 1285 self.isConfig = True
1287 1286 if self.server is None:
1288 1287 self.getData()
1289 1288 else:
1290 1289 self.getFromServer()
1291 1290
1292 1291
1293 1292 class JRODataWriter(Reader):
1294 1293
1295 1294 """
1296 1295 Esta clase permite escribir datos a archivos procesados (.r o ,pdata). La escritura
1297 1296 de los datos siempre se realiza por bloques.
1298 1297 """
1299 1298
1300 1299 setFile = None
1301 1300 profilesPerBlock = None
1302 1301 blocksPerFile = None
1303 1302 nWriteBlocks = 0
1304 1303 fileDate = None
1305 1304
1306 1305 def __init__(self, dataOut=None):
1307 1306 raise NotImplementedError
1308 1307
1309 1308 def hasAllDataInBuffer(self):
1310 1309 raise NotImplementedError
1311 1310
1312 1311 def setBlockDimension(self):
1313 1312 raise NotImplementedError
1314 1313
1315 1314 def writeBlock(self):
1316 1315 raise NotImplementedError
1317 1316
1318 1317 def putData(self):
1319 1318 raise NotImplementedError
1320 1319
1321 1320 def getDtypeWidth(self):
1322 1321
1323 1322 dtype_index = get_dtype_index(self.dtype)
1324 1323 dtype_width = get_dtype_width(dtype_index)
1325 1324
1326 1325 return dtype_width
1327 1326
1328 1327 def getProcessFlags(self):
1329 1328
1330 1329 processFlags = 0
1331 1330
1332 1331 dtype_index = get_dtype_index(self.dtype)
1333 1332 procflag_dtype = get_procflag_dtype(dtype_index)
1334 1333
1335 1334 processFlags += procflag_dtype
1336 1335
1337 1336 if self.dataOut.flagDecodeData:
1338 1337 processFlags += PROCFLAG.DECODE_DATA
1339 1338
1340 1339 if self.dataOut.flagDeflipData:
1341 1340 processFlags += PROCFLAG.DEFLIP_DATA
1342 1341
1343 1342 if self.dataOut.code is not None:
1344 1343 processFlags += PROCFLAG.DEFINE_PROCESS_CODE
1345 1344
1346 1345 if self.dataOut.nCohInt > 1:
1347 1346 processFlags += PROCFLAG.COHERENT_INTEGRATION
1348 1347
1349 1348 if self.dataOut.type == "Spectra":
1350 1349 if self.dataOut.nIncohInt > 1:
1351 1350 processFlags += PROCFLAG.INCOHERENT_INTEGRATION
1352 1351
1353 1352 if self.dataOut.data_dc is not None:
1354 1353 processFlags += PROCFLAG.SAVE_CHANNELS_DC
1355 1354
1356 1355 if self.dataOut.flagShiftFFT:
1357 1356 processFlags += PROCFLAG.SHIFT_FFT_DATA
1358 1357
1359 1358 return processFlags
1360 1359
1361 1360 def setBasicHeader(self):
1362 1361
1363 1362 self.basicHeaderObj.size = self.basicHeaderSize # bytes
1364 1363 self.basicHeaderObj.version = self.versionFile
1365 1364 self.basicHeaderObj.dataBlock = self.nTotalBlocks
1366 1365 utc = numpy.floor(self.dataOut.utctime)
1367 1366 milisecond = (self.dataOut.utctime - utc) * 1000.0
1368 1367 self.basicHeaderObj.utc = utc
1369 1368 self.basicHeaderObj.miliSecond = milisecond
1370 1369 self.basicHeaderObj.timeZone = self.dataOut.timeZone
1371 1370 self.basicHeaderObj.dstFlag = self.dataOut.dstFlag
1372 1371 self.basicHeaderObj.errorCount = self.dataOut.errorCount
1373 1372
1374 1373 def setFirstHeader(self):
1375 1374 """
1376 1375 Obtiene una copia del First Header
1377 1376
1378 1377 Affected:
1379 1378
1380 1379 self.basicHeaderObj
1381 1380 self.systemHeaderObj
1382 1381 self.radarControllerHeaderObj
1383 1382 self.processingHeaderObj self.
1384 1383
1385 1384 Return:
1386 1385 None
1387 1386 """
1388 1387
1389 1388 raise NotImplementedError
1390 1389
1391 1390 def __writeFirstHeader(self):
1392 1391 """
1393 1392 Escribe el primer header del file es decir el Basic header y el Long header (SystemHeader, RadarControllerHeader, ProcessingHeader)
1394 1393
1395 1394 Affected:
1396 1395 __dataType
1397 1396
1398 1397 Return:
1399 1398 None
1400 1399 """
1401 1400
1402 1401 # CALCULAR PARAMETROS
1403 1402
1404 1403 sizeLongHeader = self.systemHeaderObj.size + \
1405 1404 self.radarControllerHeaderObj.size + self.processingHeaderObj.size
1406 1405 self.basicHeaderObj.size = self.basicHeaderSize + sizeLongHeader
1407 1406
1408 1407 self.basicHeaderObj.write(self.fp)
1409 1408 self.systemHeaderObj.write(self.fp)
1410 1409 self.radarControllerHeaderObj.write(self.fp)
1411 1410 self.processingHeaderObj.write(self.fp)
1412 1411
1413 1412 def __setNewBlock(self):
1414 1413 """
1415 1414 Si es un nuevo file escribe el First Header caso contrario escribe solo el Basic Header
1416 1415
1417 1416 Return:
1418 1417 0 : si no pudo escribir nada
1419 1418 1 : Si escribio el Basic el First Header
1420 1419 """
1421 1420 if self.fp == None:
1422 1421 self.setNextFile()
1423 1422
1424 1423 if self.flagIsNewFile:
1425 1424 return 1
1426 1425
1427 1426 if self.blockIndex < self.processingHeaderObj.dataBlocksPerFile:
1428 1427 self.basicHeaderObj.write(self.fp)
1429 1428 return 1
1430 1429
1431 1430 if not(self.setNextFile()):
1432 1431 return 0
1433 1432
1434 1433 return 1
1435 1434
1436 1435 def writeNextBlock(self):
1437 1436 """
1438 1437 Selecciona el bloque siguiente de datos y los escribe en un file
1439 1438
1440 1439 Return:
1441 1440 0 : Si no hizo pudo escribir el bloque de datos
1442 1441 1 : Si no pudo escribir el bloque de datos
1443 1442 """
1444 1443 if not(self.__setNewBlock()):
1445 1444 return 0
1446 1445
1447 1446 self.writeBlock()
1448 1447
1449 1448 print("[Writing] Block No. %d/%d" % (self.blockIndex,
1450 1449 self.processingHeaderObj.dataBlocksPerFile))
1451 1450
1452 1451 return 1
1453 1452
1454 1453 def setNextFile(self):
1455 1454 """Determina el siguiente file que sera escrito
1456 1455
1457 1456 Affected:
1458 1457 self.filename
1459 1458 self.subfolder
1460 1459 self.fp
1461 1460 self.setFile
1462 1461 self.flagIsNewFile
1463 1462
1464 1463 Return:
1465 1464 0 : Si el archivo no puede ser escrito
1466 1465 1 : Si el archivo esta listo para ser escrito
1467 1466 """
1468 1467 ext = self.ext
1469 1468 path = self.path
1470 1469
1471 1470 if self.fp != None:
1472 1471 self.fp.close()
1473 1472
1474 1473 if not os.path.exists(path):
1475 1474 os.mkdir(path)
1476 1475
1477 1476 timeTuple = time.localtime(self.dataOut.utctime)
1478 1477 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year, timeTuple.tm_yday)
1479 1478
1480 1479 fullpath = os.path.join(path, subfolder)
1481 1480 setFile = self.setFile
1482 1481
1483 1482 if not(os.path.exists(fullpath)):
1484 1483 os.mkdir(fullpath)
1485 1484 setFile = -1 # inicializo mi contador de seteo
1486 1485 else:
1487 1486 filesList = os.listdir(fullpath)
1488 1487 if len(filesList) > 0:
1489 1488 filesList = sorted(filesList, key=str.lower)
1490 1489 filen = filesList[-1]
1491 1490 # el filename debera tener el siguiente formato
1492 1491 # 0 1234 567 89A BCDE (hex)
1493 1492 # x YYYY DDD SSS .ext
1494 1493 if isNumber(filen[8:11]):
1495 1494 # inicializo mi contador de seteo al seteo del ultimo file
1496 1495 setFile = int(filen[8:11])
1497 1496 else:
1498 1497 setFile = -1
1499 1498 else:
1500 1499 setFile = -1 # inicializo mi contador de seteo
1501 1500
1502 1501 setFile += 1
1503 1502
1504 1503 # If this is a new day it resets some values
1505 1504 if self.dataOut.datatime.date() > self.fileDate:
1506 1505 setFile = 0
1507 1506 self.nTotalBlocks = 0
1508 1507
1509 1508 filen = '{}{:04d}{:03d}{:03d}{}'.format(
1510 1509 self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
1511 1510
1512 1511 filename = os.path.join(path, subfolder, filen)
1513 1512
1514 1513 fp = open(filename, 'wb')
1515 1514
1516 1515 self.blockIndex = 0
1517 1516 self.filename = filename
1518 1517 self.subfolder = subfolder
1519 1518 self.fp = fp
1520 1519 self.setFile = setFile
1521 1520 self.flagIsNewFile = 1
1522 1521 self.fileDate = self.dataOut.datatime.date()
1523 1522 self.setFirstHeader()
1524 1523
1525 1524 print('[Writing] Opening file: %s' % self.filename)
1526 1525
1527 1526 self.__writeFirstHeader()
1528 1527
1529 1528 return 1
1530 1529
1531 1530 def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4):
1532 1531 """
1533 1532 Setea el tipo de formato en la cual sera guardada la data y escribe el First Header
1534 1533
1535 1534 Inputs:
1536 1535 path : directory where data will be saved
1537 1536 profilesPerBlock : number of profiles per block
1538 1537 set : initial file set
1539 1538 datatype : An integer number that defines data type:
1540 1539 0 : int8 (1 byte)
1541 1540 1 : int16 (2 bytes)
1542 1541 2 : int32 (4 bytes)
1543 1542 3 : int64 (8 bytes)
1544 1543 4 : float32 (4 bytes)
1545 1544 5 : double64 (8 bytes)
1546 1545
1547 1546 Return:
1548 1547 0 : Si no realizo un buen seteo
1549 1548 1 : Si realizo un buen seteo
1550 1549 """
1551 1550
1552 1551 if ext == None:
1553 1552 ext = self.ext
1554 1553
1555 1554 self.ext = ext.lower()
1556 1555
1557 1556 self.path = path
1558 1557
1559 1558 if set is None:
1560 1559 self.setFile = -1
1561 1560 else:
1562 1561 self.setFile = set - 1
1563 1562
1564 1563 self.blocksPerFile = blocksPerFile
1565 1564 self.profilesPerBlock = profilesPerBlock
1566 1565 self.dataOut = dataOut
1567 1566 self.fileDate = self.dataOut.datatime.date()
1568 1567 self.dtype = self.dataOut.dtype
1569 1568
1570 1569 if datatype is not None:
1571 1570 self.dtype = get_numpy_dtype(datatype)
1572 1571
1573 1572 if not(self.setNextFile()):
1574 1573 print("[Writing] There isn't a next file")
1575 1574 return 0
1576 1575
1577 1576 self.setBlockDimension()
1578 1577
1579 1578 return 1
1580 1579
1581 1580 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1582 1581
1583 1582 if not(self.isConfig):
1584 1583
1585 1584 self.setup(dataOut, path, blocksPerFile, profilesPerBlock=profilesPerBlock,
1586 1585 set=set, ext=ext, datatype=datatype, **kwargs)
1587 1586 self.isConfig = True
1588 1587
1589 1588 self.dataOut = dataOut
1590 1589 self.putData()
1591 1590 return self.dataOut
1592 1591
1593 1592 @MPDecorator
1594 1593 class printInfo(Operation):
1595 1594
1596 1595 def __init__(self):
1597 1596
1598 1597 Operation.__init__(self)
1599 1598 self.__printInfo = True
1600 1599
1601 1600 def run(self, dataOut, headers = ['systemHeaderObj', 'radarControllerHeaderObj', 'processingHeaderObj']):
1602 1601 if self.__printInfo == False:
1603 1602 return
1604 1603
1605 1604 for header in headers:
1606 1605 if hasattr(dataOut, header):
1607 1606 obj = getattr(dataOut, header)
1608 1607 if hasattr(obj, 'printInfo'):
1609 1608 obj.printInfo()
1610 1609 else:
1611 1610 print(obj)
1612 1611 else:
1613 1612 log.warning('Header {} Not found in object'.format(header))
1614 1613
1615 1614 self.__printInfo = False
Show file before | Show file after | Hide comments
@@ -1,630 +1,648
1 1 '''
2 2 Created on Aug 1, 2017
3 3
4 4 @author: Juan C. Espinoza
5 5 '''
6 6
7 7 import os
8 8 import sys
9 9 import time
10 10 import json
11 11 import glob
12 12 import datetime
13 13
14 14 import numpy
15 15 import h5py
16 16
17 17 import schainpy.admin
18 18 from schainpy.model.io.jroIO_base import LOCALTIME, Reader
19 19 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
20 20 from schainpy.model.data.jrodata import Parameters
21 21 from schainpy.utils import log
22 22
23 23 try:
24 24 import madrigal.cedar
25 25 except:
26 26 pass
27 27
28 28 try:
29 29 basestring
30 30 except:
31 31 basestring = str
32 32
33 33 DEF_CATALOG = {
34 34 'principleInvestigator': 'Marco Milla',
35 35 'expPurpose': '',
36 36 'cycleTime': '',
37 37 'correlativeExp': '',
38 38 'sciRemarks': '',
39 39 'instRemarks': ''
40 40 }
41 41
42 42 DEF_HEADER = {
43 43 'kindatDesc': '',
44 44 'analyst': 'Jicamarca User',
45 45 'comments': '',
46 46 'history': ''
47 47 }
48 48
49 49 MNEMONICS = {
50 50 10: 'jro',
51 51 11: 'jbr',
52 52 840: 'jul',
53 53 13: 'jas',
54 54 1000: 'pbr',
55 55 1001: 'hbr',
56 56 1002: 'obr',
57 57 400: 'clr'
58 58
59 59 }
60 60
61 61 UT1970 = datetime.datetime(1970, 1, 1) - datetime.timedelta(seconds=time.timezone)
62 62
63 63 def load_json(obj):
64 64 '''
65 65 Parse json as string instead of unicode
66 66 '''
67 67
68 68 if isinstance(obj, str):
69 69 iterable = json.loads(obj)
70 70 else:
71 71 iterable = obj
72 72
73 73 if isinstance(iterable, dict):
74 74 return {str(k): load_json(v) if isinstance(v, dict) else str(v) if isinstance(v, basestring) else v
75 75 for k, v in list(iterable.items())}
76 76 elif isinstance(iterable, (list, tuple)):
77 77 return [str(v) if isinstance(v, basestring) else v for v in iterable]
78 78
79 79 return iterable
80 80
81 81
82 82 class MADReader(Reader, ProcessingUnit):
83 83
84 84 def __init__(self):
85 85
86 86 ProcessingUnit.__init__(self)
87 87
88 88 self.dataOut = Parameters()
89 89 self.counter_records = 0
90 90 self.nrecords = None
91 91 self.flagNoMoreFiles = 0
92 92 self.filename = None
93 93 self.intervals = set()
94 94 self.datatime = datetime.datetime(1900,1,1)
95 95 self.format = None
96 96 self.filefmt = "***%Y%m%d*******"
97 97
98 98 def setup(self, **kwargs):
99 99
100 100 self.set_kwargs(**kwargs)
101 101 self.oneDDict = load_json(self.oneDDict)
102 102 self.twoDDict = load_json(self.twoDDict)
103 103 self.ind2DList = load_json(self.ind2DList)
104 104 self.independentParam = self.ind2DList[0]
105 105
106 106 if self.path is None:
107 107 raise ValueError('The path is not valid')
108 108
109 109 self.open_file = open
110 110 self.open_mode = 'rb'
111 111
112 112 if self.format is None:
113 113 raise ValueError('The format is not valid choose simple or hdf5')
114 114 elif self.format.lower() in ('simple', 'txt'):
115 115 self.ext = '.txt'
116 116 elif self.format.lower() in ('cedar',):
117 117 self.ext = '.001'
118 118 else:
119 119 self.ext = '.hdf5'
120 120 self.open_file = h5py.File
121 121 self.open_mode = 'r'
122 122
123 123 if self.online:
124 124 log.log("Searching files in online mode...", self.name)
125 125
126 126 for nTries in range(self.nTries):
127 127 fullpath = self.searchFilesOnLine(self.path, self.startDate,
128 128 self.endDate, self.expLabel, self.ext, self.walk,
129 129 self.filefmt, self.folderfmt)
130 130
131 131 try:
132 132 fullpath = next(fullpath)
133 133 except:
134 134 fullpath = None
135 135
136 136 if fullpath:
137 137 break
138 138
139 139 log.warning(
140 140 'Waiting {} sec for a valid file in {}: try {} ...'.format(
141 141 self.delay, self.path, nTries + 1),
142 142 self.name)
143 143 time.sleep(self.delay)
144 144
145 145 if not(fullpath):
146 146 raise schainpy.admin.SchainError(
147 147 'There isn\'t any valid file in {}'.format(self.path))
148 148
149 149 else:
150 150 log.log("Searching files in {}".format(self.path), self.name)
151 151 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
152 152 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
153 153
154 154 self.setNextFile()
155 155
156 156 def readFirstHeader(self):
157 157 '''Read header and data'''
158 158
159 159 self.parseHeader()
160 160 self.parseData()
161 161 self.blockIndex = 0
162 162
163 163 return
164 164
165 165 def parseHeader(self):
166 166 '''
167 167 '''
168 168
169 169 self.output = {}
170 170 self.version = '2'
171 171 s_parameters = None
172 172 if self.ext == '.txt':
173 173 self.parameters = [s.strip().lower() for s in self.fp.readline().decode().strip().split(' ') if s]
174 174 elif self.ext == '.hdf5':
175 175 self.metadata = self.fp['Metadata']
176 176 if '_record_layout' in self.metadata:
177 177 s_parameters = [s[0].lower().decode() for s in self.metadata['Independent Spatial Parameters']]
178 178 self.version = '3'
179 179 self.parameters = [s[0].lower().decode() for s in self.metadata['Data Parameters']]
180 180
181 181 log.success('Parameters found: {}'.format(self.parameters),
182 182 'MADReader')
183 183 if s_parameters:
184 184 log.success('Spatial parameters found: {}'.format(s_parameters),
185 185 'MADReader')
186 186
187 187 for param in list(self.oneDDict.keys()):
188 188 if param.lower() not in self.parameters:
189 189 log.warning(
190 190 'Parameter {} not found will be ignored'.format(
191 191 param),
192 192 'MADReader')
193 193 self.oneDDict.pop(param, None)
194 194
195 195 for param, value in list(self.twoDDict.items()):
196 196 if param.lower() not in self.parameters:
197 197 log.warning(
198 198 'Parameter {} not found, it will be ignored'.format(
199 199 param),
200 200 'MADReader')
201 201 self.twoDDict.pop(param, None)
202 202 continue
203 203 if isinstance(value, list):
204 204 if value[0] not in self.output:
205 205 self.output[value[0]] = []
206 206 self.output[value[0]].append([])
207 207
208 208 def parseData(self):
209 209 '''
210 210 '''
211 211
212 212 if self.ext == '.txt':
213 213 self.data = numpy.genfromtxt(self.fp, missing_values=('missing'))
214 214 self.nrecords = self.data.shape[0]
215 215 self.ranges = numpy.unique(self.data[:,self.parameters.index(self.independentParam.lower())])
216 216 self.counter_records = 0
217 217 elif self.ext == '.hdf5':
218 218 self.data = self.fp['Data']
219 219 self.ranges = numpy.unique(self.data['Table Layout'][self.independentParam.lower()])
220 220 self.times = numpy.unique(self.data['Table Layout']['ut1_unix'])
221 221 self.counter_records = int(self.data['Table Layout']['recno'][0])
222 222 self.nrecords = int(self.data['Table Layout']['recno'][-1])
223 223
224 224 def readNextBlock(self):
225 225
226 226 while True:
227 227 self.flagDiscontinuousBlock = 0
228 228 if self.counter_records == self.nrecords:
229 229 self.setNextFile()
230 230
231 231 self.readBlock()
232 232
233 233 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
234 234 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
235 235 log.warning(
236 236 'Reading Record No. {}/{} -> {} [Skipping]'.format(
237 237 self.counter_records,
238 238 self.nrecords,
239 239 self.datatime.ctime()),
240 240 'MADReader')
241 241 continue
242 242 break
243 243
244 244 log.log(
245 245 'Reading Record No. {}/{} -> {}'.format(
246 246 self.counter_records,
247 247 self.nrecords,
248 248 self.datatime.ctime()),
249 249 'MADReader')
250 250
251 251 return 1
252 252
253 253 def readBlock(self):
254 254 '''
255 255 '''
256 256 dum = []
257 257 if self.ext == '.txt':
258 258 dt = self.data[self.counter_records][:6].astype(int)
259 259 if datetime.datetime(dt[0], dt[1], dt[2], dt[3], dt[4], dt[5]).date() > self.datatime.date():
260 260 self.flagDiscontinuousBlock = 1
261 261 self.datatime = datetime.datetime(dt[0], dt[1], dt[2], dt[3], dt[4], dt[5])
262 262 while True:
263 263 dt = self.data[self.counter_records][:6].astype(int)
264 264 datatime = datetime.datetime(dt[0], dt[1], dt[2], dt[3], dt[4], dt[5])
265 265 if datatime == self.datatime:
266 266 dum.append(self.data[self.counter_records])
267 267 self.counter_records += 1
268 268 if self.counter_records == self.nrecords:
269 269 break
270 270 continue
271 271 self.intervals.add((datatime-self.datatime).seconds)
272 272 break
273 273 elif self.ext == '.hdf5':
274 274 datatime = datetime.datetime.utcfromtimestamp(
275 275 self.times[self.counter_records])
276 276 dum = self.data['Table Layout'][self.data['Table Layout']['recno']==self.counter_records]
277 277 self.intervals.add((datatime-self.datatime).seconds)
278 278 if datatime.date()>self.datatime.date():
279 279 self.flagDiscontinuousBlock = 1
280 280 self.datatime = datatime
281 281 self.counter_records += 1
282 282
283 283 self.buffer = numpy.array(dum)
284 284 return
285 285
286 286 def set_output(self):
287 287 '''
288 288 Storing data from buffer to dataOut object
289 289 '''
290 290
291 291 parameters = [None for __ in self.parameters]
292 292
293 293 for param, attr in list(self.oneDDict.items()):
294 294 x = self.parameters.index(param.lower())
295 295 setattr(self.dataOut, attr, self.buffer[0][x])
296 296
297 297 for param, value in list(self.twoDDict.items()):
298 298 dummy = numpy.zeros(self.ranges.shape) + numpy.nan
299 299 if self.ext == '.txt':
300 300 x = self.parameters.index(param.lower())
301 301 y = self.parameters.index(self.independentParam.lower())
302 302 ranges = self.buffer[:,y]
303 303 #if self.ranges.size == ranges.size:
304 304 # continue
305 305 index = numpy.where(numpy.in1d(self.ranges, ranges))[0]
306 306 dummy[index] = self.buffer[:,x]
307 307 else:
308 308 ranges = self.buffer[self.independentParam.lower()]
309 309 index = numpy.where(numpy.in1d(self.ranges, ranges))[0]
310 310 dummy[index] = self.buffer[param.lower()]
311 311
312 312 if isinstance(value, str):
313 313 if value not in self.independentParam:
314 314 setattr(self.dataOut, value, dummy.reshape(1,-1))
315 315 elif isinstance(value, list):
316 316 self.output[value[0]][value[1]] = dummy
317 317 parameters[value[1]] = param
318 318 for key, value in list(self.output.items()):
319 319 setattr(self.dataOut, key, numpy.array(value))
320 320
321 321 self.dataOut.parameters = [s for s in parameters if s]
322 322 self.dataOut.heightList = self.ranges
323 323 self.dataOut.utctime = (self.datatime - datetime.datetime(1970, 1, 1)).total_seconds()
324 324 self.dataOut.utctimeInit = self.dataOut.utctime
325 325 self.dataOut.paramInterval = min(self.intervals)
326 326 self.dataOut.useLocalTime = False
327 327 self.dataOut.flagNoData = False
328 328 self.dataOut.nrecords = self.nrecords
329 329 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
330 330
331 331 def getData(self):
332 332 '''
333 333 Storing data from databuffer to dataOut object
334 334 '''
335 335
336 336 if not self.readNextBlock():
337 337 self.dataOut.flagNoData = True
338 338 return 0
339 339
340 340 self.set_output()
341 341
342 342 return 1
343 343
344 344 def run(self, **kwargs):
345 345
346 346 if not(self.isConfig):
347 347 self.setup(**kwargs)
348 348 self.isConfig = True
349 349
350 350 self.getData()
351 351
352 352 return
353 353
354 354 @MPDecorator
355 355 class MADWriter(Operation):
356 356 '''Writing module for Madrigal files
357 357
358 358 type: external
359 359
360 360 Inputs:
361 361 path path where files will be created
362 362 oneDDict json of one-dimensional parameters in record where keys
363 363 are Madrigal codes (integers or mnemonics) and values the corresponding
364 364 dataOut attribute e.g: {
365 365 'gdlatr': 'lat',
366 366 'gdlonr': 'lon',
367 367 'gdlat2':'lat',
368 368 'glon2':'lon'}
369 369 ind2DList list of independent spatial two-dimensional parameters e.g:
370 370 ['heigthList']
371 371 twoDDict json of two-dimensional parameters in record where keys
372 372 are Madrigal codes (integers or mnemonics) and values the corresponding
373 373 dataOut attribute if multidimensional array specify as tupple
374 374 ('attr', pos) e.g: {
375 375 'gdalt': 'heightList',
376 376 'vn1p2': ('data_output', 0),
377 377 'vn2p2': ('data_output', 1),
378 378 'vn3': ('data_output', 2),
379 379 'snl': ('data_SNR', 'db')
380 380 }
381 381 metadata json of madrigal metadata (kinst, kindat, catalog and header)
382 382 format hdf5, cedar
383 383 blocks number of blocks per file'''
384 384
385 385 __attrs__ = ['path', 'oneDDict', 'ind2DList', 'twoDDict','metadata', 'format', 'blocks']
386 386 missing = -32767
387 currentDay = None
387 388
388 389 def __init__(self):
389 390
390 391 Operation.__init__(self)
391 392 self.dataOut = Parameters()
392 393 self.counter = 0
393 394 self.path = None
394 395 self.fp = None
395 396
396 397 def run(self, dataOut, path, oneDDict, ind2DList='[]', twoDDict='{}',
397 398 metadata='{}', format='cedar', **kwargs):
398 399
399 400
400 401 #if dataOut.AUX==1: #Modified
401 402
402 403 if not self.isConfig:
403 404 self.setup(path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs)
404 405 self.isConfig = True
405 406
406 407 self.dataOut = dataOut
407 408 self.putData()
408 409
409 410 return 1
410 411
411 412 def setup(self, path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs):
412 413 '''
413 414 Configure Operation
414 415 '''
415 416
416 417 self.path = path
417 418 self.blocks = kwargs.get('blocks', None)
418 419 self.counter = 0
419 420 self.oneDDict = load_json(oneDDict)
420 421 self.twoDDict = load_json(twoDDict)
421 422 self.ind2DList = load_json(ind2DList)
422 423 meta = load_json(metadata)
423 424 self.kinst = meta.get('kinst')
424 425 self.kindat = meta.get('kindat')
425 426 self.catalog = meta.get('catalog', DEF_CATALOG)
426 427 self.header = meta.get('header', DEF_HEADER)
427 428 if format == 'cedar':
428 429 self.ext = '.dat'
429 430 self.extra_args = {}
430 431 elif format == 'hdf5':
431 432 self.ext = '.hdf5'
432 433 self.extra_args = {'ind2DList': self.ind2DList}
433 434
434 435 self.keys = [k.lower() for k in self.twoDDict]
435 436 if 'range' in self.keys:
436 437 self.keys.remove('range')
437 438 if 'gdalt' in self.keys:
438 439 self.keys.remove('gdalt')
439 440
440 441 def setFile(self):
441 442 '''
442 443 Create new cedar file object
443 444 '''
444 445
445 446 self.mnemonic = MNEMONICS[self.kinst] #TODO get mnemonic from madrigal
446 447 date = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
447 448 #if self.dataOut.input_dat_type:
448 449 #date=datetime.datetime.fromtimestamp(self.dataOut.TimeBlockSeconds_for_dp_power)
449 450 #print("date",date)
450 451
451 452 filename = '{}{}{}'.format(self.mnemonic,
452 453 date.strftime('%Y%m%d_%H%M%S'),
453 454 self.ext)
454 455
455 456 self.fullname = os.path.join(self.path, filename)
456 457
457 458 if os.path.isfile(self.fullname) :
458 459 log.warning(
459 460 'Destination file {} already exists, previous file deleted.'.format(
460 461 self.fullname),
461 462 'MADWriter')
462 463 os.remove(self.fullname)
463 464
464 465 try:
465 466 log.success(
466 467 'Creating file: {}'.format(self.fullname),
467 468 'MADWriter')
468 469 if not os.path.exists(self.path):
469 470 os.makedirs(self.path)
470 471 self.fp = madrigal.cedar.MadrigalCedarFile(self.fullname, True)
471 472
472 473
473 474 except ValueError as e:
474 475 log.error(
475 476 'Impossible to create a cedar object with "madrigal.cedar.MadrigalCedarFile"',
476 477 'MADWriter')
477 478 return
478 479
479 480 return 1
480 481
481 482 def writeBlock(self):
482 483 '''
483 484 Add data records to cedar file taking data from oneDDict and twoDDict
484 485 attributes.
485 486 Allowed parameters in: parcodes.tab
486 487 '''
487 488 #self.dataOut.paramInterval=2
488 489 startTime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
489 490
490 491 endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
491 492
492 493 #if self.dataOut.input_dat_type:
493 494 #if self.dataOut.experiment=="DP":
494 495 #startTime=datetime.datetime.fromtimestamp(self.dataOut.TimeBlockSeconds_for_dp_power)
495 496 #endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
496 497
497 498
498 499 #print("2: ",startTime)
499 500 #print(endTime)
500 501 heights = self.dataOut.heightList
501 502 #print(heights)
502 503 #exit(1)
503 504 #print(self.blocks)
504 505 #print(startTime)
505 506 #print(endTime)
506 507 #print(heights)
507 508 #input()
508 509 if self.ext == '.dat':
509 510 for key, value in list(self.twoDDict.items()):
510 511 if isinstance(value, str):
511 512 data = getattr(self.dataOut, value)
512 513 invalid = numpy.isnan(data)
513 514 data[invalid] = self.missing
514 515 elif isinstance(value, (tuple, list)):
515 516 attr, key = value
516 517 data = getattr(self.dataOut, attr)
517 518 invalid = numpy.isnan(data)
518 519 data[invalid] = self.missing
519 520
520 521 out = {}
521 522 for key, value in list(self.twoDDict.items()):
522 523 key = key.lower()
523 524 if isinstance(value, str):
524 525 if 'db' in value.lower():
525 526 tmp = getattr(self.dataOut, value.replace('_db', ''))
526 527 SNRavg = numpy.average(tmp, axis=0)
527 528 tmp = 10*numpy.log10(SNRavg)
528 529 else:
529 530 tmp = getattr(self.dataOut, value)
530 531 out[key] = tmp.flatten()[:len(heights)]
531 532 elif isinstance(value, (tuple, list)):
532 533 attr, x = value
533 534 data = getattr(self.dataOut, attr)
534 535 #print(x)
535 536 #print(len(heights))
536 537 #print(data[int(x)][:len(heights)])
537 538 #print(numpy.shape(out))
538 539 #print(numpy.shape(data))
539 540
540 541 out[key] = data[int(x)][:len(heights)]
541 542
542 543 a = numpy.array([out[k] for k in self.keys])
543 544 #print(a)
544 545 nrows = numpy.array([numpy.isnan(a[:, x]).all() for x in range(len(heights))])
545 546 index = numpy.where(nrows == False)[0]
546 547
547 548 #print(startTime.minute)
548 549 rec = madrigal.cedar.MadrigalDataRecord(
549 550 self.kinst,
550 551 self.kindat,
551 552 startTime.year,
552 553 startTime.month,
553 554 startTime.day,
554 555 startTime.hour,
555 556 startTime.minute,
556 557 startTime.second,
557 558 startTime.microsecond/10000,
558 559 endTime.year,
559 560 endTime.month,
560 561 endTime.day,
561 562 endTime.hour,
562 563 endTime.minute,
563 564 endTime.second,
564 565 endTime.microsecond/10000,
565 566 list(self.oneDDict.keys()),
566 567 list(self.twoDDict.keys()),
567 568 len(index),
568 569 **self.extra_args
569 570 )
570 571 #print("rec",rec)
571 572 # Setting 1d values
572 573 for key in self.oneDDict:
573 574 rec.set1D(key, getattr(self.dataOut, self.oneDDict[key]))
574 575
575 576 # Setting 2d values
576 577 nrec = 0
577 578 for n in index:
578 579 for key in out:
579 580 rec.set2D(key, nrec, out[key][n])
580 581 nrec += 1
581 582
582 583 self.fp.append(rec)
583 584 if self.ext == '.hdf5' and self.counter %2 == 0 and self.counter > 0:
584 585 #print("here")
585 586 self.fp.dump()
586 587 if self.counter % 20 == 0 and self.counter > 0:
587 588 #self.fp.write()
588 589 log.log(
589 590 'Writing {} records'.format(
590 591 self.counter),
591 592 'MADWriter')
592 593
593 594 def setHeader(self):
594 595 '''
595 596 Create an add catalog and header to cedar file
596 597 '''
597 598
598 599 log.success('Closing file {}'.format(self.fullname), 'MADWriter')
599 600
600 601 if self.ext == '.dat':
601 602 self.fp.write()
602 603 else:
603 604 self.fp.dump()
604 605 self.fp.close()
605 606
606 607 header = madrigal.cedar.CatalogHeaderCreator(self.fullname)
607 608 header.createCatalog(**self.catalog)
608 609 header.createHeader(**self.header)
609 610 header.write()
610 611
612 def timeFlag(self):
613 currentTime = self.dataOut.utctime
614 timeTuple = time.localtime(currentTime)
615 dataDay = timeTuple.tm_yday
616
617 if self.currentDay is None:
618 self.currentDay = dataDay
619 return False
620
621 #Si el dia es diferente
622 if dataDay != self.currentDay:
623 self.currentDay = dataDay
624 return True
625
626 else:
627 return False
628
611 629 def putData(self):
612 630
613 631 if self.dataOut.flagNoData:
614 632 return 0
615 633
616 if self.dataOut.flagDiscontinuousBlock or self.counter == self.blocks:
634 if self.dataOut.flagDiscontinuousBlock or self.counter == self.blocks or self.timeFlag():
617 635 if self.counter > 0:
618 636 self.setHeader()
619 637 self.counter = 0
620 638
621 639 if self.counter == 0:
622 640 self.setFile()
623 641
624 642 self.writeBlock()
625 643 self.counter += 1
626 644
627 645 def close(self):
628 646
629 647 if self.counter > 0:
630 648 self.setHeader()
Show file before | Show file after | Hide comments
@@ -1,527 +1,527
1 1 '''
2 2 Created on Jul 2, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import numpy
7 7
8 8 from schainpy.model.io.jroIO_base import LOCALTIME, JRODataReader, JRODataWriter
9 9 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
10 10 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
11 11 from schainpy.model.data.jrodata import Spectra
12 12 from schainpy.utils import log
13 13
14 14
15 15 class SpectraReader(JRODataReader, ProcessingUnit):
16 16 """
17 17 Esta clase permite leer datos de espectros desde archivos procesados (.pdata). La lectura
18 18 de los datos siempre se realiza por bloques. Los datos leidos (array de 3 dimensiones)
19 19 son almacenados en tres buffer's para el Self Spectra, el Cross Spectra y el DC Channel.
20 20
21 21 paresCanalesIguales * alturas * perfiles (Self Spectra)
22 22 paresCanalesDiferentes * alturas * perfiles (Cross Spectra)
23 23 canales * alturas (DC Channels)
24 24
25 25 Esta clase contiene instancias (objetos) de las clases BasicHeader, SystemHeader,
26 26 RadarControllerHeader y Spectra. Los tres primeros se usan para almacenar informacion de la
27 27 cabecera de datos (metadata), y el cuarto (Spectra) para obtener y almacenar un bloque de
28 28 datos desde el "buffer" cada vez que se ejecute el metodo "getData".
29 29
30 30 Example:
31 31 dpath = "/home/myuser/data"
32 32
33 33 startTime = datetime.datetime(2010,1,20,0,0,0,0,0,0)
34 34
35 35 endTime = datetime.datetime(2010,1,21,23,59,59,0,0,0)
36 36
37 37 readerObj = SpectraReader()
38 38
39 39 readerObj.setup(dpath, startTime, endTime)
40 40
41 41 while(True):
42 42
43 43 readerObj.getData()
44 44
45 45 print readerObj.data_spc
46 46
47 47 print readerObj.data_cspc
48 48
49 49 print readerObj.data_dc
50 50
51 51 if readerObj.flagNoMoreFiles:
52 52 break
53 53
54 54 """
55 55
56 56 def __init__(self):#, **kwargs):
57 57 """
58 58 Inicializador de la clase SpectraReader para la lectura de datos de espectros.
59 59
60 60 Inputs:
61 61 dataOut : Objeto de la clase Spectra. Este objeto sera utilizado para
62 62 almacenar un perfil de datos cada vez que se haga un requerimiento
63 63 (getData). El perfil sera obtenido a partir del buffer de datos,
64 64 si el buffer esta vacio se hara un nuevo proceso de lectura de un
65 65 bloque de datos.
66 66 Si este parametro no es pasado se creara uno internamente.
67 67
68 68 Affected:
69 69 self.dataOut
70 70
71 71 Return : None
72 72 """
73 73
74 74 ProcessingUnit.__init__(self)
75 75
76 76 self.pts2read_SelfSpectra = 0
77 77 self.pts2read_CrossSpectra = 0
78 self.pts2read_DCchannels = 0
78 self.pts2read_DCchannels = 0
79 79 self.ext = ".pdata"
80 80 self.optchar = "P"
81 81 self.basicHeaderObj = BasicHeader(LOCALTIME)
82 82 self.systemHeaderObj = SystemHeader()
83 83 self.radarControllerHeaderObj = RadarControllerHeader()
84 84 self.processingHeaderObj = ProcessingHeader()
85 85 self.lastUTTime = 0
86 86 self.maxTimeStep = 30
87 87 self.dataOut = Spectra()
88 88 self.profileIndex = 1
89 89 self.nRdChannels = None
90 90 self.nRdPairs = None
91 91 self.rdPairList = []
92 92
93 93 def createObjByDefault(self):
94 94
95 95 dataObj = Spectra()
96 96
97 97 return dataObj
98 98
99 99 def __hasNotDataInBuffer(self):
100 100 return 1
101 101
102 102
103 103 def getBlockDimension(self):
104 104 """
105 105 Obtiene la cantidad de puntos a leer por cada bloque de datos
106 106
107 107 Affected:
108 108 self.nRdChannels
109 109 self.nRdPairs
110 110 self.pts2read_SelfSpectra
111 111 self.pts2read_CrossSpectra
112 112 self.pts2read_DCchannels
113 113 self.blocksize
114 114 self.dataOut.nChannels
115 115 self.dataOut.nPairs
116 116
117 117 Return:
118 118 None
119 119 """
120 120 self.nRdChannels = 0
121 121 self.nRdPairs = 0
122 122 self.rdPairList = []
123 123
124 124 for i in range(0, self.processingHeaderObj.totalSpectra*2, 2):
125 125 if self.processingHeaderObj.spectraComb[i] == self.processingHeaderObj.spectraComb[i+1]:
126 126 self.nRdChannels = self.nRdChannels + 1 #par de canales iguales
127 127 else:
128 128 self.nRdPairs = self.nRdPairs + 1 #par de canales diferentes
129 129 self.rdPairList.append((self.processingHeaderObj.spectraComb[i], self.processingHeaderObj.spectraComb[i+1]))
130 130
131 131 pts2read = self.processingHeaderObj.nHeights * self.processingHeaderObj.profilesPerBlock
132 132
133 133 self.pts2read_SelfSpectra = int(self.nRdChannels * pts2read)
134 134 self.blocksize = self.pts2read_SelfSpectra
135 135
136 136 if self.processingHeaderObj.flag_cspc:
137 137 self.pts2read_CrossSpectra = int(self.nRdPairs * pts2read)
138 138 self.blocksize += self.pts2read_CrossSpectra
139 139
140 140 if self.processingHeaderObj.flag_dc:
141 141 self.pts2read_DCchannels = int(self.systemHeaderObj.nChannels * self.processingHeaderObj.nHeights)
142 142 self.blocksize += self.pts2read_DCchannels
143 143
144 144 def readBlock(self):
145 145 """
146 146 Lee el bloque de datos desde la posicion actual del puntero del archivo
147 147 (self.fp) y actualiza todos los parametros relacionados al bloque de datos
148 148 (metadata + data). La data leida es almacenada en el buffer y el contador del buffer
149 149 es seteado a 0
150 150
151 151 Return: None
152 152
153 153 Variables afectadas:
154 154
155 155 self.flagIsNewFile
156 156 self.flagIsNewBlock
157 157 self.nTotalBlocks
158 158 self.data_spc
159 159 self.data_cspc
160 160 self.data_dc
161 161
162 162 Exceptions:
163 163 Si un bloque leido no es un bloque valido
164 164 """
165
165
166 166 fpointer = self.fp.tell()
167 167
168 168 spc = numpy.fromfile( self.fp, self.dtype[0], self.pts2read_SelfSpectra )
169 169 spc = spc.reshape( (self.nRdChannels, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D
170 170
171 171 if self.processingHeaderObj.flag_cspc:
172 172 cspc = numpy.fromfile( self.fp, self.dtype, self.pts2read_CrossSpectra )
173 173 cspc = cspc.reshape( (self.nRdPairs, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D
174 174
175 175 if self.processingHeaderObj.flag_dc:
176 176 dc = numpy.fromfile( self.fp, self.dtype, self.pts2read_DCchannels ) #int(self.processingHeaderObj.nHeights*self.systemHeaderObj.nChannels) )
177 177 dc = dc.reshape( (self.systemHeaderObj.nChannels, self.processingHeaderObj.nHeights) ) #transforma a un arreglo 2D
178 178
179 179 if not self.processingHeaderObj.shif_fft:
180 180 #desplaza a la derecha en el eje 2 determinadas posiciones
181 181 shift = int(self.processingHeaderObj.profilesPerBlock/2)
182 182 spc = numpy.roll( spc, shift , axis=2 )
183 183
184 184 if self.processingHeaderObj.flag_cspc:
185 185 #desplaza a la derecha en el eje 2 determinadas posiciones
186 186 cspc = numpy.roll( cspc, shift, axis=2 )
187 187
188 188 #Dimensions : nChannels, nProfiles, nSamples
189 189 spc = numpy.transpose( spc, (0,2,1) )
190 190 self.data_spc = spc
191 191
192 192 if self.processingHeaderObj.flag_cspc:
193 193 cspc = numpy.transpose( cspc, (0,2,1) )
194 194 self.data_cspc = cspc['real'] + cspc['imag']*1j
195 195 else:
196 196 self.data_cspc = None
197 197
198 198 if self.processingHeaderObj.flag_dc:
199 199 self.data_dc = dc['real'] + dc['imag']*1j
200 200 else:
201 201 self.data_dc = None
202 202
203 203 self.flagIsNewFile = 0
204 204 self.flagIsNewBlock = 1
205 205
206 206 self.nTotalBlocks += 1
207 207 self.nReadBlocks += 1
208 208
209 209 return 1
210 210
211 211 def getFirstHeader(self):
212 212
213 213 self.getBasicHeader()
214 214 self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
215 215 self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
216 216 self.dataOut.dtype = self.dtype
217 217 self.dataOut.pairsList = self.rdPairList
218 218 self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
219 219 self.dataOut.nFFTPoints = self.processingHeaderObj.profilesPerBlock
220 220 self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
221 221 self.dataOut.nIncohInt = self.processingHeaderObj.nIncohInt
222 222 xf = self.processingHeaderObj.firstHeight + self.processingHeaderObj.nHeights*self.processingHeaderObj.deltaHeight
223 223 self.dataOut.heightList = numpy.arange(self.processingHeaderObj.firstHeight, xf, self.processingHeaderObj.deltaHeight)
224 224 self.dataOut.channelList = list(range(self.systemHeaderObj.nChannels))
225 225 self.dataOut.flagShiftFFT = True #Data is always shifted
226 226 self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode #asumo q la data no esta decodificada
227 227 self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip #asumo q la data esta sin flip
228 228
229 229 def getData(self):
230 230 """
231 231 First method to execute before "RUN" is called.
232 232
233 233 Copia el buffer de lectura a la clase "Spectra",
234 234 con todos los parametros asociados a este (metadata). cuando no hay datos en el buffer de
235 235 lectura es necesario hacer una nueva lectura de los bloques de datos usando "readNextBlock"
236 236
237 237 Return:
238 238 0 : Si no hay mas archivos disponibles
239 239 1 : Si hizo una buena copia del buffer
240 240
241 241 Affected:
242 242 self.dataOut
243 243 self.flagDiscontinuousBlock
244 244 self.flagIsNewBlock
245 245 """
246 246
247 247 if self.flagNoMoreFiles:
248 248 self.dataOut.flagNoData = True
249 249 return 0
250 250
251 251 self.flagDiscontinuousBlock = 0
252 252 self.flagIsNewBlock = 0
253 253
254 254 if self.__hasNotDataInBuffer():
255 255
256 256 if not( self.readNextBlock() ):
257 257 self.dataOut.flagNoData = True
258 258 return 0
259 259
260 260 #data es un numpy array de 3 dmensiones (perfiles, alturas y canales)
261 261
262 262 if self.data_spc is None:
263 263 self.dataOut.flagNoData = True
264 264 return 0
265 265
266 266 self.getBasicHeader()
267 267 self.getFirstHeader()
268 268 self.dataOut.data_spc = self.data_spc
269 269 self.dataOut.data_cspc = self.data_cspc
270 270 self.dataOut.data_dc = self.data_dc
271 271 self.dataOut.flagNoData = False
272 272 self.dataOut.realtime = self.online
273 273
274 274 return self.dataOut.data_spc
275 275
276 276
277 277 @MPDecorator
278 278 class SpectraWriter(JRODataWriter, Operation):
279 279
280 280 """
281 281 Esta clase permite escribir datos de espectros a archivos procesados (.pdata). La escritura
282 282 de los datos siempre se realiza por bloques.
283 283 """
284 284
285 285 def __init__(self):
286 286 """
287 287 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
288 288
289 289 Affected:
290 290 self.dataOut
291 291 self.basicHeaderObj
292 292 self.systemHeaderObj
293 293 self.radarControllerHeaderObj
294 294 self.processingHeaderObj
295 295
296 296 Return: None
297 297 """
298 298
299 299 Operation.__init__(self)
300 300
301 301 self.ext = ".pdata"
302 302 self.optchar = "P"
303 303 self.shape_spc_Buffer = None
304 304 self.shape_cspc_Buffer = None
305 305 self.shape_dc_Buffer = None
306 306 self.data_spc = None
307 307 self.data_cspc = None
308 308 self.data_dc = None
309 309 self.setFile = None
310 310 self.noMoreFiles = 0
311 311 self.basicHeaderObj = BasicHeader(LOCALTIME)
312 312 self.systemHeaderObj = SystemHeader()
313 313 self.radarControllerHeaderObj = RadarControllerHeader()
314 314 self.processingHeaderObj = ProcessingHeader()
315 315
316 316 def hasAllDataInBuffer(self):
317 317 return 1
318 318
319 319
320 320 def setBlockDimension(self):
321 321 """
322 322 Obtiene las formas dimensionales del los subbloques de datos que componen un bloque
323 323
324 324 Affected:
325 325 self.shape_spc_Buffer
326 326 self.shape_cspc_Buffer
327 327 self.shape_dc_Buffer
328 328
329 329 Return: None
330 330 """
331 331 self.shape_spc_Buffer = (self.dataOut.nChannels,
332 332 self.processingHeaderObj.nHeights,
333 333 self.processingHeaderObj.profilesPerBlock)
334 334
335 335 self.shape_cspc_Buffer = (self.dataOut.nPairs,
336 336 self.processingHeaderObj.nHeights,
337 337 self.processingHeaderObj.profilesPerBlock)
338 338
339 339 self.shape_dc_Buffer = (self.dataOut.nChannels,
340 340 self.processingHeaderObj.nHeights)
341 341
342 342
343 343 def writeBlock(self):
344 344 """processingHeaderObj
345 345 Escribe el buffer en el file designado
346 346
347 347 Affected:
348 348 self.data_spc
349 349 self.data_cspc
350 350 self.data_dc
351 351 self.flagIsNewFile
352 352 self.flagIsNewBlock
353 353 self.nTotalBlocks
354 354 self.nWriteBlocks
355 355
356 356 Return: None
357 357 """
358 358
359 359 spc = numpy.transpose( self.data_spc, (0,2,1) )
360 360 if not self.processingHeaderObj.shif_fft:
361 361 spc = numpy.roll( spc, int(self.processingHeaderObj.profilesPerBlock/2), axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
362 362 data = spc.reshape((-1))
363 363 data = data.astype(self.dtype[0])
364 364 data.tofile(self.fp)
365 365
366 366 if self.data_cspc is not None:
367
367
368 368 cspc = numpy.transpose( self.data_cspc, (0,2,1) )
369 369 data = numpy.zeros( numpy.shape(cspc), self.dtype )
370 370 #print 'data.shape', self.shape_cspc_Buffer
371 371 if not self.processingHeaderObj.shif_fft:
372 372 cspc = numpy.roll( cspc, int(self.processingHeaderObj.profilesPerBlock/2), axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
373 373 data['real'] = cspc.real
374 374 data['imag'] = cspc.imag
375 375 data = data.reshape((-1))
376 376 data.tofile(self.fp)
377 377
378 378 if self.data_dc is not None:
379
379
380 380 dc = self.data_dc
381 381 data = numpy.zeros( numpy.shape(dc), self.dtype )
382 382 data['real'] = dc.real
383 383 data['imag'] = dc.imag
384 384 data = data.reshape((-1))
385 385 data.tofile(self.fp)
386 386
387 387 # self.data_spc.fill(0)
388 388 #
389 389 # if self.data_dc is not None:
390 390 # self.data_dc.fill(0)
391 391 #
392 392 # if self.data_cspc is not None:
393 393 # self.data_cspc.fill(0)
394 394
395 395 self.flagIsNewFile = 0
396 396 self.flagIsNewBlock = 1
397 397 self.nTotalBlocks += 1
398 398 self.nWriteBlocks += 1
399 399 self.blockIndex += 1
400 400
401 401 # print "[Writing] Block = %d04" %self.blockIndex
402 402
403 403 def putData(self):
404 404 """
405 405 Setea un bloque de datos y luego los escribe en un file
406 406
407 407 Affected:
408 408 self.data_spc
409 409 self.data_cspc
410 410 self.data_dc
411 411
412 412 Return:
413 413 0 : Si no hay data o no hay mas files que puedan escribirse
414 414 1 : Si se escribio la data de un bloque en un file
415 415 """
416 416
417 417 if self.dataOut.flagNoData:
418 418 return 0
419 419
420 420 self.flagIsNewBlock = 0
421 421
422 422 if self.dataOut.flagDiscontinuousBlock:
423 423 self.data_spc.fill(0)
424 424 if self.dataOut.data_cspc is not None:
425 425 self.data_cspc.fill(0)
426 426 if self.dataOut.data_dc is not None:
427 427 self.data_dc.fill(0)
428 428 self.setNextFile()
429 429
430 430 if self.flagIsNewFile == 0:
431 431 self.setBasicHeader()
432 432
433 433 self.data_spc = self.dataOut.data_spc.copy()
434 434
435 435 if self.dataOut.data_cspc is not None:
436 436 self.data_cspc = self.dataOut.data_cspc.copy()
437 437
438 438 if self.dataOut.data_dc is not None:
439 439 self.data_dc = self.dataOut.data_dc.copy()
440 440
441 441 # #self.processingHeaderObj.dataBlocksPerFile)
442 442 if self.hasAllDataInBuffer():
443 443 # self.setFirstHeader()
444 444 self.writeNextBlock()
445 445
446 446 def __getBlockSize(self):
447 447 '''
448 448 Este metodos determina el cantidad de bytes para un bloque de datos de tipo Spectra
449 449 '''
450 450
451 451 dtype_width = self.getDtypeWidth()
452 452
453 453 pts2write = self.dataOut.nHeights * self.dataOut.nFFTPoints
454 454
455 455 pts2write_SelfSpectra = int(self.dataOut.nChannels * pts2write)
456 456 blocksize = (pts2write_SelfSpectra*dtype_width)
457 457
458 458 if self.dataOut.data_cspc is not None:
459 459 pts2write_CrossSpectra = int(self.dataOut.nPairs * pts2write)
460 460 blocksize += (pts2write_CrossSpectra*dtype_width*2)
461 461
462 462 if self.dataOut.data_dc is not None:
463 463 pts2write_DCchannels = int(self.dataOut.nChannels * self.dataOut.nHeights)
464 464 blocksize += (pts2write_DCchannels*dtype_width*2)
465 465
466 466 # blocksize = blocksize #* datatypeValue * 2 #CORREGIR ESTO
467 467
468 468 return blocksize
469 469
470 470 def setFirstHeader(self):
471 471
472 472 """
473 473 Obtiene una copia del First Header
474 474
475 475 Affected:
476 476 self.systemHeaderObj
477 477 self.radarControllerHeaderObj
478 478 self.dtype
479 479
480 480 Return:
481 481 None
482 482 """
483 483
484 484 self.systemHeaderObj = self.dataOut.systemHeaderObj.copy()
485 485 self.systemHeaderObj.nChannels = self.dataOut.nChannels
486 486 self.radarControllerHeaderObj = self.dataOut.radarControllerHeaderObj.copy()
487 487
488 488 self.processingHeaderObj.dtype = 1 # Spectra
489 489 self.processingHeaderObj.blockSize = self.__getBlockSize()
490 490 self.processingHeaderObj.profilesPerBlock = self.dataOut.nFFTPoints
491 491 self.processingHeaderObj.dataBlocksPerFile = self.blocksPerFile
492 492 self.processingHeaderObj.nWindows = 1 #podria ser 1 o self.dataOut.processingHeaderObj.nWindows
493 493 self.processingHeaderObj.nCohInt = self.dataOut.nCohInt# Se requiere para determinar el valor de timeInterval
494 494 self.processingHeaderObj.nIncohInt = self.dataOut.nIncohInt
495 495 self.processingHeaderObj.totalSpectra = self.dataOut.nPairs + self.dataOut.nChannels
496 496 self.processingHeaderObj.shif_fft = self.dataOut.flagShiftFFT
497 497
498 498 if self.processingHeaderObj.totalSpectra > 0:
499 499 channelList = []
500 500 for channel in range(self.dataOut.nChannels):
501 501 channelList.append(channel)
502 502 channelList.append(channel)
503 503
504 504 pairsList = []
505 505 if self.dataOut.nPairs > 0:
506 506 for pair in self.dataOut.pairsList:
507 507 pairsList.append(pair[0])
508 508 pairsList.append(pair[1])
509 509
510 510 spectraComb = channelList + pairsList
511 511 spectraComb = numpy.array(spectraComb, dtype="u1")
512 512 self.processingHeaderObj.spectraComb = spectraComb
513 513
514 514 if self.dataOut.code is not None:
515 515 self.processingHeaderObj.code = self.dataOut.code
516 516 self.processingHeaderObj.nCode = self.dataOut.nCode
517 517 self.processingHeaderObj.nBaud = self.dataOut.nBaud
518 518
519 519 if self.processingHeaderObj.nWindows != 0:
520 520 self.processingHeaderObj.firstHeight = self.dataOut.heightList[0]
521 521 self.processingHeaderObj.deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
522 522 self.processingHeaderObj.nHeights = self.dataOut.nHeights
523 523 self.processingHeaderObj.samplesWin = self.dataOut.nHeights
524 524
525 525 self.processingHeaderObj.processFlags = self.getProcessFlags()
526 526
527 self.setBasicHeader() No newline at end of file
527 self.setBasicHeader()
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now