##// END OF EJS Templates
schain
RSS
PPI or RHI mode Detection/GeneralPlot
rflores -
r1461:773419e042d4
parent child
Show More
Side by Side Unified
Context file:
r1461:773419e042d4 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,714 +1,716
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Base class to create plot operations
6 6
7 7 """
8 8
9 9 import os
10 10 import sys
11 11 import zmq
12 12 import time
13 13 import numpy
14 14 import datetime
15 15 from collections import deque
16 16 from functools import wraps
17 17 from threading import Thread
18 18 import matplotlib
19 19
20 20 if 'BACKEND' in os.environ:
21 21 matplotlib.use(os.environ['BACKEND'])
22 22 elif 'linux' in sys.platform:
23 23 matplotlib.use("TkAgg")
24 24 elif 'darwin' in sys.platform:
25 25 matplotlib.use('MacOSX')
26 26 else:
27 27 from schainpy.utils import log
28 28 log.warning('Using default Backend="Agg"', 'INFO')
29 29 matplotlib.use('Agg')
30 30
31 31 import matplotlib.pyplot as plt
32 32 from matplotlib.patches import Polygon
33 33 from mpl_toolkits.axes_grid1 import make_axes_locatable
34 34 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
35 35
36 36 from schainpy.model.data.jrodata import PlotterData
37 37 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
38 38 from schainpy.utils import log
39 39
40 40 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
41 41 blu_values = matplotlib.pyplot.get_cmap(
42 42 'seismic_r', 20)(numpy.arange(20))[10:15]
43 43 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
44 44 'jro', numpy.vstack((blu_values, jet_values)))
45 45 matplotlib.pyplot.register_cmap(cmap=ncmap)
46 46
47 47 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
48 48 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
49 49
50 50 EARTH_RADIUS = 6.3710e3
51 51
52 52 def ll2xy(lat1, lon1, lat2, lon2):
53 53
54 54 p = 0.017453292519943295
55 55 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
56 56 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
57 57 r = 12742 * numpy.arcsin(numpy.sqrt(a))
58 58 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
59 59 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
60 60 theta = -theta + numpy.pi/2
61 61 return r*numpy.cos(theta), r*numpy.sin(theta)
62 62
63 63
64 64 def km2deg(km):
65 65 '''
66 66 Convert distance in km to degrees
67 67 '''
68 68
69 69 return numpy.rad2deg(km/EARTH_RADIUS)
70 70
71 71
72 72 def figpause(interval):
73 73 backend = plt.rcParams['backend']
74 74 if backend in matplotlib.rcsetup.interactive_bk:
75 75 figManager = matplotlib._pylab_helpers.Gcf.get_active()
76 76 if figManager is not None:
77 77 canvas = figManager.canvas
78 78 if canvas.figure.stale:
79 79 canvas.draw()
80 80 try:
81 81 canvas.start_event_loop(interval)
82 82 except:
83 83 pass
84 84 return
85 85
86 86 def popup(message):
87 87 '''
88 88 '''
89 89
90 90 fig = plt.figure(figsize=(12, 8), facecolor='r')
91 91 text = '\n'.join([s.strip() for s in message.split(':')])
92 92 fig.text(0.01, 0.5, text, ha='left', va='center',
93 93 size='20', weight='heavy', color='w')
94 94 fig.show()
95 95 figpause(1000)
96 96
97 97
98 98 class Throttle(object):
99 99 '''
100 100 Decorator that prevents a function from being called more than once every
101 101 time period.
102 102 To create a function that cannot be called more than once a minute, but
103 103 will sleep until it can be called:
104 104 @Throttle(minutes=1)
105 105 def foo():
106 106 pass
107 107
108 108 for i in range(10):
109 109 foo()
110 110 print "This function has run %s times." % i
111 111 '''
112 112
113 113 def __init__(self, seconds=0, minutes=0, hours=0):
114 114 self.throttle_period = datetime.timedelta(
115 115 seconds=seconds, minutes=minutes, hours=hours
116 116 )
117 117
118 118 self.time_of_last_call = datetime.datetime.min
119 119
120 120 def __call__(self, fn):
121 121 @wraps(fn)
122 122 def wrapper(*args, **kwargs):
123 123 coerce = kwargs.pop('coerce', None)
124 124 if coerce:
125 125 self.time_of_last_call = datetime.datetime.now()
126 126 return fn(*args, **kwargs)
127 127 else:
128 128 now = datetime.datetime.now()
129 129 time_since_last_call = now - self.time_of_last_call
130 130 time_left = self.throttle_period - time_since_last_call
131 131
132 132 if time_left > datetime.timedelta(seconds=0):
133 133 return
134 134
135 135 self.time_of_last_call = datetime.datetime.now()
136 136 return fn(*args, **kwargs)
137 137
138 138 return wrapper
139 139
140 140 def apply_throttle(value):
141 141
142 142 @Throttle(seconds=value)
143 143 def fnThrottled(fn):
144 144 fn()
145 145
146 146 return fnThrottled
147 147
148 148
149 149 @MPDecorator
150 150 class Plot(Operation):
151 151 """Base class for Schain plotting operations
152 152
153 153 This class should never be use directtly you must subclass a new operation,
154 154 children classes must be defined as follow:
155 155
156 156 ExamplePlot(Plot):
157 157
158 158 CODE = 'code'
159 159 colormap = 'jet'
160 160 plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
161 161
162 162 def setup(self):
163 163 pass
164 164
165 165 def plot(self):
166 166 pass
167 167
168 168 """
169 169
170 170 CODE = 'Figure'
171 171 colormap = 'jet'
172 172 bgcolor = 'white'
173 173 buffering = True
174 174 __missing = 1E30
175 175
176 176 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
177 177 'showprofile']
178 178
179 179 def __init__(self):
180 180
181 181 Operation.__init__(self)
182 182 self.isConfig = False
183 183 self.isPlotConfig = False
184 184 self.save_time = 0
185 185 self.sender_time = 0
186 186 self.data = None
187 187 self.firsttime = True
188 188 self.sender_queue = deque(maxlen=10)
189 189 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
190 190
191 191 def __fmtTime(self, x, pos):
192 192 '''
193 193 '''
194 194
195 195 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
196 196
197 197 def __setup(self, **kwargs):
198 198 '''
199 199 Initialize variables
200 200 '''
201 201
202 202 self.figures = []
203 203 self.axes = []
204 204 self.cb_axes = []
205 205 self.localtime = kwargs.pop('localtime', True)
206 206 self.show = kwargs.get('show', True)
207 207 self.save = kwargs.get('save', False)
208 208 self.save_period = kwargs.get('save_period', 0)
209 209 self.colormap = kwargs.get('colormap', self.colormap)
210 210 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
211 211 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
212 212 self.colormaps = kwargs.get('colormaps', None)
213 213 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
214 214 self.showprofile = kwargs.get('showprofile', False)
215 215 self.title = kwargs.get('wintitle', self.CODE.upper())
216 216 self.cb_label = kwargs.get('cb_label', None)
217 217 self.cb_labels = kwargs.get('cb_labels', None)
218 218 self.labels = kwargs.get('labels', None)
219 219 self.xaxis = kwargs.get('xaxis', 'frequency')
220 220 self.zmin = kwargs.get('zmin', None)
221 221 self.zmax = kwargs.get('zmax', None)
222 222 self.zlimits = kwargs.get('zlimits', None)
223 223 self.xmin = kwargs.get('xmin', None)
224 224 self.xmax = kwargs.get('xmax', None)
225 225 self.xrange = kwargs.get('xrange', 12)
226 226 self.xscale = kwargs.get('xscale', None)
227 227 self.ymin = kwargs.get('ymin', None)
228 228 self.ymax = kwargs.get('ymax', None)
229 229 self.yscale = kwargs.get('yscale', None)
230 230 self.xlabel = kwargs.get('xlabel', None)
231 231 self.attr_time = kwargs.get('attr_time', 'utctime')
232 232 self.attr_data = kwargs.get('attr_data', 'data_param')
233 233 self.decimation = kwargs.get('decimation', None)
234 234 self.oneFigure = kwargs.get('oneFigure', True)
235 235 self.width = kwargs.get('width', None)
236 236 self.height = kwargs.get('height', None)
237 237 self.colorbar = kwargs.get('colorbar', True)
238 238 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
239 239 self.channels = kwargs.get('channels', None)
240 240 self.titles = kwargs.get('titles', [])
241 241 self.polar = False
242 242 self.type = kwargs.get('type', 'iq')
243 243 self.grid = kwargs.get('grid', False)
244 244 self.pause = kwargs.get('pause', False)
245 245 self.save_code = kwargs.get('save_code', self.CODE)
246 246 self.throttle = kwargs.get('throttle', 0)
247 247 self.exp_code = kwargs.get('exp_code', None)
248 248 self.server = kwargs.get('server', False)
249 249 self.sender_period = kwargs.get('sender_period', 60)
250 250 self.tag = kwargs.get('tag', '')
251 251 self.height_index = kwargs.get('height_index', None)
252 252 self.__throttle_plot = apply_throttle(self.throttle)
253 253 code = self.attr_data if self.attr_data else self.CODE
254 254 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
255 255 self.ang_min = kwargs.get('ang_min', None)
256 256 self.ang_max = kwargs.get('ang_max', None)
257 257 self.mode = kwargs.get('mode', None)
258 258
259 259
260
260 261 if self.server:
261 262 if not self.server.startswith('tcp://'):
262 263 self.server = 'tcp://{}'.format(self.server)
263 264 log.success(
264 265 'Sending to server: {}'.format(self.server),
265 266 self.name
266 267 )
267 268
268 269 if isinstance(self.attr_data, str):
269 270 self.attr_data = [self.attr_data]
270 271
271 272 def __setup_plot(self):
272 273 '''
273 274 Common setup for all figures, here figures and axes are created
274 275 '''
275 276
276 277 self.setup()
277 278
278 279 self.time_label = 'LT' if self.localtime else 'UTC'
279 280
280 281 if self.width is None:
281 282 self.width = 8
282 283
283 284 self.figures = []
284 285 self.axes = []
285 286 self.cb_axes = []
286 287 self.pf_axes = []
287 288 self.cmaps = []
288 289
289 290 size = '15%' if self.ncols == 1 else '30%'
290 291 pad = '4%' if self.ncols == 1 else '8%'
291 292
292 293 if self.oneFigure:
293 294 if self.height is None:
294 295 self.height = 1.4 * self.nrows + 1
295 296 fig = plt.figure(figsize=(self.width, self.height),
296 297 edgecolor='k',
297 298 facecolor='w')
298 299 self.figures.append(fig)
299 300 for n in range(self.nplots):
300 301 ax = fig.add_subplot(self.nrows, self.ncols,
301 302 n + 1, polar=self.polar)
302 303 ax.tick_params(labelsize=8)
303 304 ax.firsttime = True
304 305 ax.index = 0
305 306 ax.press = None
306 307 self.axes.append(ax)
307 308 if self.showprofile:
308 309 cax = self.__add_axes(ax, size=size, pad=pad)
309 310 cax.tick_params(labelsize=8)
310 311 self.pf_axes.append(cax)
311 312 else:
312 313 if self.height is None:
313 314 self.height = 3
314 315 for n in range(self.nplots):
315 316 fig = plt.figure(figsize=(self.width, self.height),
316 317 edgecolor='k',
317 318 facecolor='w')
318 319 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
319 320 ax.tick_params(labelsize=8)
320 321 ax.firsttime = True
321 322 ax.index = 0
322 323 ax.press = None
323 324 self.figures.append(fig)
324 325 self.axes.append(ax)
325 326 if self.showprofile:
326 327 cax = self.__add_axes(ax, size=size, pad=pad)
327 328 cax.tick_params(labelsize=8)
328 329 self.pf_axes.append(cax)
329 330
330 331 for n in range(self.nrows):
331 332 if self.colormaps is not None:
332 333 cmap = plt.get_cmap(self.colormaps[n])
333 334 else:
334 335 cmap = plt.get_cmap(self.colormap)
335 336 cmap.set_bad(self.bgcolor, 1.)
336 337 self.cmaps.append(cmap)
337 338
338 339 def __add_axes(self, ax, size='30%', pad='8%'):
339 340 '''
340 341 Add new axes to the given figure
341 342 '''
342 343 divider = make_axes_locatable(ax)
343 344 nax = divider.new_horizontal(size=size, pad=pad)
344 345 ax.figure.add_axes(nax)
345 346 return nax
346 347
347 348 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
348 349 '''
349 350 Create a masked array for missing data
350 351 '''
351 352 if x_buffer.shape[0] < 2:
352 353 return x_buffer, y_buffer, z_buffer
353 354
354 355 deltas = x_buffer[1:] - x_buffer[0:-1]
355 356 x_median = numpy.median(deltas)
356 357
357 358 index = numpy.where(deltas > 5 * x_median)
358 359
359 360 if len(index[0]) != 0:
360 361 z_buffer[::, index[0], ::] = self.__missing
361 362 z_buffer = numpy.ma.masked_inside(z_buffer,
362 363 0.99 * self.__missing,
363 364 1.01 * self.__missing)
364 365
365 366 return x_buffer, y_buffer, z_buffer
366 367
367 368 def decimate(self):
368 369
369 370 # dx = int(len(self.x)/self.__MAXNUMX) + 1
370 371 dy = int(len(self.y) / self.decimation) + 1
371 372
372 373 # x = self.x[::dx]
373 374 x = self.x
374 375 y = self.y[::dy]
375 376 z = self.z[::, ::, ::dy]
376 377
377 378 return x, y, z
378 379
379 380 def format(self):
380 381 '''
381 382 Set min and max values, labels, ticks and titles
382 383 '''
383 384
384 385 for n, ax in enumerate(self.axes):
385 386 if ax.firsttime:
386 387 if self.xaxis != 'time':
387 388 xmin = self.xmin
388 389 xmax = self.xmax
389 390 else:
390 391 xmin = self.tmin
391 392 xmax = self.tmin + self.xrange*60*60
392 393 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
393 394 ax.xaxis.set_major_locator(LinearLocator(9))
394 395 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
395 396 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
396 397 ax.set_facecolor(self.bgcolor)
397 398 if self.xscale:
398 399 ax.xaxis.set_major_formatter(FuncFormatter(
399 400 lambda x, pos: '{0:g}'.format(x*self.xscale)))
400 401 if self.yscale:
401 402 ax.yaxis.set_major_formatter(FuncFormatter(
402 403 lambda x, pos: '{0:g}'.format(x*self.yscale)))
403 404 if self.xlabel is not None:
404 405 ax.set_xlabel(self.xlabel)
405 406 if self.ylabel is not None:
406 407 ax.set_ylabel(self.ylabel)
407 408 if self.showprofile:
408 409 self.pf_axes[n].set_ylim(ymin, ymax)
409 410 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
410 411 self.pf_axes[n].set_xlabel('dB')
411 412 self.pf_axes[n].grid(b=True, axis='x')
412 413 [tick.set_visible(False)
413 414 for tick in self.pf_axes[n].get_yticklabels()]
414 415 if self.colorbar:
415 416 ax.cbar = plt.colorbar(
416 417 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
417 418 ax.cbar.ax.tick_params(labelsize=8)
418 419 ax.cbar.ax.press = None
419 420 if self.cb_label:
420 421 ax.cbar.set_label(self.cb_label, size=8)
421 422 elif self.cb_labels:
422 423 ax.cbar.set_label(self.cb_labels[n], size=8)
423 424 else:
424 425 ax.cbar = None
425 426 ax.set_xlim(xmin, xmax)
426 427 ax.set_ylim(ymin, ymax)
427 428 ax.firsttime = False
428 429 if self.grid:
429 430 ax.grid(True)
430 431 if not self.polar:
431 432 ax.set_title('{} {} {}'.format(
432 433 self.titles[n],
433 434 self.getDateTime(self.data.max_time).strftime(
434 435 '%Y-%m-%d %H:%M:%S'),
435 436 self.time_label),
436 437 size=8)
437 438 else:
438 439 #ax.set_title('{}'.format(self.titles[n]), size=8)
439 440 ax.set_title('{} {} {}'.format(
440 441 self.titles[n],
441 442 self.getDateTime(self.data.max_time).strftime(
442 443 '%Y-%m-%d %H:%M:%S'),
443 444 self.time_label),
444 445 size=8)
445 446 ax.set_ylim(0, self.ymax)
446 447 #ax.set_yticks(numpy.arange(0, self.ymax, 20))
447 448 ax.yaxis.labelpad = 20
448 449
449 450 if self.firsttime:
450 451 for n, fig in enumerate(self.figures):
451 452 fig.subplots_adjust(**self.plots_adjust)
452 453 self.firsttime = False
453 454
454 455 def clear_figures(self):
455 456 '''
456 457 Reset axes for redraw plots
457 458 '''
458 459
459 460 for ax in self.axes+self.pf_axes+self.cb_axes:
460 461 ax.clear()
461 462 ax.firsttime = True
462 463 if hasattr(ax, 'cbar') and ax.cbar:
463 464 ax.cbar.remove()
464 465
465 466 def __plot(self):
466 467 '''
467 468 Main function to plot, format and save figures
468 469 '''
469 470
470 471 self.plot()
471 472 self.format()
472 473
473 474 for n, fig in enumerate(self.figures):
474 475 if self.nrows == 0 or self.nplots == 0:
475 476 log.warning('No data', self.name)
476 477 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
477 478 fig.canvas.manager.set_window_title(self.CODE)
478 479 continue
479 480
480 481 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
481 482 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
482 483 fig.canvas.draw()
483 484 if self.show:
484 485 fig.show()
485 486 figpause(0.01)
486 487
487 488 if self.save:
488 489 self.save_figure(n)
489 490
490 491 if self.server:
491 492 self.send_to_server()
492 493
493 494 def __update(self, dataOut, timestamp):
494 495 '''
495 496 '''
496 497
497 498 metadata = {
498 499 'yrange': dataOut.heightList,
499 500 'interval': dataOut.timeInterval,
500 501 'channels': dataOut.channelList
501 502 }
502 503
503 504 data, meta = self.update(dataOut)
504 505 metadata.update(meta)
505 506 self.data.update(data, timestamp, metadata)
506 507
507 508 def save_figure(self, n):
508 509 '''
509 510 '''
510 511 if self.oneFigure:
511 512 if (self.data.max_time - self.save_time) <= self.save_period:
512 513 return
513 514
514 515 self.save_time = self.data.max_time
515 516
516 517 fig = self.figures[n]
518 print("save_code",self.save_code)
517 519 if self.throttle == 0:
518 520 if self.oneFigure:
519 521 figname = os.path.join(
520 522 self.save,
521 523 self.save_code,
522 524 '{}_{}.png'.format(
523 525 self.save_code,
524 526 self.getDateTime(self.data.max_time).strftime(
525 527 '%Y%m%d_%H%M%S'
526 528 ),
527 529 )
528 530 )
529 531 else:
530 532 figname = os.path.join(
531 533 self.save,
532 534 self.save_code,
533 535 '{}_ch{}_{}.png'.format(
534 536 self.save_code,n,
535 537 self.getDateTime(self.data.max_time).strftime(
536 538 '%Y%m%d_%H%M%S'
537 539 ),
538 540 )
539 541 )
540 542 log.log('Saving figure: {}'.format(figname), self.name)
541 543 if not os.path.isdir(os.path.dirname(figname)):
542 544 os.makedirs(os.path.dirname(figname))
543 545 fig.savefig(figname)
544 546
545 547 figname = os.path.join(
546 548 self.save,
547 549 '{}_{}.png'.format(
548 550 self.save_code,
549 551 self.getDateTime(self.data.min_time).strftime(
550 552 '%Y%m%d'
551 553 ),
552 554 )
553 555 )
554 556
555 557 log.log('Saving figure: {}'.format(figname), self.name)
556 558 if not os.path.isdir(os.path.dirname(figname)):
557 559 os.makedirs(os.path.dirname(figname))
558 560 fig.savefig(figname)
559 561
560 562 def send_to_server(self):
561 563 '''
562 564 '''
563 565
564 566 if self.exp_code == None:
565 567 log.warning('Missing `exp_code` skipping sending to server...')
566 568
567 569 last_time = self.data.max_time
568 570 interval = last_time - self.sender_time
569 571 if interval < self.sender_period:
570 572 return
571 573
572 574 self.sender_time = last_time
573 575
574 576 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
575 577 for attr in attrs:
576 578 value = getattr(self, attr)
577 579 if value:
578 580 if isinstance(value, (numpy.float32, numpy.float64)):
579 581 value = round(float(value), 2)
580 582 self.data.meta[attr] = value
581 583 if self.colormap == 'jet':
582 584 self.data.meta['colormap'] = 'Jet'
583 585 elif 'RdBu' in self.colormap:
584 586 self.data.meta['colormap'] = 'RdBu'
585 587 else:
586 588 self.data.meta['colormap'] = 'Viridis'
587 589 self.data.meta['interval'] = int(interval)
588 590
589 591 self.sender_queue.append(last_time)
590 592
591 593 while True:
592 594 try:
593 595 tm = self.sender_queue.popleft()
594 596 except IndexError:
595 597 break
596 598 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
597 599 self.socket.send_string(msg)
598 600 socks = dict(self.poll.poll(2000))
599 601 if socks.get(self.socket) == zmq.POLLIN:
600 602 reply = self.socket.recv_string()
601 603 if reply == 'ok':
602 604 log.log("Response from server ok", self.name)
603 605 time.sleep(0.1)
604 606 continue
605 607 else:
606 608 log.warning(
607 609 "Malformed reply from server: {}".format(reply), self.name)
608 610 else:
609 611 log.warning(
610 612 "No response from server, retrying...", self.name)
611 613 self.sender_queue.appendleft(tm)
612 614 self.socket.setsockopt(zmq.LINGER, 0)
613 615 self.socket.close()
614 616 self.poll.unregister(self.socket)
615 617 self.socket = self.context.socket(zmq.REQ)
616 618 self.socket.connect(self.server)
617 619 self.poll.register(self.socket, zmq.POLLIN)
618 620 break
619 621
620 622 def setup(self):
621 623 '''
622 624 This method should be implemented in the child class, the following
623 625 attributes should be set:
624 626
625 627 self.nrows: number of rows
626 628 self.ncols: number of cols
627 629 self.nplots: number of plots (channels or pairs)
628 630 self.ylabel: label for Y axes
629 631 self.titles: list of axes title
630 632
631 633 '''
632 634 raise NotImplementedError
633 635
634 636 def plot(self):
635 637 '''
636 638 Must be defined in the child class, the actual plotting method
637 639 '''
638 640 raise NotImplementedError
639 641
640 642 def update(self, dataOut):
641 643 '''
642 644 Must be defined in the child class, update self.data with new data
643 645 '''
644 646
645 647 data = {
646 648 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
647 649 }
648 650 meta = {}
649 651
650 652 return data, meta
651 653
652 654 def run(self, dataOut, **kwargs):
653 655 '''
654 656 Main plotting routine
655 657 '''
656 658
657 659 if self.isConfig is False:
658 660 self.__setup(**kwargs)
659 661
660 662 if self.localtime:
661 663 self.getDateTime = datetime.datetime.fromtimestamp
662 664 else:
663 665 self.getDateTime = datetime.datetime.utcfromtimestamp
664 666
665 667 self.data.setup()
666 668 self.isConfig = True
667 669 if self.server:
668 670 self.context = zmq.Context()
669 671 self.socket = self.context.socket(zmq.REQ)
670 672 self.socket.connect(self.server)
671 673 self.poll = zmq.Poller()
672 674 self.poll.register(self.socket, zmq.POLLIN)
673 675
674 676 tm = getattr(dataOut, self.attr_time)
675 677
676 678 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
677 679 self.save_time = tm
678 680 self.__plot()
679 681 self.tmin += self.xrange*60*60
680 682 self.data.setup()
681 683 self.clear_figures()
682 684
683 685 self.__update(dataOut, tm)
684 686
685 687 if self.isPlotConfig is False:
686 688 self.__setup_plot()
687 689 self.isPlotConfig = True
688 690 if self.xaxis == 'time':
689 691 dt = self.getDateTime(tm)
690 692 if self.xmin is None:
691 693 self.tmin = tm
692 694 self.xmin = dt.hour
693 695 minutes = (self.xmin-int(self.xmin)) * 60
694 696 seconds = (minutes - int(minutes)) * 60
695 697 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
696 698 datetime.datetime(1970, 1, 1)).total_seconds()
697 699 if self.localtime:
698 700 self.tmin += time.timezone
699 701
700 702 if self.xmin is not None and self.xmax is not None:
701 703 self.xrange = self.xmax - self.xmin
702 704
703 705 if self.throttle == 0:
704 706 self.__plot()
705 707 else:
706 708 self.__throttle_plot(self.__plot)#, coerce=coerce)
707 709
708 710 def close(self):
709 711
710 712 if self.data and not self.data.flagNoData:
711 713 self.save_time = 0
712 714 self.__plot()
713 715 if self.data and not self.data.flagNoData and self.pause:
714 716 figpause(10)
Show file before | Show file after | Hide comments
@@ -1,1793 +1,1932
1 1 import os
2 2 import datetime
3 3 import numpy
4 4 from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
5 5
6 6 from schainpy.model.graphics.jroplot_base import Plot, plt
7 7 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
8 8 from schainpy.utils import log
9 9 # libreria wradlib
10 import wradlib as wrl
10 #import wradlib as wrl
11 11
12 12 EARTH_RADIUS = 6.3710e3
13 13
14 14
15 15 def ll2xy(lat1, lon1, lat2, lon2):
16 16
17 17 p = 0.017453292519943295
18 18 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
19 19 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
20 20 r = 12742 * numpy.arcsin(numpy.sqrt(a))
21 21 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
22 22 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
23 23 theta = -theta + numpy.pi/2
24 24 return r*numpy.cos(theta), r*numpy.sin(theta)
25 25
26 26
27 27 def km2deg(km):
28 28 '''
29 29 Convert distance in km to degrees
30 30 '''
31 31
32 32 return numpy.rad2deg(km/EARTH_RADIUS)
33 33
34 34
35 35
36 36 class SpectralMomentsPlot(SpectraPlot):
37 37 '''
38 38 Plot for Spectral Moments
39 39 '''
40 40 CODE = 'spc_moments'
41 41 # colormap = 'jet'
42 42 # plot_type = 'pcolor'
43 43
44 44 class DobleGaussianPlot(SpectraPlot):
45 45 '''
46 46 Plot for Double Gaussian Plot
47 47 '''
48 48 CODE = 'gaussian_fit'
49 49 # colormap = 'jet'
50 50 # plot_type = 'pcolor'
51 51
52 52 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
53 53 '''
54 54 Plot SpectraCut with Double Gaussian Fit
55 55 '''
56 56 CODE = 'cut_gaussian_fit'
57 57
58 58 class SnrPlot(RTIPlot):
59 59 '''
60 60 Plot for SNR Data
61 61 '''
62 62
63 63 CODE = 'snr'
64 64 colormap = 'jet'
65 65
66 66 def update(self, dataOut):
67 67
68 68 data = {
69 69 'snr': 10*numpy.log10(dataOut.data_snr)
70 70 }
71 71
72 72 return data, {}
73 73
74 74 class DopplerPlot(RTIPlot):
75 75 '''
76 76 Plot for DOPPLER Data (1st moment)
77 77 '''
78 78
79 79 CODE = 'dop'
80 80 colormap = 'jet'
81 81
82 82 def update(self, dataOut):
83 83
84 84 data = {
85 85 'dop': 10*numpy.log10(dataOut.data_dop)
86 86 }
87 87
88 88 return data, {}
89 89
90 90 class PowerPlot(RTIPlot):
91 91 '''
92 92 Plot for Power Data (0 moment)
93 93 '''
94 94
95 95 CODE = 'pow'
96 96 colormap = 'jet'
97 97
98 98 def update(self, dataOut):
99 99 data = {
100 100 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
101 101 }
102 102 return data, {}
103 103
104 104 class SpectralWidthPlot(RTIPlot):
105 105 '''
106 106 Plot for Spectral Width Data (2nd moment)
107 107 '''
108 108
109 109 CODE = 'width'
110 110 colormap = 'jet'
111 111
112 112 def update(self, dataOut):
113 113
114 114 data = {
115 115 'width': dataOut.data_width
116 116 }
117 117
118 118 return data, {}
119 119
120 120 class SkyMapPlot(Plot):
121 121 '''
122 122 Plot for meteors detection data
123 123 '''
124 124
125 125 CODE = 'param'
126 126
127 127 def setup(self):
128 128
129 129 self.ncols = 1
130 130 self.nrows = 1
131 131 self.width = 7.2
132 132 self.height = 7.2
133 133 self.nplots = 1
134 134 self.xlabel = 'Zonal Zenith Angle (deg)'
135 135 self.ylabel = 'Meridional Zenith Angle (deg)'
136 136 self.polar = True
137 137 self.ymin = -180
138 138 self.ymax = 180
139 139 self.colorbar = False
140 140
141 141 def plot(self):
142 142
143 143 arrayParameters = numpy.concatenate(self.data['param'])
144 144 error = arrayParameters[:, -1]
145 145 indValid = numpy.where(error == 0)[0]
146 146 finalMeteor = arrayParameters[indValid, :]
147 147 finalAzimuth = finalMeteor[:, 3]
148 148 finalZenith = finalMeteor[:, 4]
149 149
150 150 x = finalAzimuth * numpy.pi / 180
151 151 y = finalZenith
152 152
153 153 ax = self.axes[0]
154 154
155 155 if ax.firsttime:
156 156 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
157 157 else:
158 158 ax.plot.set_data(x, y)
159 159
160 160 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
161 161 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
162 162 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
163 163 dt2,
164 164 len(x))
165 165 self.titles[0] = title
166 166
167 167
168 168 class GenericRTIPlot(Plot):
169 169 '''
170 170 Plot for data_xxxx object
171 171 '''
172 172
173 173 CODE = 'param'
174 174 colormap = 'viridis'
175 175 plot_type = 'pcolorbuffer'
176 176
177 177 def setup(self):
178 178 self.xaxis = 'time'
179 179 self.ncols = 1
180 180 self.nrows = self.data.shape('param')[0]
181 181 self.nplots = self.nrows
182 182 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
183 183
184 184 if not self.xlabel:
185 185 self.xlabel = 'Time'
186 186
187 187 self.ylabel = 'Range [km]'
188 188 if not self.titles:
189 189 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
190 190
191 191 def update(self, dataOut):
192 192
193 193 data = {
194 194 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
195 195 }
196 196
197 197 meta = {}
198 198
199 199 return data, meta
200 200
201 201 def plot(self):
202 202 # self.data.normalize_heights()
203 203 self.x = self.data.times
204 204 self.y = self.data.yrange
205 205 self.z = self.data['param']
206 206 self.z = 10*numpy.log10(self.z)
207 207 self.z = numpy.ma.masked_invalid(self.z)
208 208
209 209 if self.decimation is None:
210 210 x, y, z = self.fill_gaps(self.x, self.y, self.z)
211 211 else:
212 212 x, y, z = self.fill_gaps(*self.decimate())
213 213
214 214 for n, ax in enumerate(self.axes):
215 215
216 216 self.zmax = self.zmax if self.zmax is not None else numpy.max(
217 217 self.z[n])
218 218 self.zmin = self.zmin if self.zmin is not None else numpy.min(
219 219 self.z[n])
220 220
221 221 if ax.firsttime:
222 222 if self.zlimits is not None:
223 223 self.zmin, self.zmax = self.zlimits[n]
224 224
225 225 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
226 226 vmin=self.zmin,
227 227 vmax=self.zmax,
228 228 cmap=self.cmaps[n]
229 229 )
230 230 else:
231 231 if self.zlimits is not None:
232 232 self.zmin, self.zmax = self.zlimits[n]
233 233 ax.collections.remove(ax.collections[0])
234 234 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
235 235 vmin=self.zmin,
236 236 vmax=self.zmax,
237 237 cmap=self.cmaps[n]
238 238 )
239 239
240 240
241 241 class PolarMapPlot(Plot):
242 242 '''
243 243 Plot for weather radar
244 244 '''
245 245
246 246 CODE = 'param'
247 247 colormap = 'seismic'
248 248
249 249 def setup(self):
250 250 self.ncols = 1
251 251 self.nrows = 1
252 252 self.width = 9
253 253 self.height = 8
254 254 self.mode = self.data.meta['mode']
255 255 if self.channels is not None:
256 256 self.nplots = len(self.channels)
257 257 self.nrows = len(self.channels)
258 258 else:
259 259 self.nplots = self.data.shape(self.CODE)[0]
260 260 self.nrows = self.nplots
261 261 self.channels = list(range(self.nplots))
262 262 if self.mode == 'E':
263 263 self.xlabel = 'Longitude'
264 264 self.ylabel = 'Latitude'
265 265 else:
266 266 self.xlabel = 'Range (km)'
267 267 self.ylabel = 'Height (km)'
268 268 self.bgcolor = 'white'
269 269 self.cb_labels = self.data.meta['units']
270 270 self.lat = self.data.meta['latitude']
271 271 self.lon = self.data.meta['longitude']
272 272 self.xmin, self.xmax = float(
273 273 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
274 274 self.ymin, self.ymax = float(
275 275 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
276 276 # self.polar = True
277 277
278 278 def plot(self):
279 279
280 280 for n, ax in enumerate(self.axes):
281 281 data = self.data['param'][self.channels[n]]
282 282
283 283 zeniths = numpy.linspace(
284 284 0, self.data.meta['max_range'], data.shape[1])
285 285 if self.mode == 'E':
286 286 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
287 287 r, theta = numpy.meshgrid(zeniths, azimuths)
288 288 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
289 289 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
290 290 x = km2deg(x) + self.lon
291 291 y = km2deg(y) + self.lat
292 292 else:
293 293 azimuths = numpy.radians(self.data.yrange)
294 294 r, theta = numpy.meshgrid(zeniths, azimuths)
295 295 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
296 296 self.y = zeniths
297 297
298 298 if ax.firsttime:
299 299 if self.zlimits is not None:
300 300 self.zmin, self.zmax = self.zlimits[n]
301 301 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
302 302 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
303 303 vmin=self.zmin,
304 304 vmax=self.zmax,
305 305 cmap=self.cmaps[n])
306 306 else:
307 307 if self.zlimits is not None:
308 308 self.zmin, self.zmax = self.zlimits[n]
309 309 ax.collections.remove(ax.collections[0])
310 310 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
311 311 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
312 312 vmin=self.zmin,
313 313 vmax=self.zmax,
314 314 cmap=self.cmaps[n])
315 315
316 316 if self.mode == 'A':
317 317 continue
318 318
319 319 # plot district names
320 320 f = open('/data/workspace/schain_scripts/distrito.csv')
321 321 for line in f:
322 322 label, lon, lat = [s.strip() for s in line.split(',') if s]
323 323 lat = float(lat)
324 324 lon = float(lon)
325 325 # ax.plot(lon, lat, '.b', ms=2)
326 326 ax.text(lon, lat, label.decode('utf8'), ha='center',
327 327 va='bottom', size='8', color='black')
328 328
329 329 # plot limites
330 330 limites = []
331 331 tmp = []
332 332 for line in open('/data/workspace/schain_scripts/lima.csv'):
333 333 if '#' in line:
334 334 if tmp:
335 335 limites.append(tmp)
336 336 tmp = []
337 337 continue
338 338 values = line.strip().split(',')
339 339 tmp.append((float(values[0]), float(values[1])))
340 340 for points in limites:
341 341 ax.add_patch(
342 342 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
343 343
344 344 # plot Cuencas
345 345 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
346 346 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
347 347 values = [line.strip().split(',') for line in f]
348 348 points = [(float(s[0]), float(s[1])) for s in values]
349 349 ax.add_patch(Polygon(points, ec='b', fc='none'))
350 350
351 351 # plot grid
352 352 for r in (15, 30, 45, 60):
353 353 ax.add_artist(plt.Circle((self.lon, self.lat),
354 354 km2deg(r), color='0.6', fill=False, lw=0.2))
355 355 ax.text(
356 356 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
357 357 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
358 358 '{}km'.format(r),
359 359 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
360 360
361 361 if self.mode == 'E':
362 362 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
363 363 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
364 364 else:
365 365 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
366 366 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
367 367
368 368 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
369 369 self.titles = ['{} {}'.format(
370 370 self.data.parameters[x], title) for x in self.channels]
371 371
372 372 class WeatherPlot(Plot):
373 373 CODE = 'weather'
374 374 plot_name = 'weather'
375 375 plot_type = 'ppistyle'
376 376 buffering = False
377 377
378 378 def setup(self):
379 379 self.ncols = 1
380 380 self.nrows = 1
381 381 self.width =8
382 382 self.height =8
383 383 self.nplots= 1
384 384 self.ylabel= 'Range [Km]'
385 385 self.titles= ['Weather']
386 386 self.colorbar=False
387 387 self.ini =0
388 388 self.len_azi =0
389 389 self.buffer_ini = None
390 390 self.buffer_azi = None
391 391 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
392 392 self.flag =0
393 393 self.indicador= 0
394 394 self.last_data_azi = None
395 395 self.val_mean = None
396 396
397 397 def update(self, dataOut):
398 398
399 399 data = {}
400 400 meta = {}
401 401 if hasattr(dataOut, 'dataPP_POWER'):
402 402 factor = 1
403 403 if hasattr(dataOut, 'nFFTPoints'):
404 404 factor = dataOut.normFactor
405 405 #print("DIME EL SHAPE PORFAVOR",dataOut.data_360.shape)
406 406 data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor))
407 407 data['azi'] = dataOut.data_azi
408 408 data['ele'] = dataOut.data_ele
409 409 return data, meta
410 410
411 411 def get2List(self,angulos):
412 412 list1=[]
413 413 list2=[]
414 414 for i in reversed(range(len(angulos))):
415 415 diff_ = angulos[i]-angulos[i-1]
416 416 if diff_ >1.5:
417 417 list1.append(i-1)
418 418 list2.append(diff_)
419 419 return list(reversed(list1)),list(reversed(list2))
420 420
421 421 def fixData360(self,list_,ang_):
422 422 if list_[0]==-1:
423 423 vec = numpy.where(ang_<ang_[0])
424 424 ang_[vec] = ang_[vec]+360
425 425 return ang_
426 426 return ang_
427 427
428 428 def fixData360HL(self,angulos):
429 429 vec = numpy.where(angulos>=360)
430 430 angulos[vec]=angulos[vec]-360
431 431 return angulos
432 432
433 433 def search_pos(self,pos,list_):
434 434 for i in range(len(list_)):
435 435 if pos == list_[i]:
436 436 return True,i
437 437 i=None
438 438 return False,i
439 439
440 440 def fixDataComp(self,ang_,list1_,list2_):
441 441 size = len(ang_)
442 442 size2 = 0
443 443 for i in range(len(list2_)):
444 444 size2=size2+round(list2_[i])-1
445 445 new_size= size+size2
446 446 ang_new = numpy.zeros(new_size)
447 447 ang_new2 = numpy.zeros(new_size)
448 448
449 449 tmp = 0
450 450 c = 0
451 451 for i in range(len(ang_)):
452 452 ang_new[tmp +c] = ang_[i]
453 453 ang_new2[tmp+c] = ang_[i]
454 454 condition , value = self.search_pos(i,list1_)
455 455 if condition:
456 456 pos = tmp + c + 1
457 457 for k in range(round(list2_[value])-1):
458 458 ang_new[pos+k] = ang_new[pos+k-1]+1
459 459 ang_new2[pos+k] = numpy.nan
460 460 tmp = pos +k
461 461 c = 0
462 462 c=c+1
463 463 return ang_new,ang_new2
464 464
465 465 def globalCheckPED(self,angulos):
466 466 l1,l2 = self.get2List(angulos)
467 467 if len(l1)>0:
468 468 angulos2 = self.fixData360(list_=l1,ang_=angulos)
469 469 l1,l2 = self.get2List(angulos2)
470 470
471 471 ang1_,ang2_ = self.fixDataComp(ang_=angulos2,list1_=l1,list2_=l2)
472 472 ang1_ = self.fixData360HL(ang1_)
473 473 ang2_ = self.fixData360HL(ang2_)
474 474 else:
475 475 ang1_= angulos
476 476 ang2_= angulos
477 477 return ang1_,ang2_
478 478
479 479 def analizeDATA(self,data_azi):
480 480 list1 = []
481 481 list2 = []
482 482 dat = data_azi
483 483 for i in reversed(range(1,len(dat))):
484 484 if dat[i]>dat[i-1]:
485 485 diff = int(dat[i])-int(dat[i-1])
486 486 else:
487 487 diff = 360+int(dat[i])-int(dat[i-1])
488 488 if diff > 1:
489 489 list1.append(i-1)
490 490 list2.append(diff-1)
491 491 return list1,list2
492 492
493 493 def fixDATANEW(self,data_azi,data_weather):
494 494 list1,list2 = self.analizeDATA(data_azi)
495 495 if len(list1)== 0:
496 496 return data_azi,data_weather
497 497 else:
498 498 resize = 0
499 499 for i in range(len(list2)):
500 500 resize= resize + list2[i]
501 501 new_data_azi = numpy.resize(data_azi,resize)
502 502 new_data_weather= numpy.resize(date_weather,resize)
503 503
504 504 for i in range(len(list2)):
505 505 j=0
506 506 position=list1[i]+1
507 507 for j in range(list2[i]):
508 508 new_data_azi[position+j]=new_data_azi[position+j-1]+1
509 509 return new_data_azi
510 510
511 511 def fixDATA(self,data_azi):
512 512 data=data_azi
513 513 for i in range(len(data)):
514 514 if numpy.isnan(data[i]):
515 515 data[i]=data[i-1]+1
516 516 return data
517 517
518 518 def replaceNAN(self,data_weather,data_azi,val):
519 519 data= data_azi
520 520 data_T= data_weather
521 521 if data.shape[0]> data_T.shape[0]:
522 522 data_N = numpy.ones( [data.shape[0],data_T.shape[1]])
523 523 c = 0
524 524 for i in range(len(data)):
525 525 if numpy.isnan(data[i]):
526 526 data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
527 527 else:
528 528 data_N[i,:]=data_T[c,:]
529 529 c=c+1
530 530 return data_N
531 531 else:
532 532 for i in range(len(data)):
533 533 if numpy.isnan(data[i]):
534 534 data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
535 535 return data_T
536 536
537 537 def const_ploteo(self,data_weather,data_azi,step,res):
538 538 if self.ini==0:
539 539 #-------
540 540 n = (360/res)-len(data_azi)
541 541 #--------------------- new -------------------------
542 542 data_azi_new ,data_azi_old= self.globalCheckPED(data_azi)
543 543 #------------------------
544 544 start = data_azi_new[-1] + res
545 545 end = data_azi_new[0] - res
546 546 #------ new
547 547 self.last_data_azi = end
548 548 if start>end:
549 549 end = end + 360
550 550 azi_vacia = numpy.linspace(start,end,int(n))
551 551 azi_vacia = numpy.where(azi_vacia>360,azi_vacia-360,azi_vacia)
552 552 data_azi = numpy.hstack((data_azi_new,azi_vacia))
553 553 # RADAR
554 554 val_mean = numpy.mean(data_weather[:,-1])
555 555 self.val_mean = val_mean
556 556 data_weather_cmp = numpy.ones([(360-data_weather.shape[0]),data_weather.shape[1]])*val_mean
557 557 data_weather = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean)
558 558 data_weather = numpy.vstack((data_weather,data_weather_cmp))
559 559 else:
560 560 # azimuth
561 561 flag=0
562 562 start_azi = self.res_azi[0]
563 563 #-----------new------------
564 564 data_azi ,data_azi_old= self.globalCheckPED(data_azi)
565 565 data_weather = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean)
566 566 #--------------------------
567 567 start = data_azi[0]
568 568 end = data_azi[-1]
569 569 self.last_data_azi= end
570 570 if start< start_azi:
571 571 start = start +360
572 572 if end <start_azi:
573 573 end = end +360
574 574
575 575 pos_ini = int((start-start_azi)/res)
576 576 len_azi = len(data_azi)
577 577 if (360-pos_ini)<len_azi:
578 578 if pos_ini+1==360:
579 579 pos_ini=0
580 580 else:
581 581 flag=1
582 582 dif= 360-pos_ini
583 583 comp= len_azi-dif
584 584 #-----------------
585 585 if flag==0:
586 586 # AZIMUTH
587 587 self.res_azi[pos_ini:pos_ini+len_azi] = data_azi
588 588 # RADAR
589 589 self.res_weather[pos_ini:pos_ini+len_azi,:] = data_weather
590 590 else:
591 591 # AZIMUTH
592 592 self.res_azi[pos_ini:pos_ini+dif] = data_azi[0:dif]
593 593 self.res_azi[0:comp] = data_azi[dif:]
594 594 # RADAR
595 595 self.res_weather[pos_ini:pos_ini+dif,:] = data_weather[0:dif,:]
596 596 self.res_weather[0:comp,:] = data_weather[dif:,:]
597 597 flag=0
598 598 data_azi = self.res_azi
599 599 data_weather = self.res_weather
600 600
601 601 return data_weather,data_azi
602 602
603 603 def plot(self):
604 604 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
605 605 data = self.data[-1]
606 606 r = self.data.yrange
607 607 delta_height = r[1]-r[0]
608 608 r_mask = numpy.where(r>=0)[0]
609 609 r = numpy.arange(len(r_mask))*delta_height
610 610 self.y = 2*r
611 611 # RADAR
612 612 #data_weather = data['weather']
613 613 # PEDESTAL
614 614 #data_azi = data['azi']
615 615 res = 1
616 616 # STEP
617 617 step = (360/(res*data['weather'].shape[0]))
618 618
619 619 self.res_weather, self.res_azi = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_azi=data['azi'],step=step,res=res)
620 620 self.res_ele = numpy.mean(data['ele'])
621 621 ################# PLOTEO ###################
622 622 for i,ax in enumerate(self.axes):
623 623 self.zmin = self.zmin if self.zmin else 20
624 624 self.zmax = self.zmax if self.zmax else 80
625 625 if ax.firsttime:
626 626 plt.clf()
627 627 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=self.zmin, vmax=self.zmax)
628 628 else:
629 629 plt.clf()
630 630 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=self.zmin, vmax=self.zmax)
631 631 caax = cgax.parasites[0]
632 632 paax = cgax.parasites[1]
633 633 cbar = plt.gcf().colorbar(pm, pad=0.075)
634 634 caax.set_xlabel('x_range [km]')
635 635 caax.set_ylabel('y_range [km]')
636 636 plt.text(1.0, 1.05, 'Azimuth '+str(thisDatetime)+" Step "+str(self.ini)+ " EL: "+str(round(self.res_ele, 1)), transform=caax.transAxes, va='bottom',ha='right')
637 637
638 638 self.ini= self.ini+1
639 639
640 640
641 641 class WeatherRHIPlot(Plot):
642 642 CODE = 'weather'
643 643 plot_name = 'weather'
644 644 plot_type = 'rhistyle'
645 645 buffering = False
646 646 data_ele_tmp = None
647 647
648 648 def setup(self):
649 649 print("********************")
650 650 print("********************")
651 651 print("********************")
652 652 print("SETUP WEATHER PLOT")
653 653 self.ncols = 1
654 654 self.nrows = 1
655 655 self.nplots= 1
656 656 self.ylabel= 'Range [Km]'
657 657 self.titles= ['Weather']
658 658 if self.channels is not None:
659 659 self.nplots = len(self.channels)
660 660 self.nrows = len(self.channels)
661 661 else:
662 662 self.nplots = self.data.shape(self.CODE)[0]
663 663 self.nrows = self.nplots
664 664 self.channels = list(range(self.nplots))
665 665 print("channels",self.channels)
666 666 print("que saldra", self.data.shape(self.CODE)[0])
667 667 self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)]
668 668 print("self.titles",self.titles)
669 669 self.colorbar=False
670 670 self.width =12
671 671 self.height =8
672 672 self.ini =0
673 673 self.len_azi =0
674 674 self.buffer_ini = None
675 675 self.buffer_ele = None
676 676 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
677 677 self.flag =0
678 678 self.indicador= 0
679 679 self.last_data_ele = None
680 680 self.val_mean = None
681 681
682 682 def update(self, dataOut):
683 683
684 684 data = {}
685 685 meta = {}
686 686 if hasattr(dataOut, 'dataPP_POWER'):
687 687 factor = 1
688 688 if hasattr(dataOut, 'nFFTPoints'):
689 689 factor = dataOut.normFactor
690 690 print("dataOut",dataOut.data_360.shape)
691 691 #
692 692 data['weather'] = 10*numpy.log10(dataOut.data_360/(factor))
693 693 #
694 694 #data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor))
695 695 data['azi'] = dataOut.data_azi
696 696 data['ele'] = dataOut.data_ele
697 697 #print("UPDATE")
698 698 #print("data[weather]",data['weather'].shape)
699 699 #print("data[azi]",data['azi'])
700 700 return data, meta
701 701
702 702 def get2List(self,angulos):
703 703 list1=[]
704 704 list2=[]
705 705 for i in reversed(range(len(angulos))):
706 706 if not i==0:#el caso de i=0 evalula el primero de la lista con el ultimo y no es relevante
707 707 diff_ = angulos[i]-angulos[i-1]
708 708 if abs(diff_) >1.5:
709 709 list1.append(i-1)
710 710 list2.append(diff_)
711 711 return list(reversed(list1)),list(reversed(list2))
712 712
713 713 def fixData90(self,list_,ang_):
714 714 if list_[0]==-1:
715 715 vec = numpy.where(ang_<ang_[0])
716 716 ang_[vec] = ang_[vec]+90
717 717 return ang_
718 718 return ang_
719 719
720 720 def fixData90HL(self,angulos):
721 721 vec = numpy.where(angulos>=90)
722 722 angulos[vec]=angulos[vec]-90
723 723 return angulos
724 724
725 725
726 726 def search_pos(self,pos,list_):
727 727 for i in range(len(list_)):
728 728 if pos == list_[i]:
729 729 return True,i
730 730 i=None
731 731 return False,i
732 732
733 733 def fixDataComp(self,ang_,list1_,list2_,tipo_case):
734 734 size = len(ang_)
735 735 size2 = 0
736 736 for i in range(len(list2_)):
737 737 size2=size2+round(abs(list2_[i]))-1
738 738 new_size= size+size2
739 739 ang_new = numpy.zeros(new_size)
740 740 ang_new2 = numpy.zeros(new_size)
741 741
742 742 tmp = 0
743 743 c = 0
744 744 for i in range(len(ang_)):
745 745 ang_new[tmp +c] = ang_[i]
746 746 ang_new2[tmp+c] = ang_[i]
747 747 condition , value = self.search_pos(i,list1_)
748 748 if condition:
749 749 pos = tmp + c + 1
750 750 for k in range(round(abs(list2_[value]))-1):
751 751 if tipo_case==0 or tipo_case==3:#subida
752 752 ang_new[pos+k] = ang_new[pos+k-1]+1
753 753 ang_new2[pos+k] = numpy.nan
754 754 elif tipo_case==1 or tipo_case==2:#bajada
755 755 ang_new[pos+k] = ang_new[pos+k-1]-1
756 756 ang_new2[pos+k] = numpy.nan
757 757
758 758 tmp = pos +k
759 759 c = 0
760 760 c=c+1
761 761 return ang_new,ang_new2
762 762
763 763 def globalCheckPED(self,angulos,tipo_case):
764 764 l1,l2 = self.get2List(angulos)
765 765 ##print("l1",l1)
766 766 ##print("l2",l2)
767 767 if len(l1)>0:
768 768 #angulos2 = self.fixData90(list_=l1,ang_=angulos)
769 769 #l1,l2 = self.get2List(angulos2)
770 770 ang1_,ang2_ = self.fixDataComp(ang_=angulos,list1_=l1,list2_=l2,tipo_case=tipo_case)
771 771 #ang1_ = self.fixData90HL(ang1_)
772 772 #ang2_ = self.fixData90HL(ang2_)
773 773 else:
774 774 ang1_= angulos
775 775 ang2_= angulos
776 776 return ang1_,ang2_
777 777
778 778
779 779 def replaceNAN(self,data_weather,data_ele,val):
780 780 data= data_ele
781 781 data_T= data_weather
782 782 if data.shape[0]> data_T.shape[0]:
783 783 data_N = numpy.ones( [data.shape[0],data_T.shape[1]])
784 784 c = 0
785 785 for i in range(len(data)):
786 786 if numpy.isnan(data[i]):
787 787 data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
788 788 else:
789 789 data_N[i,:]=data_T[c,:]
790 790 c=c+1
791 791 return data_N
792 792 else:
793 793 for i in range(len(data)):
794 794 if numpy.isnan(data[i]):
795 795 data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
796 796 return data_T
797 797
798 798 def check_case(self,data_ele,ang_max,ang_min):
799 799 start = data_ele[0]
800 800 end = data_ele[-1]
801 801 number = (end-start)
802 802 len_ang=len(data_ele)
803 803 print("start",start)
804 804 print("end",end)
805 805 print("number",number)
806 806
807 807 print("len_ang",len_ang)
808 808
809 809 #exit(1)
810 810
811 811 if start<end and (round(abs(number)+1)>=len_ang or (numpy.argmin(data_ele)==0)):#caso subida
812 812 return 0
813 813 #elif start>end and (round(abs(number)+1)>=len_ang or(numpy.argmax(data_ele)==0)):#caso bajada
814 814 # return 1
815 815 elif round(abs(number)+1)>=len_ang and (start>end or(numpy.argmax(data_ele)==0)):#caso bajada
816 816 return 1
817 817 elif round(abs(number)+1)<len_ang and data_ele[-2]>data_ele[-1]:# caso BAJADA CAMBIO ANG MAX
818 818 return 2
819 819 elif round(abs(number)+1)<len_ang and data_ele[-2]<data_ele[-1] :# caso SUBIDA CAMBIO ANG MIN
820 820 return 3
821 821
822 822
823 823 def const_ploteo(self,val_ch,data_weather,data_ele,step,res,ang_max,ang_min):
824 824 ang_max= ang_max
825 825 ang_min= ang_min
826 826 data_weather=data_weather
827 827 val_ch=val_ch
828 828 ##print("*********************DATA WEATHER**************************************")
829 829 ##print(data_weather)
830 830 if self.ini==0:
831 831 '''
832 832 print("**********************************************")
833 833 print("**********************************************")
834 834 print("***************ini**************")
835 835 print("**********************************************")
836 836 print("**********************************************")
837 837 '''
838 838 #print("data_ele",data_ele)
839 839 #----------------------------------------------------------
840 840 tipo_case = self.check_case(data_ele,ang_max,ang_min)
841 841 print("check_case",tipo_case)
842 842 #exit(1)
843 843 #--------------------- new -------------------------
844 844 data_ele_new ,data_ele_old= self.globalCheckPED(data_ele,tipo_case)
845 845
846 846 #-------------------------CAMBIOS RHI---------------------------------
847 847 start= ang_min
848 848 end = ang_max
849 849 n= (ang_max-ang_min)/res
850 850 #------ new
851 851 self.start_data_ele = data_ele_new[0]
852 852 self.end_data_ele = data_ele_new[-1]
853 853 if tipo_case==0 or tipo_case==3: # SUBIDA
854 854 n1= round(self.start_data_ele)- start
855 855 n2= end - round(self.end_data_ele)
856 856 print(self.start_data_ele)
857 857 print(self.end_data_ele)
858 858 if n1>0:
859 859 ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1)
860 860 ele1_nan= numpy.ones(n1)*numpy.nan
861 861 data_ele = numpy.hstack((ele1,data_ele_new))
862 862 print("ele1_nan",ele1_nan.shape)
863 863 print("data_ele_old",data_ele_old.shape)
864 864 data_ele_old = numpy.hstack((ele1_nan,data_ele_old))
865 865 if n2>0:
866 866 ele2= numpy.linspace(self.end_data_ele+1,end,n2)
867 867 ele2_nan= numpy.ones(n2)*numpy.nan
868 868 data_ele = numpy.hstack((data_ele,ele2))
869 869 print("ele2_nan",ele2_nan.shape)
870 870 print("data_ele_old",data_ele_old.shape)
871 871 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
872 872
873 873 if tipo_case==1 or tipo_case==2: # BAJADA
874 874 data_ele_new = data_ele_new[::-1] # reversa
875 875 data_ele_old = data_ele_old[::-1]# reversa
876 876 data_weather = data_weather[::-1,:]# reversa
877 877 vec= numpy.where(data_ele_new<ang_max)
878 878 data_ele_new = data_ele_new[vec]
879 879 data_ele_old = data_ele_old[vec]
880 880 data_weather = data_weather[vec[0]]
881 881 vec2= numpy.where(0<data_ele_new)
882 882 data_ele_new = data_ele_new[vec2]
883 883 data_ele_old = data_ele_old[vec2]
884 884 data_weather = data_weather[vec2[0]]
885 885 self.start_data_ele = data_ele_new[0]
886 886 self.end_data_ele = data_ele_new[-1]
887 887
888 888 n1= round(self.start_data_ele)- start
889 889 n2= end - round(self.end_data_ele)-1
890 890 print(self.start_data_ele)
891 891 print(self.end_data_ele)
892 892 if n1>0:
893 893 ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1)
894 894 ele1_nan= numpy.ones(n1)*numpy.nan
895 895 data_ele = numpy.hstack((ele1,data_ele_new))
896 896 data_ele_old = numpy.hstack((ele1_nan,data_ele_old))
897 897 if n2>0:
898 898 ele2= numpy.linspace(self.end_data_ele+1,end,n2)
899 899 ele2_nan= numpy.ones(n2)*numpy.nan
900 900 data_ele = numpy.hstack((data_ele,ele2))
901 901 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
902 902 # RADAR
903 903 # NOTA data_ele y data_weather es la variable que retorna
904 904 val_mean = numpy.mean(data_weather[:,-1])
905 905 self.val_mean = val_mean
906 906 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
907 907 self.data_ele_tmp[val_ch]= data_ele_old
908 908 else:
909 909 #print("**********************************************")
910 910 #print("****************VARIABLE**********************")
911 911 #-------------------------CAMBIOS RHI---------------------------------
912 912 #---------------------------------------------------------------------
913 913 ##print("INPUT data_ele",data_ele)
914 914 flag=0
915 915 start_ele = self.res_ele[0]
916 916 tipo_case = self.check_case(data_ele,ang_max,ang_min)
917 917 #print("TIPO DE DATA",tipo_case)
918 918 #-----------new------------
919 919 data_ele ,data_ele_old = self.globalCheckPED(data_ele,tipo_case)
920 920 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
921 921
922 922 #-------------------------------NEW RHI ITERATIVO-------------------------
923 923
924 924 if tipo_case==0 : # SUBIDA
925 925 vec = numpy.where(data_ele<ang_max)
926 926 data_ele = data_ele[vec]
927 927 data_ele_old = data_ele_old[vec]
928 928 data_weather = data_weather[vec[0]]
929 929
930 930 vec2 = numpy.where(0<data_ele)
931 931 data_ele= data_ele[vec2]
932 932 data_ele_old= data_ele_old[vec2]
933 933 ##print(data_ele_new)
934 934 data_weather= data_weather[vec2[0]]
935 935
936 936 new_i_ele = int(round(data_ele[0]))
937 937 new_f_ele = int(round(data_ele[-1]))
938 938 #print(new_i_ele)
939 939 #print(new_f_ele)
940 940 #print(data_ele,len(data_ele))
941 941 #print(data_ele_old,len(data_ele_old))
942 942 if new_i_ele< 2:
943 943 self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan
944 944 self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean)
945 945 self.data_ele_tmp[val_ch][new_i_ele:new_i_ele+len(data_ele)]=data_ele_old
946 946 self.res_ele[new_i_ele:new_i_ele+len(data_ele)]= data_ele
947 947 self.res_weather[val_ch][new_i_ele:new_i_ele+len(data_ele),:]= data_weather
948 948 data_ele = self.res_ele
949 949 data_weather = self.res_weather[val_ch]
950 950
951 951 elif tipo_case==1 : #BAJADA
952 952 data_ele = data_ele[::-1] # reversa
953 953 data_ele_old = data_ele_old[::-1]# reversa
954 954 data_weather = data_weather[::-1,:]# reversa
955 955 vec= numpy.where(data_ele<ang_max)
956 956 data_ele = data_ele[vec]
957 957 data_ele_old = data_ele_old[vec]
958 958 data_weather = data_weather[vec[0]]
959 959 vec2= numpy.where(0<data_ele)
960 960 data_ele = data_ele[vec2]
961 961 data_ele_old = data_ele_old[vec2]
962 962 data_weather = data_weather[vec2[0]]
963 963
964 964
965 965 new_i_ele = int(round(data_ele[0]))
966 966 new_f_ele = int(round(data_ele[-1]))
967 967 #print(data_ele)
968 968 #print(ang_max)
969 969 #print(data_ele_old)
970 970 if new_i_ele <= 1:
971 971 new_i_ele = 1
972 972 if round(data_ele[-1])>=ang_max-1:
973 973 self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan
974 974 self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean)
975 975 self.data_ele_tmp[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1]=data_ele_old
976 976 self.res_ele[new_i_ele-1:new_i_ele+len(data_ele)-1]= data_ele
977 977 self.res_weather[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1,:]= data_weather
978 978 data_ele = self.res_ele
979 979 data_weather = self.res_weather[val_ch]
980 980
981 981 elif tipo_case==2: #bajada
982 982 vec = numpy.where(data_ele<ang_max)
983 983 data_ele = data_ele[vec]
984 984 data_weather= data_weather[vec[0]]
985 985
986 986 len_vec = len(vec)
987 987 data_ele_new = data_ele[::-1] # reversa
988 988 data_weather = data_weather[::-1,:]
989 989 new_i_ele = int(data_ele_new[0])
990 990 new_f_ele = int(data_ele_new[-1])
991 991
992 992 n1= new_i_ele- ang_min
993 993 n2= ang_max - new_f_ele-1
994 994 if n1>0:
995 995 ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1)
996 996 ele1_nan= numpy.ones(n1)*numpy.nan
997 997 data_ele = numpy.hstack((ele1,data_ele_new))
998 998 data_ele_old = numpy.hstack((ele1_nan,data_ele_new))
999 999 if n2>0:
1000 1000 ele2= numpy.linspace(new_f_ele+1,ang_max,n2)
1001 1001 ele2_nan= numpy.ones(n2)*numpy.nan
1002 1002 data_ele = numpy.hstack((data_ele,ele2))
1003 1003 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1004 1004
1005 1005 self.data_ele_tmp[val_ch] = data_ele_old
1006 1006 self.res_ele = data_ele
1007 1007 self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1008 1008 data_ele = self.res_ele
1009 1009 data_weather = self.res_weather[val_ch]
1010 1010
1011 1011 elif tipo_case==3:#subida
1012 1012 vec = numpy.where(0<data_ele)
1013 1013 data_ele= data_ele[vec]
1014 1014 data_ele_new = data_ele
1015 1015 data_ele_old= data_ele_old[vec]
1016 1016 data_weather= data_weather[vec[0]]
1017 1017 pos_ini = numpy.argmin(data_ele)
1018 1018 if pos_ini>0:
1019 1019 len_vec= len(data_ele)
1020 1020 vec3 = numpy.linspace(pos_ini,len_vec-1,len_vec-pos_ini).astype(int)
1021 1021 #print(vec3)
1022 1022 data_ele= data_ele[vec3]
1023 1023 data_ele_new = data_ele
1024 1024 data_ele_old= data_ele_old[vec3]
1025 1025 data_weather= data_weather[vec3]
1026 1026
1027 1027 new_i_ele = int(data_ele_new[0])
1028 1028 new_f_ele = int(data_ele_new[-1])
1029 1029 n1= new_i_ele- ang_min
1030 1030 n2= ang_max - new_f_ele-1
1031 1031 if n1>0:
1032 1032 ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1)
1033 1033 ele1_nan= numpy.ones(n1)*numpy.nan
1034 1034 data_ele = numpy.hstack((ele1,data_ele_new))
1035 1035 data_ele_old = numpy.hstack((ele1_nan,data_ele_new))
1036 1036 if n2>0:
1037 1037 ele2= numpy.linspace(new_f_ele+1,ang_max,n2)
1038 1038 ele2_nan= numpy.ones(n2)*numpy.nan
1039 1039 data_ele = numpy.hstack((data_ele,ele2))
1040 1040 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1041 1041
1042 1042 self.data_ele_tmp[val_ch] = data_ele_old
1043 1043 self.res_ele = data_ele
1044 1044 self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1045 1045 data_ele = self.res_ele
1046 1046 data_weather = self.res_weather[val_ch]
1047 1047 #print("self.data_ele_tmp",self.data_ele_tmp)
1048 1048 return data_weather,data_ele
1049 1049
1050 1050
1051 1051 def plot(self):
1052 1052 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
1053 1053 data = self.data[-1]
1054 1054 r = self.data.yrange
1055 1055 delta_height = r[1]-r[0]
1056 1056 r_mask = numpy.where(r>=0)[0]
1057 1057 ##print("delta_height",delta_height)
1058 1058 #print("r_mask",r_mask,len(r_mask))
1059 1059 r = numpy.arange(len(r_mask))*delta_height
1060 1060 self.y = 2*r
1061 1061 res = 1
1062 1062 ###print("data['weather'].shape[0]",data['weather'].shape[0])
1063 1063 ang_max = self.ang_max
1064 1064 ang_min = self.ang_min
1065 1065 var_ang =ang_max - ang_min
1066 1066 step = (int(var_ang)/(res*data['weather'].shape[0]))
1067 1067 ###print("step",step)
1068 1068 #--------------------------------------------------------
1069 1069 ##print('weather',data['weather'].shape)
1070 1070 ##print('ele',data['ele'].shape)
1071 1071
1072 1072 ###self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min)
1073 1073 ###self.res_azi = numpy.mean(data['azi'])
1074 1074 ###print("self.res_ele",self.res_ele)
1075 1075 plt.clf()
1076 1076 subplots = [121, 122]
1077 1077 cg={'angular_spacing': 20.}
1078 1078 if self.ini==0:
1079 1079 self.data_ele_tmp = numpy.ones([self.nplots,int(var_ang)])*numpy.nan
1080 1080 self.res_weather= numpy.ones([self.nplots,int(var_ang),len(r_mask)])*numpy.nan
1081 1081 print("SHAPE",self.data_ele_tmp.shape)
1082 1082
1083 1083 for i,ax in enumerate(self.axes):
1084 1084 self.res_weather[i], self.res_ele = self.const_ploteo(val_ch=i, data_weather=data['weather'][i][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min)
1085 1085 self.res_azi = numpy.mean(data['azi'])
1086 1086 if i==0:
1087 1087 print("*****************************************************************************to plot**************************",self.res_weather[i].shape)
1088 1088 self.zmin = self.zmin if self.zmin else 20
1089 1089 self.zmax = self.zmax if self.zmax else 80
1090 1090 if ax.firsttime:
1091 1091 #plt.clf()
1092 1092 cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj=cg,vmin=self.zmin, vmax=self.zmax)
1093 1093 #fig=self.figures[0]
1094 1094 else:
1095 1095 #plt.clf()
1096 1096 if i==0:
1097 1097 print(self.res_weather[i])
1098 1098 print(self.res_ele)
1099 1099 cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj=cg,vmin=self.zmin, vmax=self.zmax)
1100 1100 caax = cgax.parasites[0]
1101 1101 paax = cgax.parasites[1]
1102 1102 cbar = plt.gcf().colorbar(pm, pad=0.075)
1103 1103 caax.set_xlabel('x_range [km]')
1104 1104 caax.set_ylabel('y_range [km]')
1105 1105 plt.text(1.0, 1.05, 'Elevacion '+str(thisDatetime)+" Step "+str(self.ini)+ " Azi: "+str(round(self.res_azi,2)), transform=caax.transAxes, va='bottom',ha='right')
1106 1106 print("***************************self.ini****************************",self.ini)
1107 1107 self.ini= self.ini+1
1108 1108
1109 1109 class Weather_vRF_Plot(Plot):
1110 1110 CODE = 'PPI'
1111 1111 plot_name = 'PPI'
1112 1112 #plot_type = 'ppistyle'
1113 1113 buffering = False
1114 1114
1115 1115 def setup(self):
1116 1116
1117 1117 self.ncols = 1
1118 1118 self.nrows = 1
1119 1119 self.width =8
1120 1120 self.height =8
1121 1121 self.nplots= 1
1122 1122 self.ylabel= 'Range [Km]'
1123 1123 self.xlabel= 'Range [Km]'
1124 1124 self.titles= ['PPI']
1125 1125 self.polar = True
1126 1126 if self.channels is not None:
1127 1127 self.nplots = len(self.channels)
1128 1128 self.nrows = len(self.channels)
1129 1129 else:
1130 1130 self.nplots = self.data.shape(self.CODE)[0]
1131 1131 self.nrows = self.nplots
1132 1132 self.channels = list(range(self.nplots))
1133 1133
1134 1134 if self.CODE == 'POWER':
1135 1135 self.cb_label = r'Power (dB)'
1136 1136 elif self.CODE == 'DOPPLER':
1137 1137 self.cb_label = r'Velocity (m/s)'
1138 1138 self.colorbar=True
1139 1139 self.width = 9
1140 1140 self.height =8
1141 1141 self.ini =0
1142 1142 self.len_azi =0
1143 1143 self.buffer_ini = None
1144 1144 self.buffer_ele = None
1145 1145 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.15, 'right': 0.9, 'bottom': 0.08})
1146 1146 self.flag =0
1147 1147 self.indicador= 0
1148 1148 self.last_data_ele = None
1149 1149 self.val_mean = None
1150 1150
1151 1151 def update(self, dataOut):
1152 1152
1153 1153 data = {}
1154 1154 meta = {}
1155 1155 if hasattr(dataOut, 'dataPP_POWER'):
1156 1156 factor = 1
1157 1157 if hasattr(dataOut, 'nFFTPoints'):
1158 1158 factor = dataOut.normFactor
1159 1159
1160 1160 if 'pow' in self.attr_data[0].lower():
1161 1161 data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))
1162 1162 else:
1163 1163 data['data'] = getattr(dataOut, self.attr_data[0])/(factor)
1164 1164
1165 1165 data['azi'] = dataOut.data_azi
1166 1166 data['ele'] = dataOut.data_ele
1167 1167
1168 1168 return data, meta
1169 1169
1170 1170 def plot(self):
1171 1171 data = self.data[-1]
1172 1172 r = self.data.yrange
1173 1173 delta_height = r[1]-r[0]
1174 1174 r_mask = numpy.where(r>=0)[0]
1175 1175 self.r_mask = r_mask
1176 1176 r = numpy.arange(len(r_mask))*delta_height
1177 1177 self.y = 2*r
1178 1178
1179 z = data['data'][self.channels[0]][:,r_mask]
1179 try:
1180 z = data['data'][self.channels[0]][:,r_mask]
1181
1182 except:
1183 z = data['data'][0][:,r_mask]
1180 1184
1181 1185 self.titles = []
1182 1186
1183 1187 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
1184 1188 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
1185 1189 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1186 1190 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1187 1191 self.ang_min = self.ang_min if self.ang_min else 0
1188 1192 self.ang_max = self.ang_max if self.ang_max else 360
1189 1193
1190 1194 r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )
1191 1195
1192 1196 for i,ax in enumerate(self.axes):
1193 1197
1194 1198 if ax.firsttime:
1195 1199 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1196 1200 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1197 1201 ax.set_theta_direction(-1)
1198 1202
1199 1203 else:
1200 1204 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1201 1205 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1202 1206 ax.set_theta_direction(-1)
1203 1207
1204 1208 ax.grid(True)
1205 1209
1206 1210 if len(self.channels) !=1:
1207 1211 self.titles = ['PPI {} at EL: {} Channel {}'.format(self.self.labels[x], str(round(numpy.mean(data['ele']),1)), x) for x in range(self.nrows)]
1208 1212 else:
1209 1213 self.titles = ['PPI {} at EL: {} Channel {}'.format(self.labels[0], str(round(numpy.mean(data['ele']),1)), self.channels[0])]
1210 1214
1211 1215 class WeatherRHI_vRF2_Plot(Plot):
1212 1216 CODE = 'weather'
1213 1217 plot_name = 'weather'
1214 1218 plot_type = 'rhistyle'
1215 1219 buffering = False
1216 1220 data_ele_tmp = None
1217 1221
1218 1222 def setup(self):
1219 1223 print("********************")
1220 1224 print("********************")
1221 1225 print("********************")
1222 1226 print("SETUP WEATHER PLOT")
1223 1227 self.ncols = 1
1224 1228 self.nrows = 1
1225 1229 self.nplots= 1
1226 1230 self.ylabel= 'Range [Km]'
1227 1231 self.titles= ['Weather']
1228 1232 if self.channels is not None:
1229 1233 self.nplots = len(self.channels)
1230 1234 self.nrows = len(self.channels)
1231 1235 else:
1232 1236 self.nplots = self.data.shape(self.CODE)[0]
1233 1237 self.nrows = self.nplots
1234 1238 self.channels = list(range(self.nplots))
1235 1239 print("channels",self.channels)
1236 1240 print("que saldra", self.data.shape(self.CODE)[0])
1237 1241 self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)]
1238 1242 print("self.titles",self.titles)
1239 1243 self.colorbar=False
1240 1244 self.width =8
1241 1245 self.height =8
1242 1246 self.ini =0
1243 1247 self.len_azi =0
1244 1248 self.buffer_ini = None
1245 1249 self.buffer_ele = None
1246 1250 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
1247 1251 self.flag =0
1248 1252 self.indicador= 0
1249 1253 self.last_data_ele = None
1250 1254 self.val_mean = None
1251 1255
1252 1256 def update(self, dataOut):
1253 1257
1254 1258 data = {}
1255 1259 meta = {}
1256 1260 if hasattr(dataOut, 'dataPP_POWER'):
1257 1261 factor = 1
1258 1262 if hasattr(dataOut, 'nFFTPoints'):
1259 1263 factor = dataOut.normFactor
1260 1264 print("dataOut",dataOut.data_360.shape)
1261 1265 #
1262 1266 data['weather'] = 10*numpy.log10(dataOut.data_360/(factor))
1263 1267 #
1264 1268 #data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor))
1265 1269 data['azi'] = dataOut.data_azi
1266 1270 data['ele'] = dataOut.data_ele
1267 1271 data['case_flag'] = dataOut.case_flag
1268 1272 #print("UPDATE")
1269 1273 #print("data[weather]",data['weather'].shape)
1270 1274 #print("data[azi]",data['azi'])
1271 1275 return data, meta
1272 1276
1273 1277 def get2List(self,angulos):
1274 1278 list1=[]
1275 1279 list2=[]
1276 1280 for i in reversed(range(len(angulos))):
1277 1281 if not i==0:#el caso de i=0 evalula el primero de la lista con el ultimo y no es relevante
1278 1282 diff_ = angulos[i]-angulos[i-1]
1279 1283 if abs(diff_) >1.5:
1280 1284 list1.append(i-1)
1281 1285 list2.append(diff_)
1282 1286 return list(reversed(list1)),list(reversed(list2))
1283 1287
1284 1288 def fixData90(self,list_,ang_):
1285 1289 if list_[0]==-1:
1286 1290 vec = numpy.where(ang_<ang_[0])
1287 1291 ang_[vec] = ang_[vec]+90
1288 1292 return ang_
1289 1293 return ang_
1290 1294
1291 1295 def fixData90HL(self,angulos):
1292 1296 vec = numpy.where(angulos>=90)
1293 1297 angulos[vec]=angulos[vec]-90
1294 1298 return angulos
1295 1299
1296 1300
1297 1301 def search_pos(self,pos,list_):
1298 1302 for i in range(len(list_)):
1299 1303 if pos == list_[i]:
1300 1304 return True,i
1301 1305 i=None
1302 1306 return False,i
1303 1307
1304 1308 def fixDataComp(self,ang_,list1_,list2_,tipo_case):
1305 1309 size = len(ang_)
1306 1310 size2 = 0
1307 1311 for i in range(len(list2_)):
1308 1312 size2=size2+round(abs(list2_[i]))-1
1309 1313 new_size= size+size2
1310 1314 ang_new = numpy.zeros(new_size)
1311 1315 ang_new2 = numpy.zeros(new_size)
1312 1316
1313 1317 tmp = 0
1314 1318 c = 0
1315 1319 for i in range(len(ang_)):
1316 1320 ang_new[tmp +c] = ang_[i]
1317 1321 ang_new2[tmp+c] = ang_[i]
1318 1322 condition , value = self.search_pos(i,list1_)
1319 1323 if condition:
1320 1324 pos = tmp + c + 1
1321 1325 for k in range(round(abs(list2_[value]))-1):
1322 1326 if tipo_case==0 or tipo_case==3:#subida
1323 1327 ang_new[pos+k] = ang_new[pos+k-1]+1
1324 1328 ang_new2[pos+k] = numpy.nan
1325 1329 elif tipo_case==1 or tipo_case==2:#bajada
1326 1330 ang_new[pos+k] = ang_new[pos+k-1]-1
1327 1331 ang_new2[pos+k] = numpy.nan
1328 1332
1329 1333 tmp = pos +k
1330 1334 c = 0
1331 1335 c=c+1
1332 1336 return ang_new,ang_new2
1333 1337
1334 1338 def globalCheckPED(self,angulos,tipo_case):
1335 1339 l1,l2 = self.get2List(angulos)
1336 1340 ##print("l1",l1)
1337 1341 ##print("l2",l2)
1338 1342 if len(l1)>0:
1339 1343 #angulos2 = self.fixData90(list_=l1,ang_=angulos)
1340 1344 #l1,l2 = self.get2List(angulos2)
1341 1345 ang1_,ang2_ = self.fixDataComp(ang_=angulos,list1_=l1,list2_=l2,tipo_case=tipo_case)
1342 1346 #ang1_ = self.fixData90HL(ang1_)
1343 1347 #ang2_ = self.fixData90HL(ang2_)
1344 1348 else:
1345 1349 ang1_= angulos
1346 1350 ang2_= angulos
1347 1351 return ang1_,ang2_
1348 1352
1349 1353
1350 1354 def replaceNAN(self,data_weather,data_ele,val):
1351 1355 data= data_ele
1352 1356 data_T= data_weather
1353 1357 if data.shape[0]> data_T.shape[0]:
1354 1358 data_N = numpy.ones( [data.shape[0],data_T.shape[1]])
1355 1359 c = 0
1356 1360 for i in range(len(data)):
1357 1361 if numpy.isnan(data[i]):
1358 1362 data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
1359 1363 else:
1360 1364 data_N[i,:]=data_T[c,:]
1361 1365 c=c+1
1362 1366 return data_N
1363 1367 else:
1364 1368 for i in range(len(data)):
1365 1369 if numpy.isnan(data[i]):
1366 1370 data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
1367 1371 return data_T
1368 1372
1369 1373 def check_case(self,data_ele,ang_max,ang_min):
1370 1374 start = data_ele[0]
1371 1375 end = data_ele[-1]
1372 1376 number = (end-start)
1373 1377 len_ang=len(data_ele)
1374 1378 print("start",start)
1375 1379 print("end",end)
1376 1380 print("number",number)
1377 1381
1378 1382 print("len_ang",len_ang)
1379 1383
1380 1384 #exit(1)
1381 1385
1382 1386 if start<end and (round(abs(number)+1)>=len_ang or (numpy.argmin(data_ele)==0)):#caso subida
1383 1387 return 0
1384 1388 #elif start>end and (round(abs(number)+1)>=len_ang or(numpy.argmax(data_ele)==0)):#caso bajada
1385 1389 # return 1
1386 1390 elif round(abs(number)+1)>=len_ang and (start>end or(numpy.argmax(data_ele)==0)):#caso bajada
1387 1391 return 1
1388 1392 elif round(abs(number)+1)<len_ang and data_ele[-2]>data_ele[-1]:# caso BAJADA CAMBIO ANG MAX
1389 1393 return 2
1390 1394 elif round(abs(number)+1)<len_ang and data_ele[-2]<data_ele[-1] :# caso SUBIDA CAMBIO ANG MIN
1391 1395 return 3
1392 1396
1393 1397
1394 1398 def const_ploteo(self,val_ch,data_weather,data_ele,step,res,ang_max,ang_min,case_flag):
1395 1399 ang_max= ang_max
1396 1400 ang_min= ang_min
1397 1401 data_weather=data_weather
1398 1402 val_ch=val_ch
1399 1403 ##print("*********************DATA WEATHER**************************************")
1400 1404 ##print(data_weather)
1401 1405 if self.ini==0:
1402 1406 '''
1403 1407 print("**********************************************")
1404 1408 print("**********************************************")
1405 1409 print("***************ini**************")
1406 1410 print("**********************************************")
1407 1411 print("**********************************************")
1408 1412 '''
1409 1413 #print("data_ele",data_ele)
1410 1414 #----------------------------------------------------------
1411 1415 tipo_case = case_flag[-1]
1412 1416 #tipo_case = self.check_case(data_ele,ang_max,ang_min)
1413 1417 print("check_case",tipo_case)
1414 1418 #exit(1)
1415 1419 #--------------------- new -------------------------
1416 1420 data_ele_new ,data_ele_old= self.globalCheckPED(data_ele,tipo_case)
1417 1421
1418 1422 #-------------------------CAMBIOS RHI---------------------------------
1419 1423 start= ang_min
1420 1424 end = ang_max
1421 1425 n= (ang_max-ang_min)/res
1422 1426 #------ new
1423 1427 self.start_data_ele = data_ele_new[0]
1424 1428 self.end_data_ele = data_ele_new[-1]
1425 1429 if tipo_case==0 or tipo_case==3: # SUBIDA
1426 1430 n1= round(self.start_data_ele)- start
1427 1431 n2= end - round(self.end_data_ele)
1428 1432 print(self.start_data_ele)
1429 1433 print(self.end_data_ele)
1430 1434 if n1>0:
1431 1435 ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1)
1432 1436 ele1_nan= numpy.ones(n1)*numpy.nan
1433 1437 data_ele = numpy.hstack((ele1,data_ele_new))
1434 1438 print("ele1_nan",ele1_nan.shape)
1435 1439 print("data_ele_old",data_ele_old.shape)
1436 1440 data_ele_old = numpy.hstack((ele1_nan,data_ele_old))
1437 1441 if n2>0:
1438 1442 ele2= numpy.linspace(self.end_data_ele+1,end,n2)
1439 1443 ele2_nan= numpy.ones(n2)*numpy.nan
1440 1444 data_ele = numpy.hstack((data_ele,ele2))
1441 1445 print("ele2_nan",ele2_nan.shape)
1442 1446 print("data_ele_old",data_ele_old.shape)
1443 1447 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1444 1448
1445 1449 if tipo_case==1 or tipo_case==2: # BAJADA
1446 1450 data_ele_new = data_ele_new[::-1] # reversa
1447 1451 data_ele_old = data_ele_old[::-1]# reversa
1448 1452 data_weather = data_weather[::-1,:]# reversa
1449 1453 vec= numpy.where(data_ele_new<ang_max)
1450 1454 data_ele_new = data_ele_new[vec]
1451 1455 data_ele_old = data_ele_old[vec]
1452 1456 data_weather = data_weather[vec[0]]
1453 1457 vec2= numpy.where(0<data_ele_new)
1454 1458 data_ele_new = data_ele_new[vec2]
1455 1459 data_ele_old = data_ele_old[vec2]
1456 1460 data_weather = data_weather[vec2[0]]
1457 1461 self.start_data_ele = data_ele_new[0]
1458 1462 self.end_data_ele = data_ele_new[-1]
1459 1463
1460 1464 n1= round(self.start_data_ele)- start
1461 1465 n2= end - round(self.end_data_ele)-1
1462 1466 print(self.start_data_ele)
1463 1467 print(self.end_data_ele)
1464 1468 if n1>0:
1465 1469 ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1)
1466 1470 ele1_nan= numpy.ones(n1)*numpy.nan
1467 1471 data_ele = numpy.hstack((ele1,data_ele_new))
1468 1472 data_ele_old = numpy.hstack((ele1_nan,data_ele_old))
1469 1473 if n2>0:
1470 1474 ele2= numpy.linspace(self.end_data_ele+1,end,n2)
1471 1475 ele2_nan= numpy.ones(n2)*numpy.nan
1472 1476 data_ele = numpy.hstack((data_ele,ele2))
1473 1477 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1474 1478 # RADAR
1475 1479 # NOTA data_ele y data_weather es la variable que retorna
1476 1480 val_mean = numpy.mean(data_weather[:,-1])
1477 1481 self.val_mean = val_mean
1478 1482 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1479 1483 print("eleold",data_ele_old)
1480 1484 print(self.data_ele_tmp[val_ch])
1481 1485 print(data_ele_old.shape[0])
1482 1486 print(self.data_ele_tmp[val_ch].shape[0])
1483 1487 if (data_ele_old.shape[0]==91 or self.data_ele_tmp[val_ch].shape[0]==91):
1484 1488 import sys
1485 1489 print("EXIT",self.ini)
1486 1490
1487 1491 sys.exit(1)
1488 1492 self.data_ele_tmp[val_ch]= data_ele_old
1489 1493 else:
1490 1494 #print("**********************************************")
1491 1495 #print("****************VARIABLE**********************")
1492 1496 #-------------------------CAMBIOS RHI---------------------------------
1493 1497 #---------------------------------------------------------------------
1494 1498 ##print("INPUT data_ele",data_ele)
1495 1499 flag=0
1496 1500 start_ele = self.res_ele[0]
1497 1501 #tipo_case = self.check_case(data_ele,ang_max,ang_min)
1498 1502 tipo_case = case_flag[-1]
1499 1503 #print("TIPO DE DATA",tipo_case)
1500 1504 #-----------new------------
1501 1505 data_ele ,data_ele_old = self.globalCheckPED(data_ele,tipo_case)
1502 1506 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1503 1507
1504 1508 #-------------------------------NEW RHI ITERATIVO-------------------------
1505 1509
1506 1510 if tipo_case==0 : # SUBIDA
1507 1511 vec = numpy.where(data_ele<ang_max)
1508 1512 data_ele = data_ele[vec]
1509 1513 data_ele_old = data_ele_old[vec]
1510 1514 data_weather = data_weather[vec[0]]
1511 1515
1512 1516 vec2 = numpy.where(0<data_ele)
1513 1517 data_ele= data_ele[vec2]
1514 1518 data_ele_old= data_ele_old[vec2]
1515 1519 ##print(data_ele_new)
1516 1520 data_weather= data_weather[vec2[0]]
1517 1521
1518 1522 new_i_ele = int(round(data_ele[0]))
1519 1523 new_f_ele = int(round(data_ele[-1]))
1520 1524 #print(new_i_ele)
1521 1525 #print(new_f_ele)
1522 1526 #print(data_ele,len(data_ele))
1523 1527 #print(data_ele_old,len(data_ele_old))
1524 1528 if new_i_ele< 2:
1525 1529 self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan
1526 1530 self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean)
1527 1531 self.data_ele_tmp[val_ch][new_i_ele:new_i_ele+len(data_ele)]=data_ele_old
1528 1532 self.res_ele[new_i_ele:new_i_ele+len(data_ele)]= data_ele
1529 1533 self.res_weather[val_ch][new_i_ele:new_i_ele+len(data_ele),:]= data_weather
1530 1534 data_ele = self.res_ele
1531 1535 data_weather = self.res_weather[val_ch]
1532 1536
1533 1537 elif tipo_case==1 : #BAJADA
1534 1538 data_ele = data_ele[::-1] # reversa
1535 1539 data_ele_old = data_ele_old[::-1]# reversa
1536 1540 data_weather = data_weather[::-1,:]# reversa
1537 1541 vec= numpy.where(data_ele<ang_max)
1538 1542 data_ele = data_ele[vec]
1539 1543 data_ele_old = data_ele_old[vec]
1540 1544 data_weather = data_weather[vec[0]]
1541 1545 vec2= numpy.where(0<data_ele)
1542 1546 data_ele = data_ele[vec2]
1543 1547 data_ele_old = data_ele_old[vec2]
1544 1548 data_weather = data_weather[vec2[0]]
1545 1549
1546 1550
1547 1551 new_i_ele = int(round(data_ele[0]))
1548 1552 new_f_ele = int(round(data_ele[-1]))
1549 1553 #print(data_ele)
1550 1554 #print(ang_max)
1551 1555 #print(data_ele_old)
1552 1556 if new_i_ele <= 1:
1553 1557 new_i_ele = 1
1554 1558 if round(data_ele[-1])>=ang_max-1:
1555 1559 self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan
1556 1560 self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean)
1557 1561 self.data_ele_tmp[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1]=data_ele_old
1558 1562 self.res_ele[new_i_ele-1:new_i_ele+len(data_ele)-1]= data_ele
1559 1563 self.res_weather[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1,:]= data_weather
1560 1564 data_ele = self.res_ele
1561 1565 data_weather = self.res_weather[val_ch]
1562 1566
1563 1567 elif tipo_case==2: #bajada
1564 1568 vec = numpy.where(data_ele<ang_max)
1565 1569 data_ele = data_ele[vec]
1566 1570 data_weather= data_weather[vec[0]]
1567 1571
1568 1572 len_vec = len(vec)
1569 1573 data_ele_new = data_ele[::-1] # reversa
1570 1574 data_weather = data_weather[::-1,:]
1571 1575 new_i_ele = int(data_ele_new[0])
1572 1576 new_f_ele = int(data_ele_new[-1])
1573 1577
1574 1578 n1= new_i_ele- ang_min
1575 1579 n2= ang_max - new_f_ele-1
1576 1580 if n1>0:
1577 1581 ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1)
1578 1582 ele1_nan= numpy.ones(n1)*numpy.nan
1579 1583 data_ele = numpy.hstack((ele1,data_ele_new))
1580 1584 data_ele_old = numpy.hstack((ele1_nan,data_ele_new))
1581 1585 if n2>0:
1582 1586 ele2= numpy.linspace(new_f_ele+1,ang_max,n2)
1583 1587 ele2_nan= numpy.ones(n2)*numpy.nan
1584 1588 data_ele = numpy.hstack((data_ele,ele2))
1585 1589 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1586 1590
1587 1591 self.data_ele_tmp[val_ch] = data_ele_old
1588 1592 self.res_ele = data_ele
1589 1593 self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1590 1594 data_ele = self.res_ele
1591 1595 data_weather = self.res_weather[val_ch]
1592 1596
1593 1597 elif tipo_case==3:#subida
1594 1598 vec = numpy.where(0<data_ele)
1595 1599 data_ele= data_ele[vec]
1596 1600 data_ele_new = data_ele
1597 1601 data_ele_old= data_ele_old[vec]
1598 1602 data_weather= data_weather[vec[0]]
1599 1603 pos_ini = numpy.argmin(data_ele)
1600 1604 if pos_ini>0:
1601 1605 len_vec= len(data_ele)
1602 1606 vec3 = numpy.linspace(pos_ini,len_vec-1,len_vec-pos_ini).astype(int)
1603 1607 #print(vec3)
1604 1608 data_ele= data_ele[vec3]
1605 1609 data_ele_new = data_ele
1606 1610 data_ele_old= data_ele_old[vec3]
1607 1611 data_weather= data_weather[vec3]
1608 1612
1609 1613 new_i_ele = int(data_ele_new[0])
1610 1614 new_f_ele = int(data_ele_new[-1])
1611 1615 n1= new_i_ele- ang_min
1612 1616 n2= ang_max - new_f_ele-1
1613 1617 if n1>0:
1614 1618 ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1)
1615 1619 ele1_nan= numpy.ones(n1)*numpy.nan
1616 1620 data_ele = numpy.hstack((ele1,data_ele_new))
1617 1621 data_ele_old = numpy.hstack((ele1_nan,data_ele_new))
1618 1622 if n2>0:
1619 1623 ele2= numpy.linspace(new_f_ele+1,ang_max,n2)
1620 1624 ele2_nan= numpy.ones(n2)*numpy.nan
1621 1625 data_ele = numpy.hstack((data_ele,ele2))
1622 1626 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1623 1627
1624 1628 self.data_ele_tmp[val_ch] = data_ele_old
1625 1629 self.res_ele = data_ele
1626 1630 self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1627 1631 data_ele = self.res_ele
1628 1632 data_weather = self.res_weather[val_ch]
1629 1633 #print("self.data_ele_tmp",self.data_ele_tmp)
1630 1634 return data_weather,data_ele
1631 1635
1632 1636
1633 1637 def plot(self):
1634 1638 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
1635 1639 data = self.data[-1]
1636 1640 r = self.data.yrange
1637 1641 delta_height = r[1]-r[0]
1638 1642 r_mask = numpy.where(r>=0)[0]
1639 1643 ##print("delta_height",delta_height)
1640 1644 #print("r_mask",r_mask,len(r_mask))
1641 1645 r = numpy.arange(len(r_mask))*delta_height
1642 1646 self.y = 2*r
1643 1647 res = 1
1644 1648 ###print("data['weather'].shape[0]",data['weather'].shape[0])
1645 1649 ang_max = self.ang_max
1646 1650 ang_min = self.ang_min
1647 1651 var_ang =ang_max - ang_min
1648 1652 step = (int(var_ang)/(res*data['weather'].shape[0]))
1649 1653 ###print("step",step)
1650 1654 #--------------------------------------------------------
1651 1655 ##print('weather',data['weather'].shape)
1652 1656 ##print('ele',data['ele'].shape)
1653 1657
1654 1658 ###self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min)
1655 1659 ###self.res_azi = numpy.mean(data['azi'])
1656 1660 ###print("self.res_ele",self.res_ele)
1657 1661 plt.clf()
1658 1662 subplots = [121, 122]
1659 1663 try:
1660 1664 if self.data[-2]['ele'].max()<data['ele'].max():
1661 1665 self.ini=0
1662 1666 except:
1663 1667 pass
1664 1668 if self.ini==0:
1665 1669 self.data_ele_tmp = numpy.ones([self.nplots,int(var_ang)])*numpy.nan
1666 1670 self.res_weather= numpy.ones([self.nplots,int(var_ang),len(r_mask)])*numpy.nan
1667 1671 print("SHAPE",self.data_ele_tmp.shape)
1668 1672
1669 1673 for i,ax in enumerate(self.axes):
1670 1674 self.res_weather[i], self.res_ele = self.const_ploteo(val_ch=i, data_weather=data['weather'][i][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min,case_flag=self.data['case_flag'])
1671 1675 self.res_azi = numpy.mean(data['azi'])
1672 1676
1673 1677 if ax.firsttime:
1674 1678 #plt.clf()
1675 1679 print("Frist Plot")
1676 1680 cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj='cg',vmin=20, vmax=80)
1677 1681 #fig=self.figures[0]
1678 1682 else:
1679 1683 #plt.clf()
1680 1684 print("ELSE PLOT")
1681 1685 cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj='cg',vmin=20, vmax=80)
1682 1686 caax = cgax.parasites[0]
1683 1687 paax = cgax.parasites[1]
1684 1688 cbar = plt.gcf().colorbar(pm, pad=0.075)
1685 1689 caax.set_xlabel('x_range [km]')
1686 1690 caax.set_ylabel('y_range [km]')
1687 1691 plt.text(1.0, 1.05, 'Elevacion '+str(thisDatetime)+" Step "+str(self.ini)+ " Azi: "+str(round(self.res_azi,2)), transform=caax.transAxes, va='bottom',ha='right')
1688 1692 print("***************************self.ini****************************",self.ini)
1689 1693 self.ini= self.ini+1
1690 1694
1691 1695
1692 1696
1693 1697
1694 1698
1695 1699 class WeatherRHI_vRF4_Plot(Plot):
1696 1700 CODE = 'RHI'
1697 1701 plot_name = 'RHI'
1698 1702 #plot_type = 'rhistyle'
1699 1703 buffering = False
1700 1704
1701 1705 def setup(self):
1702 1706
1703 1707 self.ncols = 1
1704 1708 self.nrows = 1
1705 1709 self.nplots= 1
1706 1710 self.ylabel= 'Range [Km]'
1707 1711 self.xlabel= 'Range [Km]'
1708 1712 self.titles= ['RHI']
1709 1713 self.polar = True
1710 1714 self.grid = True
1711 1715 if self.channels is not None:
1712 1716 self.nplots = len(self.channels)
1713 1717 self.nrows = len(self.channels)
1714 1718 else:
1715 1719 self.nplots = self.data.shape(self.CODE)[0]
1716 1720 self.nrows = self.nplots
1717 1721 self.channels = list(range(self.nplots))
1718 1722
1719 1723 if self.CODE == 'Power':
1720 1724 self.cb_label = r'Power (dB)'
1721 1725 elif self.CODE == 'Doppler':
1722 1726 self.cb_label = r'Velocity (m/s)'
1723 1727 self.colorbar=True
1724 1728 self.width =8
1725 1729 self.height =8
1726 1730 self.ini =0
1727 1731 self.len_azi =0
1728 1732 self.buffer_ini = None
1729 1733 self.buffer_ele = None
1730 1734 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
1731 1735 self.flag =0
1732 1736 self.indicador= 0
1733 1737 self.last_data_ele = None
1734 1738 self.val_mean = None
1735 1739
1736 1740 def update(self, dataOut):
1737 1741
1738 1742 data = {}
1739 1743 meta = {}
1740 1744 if hasattr(dataOut, 'dataPP_POWER'):
1741 1745 factor = 1
1742 1746 if hasattr(dataOut, 'nFFTPoints'):
1743 1747 factor = dataOut.normFactor
1744 1748
1745 1749 if 'pow' in self.attr_data[0].lower():
1746 1750 data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))
1747 1751 else:
1748 1752 data['data'] = getattr(dataOut, self.attr_data[0])/(factor)
1749 1753
1750 1754 data['azi'] = dataOut.data_azi
1751 1755 data['ele'] = dataOut.data_ele
1752 1756
1753 1757 return data, meta
1754 1758
1755 1759 def plot(self):
1756 1760 data = self.data[-1]
1757 1761 r = self.data.yrange
1758 1762 delta_height = r[1]-r[0]
1759 1763 r_mask = numpy.where(r>=0)[0]
1760 1764 self.r_mask =r_mask
1761 1765 r = numpy.arange(len(r_mask))*delta_height
1762 1766 self.y = 2*r
1763 1767
1764 1768 try:
1765 1769 z = data['data'][self.channels[0]][:,r_mask]
1766 1770 except:
1767 1771 z = data['data'][0][:,r_mask]
1768 1772
1769 1773 self.titles = []
1770 1774
1771 1775 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
1772 1776 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
1773 1777 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1774 1778 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1775 1779 self.ang_min = self.ang_min if self.ang_min else 0
1776 1780 self.ang_max = self.ang_max if self.ang_max else 90
1777 1781
1778 1782 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )
1779 1783
1780 1784 for i,ax in enumerate(self.axes):
1781 1785
1782 1786 if ax.firsttime:
1783 1787 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1784 1788 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1785 1789
1786 1790 else:
1787 1791 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1788 1792 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1789 1793 ax.grid(True)
1790 1794 if len(self.channels) !=1:
1791 1795 self.titles = ['RHI {} at AZ: {} Channel {}'.format(self.labels[x], str(round(numpy.mean(data['azi']),1)), x) for x in range(self.nrows)]
1792 1796 else:
1793 1797 self.titles = ['RHI {} at AZ: {} Channel {}'.format(self.labels[0], str(round(numpy.mean(data['azi']),1)), self.channels[0])]
1798
1799 class WeatherParamsPlot(Plot):
1800 #CODE = 'RHI'
1801 #plot_name = 'RHI'
1802 #plot_type = 'rhistyle'
1803 buffering = False
1804
1805 def setup(self):
1806
1807 self.ncols = 1
1808 self.nrows = 1
1809 self.nplots= 1
1810 self.ylabel= 'Range [Km]'
1811 self.xlabel= 'Range [Km]'
1812 self.polar = True
1813 self.grid = True
1814 if self.channels is not None:
1815 self.nplots = len(self.channels)
1816 self.nrows = len(self.channels)
1817 else:
1818 self.nplots = self.data.shape(self.CODE)[0]
1819 self.nrows = self.nplots
1820 self.channels = list(range(self.nplots))
1821
1822 self.colorbar=True
1823 self.width =8
1824 self.height =8
1825 self.ini =0
1826 self.len_azi =0
1827 self.buffer_ini = None
1828 self.buffer_ele = None
1829 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
1830 self.flag =0
1831 self.indicador= 0
1832 self.last_data_ele = None
1833 self.val_mean = None
1834
1835 def update(self, dataOut):
1836
1837 data = {}
1838 meta = {}
1839 if hasattr(dataOut, 'dataPP_POWER'):
1840 factor = 1
1841 if hasattr(dataOut, 'nFFTPoints'):
1842 factor = dataOut.normFactor
1843
1844 if 'pow' in self.attr_data[0].lower():
1845 data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))
1846 else:
1847 data['data'] = getattr(dataOut, self.attr_data[0])/(factor)
1848
1849 if dataOut.mode_op == 'PPI':
1850 self.CODE = 'PPI'
1851 self.title = self.CODE
1852 elif dataOut.mode_op == 'RHI':
1853 self.CODE = 'RHI'
1854 self.title = self.CODE
1855
1856 data['azi'] = dataOut.data_azi
1857 data['ele'] = dataOut.data_ele
1858 data['mode_op'] = dataOut.mode_op
1859
1860 return data, meta
1861
1862 def plot(self):
1863 data = self.data[-1]
1864 r = self.data.yrange
1865 delta_height = r[1]-r[0]
1866 r_mask = numpy.where(r>=0)[0]
1867 self.r_mask =r_mask
1868 r = numpy.arange(len(r_mask))*delta_height
1869 self.y = 2*r
1870
1871 try:
1872 z = data['data'][self.channels[0]][:,r_mask]
1873 except:
1874 z = data['data'][0][:,r_mask]
1875
1876 self.titles = []
1877
1878 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
1879 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
1880 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1881 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1882 print("mode inside plot",self.data['mode_op'],data['mode_op'])
1883 if data['mode_op'] == 'RHI':
1884 try:
1885 if self.data['mode_op'][-2] == 'PPI':
1886 self.ang_min = None
1887 self.ang_max = None
1888 except:
1889 pass
1890 self.ang_min = self.ang_min if self.ang_min else 0
1891 self.ang_max = self.ang_max if self.ang_max else 90
1892 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )
1893 elif data['mode_op'] == 'PPI':
1894 try:
1895 if self.data['mode_op'][-2] == 'RHI':
1896 self.ang_min = None
1897 self.ang_max = None
1898 except:
1899 pass
1900 self.ang_min = self.ang_min if self.ang_min else 0
1901 self.ang_max = self.ang_max if self.ang_max else 360
1902 r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )
1903
1904 self.clear_figures()
1905
1906 for i,ax in enumerate(self.axes):
1907
1908 if ax.firsttime:
1909 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1910 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1911 if data['mode_op'] == 'PPI':
1912 ax.set_theta_direction(-1)
1913 else:
1914 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1915 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1916 if data['mode_op'] == 'PPI':
1917 ax.set_theta_direction(-1)
1918 ax.grid(True)
1919 if data['mode_op'] == 'RHI':
1920 len_aux = int(data['azi'].shape[0]/4)
1921 mean = numpy.mean(data['azi'][len_aux:-len_aux])
1922 if len(self.channels) !=1:
1923 self.titles = ['RHI {} at AZ: {} Channel {}'.format(self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
1924 else:
1925 self.titles = ['RHI {} at AZ: {} Channel {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
1926 elif data['mode_op'] == 'PPI':
1927 len_aux = int(data['ele'].shape[0]/4)
1928 mean = numpy.mean(data['ele'][len_aux:-len_aux])
1929 if len(self.channels) !=1:
1930 self.titles = ['PPI {} at EL: {} Channel {}'.format(self.self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
1931 else:
1932 self.titles = ['PPI {} at EL: {} Channel {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
Show file before | Show file after | Hide comments
@@ -1,838 +1,834
1 1 '''
2 2 Created on Jul 3, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 # SUBCHANNELS EN VEZ DE CHANNELS
7 7 # BENCHMARKS -> PROBLEMAS CON ARCHIVOS GRANDES -> INCONSTANTE EN EL TIEMPO
8 8 # ACTUALIZACION DE VERSION
9 9 # HEADERS
10 10 # MODULO DE ESCRITURA
11 11 # METADATA
12 12
13 13 import os
14 14 import time
15 15 import datetime
16 16 import numpy
17 17 import timeit
18 18 from fractions import Fraction
19 19 from time import time
20 20 from time import sleep
21 21
22 22 import schainpy.admin
23 23 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
24 24 from schainpy.model.data.jrodata import Voltage
25 25 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
26 26
27 27 import pickle
28 28 try:
29 29 import digital_rf
30 30 except:
31 31 pass
32 32
33 33
34 34 class DigitalRFReader(ProcessingUnit):
35 35 '''
36 36 classdocs
37 37 '''
38 38
39 39 def __init__(self):
40 40 '''
41 41 Constructor
42 42 '''
43 43
44 44 ProcessingUnit.__init__(self)
45 45
46 46 self.dataOut = Voltage()
47 47 self.__printInfo = True
48 48 self.__flagDiscontinuousBlock = False
49 49 self.__bufferIndex = 9999999
50 50 self.__codeType = 0
51 51 self.__ippKm = None
52 52 self.__nCode = None
53 53 self.__nBaud = None
54 54 self.__code = None
55 55 self.dtype = None
56 56 self.oldAverage = None
57 57 self.path = None
58 58
59 59 def close(self):
60 60 print('Average of writing to digital rf format is ', self.oldAverage * 1000)
61 61 return
62 62
63 63 def __getCurrentSecond(self):
64 64
65 65 return self.__thisUnixSample / self.__sample_rate
66 66
67 67 thisSecond = property(__getCurrentSecond, "I'm the 'thisSecond' property.")
68 68
69 69 def __setFileHeader(self):
70 70 '''
71 71 In this method will be initialized every parameter of dataOut object (header, no data)
72 72 '''
73 73 ippSeconds = 1.0 * self.__nSamples / self.__sample_rate
74 74 if not self.getByBlock:
75 75 nProfiles = 1.0 / ippSeconds # Number of profiles in one second
76 76 else:
77 77 nProfiles = self.nProfileBlocks # Number of profiles in one block
78 78
79 79 try:
80 80 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
81 81 self.__radarControllerHeader)
82 82 except:
83 83 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
84 84 txA=0,
85 85 txB=0,
86 86 nWindows=1,
87 87 nHeights=self.__nSamples,
88 88 firstHeight=self.__firstHeigth,
89 89 deltaHeight=self.__deltaHeigth,
90 90 codeType=self.__codeType,
91 91 nCode=self.__nCode, nBaud=self.__nBaud,
92 92 code=self.__code)
93 93
94 94 try:
95 95 self.dataOut.systemHeaderObj = SystemHeader(self.__systemHeader)
96 96 except:
97 97 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
98 98 nProfiles=nProfiles,
99 99 nChannels=len(
100 100 self.__channelList),
101 101 adcResolution=14)
102 102 self.dataOut.type = "Voltage"
103 103
104 104 self.dataOut.data = None
105 105
106 106 self.dataOut.dtype = self.dtype
107 107
108 108 # self.dataOut.nChannels = 0
109 109
110 110 # self.dataOut.nHeights = 0
111 111
112 112 self.dataOut.nProfiles = int(nProfiles)
113 113
114 114 self.dataOut.heightList = self.__firstHeigth + \
115 115 numpy.arange(self.__nSamples, dtype=numpy.float) * \
116 116 self.__deltaHeigth
117 117
118 118 #self.dataOut.channelList = list(range(self.__num_subchannels))
119 119 self.dataOut.channelList = list(range(len(self.__channelList)))
120 120 if not self.getByBlock:
121 121
122 122 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
123 123 else:
124 124 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights*self.nProfileBlocks
125 125
126 126 # self.dataOut.channelIndexList = None
127 127
128 128 self.dataOut.flagNoData = True
129 129 if not self.getByBlock:
130 130 self.dataOut.flagDataAsBlock = False
131 131 else:
132 132 self.dataOut.flagDataAsBlock = True
133 133 # Set to TRUE if the data is discontinuous
134 134 self.dataOut.flagDiscontinuousBlock = False
135 135
136 136 self.dataOut.utctime = None
137 137
138 138 # timezone like jroheader, difference in minutes between UTC and localtime
139 139 self.dataOut.timeZone = self.__timezone / 60
140 140
141 141 self.dataOut.dstFlag = 0
142 142
143 143 self.dataOut.errorCount = 0
144 144
145 145 try:
146 146 self.dataOut.nCohInt = self.fixed_metadata_dict.get(
147 147 'nCohInt', self.nCohInt)
148 148
149 149 # asumo que la data esta decodificada
150 150 self.dataOut.flagDecodeData = self.fixed_metadata_dict.get(
151 151 'flagDecodeData', self.flagDecodeData)
152 152
153 153 # asumo que la data esta sin flip
154 154 self.dataOut.flagDeflipData = self.fixed_metadata_dict['flagDeflipData']
155 155
156 156 self.dataOut.flagShiftFFT = self.fixed_metadata_dict['flagShiftFFT']
157 157
158 158 self.dataOut.useLocalTime = self.fixed_metadata_dict['useLocalTime']
159 159 except:
160 160 pass
161 161
162 162 self.dataOut.ippSeconds = ippSeconds
163 163
164 164 # Time interval between profiles
165 165 # self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
166 166
167 167 self.dataOut.frequency = self.__frequency
168 168
169 169 self.dataOut.realtime = self.__online
170 170
171 171 def findDatafiles(self, path, startDate=None, endDate=None):
172 172
173 173 if not os.path.isdir(path):
174 174 return []
175 175
176 176 try:
177 177 digitalReadObj = digital_rf.DigitalRFReader(
178 178 path, load_all_metadata=True)
179 179 except:
180 180 digitalReadObj = digital_rf.DigitalRFReader(path)
181 181
182 182 channelNameList = digitalReadObj.get_channels()
183 183
184 184 if not channelNameList:
185 185 return []
186 186
187 187 metadata_dict = digitalReadObj.get_rf_file_metadata(channelNameList[0])
188 188
189 189 sample_rate = metadata_dict['sample_rate'][0]
190 190
191 191 this_metadata_file = digitalReadObj.get_metadata(channelNameList[0])
192 192
193 193 try:
194 194 timezone = this_metadata_file['timezone'].value
195 195 except:
196 196 timezone = 0
197 197
198 198 startUTCSecond, endUTCSecond = digitalReadObj.get_bounds(
199 199 channelNameList[0]) / sample_rate - timezone
200 200
201 201 startDatetime = datetime.datetime.utcfromtimestamp(startUTCSecond)
202 202 endDatatime = datetime.datetime.utcfromtimestamp(endUTCSecond)
203 203
204 204 if not startDate:
205 205 startDate = startDatetime.date()
206 206
207 207 if not endDate:
208 208 endDate = endDatatime.date()
209 209
210 210 dateList = []
211 211
212 212 thisDatetime = startDatetime
213 213
214 214 while(thisDatetime <= endDatatime):
215 215
216 216 thisDate = thisDatetime.date()
217 217
218 218 if thisDate < startDate:
219 219 continue
220 220
221 221 if thisDate > endDate:
222 222 break
223 223
224 224 dateList.append(thisDate)
225 225 thisDatetime += datetime.timedelta(1)
226 226
227 227 return dateList
228 228
229 229 def setup(self, path=None,
230 230 startDate=None,
231 231 endDate=None,
232 232 startTime=datetime.time(0, 0, 0),
233 233 endTime=datetime.time(23, 59, 59),
234 234 channelList=None,
235 235 nSamples=None,
236 236 online=False,
237 237 delay=60,
238 238 buffer_size=1024,
239 239 ippKm=None,
240 240 nCohInt=1,
241 241 nCode=1,
242 242 nBaud=1,
243 243 flagDecodeData=False,
244 244 code=numpy.ones((1, 1), dtype=numpy.int),
245 245 getByBlock=0,
246 246 nProfileBlocks=1,
247 247 **kwargs):
248 248 '''
249 249 In this method we should set all initial parameters.
250 250
251 251 Inputs:
252 252 path
253 253 startDate
254 254 endDate
255 255 startTime
256 256 endTime
257 257 set
258 258 expLabel
259 259 ext
260 260 online
261 261 delay
262 262 '''
263 263 self.path = path
264 264 self.nCohInt = nCohInt
265 265 self.flagDecodeData = flagDecodeData
266 266 self.i = 0
267 267
268 268 self.getByBlock = getByBlock
269 269 self.nProfileBlocks = nProfileBlocks
270 270 if not os.path.isdir(path):
271 271 raise ValueError("[Reading] Directory %s does not exist" % path)
272 272
273 273 try:
274 274 self.digitalReadObj = digital_rf.DigitalRFReader(
275 275 path, load_all_metadata=True)
276 276 except:
277 277 self.digitalReadObj = digital_rf.DigitalRFReader(path)
278 278
279 279 channelNameList = self.digitalReadObj.get_channels()
280 280
281 281 if not channelNameList:
282 282 raise ValueError("[Reading] Directory %s does not have any files" % path)
283 283
284 284 if not channelList:
285 285 channelList = list(range(len(channelNameList)))
286 286
287 287 ########## Reading metadata ######################
288 288
289 289 top_properties = self.digitalReadObj.get_properties(
290 290 channelNameList[channelList[0]])
291 291
292 292 self.__num_subchannels = top_properties['num_subchannels']
293 293 self.__sample_rate = 1.0 * \
294 294 top_properties['sample_rate_numerator'] / \
295 295 top_properties['sample_rate_denominator']
296 296 # self.__samples_per_file = top_properties['samples_per_file'][0]
297 297 self.__deltaHeigth = 1e6 * 0.15 / self.__sample_rate # why 0.15?
298 298
299 299 this_metadata_file = self.digitalReadObj.get_digital_metadata(
300 300 channelNameList[channelList[0]])
301 301 metadata_bounds = this_metadata_file.get_bounds()
302 302 self.fixed_metadata_dict = this_metadata_file.read(
303 303 metadata_bounds[0])[metadata_bounds[0]] # GET FIRST HEADER
304 304
305 305 try:
306 306 self.__processingHeader = self.fixed_metadata_dict['processingHeader']
307 307 self.__radarControllerHeader = self.fixed_metadata_dict['radarControllerHeader']
308 308 self.__systemHeader = self.fixed_metadata_dict['systemHeader']
309 309 self.dtype = pickle.loads(self.fixed_metadata_dict['dtype'])
310 310 except:
311 311 pass
312 312
313 313 self.__frequency = None
314 314
315 315 self.__frequency = self.fixed_metadata_dict.get('frequency', 1)
316 316
317 317 self.__timezone = self.fixed_metadata_dict.get('timezone', 18000)
318 318
319 319 try:
320 320 nSamples = self.fixed_metadata_dict['nSamples']
321 321 except:
322 322 nSamples = None
323 323
324 324 self.__firstHeigth = 0
325 325
326 326 try:
327 327 codeType = self.__radarControllerHeader['codeType']
328 328 except:
329 329 codeType = 0
330 330
331 331 try:
332 332 if codeType:
333 333 nCode = self.__radarControllerHeader['nCode']
334 334 nBaud = self.__radarControllerHeader['nBaud']
335 335 code = self.__radarControllerHeader['code']
336 336 except:
337 337 pass
338 338
339 339 if not ippKm:
340 340 try:
341 341 # seconds to km
342 342 ippKm = self.__radarControllerHeader['ipp']
343 343 except:
344 344 ippKm = None
345 345 ####################################################
346 346 self.__ippKm = ippKm
347 347 startUTCSecond = None
348 348 endUTCSecond = None
349 349
350 350 if startDate:
351 351 startDatetime = datetime.datetime.combine(startDate, startTime)
352 352 startUTCSecond = (
353 353 startDatetime - datetime.datetime(1970, 1, 1)).total_seconds() + self.__timezone
354 354
355 355 if endDate:
356 356 endDatetime = datetime.datetime.combine(endDate, endTime)
357 357 endUTCSecond = (endDatetime - datetime.datetime(1970,
358 358 1, 1)).total_seconds() + self.__timezone
359 359
360 360
361 361 #print(startUTCSecond,endUTCSecond)
362 362 start_index, end_index = self.digitalReadObj.get_bounds(
363 363 channelNameList[channelList[0]])
364 364
365 365 #print("*****",start_index,end_index)
366 366 if not startUTCSecond:
367 367 startUTCSecond = start_index / self.__sample_rate
368 368
369 369 if start_index > startUTCSecond * self.__sample_rate:
370 370 startUTCSecond = start_index / self.__sample_rate
371 371
372 372 if not endUTCSecond:
373 373 endUTCSecond = end_index / self.__sample_rate
374 374 if end_index < endUTCSecond * self.__sample_rate:
375 375 endUTCSecond = end_index / self.__sample_rate #Check UTC and LT time
376 376 if not nSamples:
377 377 if not ippKm:
378 378 raise ValueError("[Reading] nSamples or ippKm should be defined")
379 379 nSamples = int(ippKm / (1e6 * 0.15 / self.__sample_rate))
380 380
381 381 channelBoundList = []
382 382 channelNameListFiltered = []
383 383
384 384 for thisIndexChannel in channelList:
385 385 thisChannelName = channelNameList[thisIndexChannel]
386 386 start_index, end_index = self.digitalReadObj.get_bounds(
387 387 thisChannelName)
388 388 channelBoundList.append((start_index, end_index))
389 389 channelNameListFiltered.append(thisChannelName)
390 390
391 391 self.profileIndex = 0
392 392 self.i = 0
393 393 self.__delay = delay
394 394
395 395 self.__codeType = codeType
396 396 self.__nCode = nCode
397 397 self.__nBaud = nBaud
398 398 self.__code = code
399 399
400 400 self.__datapath = path
401 401 self.__online = online
402 402 self.__channelList = channelList
403 403 self.__channelNameList = channelNameListFiltered
404 404 self.__channelBoundList = channelBoundList
405 405 self.__nSamples = nSamples
406 406 if self.getByBlock:
407 407 nSamples = nSamples*nProfileBlocks
408 408
409 409
410 410 self.__samples_to_read = int(nSamples) # FIJO: AHORA 40
411 #self.__samples_to_read = int(1000000) # FIJO: AHORA 40
412 411 self.__nChannels = len(self.__channelList)
413 412 #print("------------------------------------------")
414 413 #print("self.__samples_to_read",self.__samples_to_read)
415 414 #print("self.__nSamples",self.__nSamples)
416 415 # son iguales y el buffer_index da 0
417 416 self.__startUTCSecond = startUTCSecond
418 417 self.__endUTCSecond = endUTCSecond
419 418
420 419 self.__timeInterval = 1.0 * self.__samples_to_read / \
421 420 self.__sample_rate # Time interval
422 421
423 422 if online:
424 423 # self.__thisUnixSample = int(endUTCSecond*self.__sample_rate - 4*self.__samples_to_read)
425 424 startUTCSecond = numpy.floor(endUTCSecond)
426 425
427 426 # por que en el otro metodo lo primero q se hace es sumar samplestoread
428 427 self.__thisUnixSample = int(startUTCSecond * self.__sample_rate) - self.__samples_to_read
429 428
430 429 #self.__data_buffer = numpy.zeros(
431 430 # (self.__num_subchannels, self.__samples_to_read), dtype=numpy.complex)
432 431 self.__data_buffer = numpy.zeros((int(len(channelList)), self.__samples_to_read), dtype=numpy.complex)
433 432
434 433
435 434 self.__setFileHeader()
436 435 self.isConfig = True
437 436
438 437 print("[Reading] Digital RF Data was found from %s to %s " % (
439 438 datetime.datetime.utcfromtimestamp(
440 439 self.__startUTCSecond - self.__timezone),
441 440 datetime.datetime.utcfromtimestamp(
442 441 self.__endUTCSecond - self.__timezone)
443 442 ))
444 443
445 444 print("[Reading] Starting process from %s to %s" % (datetime.datetime.utcfromtimestamp(startUTCSecond - self.__timezone),
446 445 datetime.datetime.utcfromtimestamp(
447 446 endUTCSecond - self.__timezone)
448 447 ))
449 448 self.oldAverage = None
450 449 self.count = 0
451 450 self.executionTime = 0
452 451
453 452 def __reload(self):
454 453 # print
455 454 # print "%s not in range [%s, %s]" %(
456 455 # datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
457 456 # datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone),
458 457 # datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone)
459 458 # )
460 459 print("[Reading] reloading metadata ...")
461 460
462 461 try:
463 462 self.digitalReadObj.reload(complete_update=True)
464 463 except:
465 464 self.digitalReadObj = digital_rf.DigitalRFReader(self.path)
466 465
467 466 start_index, end_index = self.digitalReadObj.get_bounds(
468 467 self.__channelNameList[self.__channelList[0]])
469 468
470 469 if start_index > self.__startUTCSecond * self.__sample_rate:
471 470 self.__startUTCSecond = 1.0 * start_index / self.__sample_rate
472 471
473 472 if end_index > self.__endUTCSecond * self.__sample_rate:
474 473 self.__endUTCSecond = 1.0 * end_index / self.__sample_rate
475 474 print()
476 475 print("[Reading] New timerange found [%s, %s] " % (
477 476 datetime.datetime.utcfromtimestamp(
478 477 self.__startUTCSecond - self.__timezone),
479 478 datetime.datetime.utcfromtimestamp(
480 479 self.__endUTCSecond - self.__timezone)
481 480 ))
482 481
483 482 return True
484 483
485 484 return False
486 485
487 486 def timeit(self, toExecute):
488 487 t0 = time.time()
489 488 toExecute()
490 489 self.executionTime = time.time() - t0
491 490 if self.oldAverage is None:
492 491 self.oldAverage = self.executionTime
493 492 self.oldAverage = (self.executionTime + self.count *
494 493 self.oldAverage) / (self.count + 1.0)
495 494 self.count = self.count + 1.0
496 495 return
497 496
498 497 def __readNextBlock(self, seconds=30, volt_scale=1):
499 498 '''
500 499 '''
501 500
502 501 # Set the next data
503 502 self.__flagDiscontinuousBlock = False
504 503 self.__thisUnixSample += self.__samples_to_read
505 504
506 505 if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
507 506 print ("[Reading] There are no more data into selected time-range")
508 507 if self.__online:
509 508 sleep(3)
510 509 self.__reload()
511 510 else:
512 511 return False
513 512
514 513 if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
515 514 return False
516 515 self.__thisUnixSample -= self.__samples_to_read
517 516
518 517 indexChannel = 0
519 518
520 519 dataOk = False
521 520
522 521 for thisChannelName in self.__channelNameList: # TODO VARIOS CHANNELS?
523 522 for indexSubchannel in range(self.__num_subchannels):
524 523 try:
525 524 t0 = time()
526 #print("Unitindex",self.__thisUnixSample)
527 #print("__samples_to_read",self.__samples_to_read)
528 525 result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample,
529 526 self.__samples_to_read,
530 527 thisChannelName, sub_channel=indexSubchannel)
531 528 self.executionTime = time() - t0
532 529 if self.oldAverage is None:
533 530 self.oldAverage = self.executionTime
534 531 self.oldAverage = (
535 532 self.executionTime + self.count * self.oldAverage) / (self.count + 1.0)
536 533 self.count = self.count + 1.0
537 534
538 535 except IOError as e:
539 536 # read next profile
540 537 self.__flagDiscontinuousBlock = True
541 538 print("[Reading] %s" % datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), e)
542 539 break
543 540
544 541 if result.shape[0] != self.__samples_to_read:
545 542 self.__flagDiscontinuousBlock = True
546 543 print("[Reading] %s: Too few samples were found, just %d/%d samples" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
547 544 result.shape[0],
548 545 self.__samples_to_read))
549 546 break
550 547
551 548 self.__data_buffer[indexChannel, :] = result * volt_scale
552 549 indexChannel+=1
553 550
554 551 dataOk = True
555 552
556 553 self.__utctime = self.__thisUnixSample / self.__sample_rate
557 554
558 555 if not dataOk:
559 556 return False
560 557
561 558 print("[Reading] %s: %d samples <> %f sec" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
562 559 self.__samples_to_read,
563 560 self.__timeInterval))
564 561
565 562 self.__bufferIndex = 0
566 563
567 564 return True
568 565
569 566 def __isBufferEmpty(self):
570 567
571 568 return self.__bufferIndex > self.__samples_to_read - self.__nSamples # 40960 - 40
572 569
573 570 def getData(self, seconds=30, nTries=5):
574 571 '''
575 572 This method gets the data from files and put the data into the dataOut object
576 573
577 574 In addition, increase el the buffer counter in one.
578 575
579 576 Return:
580 577 data : retorna un perfil de voltages (alturas * canales) copiados desde el
581 578 buffer. Si no hay mas archivos a leer retorna None.
582 579
583 580 Affected:
584 581 self.dataOut
585 582 self.profileIndex
586 583 self.flagDiscontinuousBlock
587 584 self.flagIsNewBlock
588 585 '''
589 586 #print("getdata")
590 587 err_counter = 0
591 588 self.dataOut.flagNoData = True
592 589
593 590
594 591 if self.__isBufferEmpty():
595 592 #print("hi")
596 593 self.__flagDiscontinuousBlock = False
597 594
598 595 while True:
599 596 if self.__readNextBlock():
600 597 break
601 598 if self.__thisUnixSample > self.__endUTCSecond * self.__sample_rate:
602 599 raise schainpy.admin.SchainError('Error')
603 600 return
604 601
605 602 if self.__flagDiscontinuousBlock:
606 603 raise schainpy.admin.SchainError('discontinuous block found')
607 604 return
608 605
609 606 if not self.__online:
610 607 raise schainpy.admin.SchainError('Online?')
611 608 return
612 609
613 610 err_counter += 1
614 611 if err_counter > nTries:
615 612 raise schainpy.admin.SchainError('Max retrys reach')
616 613 return
617 614
618 615 print('[Reading] waiting %d seconds to read a new block' % seconds)
619 616 sleep(seconds)
620 617
621 618
622 619 if not self.getByBlock:
623 620
624 621 #print("self.__bufferIndex",self.__bufferIndex)# este valor siempre es cero aparentemente
625 622 self.dataOut.data = self.__data_buffer[:, self.__bufferIndex:self.__bufferIndex + self.__nSamples]
626 623 self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
627 624 self.dataOut.flagNoData = False
628 625 self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
629 626 self.dataOut.profileIndex = self.profileIndex
630 627
631 628 self.__bufferIndex += self.__nSamples
632 629 self.profileIndex += 1
633 630
634 631 if self.profileIndex == self.dataOut.nProfiles:
635 632 self.profileIndex = 0
636 633 else:
637 634 # ojo debo anadir el readNextBLock y el __isBufferEmpty(
638 635 self.dataOut.flagNoData = False
639 636 buffer = self.__data_buffer[:,self.__bufferIndex:self.__bufferIndex + self.__samples_to_read]
640 #print("test",self.__bufferIndex)
641 637 buffer = buffer.reshape((self.__nChannels, self.nProfileBlocks, int(self.__samples_to_read/self.nProfileBlocks)))
642 638 self.dataOut.nProfileBlocks = self.nProfileBlocks
643 639 self.dataOut.data = buffer
644 640 self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
645 641 self.profileIndex += self.__samples_to_read
646 642 self.__bufferIndex += self.__samples_to_read
647 643 self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
648 644 return True
649 645
650 646
651 647 def printInfo(self):
652 648 '''
653 649 '''
654 650 if self.__printInfo == False:
655 651 return
656 652
657 653 # self.systemHeaderObj.printInfo()
658 654 # self.radarControllerHeaderObj.printInfo()
659 655
660 656 self.__printInfo = False
661 657
662 658 def printNumberOfBlock(self):
663 659 '''
664 660 '''
665 661 return
666 662 # print self.profileIndex
667 663
668 664 def run(self, **kwargs):
669 665 '''
670 666 This method will be called many times so here you should put all your code
671 667 '''
672 668
673 669 if not self.isConfig:
674 670 self.setup(**kwargs)
675 671
676 672 self.getData(seconds=self.__delay)
677 673
678 674 return
679 675
680 676 @MPDecorator
681 677 class DigitalRFWriter(Operation):
682 678 '''
683 679 classdocs
684 680 '''
685 681
686 682 def __init__(self, **kwargs):
687 683 '''
688 684 Constructor
689 685 '''
690 686 Operation.__init__(self, **kwargs)
691 687 self.metadata_dict = {}
692 688 self.dataOut = None
693 689 self.dtype = None
694 690 self.oldAverage = 0
695 691
696 692 def setHeader(self):
697 693
698 694 self.metadata_dict['frequency'] = self.dataOut.frequency
699 695 self.metadata_dict['timezone'] = self.dataOut.timeZone
700 696 self.metadata_dict['dtype'] = pickle.dumps(self.dataOut.dtype)
701 697 self.metadata_dict['nProfiles'] = self.dataOut.nProfiles
702 698 self.metadata_dict['heightList'] = self.dataOut.heightList
703 699 self.metadata_dict['channelList'] = self.dataOut.channelList
704 700 self.metadata_dict['flagDecodeData'] = self.dataOut.flagDecodeData
705 701 self.metadata_dict['flagDeflipData'] = self.dataOut.flagDeflipData
706 702 self.metadata_dict['flagShiftFFT'] = self.dataOut.flagShiftFFT
707 703 self.metadata_dict['useLocalTime'] = self.dataOut.useLocalTime
708 704 self.metadata_dict['nCohInt'] = self.dataOut.nCohInt
709 705 self.metadata_dict['type'] = self.dataOut.type
710 706 self.metadata_dict['flagDataAsBlock']= getattr(
711 707 self.dataOut, 'flagDataAsBlock', None) # chequear
712 708
713 709 def setup(self, dataOut, path, frequency, fileCadence, dirCadence, metadataCadence, set=0, metadataFile='metadata', ext='.h5'):
714 710 '''
715 711 In this method we should set all initial parameters.
716 712 Input:
717 713 dataOut: Input data will also be outputa data
718 714 '''
719 715 self.setHeader()
720 716 self.__ippSeconds = dataOut.ippSeconds
721 717 self.__deltaH = dataOut.getDeltaH()
722 718 self.__sample_rate = 1e6 * 0.15 / self.__deltaH
723 719 self.__dtype = dataOut.dtype
724 720 if len(dataOut.dtype) == 2:
725 721 self.__dtype = dataOut.dtype[0]
726 722 self.__nSamples = dataOut.systemHeaderObj.nSamples
727 723 self.__nProfiles = dataOut.nProfiles
728 724
729 725 if self.dataOut.type != 'Voltage':
730 726 raise 'Digital RF cannot be used with this data type'
731 727 self.arr_data = numpy.ones((1, dataOut.nFFTPoints * len(
732 728 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
733 729 else:
734 730 self.arr_data = numpy.ones((self.__nSamples, len(
735 731 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
736 732
737 733 file_cadence_millisecs = 1000
738 734
739 735 sample_rate_fraction = Fraction(self.__sample_rate).limit_denominator()
740 736 sample_rate_numerator = int(sample_rate_fraction.numerator)
741 737 sample_rate_denominator = int(sample_rate_fraction.denominator)
742 738 start_global_index = dataOut.utctime * self.__sample_rate
743 739
744 740 uuid = 'prueba'
745 741 compression_level = 0
746 742 checksum = False
747 743 is_complex = True
748 744 num_subchannels = len(dataOut.channelList)
749 745 is_continuous = True
750 746 marching_periods = False
751 747
752 748 self.digitalWriteObj = digital_rf.DigitalRFWriter(path, self.__dtype, dirCadence,
753 749 fileCadence, start_global_index,
754 750 sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum,
755 751 is_complex, num_subchannels, is_continuous, marching_periods)
756 752 metadata_dir = os.path.join(path, 'metadata')
757 753 os.system('mkdir %s' % (metadata_dir))
758 754 self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, # 236, file_cadence_millisecs / 1000
759 755 sample_rate_numerator, sample_rate_denominator,
760 756 metadataFile)
761 757 self.isConfig = True
762 758 self.currentSample = 0
763 759 self.oldAverage = 0
764 760 self.count = 0
765 761 return
766 762
767 763 def writeMetadata(self):
768 764 start_idx = self.__sample_rate * self.dataOut.utctime
769 765
770 766 self.metadata_dict['processingHeader'] = self.dataOut.processingHeaderObj.getAsDict(
771 767 )
772 768 self.metadata_dict['radarControllerHeader'] = self.dataOut.radarControllerHeaderObj.getAsDict(
773 769 )
774 770 self.metadata_dict['systemHeader'] = self.dataOut.systemHeaderObj.getAsDict(
775 771 )
776 772 self.digitalMetadataWriteObj.write(start_idx, self.metadata_dict)
777 773 return
778 774
779 775 def timeit(self, toExecute):
780 776 t0 = time()
781 777 toExecute()
782 778 self.executionTime = time() - t0
783 779 if self.oldAverage is None:
784 780 self.oldAverage = self.executionTime
785 781 self.oldAverage = (self.executionTime + self.count *
786 782 self.oldAverage) / (self.count + 1.0)
787 783 self.count = self.count + 1.0
788 784 return
789 785
790 786 def writeData(self):
791 787 if self.dataOut.type != 'Voltage':
792 788 raise 'Digital RF cannot be used with this data type'
793 789 for channel in self.dataOut.channelList:
794 790 for i in range(self.dataOut.nFFTPoints):
795 791 self.arr_data[1][channel * self.dataOut.nFFTPoints +
796 792 i]['r'] = self.dataOut.data[channel][i].real
797 793 self.arr_data[1][channel * self.dataOut.nFFTPoints +
798 794 i]['i'] = self.dataOut.data[channel][i].imag
799 795 else:
800 796 for i in range(self.dataOut.systemHeaderObj.nSamples):
801 797 for channel in self.dataOut.channelList:
802 798 self.arr_data[i][channel]['r'] = self.dataOut.data[channel][i].real
803 799 self.arr_data[i][channel]['i'] = self.dataOut.data[channel][i].imag
804 800
805 801 def f(): return self.digitalWriteObj.rf_write(self.arr_data)
806 802 self.timeit(f)
807 803
808 804 return
809 805
810 806 def run(self, dataOut, frequency=49.92e6, path=None, fileCadence=1000, dirCadence=36000, metadataCadence=1, **kwargs):
811 807 '''
812 808 This method will be called many times so here you should put all your code
813 809 Inputs:
814 810 dataOut: object with the data
815 811 '''
816 812 # print dataOut.__dict__
817 813 self.dataOut = dataOut
818 814 if not self.isConfig:
819 815 self.setup(dataOut, path, frequency, fileCadence,
820 816 dirCadence, metadataCadence, **kwargs)
821 817 self.writeMetadata()
822 818
823 819 self.writeData()
824 820
825 821 ## self.currentSample += 1
826 822 # if self.dataOut.flagDataAsBlock or self.currentSample == 1:
827 823 # self.writeMetadata()
828 824 ## if self.currentSample == self.__nProfiles: self.currentSample = 0
829 825
830 826 return dataOut# en la version 2.7 no aparece este return
831 827
832 828 def close(self):
833 829 print('[Writing] - Closing files ')
834 830 print('Average of writing to digital rf format is ', self.oldAverage * 1000)
835 831 try:
836 832 self.digitalWriteObj.close()
837 833 except:
838 834 pass
Show file before | Show file after | Hide comments
@@ -1,731 +1,735
1 1 import os
2 2 import time
3 3 import datetime
4 4
5 5 import numpy
6 6 import h5py
7 7
8 8 import schainpy.admin
9 9 from schainpy.model.data.jrodata import *
10 10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 11 from schainpy.model.io.jroIO_base import *
12 12 from schainpy.utils import log
13 13
14 14
15 15 class HDFReader(Reader, ProcessingUnit):
16 16 """Processing unit to read HDF5 format files
17 17
18 18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 20 attributes.
21 21 It is possible to read any HDF5 file by given the structure in the `description`
22 22 parameter, also you can add extra values to metadata with the parameter `extras`.
23 23
24 24 Parameters:
25 25 -----------
26 26 path : str
27 27 Path where files are located.
28 28 startDate : date
29 29 Start date of the files
30 30 endDate : list
31 31 End date of the files
32 32 startTime : time
33 33 Start time of the files
34 34 endTime : time
35 35 End time of the files
36 36 description : dict, optional
37 37 Dictionary with the description of the HDF5 file
38 38 extras : dict, optional
39 39 Dictionary with extra metadata to be be added to `dataOut`
40 40
41 41 Examples
42 42 --------
43 43
44 44 desc = {
45 45 'Data': {
46 46 'data_output': ['u', 'v', 'w'],
47 47 'utctime': 'timestamps',
48 48 } ,
49 49 'Metadata': {
50 50 'heightList': 'heights'
51 51 }
52 52 }
53 53
54 54 desc = {
55 55 'Data': {
56 56 'data_output': 'winds',
57 57 'utctime': 'timestamps'
58 58 },
59 59 'Metadata': {
60 60 'heightList': 'heights'
61 61 }
62 62 }
63 63
64 64 extras = {
65 65 'timeZone': 300
66 66 }
67 67
68 68 reader = project.addReadUnit(
69 69 name='HDFReader',
70 70 path='/path/to/files',
71 71 startDate='2019/01/01',
72 72 endDate='2019/01/31',
73 73 startTime='00:00:00',
74 74 endTime='23:59:59',
75 75 # description=json.dumps(desc),
76 76 # extras=json.dumps(extras),
77 77 )
78 78
79 79 """
80 80
81 81 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
82 82
83 83 def __init__(self):
84 84 ProcessingUnit.__init__(self)
85 85 self.dataOut = Parameters()
86 86 self.ext = ".hdf5"
87 87 self.optchar = "D"
88 88 self.meta = {}
89 89 self.data = {}
90 90 self.open_file = h5py.File
91 91 self.open_mode = 'r'
92 92 self.description = {}
93 93 self.extras = {}
94 94 self.filefmt = "*%Y%j***"
95 95 self.folderfmt = "*%Y%j"
96 96 self.utcoffset = 0
97 97
98 98 def setup(self, **kwargs):
99 99
100 100 self.set_kwargs(**kwargs)
101 101 if not self.ext.startswith('.'):
102 102 self.ext = '.{}'.format(self.ext)
103 103
104 104 if self.online:
105 105 log.log("Searching files in online mode...", self.name)
106 106
107 107 for nTries in range(self.nTries):
108 108 fullpath = self.searchFilesOnLine(self.path, self.startDate,
109 109 self.endDate, self.expLabel, self.ext, self.walk,
110 110 self.filefmt, self.folderfmt)
111 111 try:
112 112 fullpath = next(fullpath)
113 113 except:
114 114 fullpath = None
115 115
116 116 if fullpath:
117 117 break
118 118
119 119 log.warning(
120 120 'Waiting {} sec for a valid file in {}: try {} ...'.format(
121 121 self.delay, self.path, nTries + 1),
122 122 self.name)
123 123 time.sleep(self.delay)
124 124
125 125 if not(fullpath):
126 126 raise schainpy.admin.SchainError(
127 127 'There isn\'t any valid file in {}'.format(self.path))
128 128
129 129 pathname, filename = os.path.split(fullpath)
130 130 self.year = int(filename[1:5])
131 131 self.doy = int(filename[5:8])
132 132 self.set = int(filename[8:11]) - 1
133 133 else:
134 134 log.log("Searching files in {}".format(self.path), self.name)
135 135 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
136 136 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
137 137
138 138 self.setNextFile()
139 139
140 140 return
141 141
142 142 def readFirstHeader(self):
143 143 '''Read metadata and data'''
144 144
145 145 self.__readMetadata()
146 146 self.__readData()
147 147 self.__setBlockList()
148 148
149 149 if 'type' in self.meta:
150 150 self.dataOut = eval(self.meta['type'])()
151 151
152 152 for attr in self.meta:
153 153 setattr(self.dataOut, attr, self.meta[attr])
154 154
155 155 self.blockIndex = 0
156 156
157 157 return
158 158
159 159 def __setBlockList(self):
160 160 '''
161 161 Selects the data within the times defined
162 162
163 163 self.fp
164 164 self.startTime
165 165 self.endTime
166 166 self.blockList
167 167 self.blocksPerFile
168 168
169 169 '''
170 170
171 171 startTime = self.startTime
172 172 endTime = self.endTime
173 173 thisUtcTime = self.data['utctime'] + self.utcoffset
174 174 self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
175 175 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
176 176
177 177 thisDate = thisDatetime.date()
178 178 thisTime = thisDatetime.time()
179 179
180 180 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
181 181 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
182 182
183 183 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
184 184
185 185 self.blockList = ind
186 186 self.blocksPerFile = len(ind)
187 187 return
188 188
189 189 def __readMetadata(self):
190 190 '''
191 191 Reads Metadata
192 192 '''
193 193
194 194 meta = {}
195 195
196 196 if self.description:
197 197 for key, value in self.description['Metadata'].items():
198 198 meta[key] = self.fp[value][()]
199 199 else:
200 200 grp = self.fp['Metadata']
201 201 for name in grp:
202 202 meta[name] = grp[name][()]
203 203
204 204 if self.extras:
205 205 for key, value in self.extras.items():
206 206 meta[key] = value
207 207 self.meta = meta
208 208
209 209 return
210 210
211 211 def __readData(self):
212 212
213 213 data = {}
214 214
215 215 if self.description:
216 216 for key, value in self.description['Data'].items():
217 217 if isinstance(value, str):
218 218 if isinstance(self.fp[value], h5py.Dataset):
219 219 data[key] = self.fp[value][()]
220 220 elif isinstance(self.fp[value], h5py.Group):
221 221 array = []
222 222 for ch in self.fp[value]:
223 223 array.append(self.fp[value][ch][()])
224 224 data[key] = numpy.array(array)
225 225 elif isinstance(value, list):
226 226 array = []
227 227 for ch in value:
228 228 array.append(self.fp[ch][()])
229 229 data[key] = numpy.array(array)
230 230 else:
231 231 grp = self.fp['Data']
232 232 for name in grp:
233 233 if isinstance(grp[name], h5py.Dataset):
234 234 array = grp[name][()]
235 235 elif isinstance(grp[name], h5py.Group):
236 236 array = []
237 237 for ch in grp[name]:
238 238 array.append(grp[name][ch][()])
239 239 array = numpy.array(array)
240 240 else:
241 241 log.warning('Unknown type: {}'.format(name))
242 242
243 243 if name in self.description:
244 244 key = self.description[name]
245 245 else:
246 246 key = name
247 247 data[key] = array
248 248
249 249 self.data = data
250 250 return
251 251
252 252 def getData(self):
253 253
254 254 for attr in self.data:
255 255 if self.data[attr].ndim == 1:
256 256 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
257 257 else:
258 258 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
259 259
260 260 self.dataOut.flagNoData = False
261 261 self.blockIndex += 1
262 262
263 263 log.log("Block No. {}/{} -> {}".format(
264 264 self.blockIndex,
265 265 self.blocksPerFile,
266 266 self.dataOut.datatime.ctime()), self.name)
267 267
268 268 return
269 269
270 270 def run(self, **kwargs):
271 271
272 272 if not(self.isConfig):
273 273 self.setup(**kwargs)
274 274 self.isConfig = True
275 275
276 276 if self.blockIndex == self.blocksPerFile:
277 277 self.setNextFile()
278 278
279 279 self.getData()
280 280
281 281 return
282 282
283 283 @MPDecorator
284 284 class HDFWriter(Operation):
285 285 """Operation to write HDF5 files.
286 286
287 287 The HDF5 file contains by default two groups Data and Metadata where
288 288 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
289 289 parameters, data attributes are normaly time dependent where the metadata
290 290 are not.
291 291 It is possible to customize the structure of the HDF5 file with the
292 292 optional description parameter see the examples.
293 293
294 294 Parameters:
295 295 -----------
296 296 path : str
297 297 Path where files will be saved.
298 298 blocksPerFile : int
299 299 Number of blocks per file
300 300 metadataList : list
301 301 List of the dataOut attributes that will be saved as metadata
302 302 dataList : int
303 303 List of the dataOut attributes that will be saved as data
304 304 setType : bool
305 305 If True the name of the files corresponds to the timestamp of the data
306 306 description : dict, optional
307 307 Dictionary with the desired description of the HDF5 file
308 308
309 309 Examples
310 310 --------
311 311
312 312 desc = {
313 313 'data_output': {'winds': ['z', 'w', 'v']},
314 314 'utctime': 'timestamps',
315 315 'heightList': 'heights'
316 316 }
317 317 desc = {
318 318 'data_output': ['z', 'w', 'v'],
319 319 'utctime': 'timestamps',
320 320 'heightList': 'heights'
321 321 }
322 322 desc = {
323 323 'Data': {
324 324 'data_output': 'winds',
325 325 'utctime': 'timestamps'
326 326 },
327 327 'Metadata': {
328 328 'heightList': 'heights'
329 329 }
330 330 }
331 331
332 332 writer = proc_unit.addOperation(name='HDFWriter')
333 333 writer.addParameter(name='path', value='/path/to/file')
334 334 writer.addParameter(name='blocksPerFile', value='32')
335 335 writer.addParameter(name='metadataList', value='heightList,timeZone')
336 336 writer.addParameter(name='dataList',value='data_output,utctime')
337 337 # writer.addParameter(name='description',value=json.dumps(desc))
338 338
339 339 """
340 340
341 341 ext = ".hdf5"
342 342 optchar = "D"
343 343 filename = None
344 344 path = None
345 345 setFile = None
346 346 fp = None
347 347 firsttime = True
348 348 #Configurations
349 349 blocksPerFile = None
350 350 blockIndex = None
351 351 dataOut = None
352 352 #Data Arrays
353 353 dataList = None
354 354 metadataList = None
355 355 currentDay = None
356 356 lastTime = None
357 357 last_Azipos = None
358 358 last_Elepos = None
359 359 mode = None
360 360 #-----------------------
361 361 Typename = None
362 362
363 363
364 364
365 365 def __init__(self):
366 366
367 367 Operation.__init__(self)
368 368 return
369 369
370 370
371 371 def set_kwargs(self, **kwargs):
372 372
373 373 for key, value in kwargs.items():
374 374 setattr(self, key, value)
375 375
376 376 def set_kwargs_obj(self,obj, **kwargs):
377 377
378 378 for key, value in kwargs.items():
379 379 setattr(obj, key, value)
380 380
381 381 def generalFlag(self):
382 382 ####rint("GENERALFLAG")
383 383 if self.mode== "weather":
384 384 if self.last_Azipos == None:
385 385 tmp = self.dataOut.azimuth
386 386 ####print("ang azimuth writer",tmp)
387 387 self.last_Azipos = tmp
388 388 flag = False
389 389 return flag
390 390 ####print("ang_azimuth writer",self.dataOut.azimuth)
391 391 result = self.dataOut.azimuth - self.last_Azipos
392 392 self.last_Azipos = self.dataOut.azimuth
393 393 if result<0:
394 394 flag = True
395 395 return flag
396 396
397 397 def generalFlag_vRF(self):
398 398 ####rint("GENERALFLAG")
399 399
400 400 try:
401 401 self.dataOut.flagBlock360Done
402 402 return self.dataOut.flagBlock360Done
403 403 except:
404 404 return 0
405 405
406 406
407 407 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None,type_data=None,**kwargs):
408 408 self.path = path
409 409 self.blocksPerFile = blocksPerFile
410 410 self.metadataList = metadataList
411 411 self.dataList = [s.strip() for s in dataList]
412 412 self.setType = setType
413 413 if self.mode == "weather":
414 414 self.setType = "weather"
415 415 self.set_kwargs(**kwargs)
416 416 self.set_kwargs_obj(self.dataOut,**kwargs)
417 417
418 418
419 419 self.description = description
420 420 self.type_data=type_data
421 421
422 422 if self.metadataList is None:
423 423 self.metadataList = self.dataOut.metadata_list
424 424
425 425 tableList = []
426 426 dsList = []
427 427
428 428 for i in range(len(self.dataList)):
429 429 dsDict = {}
430 430 if hasattr(self.dataOut, self.dataList[i]):
431 431 dataAux = getattr(self.dataOut, self.dataList[i])
432 432 dsDict['variable'] = self.dataList[i]
433 433 else:
434 434 log.warning('Attribute {} not found in dataOut', self.name)
435 435 continue
436 436
437 437 if dataAux is None:
438 438 continue
439 439 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
440 440 dsDict['nDim'] = 0
441 441 else:
442 442 dsDict['nDim'] = len(dataAux.shape)
443 443 dsDict['shape'] = dataAux.shape
444 444 dsDict['dsNumber'] = dataAux.shape[0]
445 445 dsDict['dtype'] = dataAux.dtype
446 446 dsList.append(dsDict)
447 447
448 448 self.dsList = dsList
449 449 self.currentDay = self.dataOut.datatime.date()
450 450
451 451 def timeFlag(self):
452 452 currentTime = self.dataOut.utctime
453 453 timeTuple = time.localtime(currentTime)
454 454 dataDay = timeTuple.tm_yday
455 455
456 456 if self.lastTime is None:
457 457 self.lastTime = currentTime
458 458 self.currentDay = dataDay
459 459 return False
460 460
461 461 timeDiff = currentTime - self.lastTime
462 462
463 463 #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
464 464 if dataDay != self.currentDay:
465 465 self.currentDay = dataDay
466 466 return True
467 467 elif timeDiff > 3*60*60:
468 468 self.lastTime = currentTime
469 469 return True
470 470 else:
471 471 self.lastTime = currentTime
472 472 return False
473 473
474 474 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
475 475 dataList=[], setType=None, description={},mode= None,type_data=None,Reset = False,**kwargs):
476 476
477 477 if Reset:
478 478 self.isConfig = False
479 479 self.closeFile()
480 480 self.lastTime = None
481 481 self.blockIndex = 0
482 482
483 483 self.dataOut = dataOut
484 484 self.mode = mode
485 485 self.var = dataList[0]
486 486
487 487 if not(self.isConfig):
488 488 self.setup(path=path, blocksPerFile=blocksPerFile,
489 489 metadataList=metadataList, dataList=dataList,
490 490 setType=setType, description=description,type_data=type_data,**kwargs)
491 491
492 492 self.isConfig = True
493 493 self.setNextFile()
494 494
495 495 self.putData()
496 496 return
497 497
498 498 def setNextFile(self):
499 499 ###print("HELLO WORLD--------------------------------")
500 500 ext = self.ext
501 501 path = self.path
502 502 setFile = self.setFile
503 503 type_data = self.type_data
504 504
505 505 timeTuple = time.localtime(self.dataOut.utctime)
506 506 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
507 507 fullpath = os.path.join(path, subfolder)
508 508
509 509 if os.path.exists(fullpath):
510 510 filesList = os.listdir(fullpath)
511 511 filesList = [k for k in filesList if k.startswith(self.optchar)]
512 512 if len( filesList ) > 0:
513 513 filesList = sorted(filesList, key=str.lower)
514 514 filen = filesList[-1]
515 515 # el filename debera tener el siguiente formato
516 516 # 0 1234 567 89A BCDE (hex)
517 517 # x YYYY DDD SSS .ext
518 518 if isNumber(filen[8:11]):
519 519 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
520 520 else:
521 521 setFile = -1
522 522 else:
523 523 setFile = -1 #inicializo mi contador de seteo
524 524 else:
525 525 os.makedirs(fullpath)
526 526 setFile = -1 #inicializo mi contador de seteo
527 527
528 528 ###print("**************************",self.setType)
529 529 if self.setType is None:
530 530 setFile += 1
531 531 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
532 532 timeTuple.tm_year,
533 533 timeTuple.tm_yday,
534 534 setFile,
535 535 ext )
536 536 elif self.setType == "weather":
537 537
538 538 if self.var.lower() == 'Zdb'.lower():
539 539 wr_type = 'Z'
540 540 elif self.var.lower() == 'Zdb_D'.lower():
541 541 wr_type = 'D'
542 542 elif self.var.lower() == 'PhiD_P'.lower():
543 543 wr_type = 'P'
544 544 elif self.var.lower() == 'RhoHV_R'.lower():
545 545 wr_type = 'R'
546 546 elif self.var.lower() == 'velRadial_V'.lower():
547 547 wr_type = 'V'
548 548 elif self.var.lower() == 'Sigmav_W'.lower():
549 549 wr_type = 'S'
550 550 elif self.var.lower() == 'dataPP_POWER'.lower():
551 551 wr_type = 'Pow'
552 552 elif self.var.lower() == 'dataPP_DOP'.lower():
553 553 wr_type = 'Dop'
554 554
555 555
556 556 #Z_SOPHy_El10.0_20200505_14:02:15.h5
557 557 #Z_SOPHy_Az40.0_20200505_14:02:15.h5
558 558 if self.dataOut.flagMode == 1: #'AZI' #PPI
559 559 ang_type = 'El'
560 ang_ = round(numpy.mean(self.dataOut.data_ele),1)
560 len_aux = int(self.dataOut.data_ele.shape[0]/4)
561 mean = numpy.mean(self.dataOut.data_ele[len_aux:-len:aux])
562 ang_ = round(mean,1)
561 563 elif self.dataOut.flagMode == 0: #'ELE' #RHI
562 564 ang_type = 'Az'
563 ang_ = round(numpy.mean(self.dataOut.data_azi),1)
565 len_aux = int(self.dataOut.data_azi.shape[0]/4)
566 mean = numpy.mean(self.dataOut.data_azi[len_aux:-len:aux])
567 ang_ = round(mean,1)
564 568
565 569 file = '%s%s%s%2.1f%s%2.2d%2.2d%2.2d%s%2.2d%2.2d%2.2d%s' % (wr_type,
566 570 '_SOPHy_',
567 571 ang_type,
568 572 ang_,
569 573 '_',
570 574 timeTuple.tm_year,
571 575 timeTuple.tm_mon,
572 576 timeTuple.tm_mday,
573 577 '_',
574 578 timeTuple.tm_hour,
575 579 timeTuple.tm_min,
576 580 timeTuple.tm_sec,
577 581 ext )
578 582
579 583 else:
580 584 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
581 585 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
582 586 timeTuple.tm_year,
583 587 timeTuple.tm_yday,
584 588 setFile,
585 589 ext )
586 590
587 591 self.filename = os.path.join( path, subfolder, file )
588 592
589 593 #Setting HDF5 File
590 print("filename",self.filename)
594 #print("filename",self.filename)
591 595 self.fp = h5py.File(self.filename, 'w')
592 596 #write metadata
593 597 self.writeMetadata(self.fp)
594 598 #Write data
595 599 self.writeData(self.fp)
596 600
597 601 def getLabel(self, name, x=None):
598 602
599 603 if x is None:
600 604 if 'Data' in self.description:
601 605 data = self.description['Data']
602 606 if 'Metadata' in self.description:
603 607 data.update(self.description['Metadata'])
604 608 else:
605 609 data = self.description
606 610 if name in data:
607 611 if isinstance(data[name], str):
608 612 return data[name]
609 613 elif isinstance(data[name], list):
610 614 return None
611 615 elif isinstance(data[name], dict):
612 616 for key, value in data[name].items():
613 617 return key
614 618 return name
615 619 else:
616 620 if 'Metadata' in self.description:
617 621 meta = self.description['Metadata']
618 622 else:
619 623 meta = self.description
620 624 if name in meta:
621 625 if isinstance(meta[name], list):
622 626 return meta[name][x]
623 627 elif isinstance(meta[name], dict):
624 628 for key, value in meta[name].items():
625 629 return value[x]
626 630 if 'cspc' in name:
627 631 return 'pair{:02d}'.format(x)
628 632 else:
629 633 return 'channel{:02d}'.format(x)
630 634
631 635 def writeMetadata(self, fp):
632 636
633 637 if self.description:
634 638 if 'Metadata' in self.description:
635 639 grp = fp.create_group('Metadata')
636 640 else:
637 641 grp = fp
638 642 else:
639 643 grp = fp.create_group('Metadata')
640 644
641 645 for i in range(len(self.metadataList)):
642 646 if not hasattr(self.dataOut, self.metadataList[i]):
643 647 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
644 648 continue
645 649 value = getattr(self.dataOut, self.metadataList[i])
646 650 if isinstance(value, bool):
647 651 if value is True:
648 652 value = 1
649 653 else:
650 654 value = 0
651 655 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
652 656 return
653 657
654 658 def writeData(self, fp):
655 659
656 660 if self.description:
657 661 if 'Data' in self.description:
658 662 grp = fp.create_group('Data')
659 663 else:
660 664 grp = fp
661 665 else:
662 666 grp = fp.create_group('Data')
663 667
664 668 dtsets = []
665 669 data = []
666 670
667 671 for dsInfo in self.dsList:
668 672 if dsInfo['nDim'] == 0:
669 673 ds = grp.create_dataset(
670 674 self.getLabel(dsInfo['variable']),
671 675 (self.blocksPerFile, ),
672 676 chunks=True,
673 677 dtype=numpy.float64)
674 678 dtsets.append(ds)
675 679 data.append((dsInfo['variable'], -1))
676 680 else:
677 681 label = self.getLabel(dsInfo['variable'])
678 682 if label is not None:
679 683 sgrp = grp.create_group(label)
680 684 else:
681 685 sgrp = grp
682 686 for i in range(dsInfo['dsNumber']):
683 687 ds = sgrp.create_dataset(
684 688 self.getLabel(dsInfo['variable'], i),
685 689 (self.blocksPerFile, ) + dsInfo['shape'][1:],
686 690 chunks=True,
687 691 dtype=dsInfo['dtype'])
688 692 dtsets.append(ds)
689 693 data.append((dsInfo['variable'], i))
690 694 fp.flush()
691 695
692 696 log.log('Creating file: {}'.format(fp.filename), self.name)
693 697
694 698 self.ds = dtsets
695 699 self.data = data
696 700 self.firsttime = True
697 701 self.blockIndex = 0
698 702 return
699 703
700 704 def putData(self):
701 705
702 706 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():# or self.generalFlag_vRF():
703 707 self.closeFile()
704 708 self.setNextFile()
705 709
706 710 for i, ds in enumerate(self.ds):
707 711 attr, ch = self.data[i]
708 712 if ch == -1:
709 713 ds[self.blockIndex] = getattr(self.dataOut, attr)
710 714 else:
711 715 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
712 716
713 717 self.fp.flush()
714 718 self.blockIndex += 1
715 719 log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
716 720
717 721 return
718 722
719 723 def closeFile(self):
720 724
721 725 if self.blockIndex != self.blocksPerFile:
722 726 for ds in self.ds:
723 727 ds.resize(self.blockIndex, axis=0)
724 728
725 729 if self.fp:
726 730 self.fp.flush()
727 731 self.fp.close()
728 732
729 733 def close(self):
730 734
731 735 self.closeFile()
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file chmod 100755 => 100644
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,164 +1,164
1 1 # SOPHY PROC script
2 2 import os, sys, json, argparse
3 3 import datetime
4 4 import time
5 5
6 6 PATH = '/DATA_RM/DATA'
7 7 # PATH = '/Users/jespinoza/workspace/data/'
8 8
9 9 PARAM = {
10 'P': {'name': 'dataPP_POWER', 'zmin': 35, 'zmax': 60, 'colormap': 'viridis', 'label': 'Power', 'cb_label': 'dB'},
10 'P': {'name': 'dataPP_POWER', 'zmin': 35, 'zmax': 60, 'colormap': 'jet', 'label': 'Power', 'cb_label': 'dB'},
11 11 'V': {'name': 'dataPP_DOP', 'zmin': -20, 'zmax': 20, 'colormap': 'seismic', 'label': 'Velocity', 'cb_label': 'm/s'},
12 12 'RH': {'name': 'RhoHV_R', 'zmin': 0, 'zmax': 1, 'colormap': 'jet', 'label': 'CoeficienteCorrelacion', 'cb_label': '*'},
13 13 'FD': {'name': 'PhiD_P', 'zmin': -180, 'zmax': 180, 'colormap': 'RdBu_r', 'label': 'Fase Diferencial', 'cb_label': 'º'},
14 14 'ZD': {'name': 'Zdb_D', 'zmin': -20, 'zmax': 80, 'colormap': 'viridis', 'label': 'ReflectividadDiferencial', 'cb_label': 'dB'},
15 'Z': {'name': 'Zdb', 'zmin': 100, 'zmax': 200, 'colormap': 'viridis', 'label': 'Reflectividad', 'cb_label': 'dB'},
15 'Z': {'name': 'Zdb', 'zmin': -20, 'zmax': 60, 'colormap': 'viridis', 'label': 'Reflectividad', 'cb_label': 'dB'},
16 16 'W': {'name': 'Sigmav_W', 'zmin': -20, 'zmax': 60, 'colormap': 'viridis', 'label': 'AnchoEspectral', 'cb_label': 'hz'}
17 17 }
18 18 #Z,ZD 'mm^6/m^3'
19 19
20 20 PATH = '/home/soporte/Downloads/data_WR_RHI'
21 21
22 22
23 23
24 24 def main(args):
25 25
26 26 experiment = args.experiment
27 27 fp = open(os.path.join(PATH, experiment, 'experiment.conf'))
28 28 conf = json.loads(fp.read())
29 29
30 30 ipp_km = conf['usrp_tx']['ipp']
31 31 ipp = ipp_km * 2 /300000
32 32 samp_rate = conf['usrp_rx']['sample_rate']
33 33
34 34 #axis = ['0' if x=='elevation' else '1' for x in conf['pedestal']['speed']] # AZIMUTH 1 ELEVACION 0
35 35 axis = ['0' if x=='elevation' else '1' for x in conf['pedestal']['axis']] # AZIMUTH 1 ELEVACION 0
36 36 speed_axis = conf['pedestal']['speed']
37 37 steeps = conf['pedestal']['table']
38 38 time_offset = args.time_offset
39 39 parameters = args.parameters
40 40 #start_date = experiment.split('@')[1].split('T')[0].replace('-', '/')
41 41 start_date = '2022/04/22'
42 42 end_date = start_date
43 43 #start_time = experiment.split('@')[1].split('T')[1]
44 44 start_time = '00:00:01'
45 45 end_time = '23:59:59'
46 46 max_index = int(samp_rate*ipp*1e6 * args.range / 60) + int(samp_rate*ipp*1e6 * 1.2 / 60)
47 47 N = int(1/(speed_axis[0]*ipp)) # 1 GRADO DE RESOLUCION
48 48 path = os.path.join(PATH, experiment, 'rawdata')
49 49 path_ped = os.path.join(PATH, experiment, 'position')
50 50 path_plots = os.path.join(PATH, experiment, 'plots')
51 51 path_save = os.path.join(PATH, experiment, 'param')
52 52
53 53 dBmin = 35
54 54 dBmax = 60
55 55 Vmin = -20
56 56 Vmax = 20
57 57
58 58 from schainpy.controller import Project
59 59
60 60 project = Project()
61 61 project.setup(id='1', name='Sophy', description='sophy proc')
62 62
63 63 reader = project.addReadUnit(datatype='DigitalRFReader',
64 64 path=path,
65 65 startDate=start_date,
66 66 endDate=end_date,
67 67 startTime=start_time,
68 68 endTime=end_time,
69 69 delay=0,
70 70 online=0,
71 71 walk=1,
72 72 ippKm = ipp_km,
73 73 getByBlock = 1,
74 74 nProfileBlocks = N,
75 75 )
76 76
77 77 voltage = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
78 78 op = voltage.addOperation(name='setH0')
79 79 op.addParameter(name='h0', value='-1.2')
80 80
81 81 if args.range > 0:
82 82 op = voltage.addOperation(name='selectHeights')
83 83 op.addParameter(name='minIndex', value='0', format='int')
84 84 op.addParameter(name='maxIndex', value=max_index, format='int')
85 85
86 86 op = voltage.addOperation(name='PulsePair_vRF', optype='other')
87 87 op.addParameter(name='n', value=int(N), format='int')
88 88
89 89 proc = project.addProcUnit(datatype='ParametersProc', inputId=voltage.getId())
90 90
91 91 #-----------------------new--------- variables polarimetricas---------------
92 92 opObj10 = proc.addOperation(name="WeatherRadar")
93 93 opObj10.addParameter(name='variableList',value='Reflectividad,ReflectividadDiferencial,CoeficienteCorrelacion,FaseDiferencial,VelocidadRadial,AnchoEspectral')
94 94
95 95 #---------------------------------------------------------------------------
96 96 op = proc.addOperation(name='PedestalInformation')
97 97 op.addParameter(name='path', value=path_ped, format='str')
98 98 op.addParameter(name='interval', value='0.04', format='float')
99 99 op.addParameter(name='time_offset', value=time_offset)
100 100
101 101 for param in parameters:
102 102 op = proc.addOperation(name='Block360_vRF4')
103 103 op.addParameter(name='axis', value=','.join(axis))
104 104 op.addParameter(name='attr_data', value=PARAM[param]['name'])
105 105
106 106 if axis[0] == '1':
107 107 path_fig = '/PPI-{}km'.format(args.range)
108 108 op= proc.addOperation(name='Weather_vRF_Plot')
109 109 op.addParameter(name='save', value=path_plots+path_fig, format='str')
110 110 op.addParameter(name='save_period', value=-1)
111 111 op.addParameter(name='show', value=args.show)
112 112 op.addParameter(name='channels', value='1,')
113 113 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
114 114 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
115 115 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
116 116 op.addParameter(name='labels', value=[PARAM[param]['label']])
117 117 op.addParameter(name='save_code', value=param)
118 118 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
119 119 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
120 120 if axis[0] == '0':
121 121 path_fig = '/RHI{}km'.format(args.range)
122 122 op= proc.addOperation(name='WeatherRHI_vRF4_Plot')
123 123 op.addParameter(name='save', value=path_plots+path_fig, format='str')
124 124 op.addParameter(name='save_period', value=-1)
125 125 op.addParameter(name='show', value=args.show)
126 126 op.addParameter(name='channels', value='(1,)')
127 127 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
128 128 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
129 129 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
130 130 op.addParameter(name='labels', value=[PARAM[param]['label']])
131 131 op.addParameter(name='save_code', value=param)
132 132 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
133 133 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
134 134
135 135 if args.save:
136 136 opObj10 = proc.addOperation(name='HDFWriter')
137 137 opObj10.addParameter(name='path',value=path_save, format='str')
138 138 opObj10.addParameter(name='Reset',value=True)
139 139 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
140 140 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
141 141 opObj10.addParameter(name='dataList',value='dataPP_POWER,utctime')
142 142
143 143
144 144 project.start()
145 145
146 146 if __name__ == '__main__':
147 147
148 148 parser = argparse.ArgumentParser(description='Script to process SOPHy data.')
149 149 parser.add_argument('experiment',
150 150 help='Experiment name')
151 151 parser.add_argument('--parameters', nargs='*', default=['P'],
152 152 help='Variables to process: P, Z, V')
153 153 parser.add_argument('--time_offset', default=0,
154 154 help='Fix time offset')
155 155 parser.add_argument('--range', default=0, type=int,
156 156 help='Max range to plot')
157 157 parser.add_argument('--save', action='store_true',
158 158 help='Create output files')
159 159 parser.add_argument('--show', action='store_true',
160 160 help='Show matplotlib plot.')
161 161
162 162 args = parser.parse_args()
163 163 print (args)
164 164 main(args)
General Comments 0
You need to be logged in to leave comments. Login now