##// END OF EJS Templates
schain
RSS
Update ppi and rhi plot fix processing ppi+rhi
jespinoza -
r1562:dab94d864b0c
parent child
Show More
Side by Side Unified
Context file:
r1562:dab94d864b0c Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,723 +1,749
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,re
19 19
20 20 if 'BACKEND' in os.environ:
21 21 matplotlib.use(os.environ['BACKEND'])
22 22 elif 'linux' in sys.platform:
23 matplotlib.use("TkAgg")
23 matplotlib.use("Agg")
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 .plotting_codes import *
37 37
38 38 from schainpy.model.data.jrodata import PlotterData
39 39 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
40 40 from schainpy.utils import log
41 41
42 42 for name, cb_table in sophy_cb_tables:
43 43 ncmap = matplotlib.colors.ListedColormap(cb_table, name=name)
44 44 matplotlib.pyplot.register_cmap(cmap=ncmap)
45 45
46 46 EARTH_RADIUS = 6.3710e3
47 47
48 48 def ll2xy(lat1, lon1, lat2, lon2):
49 49
50 50 p = 0.017453292519943295
51 51 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
52 52 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
53 53 r = 12742 * numpy.arcsin(numpy.sqrt(a))
54 54 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
55 55 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
56 56 theta = -theta + numpy.pi/2
57 57 return r*numpy.cos(theta), r*numpy.sin(theta)
58 58
59 59
60 60 def km2deg(km):
61 61 '''
62 62 Convert distance in km to degrees
63 63 '''
64 64
65 65 return numpy.rad2deg(km/EARTH_RADIUS)
66 66
67 67
68 68 def figpause(interval):
69 69 backend = plt.rcParams['backend']
70 70 if backend in matplotlib.rcsetup.interactive_bk:
71 71 figManager = matplotlib._pylab_helpers.Gcf.get_active()
72 72 if figManager is not None:
73 73 canvas = figManager.canvas
74 74 if canvas.figure.stale:
75 75 canvas.draw()
76 76 try:
77 77 canvas.start_event_loop(interval)
78 78 except:
79 79 pass
80 80 return
81 81
82 82 def popup(message):
83 83 '''
84 84 '''
85 85
86 86 fig = plt.figure(figsize=(12, 8), facecolor='r')
87 87 text = '\n'.join([s.strip() for s in message.split(':')])
88 88 fig.text(0.01, 0.5, text, ha='left', va='center',
89 89 size='20', weight='heavy', color='w')
90 90 fig.show()
91 91 figpause(1000)
92 92
93 93
94 94 class Throttle(object):
95 95 '''
96 96 Decorator that prevents a function from being called more than once every
97 97 time period.
98 98 To create a function that cannot be called more than once a minute, but
99 99 will sleep until it can be called:
100 100 @Throttle(minutes=1)
101 101 def foo():
102 102 pass
103 103
104 104 for i in range(10):
105 105 foo()
106 106 print "This function has run %s times." % i
107 107 '''
108 108
109 109 def __init__(self, seconds=0, minutes=0, hours=0):
110 110 self.throttle_period = datetime.timedelta(
111 111 seconds=seconds, minutes=minutes, hours=hours
112 112 )
113 113
114 114 self.time_of_last_call = datetime.datetime.min
115 115
116 116 def __call__(self, fn):
117 117 @wraps(fn)
118 118 def wrapper(*args, **kwargs):
119 119 coerce = kwargs.pop('coerce', None)
120 120 if coerce:
121 121 self.time_of_last_call = datetime.datetime.now()
122 122 return fn(*args, **kwargs)
123 123 else:
124 124 now = datetime.datetime.now()
125 125 time_since_last_call = now - self.time_of_last_call
126 126 time_left = self.throttle_period - time_since_last_call
127 127
128 128 if time_left > datetime.timedelta(seconds=0):
129 129 return
130 130
131 131 self.time_of_last_call = datetime.datetime.now()
132 132 return fn(*args, **kwargs)
133 133
134 134 return wrapper
135 135
136 136 def apply_throttle(value):
137 137
138 138 @Throttle(seconds=value)
139 139 def fnThrottled(fn):
140 140 fn()
141 141
142 142 return fnThrottled
143 143
144 144
145 145 @MPDecorator
146 146 class Plot(Operation):
147 147 """Base class for Schain plotting operations
148 148
149 149 This class should never be use directtly you must subclass a new operation,
150 150 children classes must be defined as follow:
151 151
152 152 ExamplePlot(Plot):
153 153
154 154 CODE = 'code'
155 155 colormap = 'jet'
156 156 plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
157 157
158 158 def setup(self):
159 159 pass
160 160
161 161 def plot(self):
162 162 pass
163 163
164 164 """
165 165
166 166 CODE = 'Figure'
167 167 colormap = 'jet'
168 168 bgcolor = 'white'
169 169 buffering = True
170 170 __missing = 1E30
171 projection = None
171 172
172 173 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
173 174 'showprofile']
174 175
175 176 def __init__(self):
176 177
177 178 Operation.__init__(self)
178 179 self.isConfig = False
179 180 self.isPlotConfig = False
180 181 self.save_time = 0
181 182 self.sender_time = 0
182 183 self.data = None
183 184 self.firsttime = True
184 185 self.sender_queue = deque(maxlen=10)
185 186 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
186 187
187 188 def __fmtTime(self, x, pos):
188 189 '''
189 190 '''
190 191
191 192 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
192 193
193 194 def __setup(self, **kwargs):
194 195 '''
195 196 Initialize variables
196 197 '''
197 198
198 199 self.figures = []
199 200 self.axes = []
200 201 self.cb_axes = []
201 202 self.localtime = kwargs.pop('localtime', True)
202 203 self.show = kwargs.get('show', True)
203 204 self.save = kwargs.get('save', False)
204 205 self.save_period = kwargs.get('save_period', 0)
205 206 self.colormap = kwargs.get('colormap', self.colormap)
206 207 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
207 208 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
208 209 self.colormaps = kwargs.get('colormaps', None)
209 210 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
210 211 self.showprofile = kwargs.get('showprofile', False)
211 212 self.title = kwargs.get('wintitle', self.CODE.upper())
212 213 self.cb_label = kwargs.get('cb_label', None)
213 214 self.cb_labels = kwargs.get('cb_labels', None)
214 215 self.labels = kwargs.get('labels', None)
215 216 self.xaxis = kwargs.get('xaxis', 'frequency')
216 217 self.zmin = kwargs.get('zmin', None)
217 218 self.zmax = kwargs.get('zmax', None)
218 219 self.zlimits = kwargs.get('zlimits', None)
219 220 self.xmin = kwargs.get('xmin', None)
220 221 self.xmax = kwargs.get('xmax', None)
221 222 self.xrange = kwargs.get('xrange', 12)
222 223 self.xscale = kwargs.get('xscale', None)
223 224 self.ymin = kwargs.get('ymin', None)
224 225 self.ymax = kwargs.get('ymax', None)
225 226 self.yscale = kwargs.get('yscale', None)
226 227 self.xlabel = kwargs.get('xlabel', None)
227 228 self.attr_time = kwargs.get('attr_time', 'utctime')
228 229 self.attr_data = kwargs.get('attr_data', 'data_param')
229 230 self.decimation = kwargs.get('decimation', None)
230 231 self.oneFigure = kwargs.get('oneFigure', True)
231 232 self.width = kwargs.get('width', None)
232 233 self.height = kwargs.get('height', None)
233 234 self.colorbar = kwargs.get('colorbar', True)
234 235 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
235 236 self.channels = kwargs.get('channels', None)
236 237 self.titles = kwargs.get('titles', [])
237 238 self.polar = False
238 239 self.type = kwargs.get('type', 'iq')
239 240 self.grid = kwargs.get('grid', False)
240 241 self.pause = kwargs.get('pause', False)
241 242 self.save_code = kwargs.get('save_code', self.CODE)
242 243 self.throttle = kwargs.get('throttle', 0)
243 244 self.exp_code = kwargs.get('exp_code', None)
244 245 self.server = kwargs.get('server', False)
245 246 self.sender_period = kwargs.get('sender_period', 60)
246 247 self.tag = kwargs.get('tag', '')
247 248 self.height_index = kwargs.get('height_index', None)
248 249 self.__throttle_plot = apply_throttle(self.throttle)
249 250 code = self.attr_data if self.attr_data else self.CODE
250 251 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
251 252 self.ang_min = kwargs.get('ang_min', None)
252 253 self.ang_max = kwargs.get('ang_max', None)
253 254 self.mode = kwargs.get('mode', None)
254 255 self.mask = kwargs.get('mask', False)
255
256 self.shapes = kwargs.get('shapes', './')
256 257
257 258 if self.server:
258 259 if not self.server.startswith('tcp://'):
259 260 self.server = 'tcp://{}'.format(self.server)
260 261 log.success(
261 262 'Sending to server: {}'.format(self.server),
262 263 self.name
263 264 )
264 265
265 266 if isinstance(self.attr_data, str):
266 267 self.attr_data = [self.attr_data]
267 268
268 269 def __setup_plot(self):
269 270 '''
270 271 Common setup for all figures, here figures and axes are created
271 272 '''
272 273
273 274 self.setup()
274 275
275 276 self.time_label = 'LT' if self.localtime else 'UTC'
276 277
277 278 if self.width is None:
278 279 self.width = 8
279 280
280 self.figures = []
281 self.axes = []
281 self.figures = {'PPI':[], 'RHI':[]}
282 self.axes = {'PPI':[], 'RHI':[]}
282 283 self.cb_axes = []
283 284 self.pf_axes = []
284 285 self.cmaps = []
285 286
286 287 size = '15%' if self.ncols == 1 else '30%'
287 288 pad = '4%' if self.ncols == 1 else '8%'
288 289
289 290 if self.oneFigure:
290 291 if self.height is None:
291 292 self.height = 1.4 * self.nrows + 1
292 fig = plt.figure(figsize=(self.width, self.height),
293 fig_p = plt.figure(figsize=(self.width, self.height),
293 294 edgecolor='k',
294 295 facecolor='w')
295 self.figures.append(fig)
296 for n in range(self.nplots):
297 ax = fig.add_subplot(self.nrows, self.ncols,
298 n + 1, polar=self.polar)
299 ax.tick_params(labelsize=8)
300 ax.firsttime = True
301 ax.index = 0
302 ax.press = None
303 self.axes.append(ax)
296 fig_r = plt.figure(figsize=(self.width, self.height),
297 edgecolor='k',
298 facecolor='w')
299 self.figures['PPI'].append(fig_p)
300 self.figures['RHI'].append(fig_r)
301 for n in range(self.nplots):
302 ax_p = fig_p.add_subplot(self.nrows, self.ncols, n+1, polar=self.polar, projection=self.projection)
303 ax_r = fig_r.add_subplot(self.nrows, self.ncols, n+1, polar=self.polar)
304 ax_p.tick_params(labelsize=8)
305 ax_p.firsttime = True
306 ax_p.index = 0
307 ax_p.press = None
308 ax_r.tick_params(labelsize=8)
309 ax_r.firsttime = True
310 ax_r.index = 0
311 ax_r.press = None
312
313 self.axes['PPI'].append(ax_p)
314 self.axes['RHI'].append(ax_r)
315
304 316 if self.showprofile:
305 317 cax = self.__add_axes(ax, size=size, pad=pad)
306 318 cax.tick_params(labelsize=8)
307 319 self.pf_axes.append(cax)
308 320 else:
309 321 if self.height is None:
310 322 self.height = 3
311 323 for n in range(self.nplots):
312 324 fig = plt.figure(figsize=(self.width, self.height),
313 325 edgecolor='k',
314 326 facecolor='w')
315 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
316 ax.tick_params(labelsize=8)
317 ax.firsttime = True
318 ax.index = 0
319 ax.press = None
327 ax_p = fig.add_subplot(1, 1, 1, polar=self.polar, projection=self.projection)
328 ax_r = fig.add_subplot(1, 1, 1, polar=self.polar)
329 ax_p.tick_params(labelsize=8)
330 ax_p.firsttime = True
331 ax_p.index = 0
332 ax_p.press = None
333 ax_r.tick_params(labelsize=8)
334 ax_r.firsttime = True
335 ax_r.index = 0
336 ax_r.press = None
320 337 self.figures.append(fig)
321 self.axes.append(ax)
338 self.axes['PPI'].append(ax_p)
339 self.axes['RHI'].append(ax_r)
322 340 if self.showprofile:
323 341 cax = self.__add_axes(ax, size=size, pad=pad)
324 342 cax.tick_params(labelsize=8)
325 343 self.pf_axes.append(cax)
326 344
327 345 for n in range(self.nrows):
328 346 if self.colormaps is not None:
329 347 cmap = plt.get_cmap(self.colormaps[n])
330 348 else:
331 349 cmap = plt.get_cmap(self.colormap)
332 350 cmap.set_bad(self.bgcolor, 1.)
333 351 self.cmaps.append(cmap)
334 352
335 353 def __add_axes(self, ax, size='30%', pad='8%'):
336 354 '''
337 355 Add new axes to the given figure
338 356 '''
339 357 divider = make_axes_locatable(ax)
340 358 nax = divider.new_horizontal(size=size, pad=pad)
341 359 ax.figure.add_axes(nax)
342 360 return nax
343 361
344 362 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
345 363 '''
346 364 Create a masked array for missing data
347 365 '''
348 366 if x_buffer.shape[0] < 2:
349 367 return x_buffer, y_buffer, z_buffer
350 368
351 369 deltas = x_buffer[1:] - x_buffer[0:-1]
352 370 x_median = numpy.median(deltas)
353 371
354 372 index = numpy.where(deltas > 5 * x_median)
355 373
356 374 if len(index[0]) != 0:
357 375 z_buffer[::, index[0], ::] = self.__missing
358 376 z_buffer = numpy.ma.masked_inside(z_buffer,
359 377 0.99 * self.__missing,
360 378 1.01 * self.__missing)
361 379
362 380 return x_buffer, y_buffer, z_buffer
363 381
364 382 def decimate(self):
365 383
366 384 # dx = int(len(self.x)/self.__MAXNUMX) + 1
367 385 dy = int(len(self.y) / self.decimation) + 1
368 386
369 387 # x = self.x[::dx]
370 388 x = self.x
371 389 y = self.y[::dy]
372 390 z = self.z[::, ::, ::dy]
373 391
374 392 return x, y, z
375 393
376 394 def format(self):
377 395 '''
378 396 Set min and max values, labels, ticks and titles
379 397 '''
380 398
381 for n, ax in enumerate(self.axes):
399 for n, ax in enumerate(self.axes[self.mode]):
382 400 if ax.firsttime:
383 401 if self.xaxis != 'time':
384 402 xmin = self.xmin
385 403 xmax = self.xmax
386 404 else:
387 405 xmin = self.tmin
388 406 xmax = self.tmin + self.xrange*60*60
389 407 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
390 408 ax.xaxis.set_major_locator(LinearLocator(9))
391 409 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
392 410 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
411
393 412 ax.set_facecolor(self.bgcolor)
413
394 414 if self.xscale:
395 415 ax.xaxis.set_major_formatter(FuncFormatter(
396 416 lambda x, pos: '{0:g}'.format(x*self.xscale)))
397 417 if self.yscale:
398 418 ax.yaxis.set_major_formatter(FuncFormatter(
399 419 lambda x, pos: '{0:g}'.format(x*self.yscale)))
400 420 if self.xlabel is not None:
401 421 ax.set_xlabel(self.xlabel)
402 422 if self.ylabel is not None:
403 423 ax.set_ylabel(self.ylabel)
404 424 if self.showprofile:
405 425 self.pf_axes[n].set_ylim(ymin, ymax)
406 426 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
407 427 self.pf_axes[n].set_xlabel('dB')
408 428 self.pf_axes[n].grid(b=True, axis='x')
409 429 [tick.set_visible(False)
410 for tick in self.pf_axes[n].get_yticklabels()]
430 for tick in self.pf_axes[n].get_yticklabels()]
411 431 if self.colorbar:
412 432 ax.cbar = plt.colorbar(
413 433 ax.plt, ax=ax, fraction=0.05, pad=0.06, aspect=10)
414 434 ax.cbar.ax.tick_params(labelsize=8)
415 435 ax.cbar.ax.press = None
416 436 if self.cb_label:
417 437 ax.cbar.set_label(self.cb_label, size=8)
418 438 elif self.cb_labels:
419 439 ax.cbar.set_label(self.cb_labels[n], size=8)
420 440 else:
421 441 ax.cbar = None
422 ax.set_xlim(xmin, xmax)
423 ax.set_ylim(ymin, ymax)
442 if self.mode == 'RHI':
443 ax.set_xlim(xmin, xmax)
444 ax.set_ylim(ymin, ymax)
424 445 ax.firsttime = False
425 446 if self.grid:
426 447 ax.grid(True)
427 448 if not self.polar:
428 449 ax.set_title('{} {} {}'.format(
429 450 self.titles[n],
430 451 self.getDateTime(self.data.max_time).strftime(
431 452 '%Y-%m-%d %H:%M:%S'),
432 453 self.time_label),
433 454 size=8)
434 455 else:
435 456 #ax.set_title('{}'.format(self.titles[n]), size=8)
436 457 ax.set_title('{} {} {}'.format(
437 458 self.titles[n],
438 459 self.getDateTime(self.data.max_time).strftime(
439 460 '%Y-%m-%d %H:%M:%S'),
440 461 self.time_label),
441 462 size=8)
442 463 ax.set_ylim(0, self.ymax)
443 464 if self.mode == 'PPI':
444 465 ax.set_yticks(ax.get_yticks(), labels=ax.get_yticks(), color='white')
445 466 ax.yaxis.labelpad = 28
446 467 elif self.mode == 'RHI':
447 468 ax.xaxis.labelpad = 16
448 469
449 470 if self.firsttime:
450 for n, fig in enumerate(self.figures):
471 for fig in self.figures['PPI'] + self.figures['RHI']:
451 472 fig.subplots_adjust(**self.plots_adjust)
452 473 self.firsttime = False
453 474
454 475 def clear_figures(self):
455 476 '''
456 477 Reset axes for redraw plots
457 478 '''
458 479
459 for ax in self.axes+self.pf_axes+self.cb_axes:
480 axes = self.pf_axes + self.cb_axes + self.axes[self.mode]
481
482 for ax in axes:
460 483 ax.clear()
461 484 ax.firsttime = True
462 485 if hasattr(ax, 'cbar') and ax.cbar:
463 486 ax.cbar.remove()
464 487
465 488 def __plot(self):
466 489 '''
467 490 Main function to plot, format and save figures
468 491 '''
469 492
470 493 self.plot()
471 494 self.format()
472
473 for n, fig in enumerate(self.figures):
495 figures = self.figures[self.mode]
496 for n, fig in enumerate(figures):
474 497 if self.nrows == 0 or self.nplots == 0:
475 498 log.warning('No data', self.name)
476 499 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
477 500 fig.canvas.manager.set_window_title(self.CODE)
478 501 continue
479 502
480 503 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
481 504 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
482 505 fig.canvas.draw()
483 506 if self.show:
484 507 fig.show()
485 508 figpause(0.01)
486 509
487 510 if self.save:
488 511 self.save_figure(n)
489 512
490 513 if self.server:
491 514 if self.mode and self.mode == 'RHI':
492 515 return
493 516 self.send_to_server()
494 517
495 518 def __update(self, dataOut, timestamp):
496 519 '''
497 520 '''
498 521
499 522 metadata = {
500 523 'yrange': dataOut.heightList,
501 524 'interval': dataOut.timeInterval,
502 525 'channels': dataOut.channelList
503 526 }
504 527
505 528 data, meta = self.update(dataOut)
506 529 metadata.update(meta)
507 530 self.data.update(data, timestamp, metadata)
508 531
509 532 def save_figure(self, n):
510 533 '''
511 534 '''
512 535 if self.mode is not None:
513 536 ang = 'AZ' if self.mode == 'RHI' else 'EL'
514 label = '_{}_{}_{}'.format(self.mode, ang, self.mode_value)
537 folder = '_{}_{}_{}'.format(self.mode, ang, self.mode_value)
538 label = '{}{}_{}'.format(ang[0], self.mode_value, self.save_code)
515 539 else:
540 folder = ''
516 541 label = ''
517 542
518 543 if self.oneFigure:
519 544 if (self.data.max_time - self.save_time) <= self.save_period:
520 545 return
521 546
522 547 self.save_time = self.data.max_time
523 548
524 fig = self.figures[n]
549 fig = self.figures[self.mode][n]
525 550
526 551 if self.throttle == 0:
527 552 if self.oneFigure:
528 553 figname = os.path.join(
529 554 self.save,
530 self.save_code + label,
531 '{}_{}.png'.format(
532 self.save_code + label,
555 self.save_code + folder,
556 '{}_{}_{}.png'.format(
557 'SOPHY',
533 558 self.getDateTime(self.data.max_time).strftime(
534 559 '%Y%m%d_%H%M%S'
535 560 ),
561 label
536 562 )
537 563 )
538 564 else:
539 565 figname = os.path.join(
540 566 self.save,
541 567 self.save_code,
542 568 '{}_ch{}_{}.png'.format(
543 569 self.save_code, n,
544 570 self.getDateTime(self.data.max_time).strftime(
545 571 '%Y%m%d_%H%M%S'
546 572 ),
547 573 )
548 574 )
549 575 log.log('Saving figure: {}'.format(figname), self.name)
550 576 if not os.path.isdir(os.path.dirname(figname)):
551 577 os.makedirs(os.path.dirname(figname))
552 578 fig.savefig(figname)
553 579
554 580 figname = os.path.join(
555 581 self.save,
556 582 '{}_{}.png'.format(
557 583 self.save_code,
558 584 self.getDateTime(self.data.min_time).strftime(
559 585 '%Y%m%d'
560 586 ),
561 587 )
562 588 )
563 589
564 590 log.log('Saving figure: {}'.format(figname), self.name)
565 591 if not os.path.isdir(os.path.dirname(figname)):
566 592 os.makedirs(os.path.dirname(figname))
567 593 fig.savefig(figname)
568 594
569 595 def send_to_server(self):
570 596 '''
571 597 '''
572 598
573 599 if self.exp_code == None:
574 600 log.warning('Missing `exp_code` skipping sending to server...')
575 601
576 602 last_time = self.data.max_time
577 603 interval = last_time - self.sender_time
578 604 if interval < self.sender_period:
579 605 return
580 606
581 607 self.sender_time = last_time
582 608
583 609 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
584 610 for attr in attrs:
585 611 value = getattr(self, attr)
586 612 if value:
587 613 if isinstance(value, (numpy.float32, numpy.float64)):
588 614 value = round(float(value), 2)
589 615 self.data.meta[attr] = value
590 616 if self.colormap == 'jet' or self.colormap == 'sophy_w':
591 617 self.data.meta['colormap'] = 'Jet'
592 618 elif 'sophy_v' in self.colormap:
593 619 self.data.meta['colormap'] = 'RdBu'
594 620 else:
595 621 self.data.meta['colormap'] = 'Viridis'
596 622 self.data.meta['interval'] = int(interval)
597 623
598 624 self.sender_queue.append(last_time)
599 625
600 626 while True:
601 627 try:
602 628 tm = self.sender_queue.popleft()
603 629 except IndexError:
604 630 break
605 631 msg = self.data.jsonify(tm, self.save_code, self.plot_type, key='var')
606 632 self.socket.send_string(msg)
607 633 socks = dict(self.poll.poll(2000))
608 634 if socks.get(self.socket) == zmq.POLLIN:
609 635 reply = self.socket.recv_string()
610 636 if reply == 'ok':
611 637 log.log("Response from server ok", self.name)
612 638 time.sleep(0.1)
613 639 continue
614 640 else:
615 641 log.warning(
616 642 "Malformed reply from server: {}".format(reply), self.name)
617 643 else:
618 644 log.warning(
619 645 "No response from server, retrying...", self.name)
620 646 self.sender_queue.appendleft(tm)
621 647 self.socket.setsockopt(zmq.LINGER, 0)
622 648 self.socket.close()
623 649 self.poll.unregister(self.socket)
624 650 self.socket = self.context.socket(zmq.REQ)
625 651 self.socket.connect(self.server)
626 652 self.poll.register(self.socket, zmq.POLLIN)
627 653 break
628 654
629 655 def setup(self):
630 656 '''
631 657 This method should be implemented in the child class, the following
632 658 attributes should be set:
633 659
634 660 self.nrows: number of rows
635 661 self.ncols: number of cols
636 662 self.nplots: number of plots (channels or pairs)
637 663 self.ylabel: label for Y axes
638 664 self.titles: list of axes title
639 665
640 666 '''
641 667 raise NotImplementedError
642 668
643 669 def plot(self):
644 670 '''
645 671 Must be defined in the child class, the actual plotting method
646 672 '''
647 673 raise NotImplementedError
648 674
649 675 def update(self, dataOut):
650 676 '''
651 677 Must be defined in the child class, update self.data with new data
652 678 '''
653 679
654 680 data = {
655 681 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
656 682 }
657 683 meta = {}
658 684
659 685 return data, meta
660 686
661 687 def run(self, dataOut, **kwargs):
662 688 '''
663 689 Main plotting routine
664 690 '''
665 691
666 692 if self.isConfig is False:
667 693 self.__setup(**kwargs)
668 694
669 695 if self.localtime:
670 696 self.getDateTime = datetime.datetime.fromtimestamp
671 697 else:
672 698 self.getDateTime = datetime.datetime.utcfromtimestamp
673 699
674 700 self.data.setup()
675 701 self.isConfig = True
676 702 if self.server:
677 703 self.context = zmq.Context()
678 704 self.socket = self.context.socket(zmq.REQ)
679 705 self.socket.connect(self.server)
680 706 self.poll = zmq.Poller()
681 707 self.poll.register(self.socket, zmq.POLLIN)
682 708
683 709 tm = getattr(dataOut, self.attr_time)
684 710
685 711 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
686 712 self.save_time = tm
687 713 self.__plot()
688 714 self.tmin += self.xrange*60*60
689 715 self.data.setup()
690 716 self.clear_figures()
691 717
692 718 self.__update(dataOut, tm)
693 719
694 720 if self.isPlotConfig is False:
695 721 self.__setup_plot()
696 722 self.isPlotConfig = True
697 723 if self.xaxis == 'time':
698 724 dt = self.getDateTime(tm)
699 725 if self.xmin is None:
700 726 self.tmin = tm
701 727 self.xmin = dt.hour
702 728 minutes = (self.xmin-int(self.xmin)) * 60
703 729 seconds = (minutes - int(minutes)) * 60
704 730 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
705 731 datetime.datetime(1970, 1, 1)).total_seconds()
706 732 if self.localtime:
707 733 self.tmin += time.timezone
708 734
709 735 if self.xmin is not None and self.xmax is not None:
710 736 self.xrange = self.xmax - self.xmin
711 737
712 738 if self.throttle == 0:
713 739 self.__plot()
714 740 else:
715 741 self.__throttle_plot(self.__plot)#, coerce=coerce)
716 742
717 743 def close(self):
718 744
719 745 if self.data and not self.data.flagNoData:
720 746 self.save_time = 0
721 747 self.__plot()
722 748 if self.data and not self.data.flagNoData and self.pause:
723 749 figpause(10)
Show file before | Show file after | Hide comments
@@ -1,697 +1,716
1 1 import os
2 2 import datetime
3 3 import warnings
4 4 import numpy
5 5 from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
6 from matplotlib.patches import Circle
7 import cartopy.crs as ccrs
8 from cartopy.feature import ShapelyFeature
9 import cartopy.io.shapereader as shpreader
6 10
7 11 from schainpy.model.graphics.jroplot_base import Plot, plt
8 12 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
9 13 from schainpy.utils import log
10 14
11 15
12 16 EARTH_RADIUS = 6.3710e3
13 17
14 18
15 19 def antenna_to_cartesian(ranges, azimuths, elevations):
16 20 """
17 21 Return Cartesian coordinates from antenna coordinates.
18 22
19 23 Parameters
20 24 ----------
21 25 ranges : array
22 26 Distances to the center of the radar gates (bins) in kilometers.
23 27 azimuths : array
24 28 Azimuth angle of the radar in degrees.
25 29 elevations : array
26 30 Elevation angle of the radar in degrees.
27 31
28 32 Returns
29 33 -------
30 34 x, y, z : array
31 35 Cartesian coordinates in meters from the radar.
32 36
33 37 Notes
34 38 -----
35 39 The calculation for Cartesian coordinate is adapted from equations
36 40 2.28(b) and 2.28(c) of Doviak and Zrnic [1]_ assuming a
37 41 standard atmosphere (4/3 Earth's radius model).
38 42
39 43 .. math::
40 44
41 45 z = \\sqrt{r^2+R^2+2*r*R*sin(\\theta_e)} - R
42 46
43 47 s = R * arcsin(\\frac{r*cos(\\theta_e)}{R+z})
44 48
45 49 x = s * sin(\\theta_a)
46 50
47 51 y = s * cos(\\theta_a)
48 52
49 53 Where r is the distance from the radar to the center of the gate,
50 54 :math:`\\theta_a` is the azimuth angle, :math:`\\theta_e` is the
51 55 elevation angle, s is the arc length, and R is the effective radius
52 56 of the earth, taken to be 4/3 the mean radius of earth (6371 km).
53 57
54 58 References
55 59 ----------
56 60 .. [1] Doviak and Zrnic, Doppler Radar and Weather Observations, Second
57 61 Edition, 1993, p. 21.
58 62
59 63 """
60 64 theta_e = numpy.deg2rad(elevations) # elevation angle in radians.
61 65 theta_a = numpy.deg2rad(azimuths) # azimuth angle in radians.
62 66 R = 6371.0 * 1000.0 * 4.0 / 3.0 # effective radius of earth in meters.
63 67 r = ranges * 1000.0 # distances to gates in meters.
64 68
65 69 z = (r ** 2 + R ** 2 + 2.0 * r * R * numpy.sin(theta_e)) ** 0.5 - R
66 70 s = R * numpy.arcsin(r * numpy.cos(theta_e) / (R + z)) # arc length in m.
67 71 x = s * numpy.sin(theta_a)
68 72 y = s * numpy.cos(theta_a)
69 73 return x, y, z
70 74
71 75 def cartesian_to_geographic_aeqd(x, y, lon_0, lat_0, R=EARTH_RADIUS):
72 76 """
73 77 Azimuthal equidistant Cartesian to geographic coordinate transform.
74 78
75 79 Transform a set of Cartesian/Cartographic coordinates (x, y) to
76 80 geographic coordinate system (lat, lon) using a azimuthal equidistant
77 81 map projection [1]_.
78 82
79 83 .. math::
80 84
81 85 lat = \\arcsin(\\cos(c) * \\sin(lat_0) +
82 86 (y * \\sin(c) * \\cos(lat_0) / \\rho))
83 87
84 88 lon = lon_0 + \\arctan2(
85 89 x * \\sin(c),
86 90 \\rho * \\cos(lat_0) * \\cos(c) - y * \\sin(lat_0) * \\sin(c))
87 91
88 92 \\rho = \\sqrt(x^2 + y^2)
89 93
90 94 c = \\rho / R
91 95
92 96 Where x, y are the Cartesian position from the center of projection;
93 97 lat, lon the corresponding latitude and longitude; lat_0, lon_0 are the
94 98 latitude and longitude of the center of the projection; R is the radius of
95 99 the earth (defaults to ~6371 km). lon is adjusted to be between -180 and
96 100 180.
97 101
98 102 Parameters
99 103 ----------
100 104 x, y : array-like
101 105 Cartesian coordinates in the same units as R, typically meters.
102 106 lon_0, lat_0 : float
103 107 Longitude and latitude, in degrees, of the center of the projection.
104 108 R : float, optional
105 109 Earth radius in the same units as x and y. The default value is in
106 110 units of meters.
107 111
108 112 Returns
109 113 -------
110 114 lon, lat : array
111 115 Longitude and latitude of Cartesian coordinates in degrees.
112 116
113 117 References
114 118 ----------
115 119 .. [1] Snyder, J. P. Map Projections--A Working Manual. U. S. Geological
116 120 Survey Professional Paper 1395, 1987, pp. 191-202.
117 121
118 122 """
119 123 x = numpy.atleast_1d(numpy.asarray(x))
120 124 y = numpy.atleast_1d(numpy.asarray(y))
121 125
122 126 lat_0_rad = numpy.deg2rad(lat_0)
123 127 lon_0_rad = numpy.deg2rad(lon_0)
124 128
125 129 rho = numpy.sqrt(x*x + y*y)
126 130 c = rho / R
127 131
128 132 with warnings.catch_warnings():
129 133 # division by zero may occur here but is properly addressed below so
130 134 # the warnings can be ignored
131 135 warnings.simplefilter("ignore", RuntimeWarning)
132 136 lat_rad = numpy.arcsin(numpy.cos(c) * numpy.sin(lat_0_rad) +
133 137 y * numpy.sin(c) * numpy.cos(lat_0_rad) / rho)
134 138 lat_deg = numpy.rad2deg(lat_rad)
135 139 # fix cases where the distance from the center of the projection is zero
136 140 lat_deg[rho == 0] = lat_0
137 141
138 142 x1 = x * numpy.sin(c)
139 143 x2 = rho*numpy.cos(lat_0_rad)*numpy.cos(c) - y*numpy.sin(lat_0_rad)*numpy.sin(c)
140 144 lon_rad = lon_0_rad + numpy.arctan2(x1, x2)
141 145 lon_deg = numpy.rad2deg(lon_rad)
142 146 # Longitudes should be from -180 to 180 degrees
143 147 lon_deg[lon_deg > 180] -= 360.
144 148 lon_deg[lon_deg < -180] += 360.
145 149
146 150 return lon_deg, lat_deg
147 151
148 152 def antenna_to_geographic(ranges, azimuths, elevations, site):
149 153
150 154 x, y, z = antenna_to_cartesian(numpy.array(ranges), numpy.array(azimuths), numpy.array(elevations))
151 155 lon, lat = cartesian_to_geographic_aeqd(x, y, site[0], site[1], R=6370997.)
152 156
153 157 return lon, lat
154 158
155 159 def ll2xy(lat1, lon1, lat2, lon2):
156 160
157 161 p = 0.017453292519943295
158 162 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
159 163 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
160 164 r = 12742 * numpy.arcsin(numpy.sqrt(a))
161 165 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
162 166 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
163 167 theta = -theta + numpy.pi/2
164 168 return r*numpy.cos(theta), r*numpy.sin(theta)
165 169
166 170
167 171 def km2deg(km):
168 172 '''
169 173 Convert distance in km to degrees
170 174 '''
171 175
172 176 return numpy.rad2deg(km/EARTH_RADIUS)
173 177
174 178
175 179
176 180 class SpectralMomentsPlot(SpectraPlot):
177 181 '''
178 182 Plot for Spectral Moments
179 183 '''
180 184 CODE = 'spc_moments'
181 185 # colormap = 'jet'
182 186 # plot_type = 'pcolor'
183 187
184 188 class DobleGaussianPlot(SpectraPlot):
185 189 '''
186 190 Plot for Double Gaussian Plot
187 191 '''
188 192 CODE = 'gaussian_fit'
189 193 # colormap = 'jet'
190 194 # plot_type = 'pcolor'
191 195
192 196 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
193 197 '''
194 198 Plot SpectraCut with Double Gaussian Fit
195 199 '''
196 200 CODE = 'cut_gaussian_fit'
197 201
198 202 class SnrPlot(RTIPlot):
199 203 '''
200 204 Plot for SNR Data
201 205 '''
202 206
203 207 CODE = 'snr'
204 208 colormap = 'jet'
205 209
206 210 def update(self, dataOut):
207 211
208 212 data = {
209 213 'snr': 10*numpy.log10(dataOut.data_snr)
210 214 }
211 215
212 216 return data, {}
213 217
214 218 class DopplerPlot(RTIPlot):
215 219 '''
216 220 Plot for DOPPLER Data (1st moment)
217 221 '''
218 222
219 223 CODE = 'dop'
220 224 colormap = 'jet'
221 225
222 226 def update(self, dataOut):
223 227
224 228 data = {
225 229 'dop': 10*numpy.log10(dataOut.data_dop)
226 230 }
227 231
228 232 return data, {}
229 233
230 234 class PowerPlot(RTIPlot):
231 235 '''
232 236 Plot for Power Data (0 moment)
233 237 '''
234 238
235 239 CODE = 'pow'
236 240 colormap = 'jet'
237 241
238 242 def update(self, dataOut):
239 243 data = {
240 244 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
241 245 }
242 246 return data, {}
243 247
244 248 class SpectralWidthPlot(RTIPlot):
245 249 '''
246 250 Plot for Spectral Width Data (2nd moment)
247 251 '''
248 252
249 253 CODE = 'width'
250 254 colormap = 'jet'
251 255
252 256 def update(self, dataOut):
253 257
254 258 data = {
255 259 'width': dataOut.data_width
256 260 }
257 261
258 262 return data, {}
259 263
260 264 class SkyMapPlot(Plot):
261 265 '''
262 266 Plot for meteors detection data
263 267 '''
264 268
265 269 CODE = 'param'
266 270
267 271 def setup(self):
268 272
269 273 self.ncols = 1
270 274 self.nrows = 1
271 275 self.width = 7.2
272 276 self.height = 7.2
273 277 self.nplots = 1
274 278 self.xlabel = 'Zonal Zenith Angle (deg)'
275 279 self.ylabel = 'Meridional Zenith Angle (deg)'
276 280 self.polar = True
277 281 self.ymin = -180
278 282 self.ymax = 180
279 283 self.colorbar = False
280 284
281 285 def plot(self):
282 286
283 287 arrayParameters = numpy.concatenate(self.data['param'])
284 288 error = arrayParameters[:, -1]
285 289 indValid = numpy.where(error == 0)[0]
286 290 finalMeteor = arrayParameters[indValid, :]
287 291 finalAzimuth = finalMeteor[:, 3]
288 292 finalZenith = finalMeteor[:, 4]
289 293
290 294 x = finalAzimuth * numpy.pi / 180
291 295 y = finalZenith
292 296
293 297 ax = self.axes[0]
294 298
295 299 if ax.firsttime:
296 300 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
297 301 else:
298 302 ax.plot.set_data(x, y)
299 303
300 304 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
301 305 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
302 306 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
303 307 dt2,
304 308 len(x))
305 309 self.titles[0] = title
306 310
307 311
308 312 class GenericRTIPlot(Plot):
309 313 '''
310 314 Plot for data_xxxx object
311 315 '''
312 316
313 317 CODE = 'param'
314 318 colormap = 'viridis'
315 319 plot_type = 'pcolorbuffer'
316 320
317 321 def setup(self):
318 322 self.xaxis = 'time'
319 323 self.ncols = 1
320 324 self.nrows = self.data.shape('param')[0]
321 325 self.nplots = self.nrows
322 326 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
323 327
324 328 if not self.xlabel:
325 329 self.xlabel = 'Time'
326 330
327 331 self.ylabel = 'Range [km]'
328 332 if not self.titles:
329 333 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
330 334
331 335 def update(self, dataOut):
332 336
333 337 data = {
334 338 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
335 339 }
336 340
337 341 meta = {}
338 342
339 343 return data, meta
340 344
341 345 def plot(self):
342 346 # self.data.normalize_heights()
343 347 self.x = self.data.times
344 348 self.y = self.data.yrange
345 349 self.z = self.data['param']
346 350 self.z = 10*numpy.log10(self.z)
347 351 self.z = numpy.ma.masked_invalid(self.z)
348 352
349 353 if self.decimation is None:
350 354 x, y, z = self.fill_gaps(self.x, self.y, self.z)
351 355 else:
352 356 x, y, z = self.fill_gaps(*self.decimate())
353 357
354 358 for n, ax in enumerate(self.axes):
355 359
356 360 self.zmax = self.zmax if self.zmax is not None else numpy.max(
357 361 self.z[n])
358 362 self.zmin = self.zmin if self.zmin is not None else numpy.min(
359 363 self.z[n])
360 364
361 365 if ax.firsttime:
362 366 if self.zlimits is not None:
363 367 self.zmin, self.zmax = self.zlimits[n]
364 368
365 369 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
366 370 vmin=self.zmin,
367 371 vmax=self.zmax,
368 372 cmap=self.cmaps[n]
369 373 )
370 374 else:
371 375 if self.zlimits is not None:
372 376 self.zmin, self.zmax = self.zlimits[n]
373 377 ax.collections.remove(ax.collections[0])
374 378 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
375 379 vmin=self.zmin,
376 380 vmax=self.zmax,
377 381 cmap=self.cmaps[n]
378 382 )
379 383
380 384
381 385 class PolarMapPlot(Plot):
382 386 '''
383 387 Plot for weather radar
384 388 '''
385 389
386 390 CODE = 'param'
387 391 colormap = 'seismic'
388 392
389 393 def setup(self):
390 394 self.ncols = 1
391 395 self.nrows = 1
392 396 self.width = 9
393 397 self.height = 8
394 398 self.mode = self.data.meta['mode']
395 399 if self.channels is not None:
396 400 self.nplots = len(self.channels)
397 401 self.nrows = len(self.channels)
398 402 else:
399 403 self.nplots = self.data.shape(self.CODE)[0]
400 404 self.nrows = self.nplots
401 405 self.channels = list(range(self.nplots))
402 406 if self.mode == 'E':
403 407 self.xlabel = 'Longitude'
404 408 self.ylabel = 'Latitude'
405 409 else:
406 410 self.xlabel = 'Range (km)'
407 411 self.ylabel = 'Height (km)'
408 412 self.bgcolor = 'white'
409 413 self.cb_labels = self.data.meta['units']
410 414 self.lat = self.data.meta['latitude']
411 415 self.lon = self.data.meta['longitude']
412 416 self.xmin, self.xmax = float(
413 417 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
414 418 self.ymin, self.ymax = float(
415 419 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
416 420 # self.polar = True
417 421
418 422 def plot(self):
419 423
420 424 for n, ax in enumerate(self.axes):
421 425 data = self.data['param'][self.channels[n]]
422 426
423 427 zeniths = numpy.linspace(
424 428 0, self.data.meta['max_range'], data.shape[1])
425 429 if self.mode == 'E':
426 430 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
427 431 r, theta = numpy.meshgrid(zeniths, azimuths)
428 432 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
429 433 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
430 434 x = km2deg(x) + self.lon
431 435 y = km2deg(y) + self.lat
432 436 else:
433 437 azimuths = numpy.radians(self.data.yrange)
434 438 r, theta = numpy.meshgrid(zeniths, azimuths)
435 439 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
436 440 self.y = zeniths
437 441
438 442 if ax.firsttime:
439 443 if self.zlimits is not None:
440 444 self.zmin, self.zmax = self.zlimits[n]
441 445 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
442 446 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
443 447 vmin=self.zmin,
444 448 vmax=self.zmax,
445 449 cmap=self.cmaps[n])
446 450 else:
447 451 if self.zlimits is not None:
448 452 self.zmin, self.zmax = self.zlimits[n]
449 453 ax.collections.remove(ax.collections[0])
450 454 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
451 455 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
452 456 vmin=self.zmin,
453 457 vmax=self.zmax,
454 458 cmap=self.cmaps[n])
455 459
456 460 if self.mode == 'A':
457 461 continue
458 462
459 463 # plot district names
460 464 f = open('/data/workspace/schain_scripts/distrito.csv')
461 465 for line in f:
462 466 label, lon, lat = [s.strip() for s in line.split(',') if s]
463 467 lat = float(lat)
464 468 lon = float(lon)
465 469 # ax.plot(lon, lat, '.b', ms=2)
466 470 ax.text(lon, lat, label.decode('utf8'), ha='center',
467 471 va='bottom', size='8', color='black')
468 472
469 473 # plot limites
470 474 limites = []
471 475 tmp = []
472 476 for line in open('/data/workspace/schain_scripts/lima.csv'):
473 477 if '#' in line:
474 478 if tmp:
475 479 limites.append(tmp)
476 480 tmp = []
477 481 continue
478 482 values = line.strip().split(',')
479 483 tmp.append((float(values[0]), float(values[1])))
480 484 for points in limites:
481 485 ax.add_patch(
482 486 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
483 487
484 488 # plot Cuencas
485 489 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
486 490 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
487 491 values = [line.strip().split(',') for line in f]
488 492 points = [(float(s[0]), float(s[1])) for s in values]
489 493 ax.add_patch(Polygon(points, ec='b', fc='none'))
490 494
491 495 # plot grid
492 496 for r in (15, 30, 45, 60):
493 497 ax.add_artist(plt.Circle((self.lon, self.lat),
494 498 km2deg(r), color='0.6', fill=False, lw=0.2))
495 499 ax.text(
496 500 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
497 501 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
498 502 '{}km'.format(r),
499 503 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
500 504
501 505 if self.mode == 'E':
502 506 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
503 507 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
504 508 else:
505 509 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
506 510 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
507 511
508 512 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
509 513 self.titles = ['{} {}'.format(
510 514 self.data.parameters[x], title) for x in self.channels]
511 515
512 516 class WeatherParamsPlot(Plot):
513 517 #CODE = 'RHI'
514 518 #plot_name = 'RHI'
515 519 plot_type = 'scattermap'
516 520 buffering = False
521 projection = ccrs.PlateCarree()
517 522
518 523 def setup(self):
519 524
520 525 self.ncols = 1
521 526 self.nrows = 1
522 527 self.nplots= 1
523 528 self.ylabel= 'Range [km]'
524 529 self.xlabel= 'Range [km]'
525 self.polar = True
526 self.grid = True
530
527 531 if self.channels is not None:
528 532 self.nplots = len(self.channels)
529 533 self.ncols = len(self.channels)
530 534 else:
531 535 self.nplots = self.data.shape(self.CODE)[0]
532 536 self.ncols = self.nplots
533 537 self.channels = list(range(self.nplots))
534 538
535 539 self.colorbar=True
536 540 if len(self.channels)>1:
537 541 self.width = 12
538 542 else:
539 543 self.width =8
540 self.height =8
544 self.height =7
541 545 self.ini =0
542 546 self.len_azi =0
543 547 self.buffer_ini = None
544 548 self.buffer_ele = None
545 549 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
546 550 self.flag =0
547 551 self.indicador= 0
548 552 self.last_data_ele = None
549 553 self.val_mean = None
550 554
551 555 def update(self, dataOut):
552 556
553 557 vars = {
554 558 'S' : 0,
555 559 'V' : 1,
556 560 'W' : 2,
557 561 'SNR' : 3,
558 562 'Z' : 4,
559 563 'D' : 5,
560 564 'P' : 6,
561 565 'R' : 7,
562 566 }
563 567
564 568 data = {}
565 569 meta = {}
566 570
567 571 if hasattr(dataOut, 'nFFTPoints'):
568 572 factor = dataOut.normFactor
569 573 else:
570 574 factor = 1
571 575
572 576 if hasattr(dataOut, 'dparam'):
573 577 tmp = getattr(dataOut, 'data_param')
574 578 else:
575 579
576 580 if 'S' in self.attr_data[0]:
577 581 tmp = 10*numpy.log10(10.0*getattr(dataOut, 'data_param')[:,0,:]/(factor))
578 582 else:
579 583 tmp = getattr(dataOut, 'data_param')[:,vars[self.attr_data[0]],:]
580 584
581 585 if self.mask:
582 586 mask = dataOut.data_param[:,3,:] < self.mask
583 587 tmp = numpy.ma.masked_array(tmp, mask=mask)
584 588
585 589 r = dataOut.heightList
586 590 delta_height = r[1]-r[0]
587 591 valid = numpy.where(r>=0)[0]
588 592 data['r'] = numpy.arange(len(valid))*delta_height
589 593
590 594 data['data'] = [0, 0]
591 595
592 #try:
593 data['data'][0] = tmp[0][:,valid]
594 data['data'][1] = tmp[1][:,valid]
595 #except:
596 # data['data'] = tmp[0][:,valid]
596 try:
597 data['data'][0] = tmp[0][:,valid]
598 data['data'][1] = tmp[1][:,valid]
599 except:
600 data['data'][0] = tmp[0][:,valid]
601 data['data'][1] = tmp[0][:,valid]
597 602
598 603 if dataOut.mode_op == 'PPI':
599 604 self.CODE = 'PPI'
600 605 self.title = self.CODE
601 606 elif dataOut.mode_op == 'RHI':
602 607 self.CODE = 'RHI'
603 608 self.title = self.CODE
604 609
605 610 data['azi'] = dataOut.data_azi
606 611 data['ele'] = dataOut.data_ele
607 612 data['mode_op'] = dataOut.mode_op
608 613 self.mode = dataOut.mode_op
609 var = data['data'][0].flatten()
610 r = numpy.tile(data['r'], data['data'][0].shape[0])
611 az = numpy.repeat(data['azi'], data['data'][0].shape[1])
612 el = numpy.repeat(data['ele'], data['data'][0].shape[1])
613
614 # lla = georef.spherical_to_proj(r, data['azi'], data['ele'], (-75.295893, -12.040436, 3379.2147))
615
616 latlon = antenna_to_geographic(r, az, el, (-75.295893, -12.040436))
617
618 if self.mask:
619 meta['lat'] = latlon[1][var.mask==False]
620 meta['lon'] = latlon[0][var.mask==False]
621 data['var'] = numpy.array([var[var.mask==False]])
622 else:
623 meta['lat'] = latlon[1]
624 meta['lon'] = latlon[0]
625 data['var'] = numpy.array([var])
626 614
627 615 return data, meta
628 616
629 617 def plot(self):
630 618 data = self.data[-1]
631 619 z = data['data']
632 620 r = data['r']
633 621 self.titles = []
634 622
635 623 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
636 624 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
637 625 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
638 626 self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z)
639 627
640 628 if isinstance(data['mode_op'], bytes):
641 629 data['mode_op'] = data['mode_op'].decode()
642 630
643 631 if data['mode_op'] == 'RHI':
644 try:
645 if self.data['mode_op'][-2] == 'PPI':
646 self.ang_min = None
647 self.ang_max = None
648 except:
649 pass
650 self.ang_min = self.ang_min if self.ang_min else 0
651 self.ang_max = self.ang_max if self.ang_max else 90
652 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )
653 elif data['mode_op'] == 'PPI':
654 try:
655 if self.data['mode_op'][-2] == 'RHI':
656 self.ang_min = None
657 self.ang_max = None
658 except:
659 pass
660 self.ang_min = self.ang_min if self.ang_min else 0
661 self.ang_max = self.ang_max if self.ang_max else 360
662 r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )
632 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']))
633 len_aux = int(data['azi'].shape[0]/4)
634 mean = numpy.mean(data['azi'][len_aux:-len_aux])
635 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
636 elif data['mode_op'] == 'PPI':
637 r, theta = numpy.meshgrid(r, -numpy.radians(data['azi'])+numpy.pi/2)
638 len_aux = int(data['ele'].shape[0]/4)
639 mean = numpy.mean(data['ele'][len_aux:-len_aux])
640 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(mean)), r*numpy.sin(
641 theta)*numpy.cos(numpy.radians(mean))
642 x = km2deg(x) + -75.295893
643 y = km2deg(y) + -12.040436
663 644
664 645 self.clear_figures()
665 646
666 for i,ax in enumerate(self.axes):
667
668 if ax.firsttime:
669 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
670 ax.plt = ax.pcolormesh(theta, r, z[i], cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
671 if data['mode_op'] == 'PPI':
672 ax.set_theta_direction(-1)
673 ax.set_theta_offset(numpy.pi/2)
647 if data['mode_op'] == 'PPI':
648 axes = self.axes['PPI']
649 else:
650 axes = self.axes['RHI']
674 651
675 else:
676 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
677 ax.plt = ax.pcolormesh(theta, r, z[i], cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
678 if data['mode_op'] == 'PPI':
679 ax.set_theta_direction(-1)
680 ax.set_theta_offset(numpy.pi/2)
652 for i, ax in enumerate(axes):
653 if data['mode_op'] == 'PPI':
654 ax.set_extent([-75.745893, -74.845893, -12.490436, -11.590436])
655
656 ax.plt = ax.pcolormesh(x, y, z[i], cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
681 657
682 ax.grid(True)
683 658 if data['mode_op'] == 'RHI':
684 659 len_aux = int(data['azi'].shape[0]/4)
685 660 mean = numpy.mean(data['azi'][len_aux:-len_aux])
686 661 if len(self.channels) !=1:
687 662 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels]
688 663 else:
689 664 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
690 665 elif data['mode_op'] == 'PPI':
691 666 len_aux = int(data['ele'].shape[0]/4)
692 667 mean = numpy.mean(data['ele'][len_aux:-len_aux])
693 668 if len(self.channels) !=1:
694 669 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels]
695 670 else:
696 671 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
697 self.mode_value = round(mean,1) No newline at end of file
672 self.mode_value = round(mean,1)
673
674 if data['mode_op'] == 'PPI':
675 gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
676 linewidth=1, color='gray', alpha=0.5, linestyle='--')
677 gl.xlabel_style = {'size': 8}
678 gl.ylabel_style = {'size': 8}
679 gl.xlabels_top = False
680 gl.ylabels_right = False
681 shape_p = os.path.join(self.shapes,'PER_ADM2/PER_ADM2.shp')
682 shape_d = os.path.join(self.shapes,'PER_ADM1/PER_ADM1.shp')
683 capitales = os.path.join(self.shapes,'CAPITALES/cap_provincia.shp')
684 vias = os.path.join(self.shapes,'Carreteras/VIAS_NACIONAL_250000.shp')
685 reader_d = shpreader.BasicReader(shape_p, encoding='latin1')
686 reader_p = shpreader.BasicReader(shape_d, encoding='latin1')
687 reader_c = shpreader.BasicReader(capitales, encoding='latin1')
688 reader_v = shpreader.BasicReader(vias, encoding='latin1')
689 caps = [x for x in reader_c.records() if x.attributes["Departa"] in ("JUNIN", "LIMA", "AYACUCHO", "HUANCAVELICA")]
690 districts = [x for x in reader_d.records() if x.attributes["Name"] in ("JUNÍN", "CHANCHAMAYO", "CHUPACA", "CONCEPCIÓN", "HUANCAYO", "JAUJA", "SATIPO", "TARMA", "YAUYOS", "HUAROCHIRÍ", "CANTA", "HUANTA", "TAYACAJA")]
691 provs = [x for x in reader_p.records() if x.attributes["NAME"] in ("Junín", "Lima")]
692 vias = [x for x in reader_v.records() if x.attributes["DEP"] in ("JUNIN", "LIMA")]
693
694 # Display Kenya's shape
695 shape_feature = ShapelyFeature([x.geometry for x in districts], ccrs.PlateCarree(), facecolor="none", edgecolor='grey', lw=0.5)
696 ax.add_feature(shape_feature)
697 shape_feature = ShapelyFeature([x.geometry for x in provs], ccrs.PlateCarree(), facecolor="none", edgecolor='white', lw=1)
698 ax.add_feature(shape_feature)
699 shape_feature = ShapelyFeature([x.geometry for x in vias], ccrs.PlateCarree(), facecolor="none", edgecolor='yellow', lw=1)
700 ax.add_feature(shape_feature)
701
702 for cap in caps:
703 if cap.attributes['Nombre'] in ("LA OROYA", "CONCEPCIÓN", "HUANCAYO", "JAUJA", "CHUPACA", "YAUYOS", "HUANTA", "PAMPAS"):
704 ax.text(cap.attributes['X'], cap.attributes['Y'], cap.attributes['Nombre'].title(), size=7, color='white')
705 ax.text(-75.052003, -11.915552, 'Huaytapallana', size=7, color='cyan')
706 ax.plot(-75.052003, -11.915552, '*')
707
708 for R in (10, 20, 30 , 40, 50):
709 circle = Circle((-75.295893, -12.040436), km2deg(R), facecolor='none',
710 edgecolor='skyblue', linewidth=1, alpha=0.5)
711 ax.add_patch(circle)
712 ax.text(km2deg(R)*numpy.cos(numpy.radians(45))-75.295893,
713 km2deg(R)*numpy.sin(numpy.radians(45))-12.040436,
714 '{}km'.format(R), color='skyblue', size=7)
715 elif data['mode_op'] == 'RHI':
716 ax.grid(color='grey', alpha=0.5, linestyle='--', linewidth=1)
Show file before | Show file after | Hide comments
@@ -1,733 +1,739
1 1 from email.utils import localtime
2 2 import os
3 3 import time
4 4 import datetime
5 5
6 6 import numpy
7 7 import h5py
8 8
9 9 import schainpy.admin
10 10 from schainpy.model.data.jrodata import *
11 11 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
12 12 from schainpy.model.io.jroIO_base import *
13 13 from schainpy.utils import log
14 14
15 15
16 16 class HDFReader(Reader, ProcessingUnit):
17 17 """Processing unit to read HDF5 format files
18 18
19 19 This unit reads HDF5 files created with `HDFWriter` operation contains
20 20 by default two groups Data and Metadata all variables would be saved as `dataOut`
21 21 attributes.
22 22 It is possible to read any HDF5 file by given the structure in the `description`
23 23 parameter, also you can add extra values to metadata with the parameter `extras`.
24 24
25 25 Parameters:
26 26 -----------
27 27 path : str
28 28 Path where files are located.
29 29 startDate : date
30 30 Start date of the files
31 31 endDate : list
32 32 End date of the files
33 33 startTime : time
34 34 Start time of the files
35 35 endTime : time
36 36 End time of the files
37 37 description : dict, optional
38 38 Dictionary with the description of the HDF5 file
39 39 extras : dict, optional
40 40 Dictionary with extra metadata to be be added to `dataOut`
41 41
42 42 Examples
43 43 --------
44 44
45 45 desc = {
46 46 'Data': {
47 47 'data_output': ['u', 'v', 'w'],
48 48 'utctime': 'timestamps',
49 49 } ,
50 50 'Metadata': {
51 51 'heightList': 'heights'
52 52 }
53 53 }
54 54
55 55 desc = {
56 56 'Data': {
57 57 'data_output': 'winds',
58 58 'utctime': 'timestamps'
59 59 },
60 60 'Metadata': {
61 61 'heightList': 'heights'
62 62 }
63 63 }
64 64
65 65 extras = {
66 66 'timeZone': 300
67 67 }
68 68
69 69 reader = project.addReadUnit(
70 70 name='HDFReader',
71 71 path='/path/to/files',
72 72 startDate='2019/01/01',
73 73 endDate='2019/01/31',
74 74 startTime='00:00:00',
75 75 endTime='23:59:59',
76 76 # description=json.dumps(desc),
77 77 # extras=json.dumps(extras),
78 78 )
79 79
80 80 """
81 81
82 82 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
83 83
84 84 def __init__(self):
85 85 ProcessingUnit.__init__(self)
86 86 self.dataOut = Parameters()
87 87 self.ext = ".hdf5"
88 88 self.optchar = "D"
89 89 self.meta = {}
90 90 self.data = {}
91 91 self.open_file = h5py.File
92 92 self.open_mode = 'r'
93 93 self.description = {}
94 94 self.extras = {}
95 95 self.filefmt = "*%Y%j***"
96 96 self.folderfmt = "*%Y%j"
97 97 self.utcoffset = 0
98 98 self.filter = None
99 99 self.dparam = None
100 100
101 101 def setup(self, **kwargs):
102 102
103 103 self.set_kwargs(**kwargs)
104 104 if not self.ext.startswith('.'):
105 105 self.ext = '.{}'.format(self.ext)
106 106
107 107 if self.online:
108 108 log.log("Searching files in online mode...", self.name)
109 109
110 110 for nTries in range(self.nTries):
111 111 fullpath = self.searchFilesOnLine(self.path, self.startDate,
112 112 self.endDate, self.expLabel, self.ext, self.walk,
113 113 self.filefmt, self.folderfmt,self.filter)
114 114 try:
115 115 fullpath = next(fullpath)
116 116 except:
117 117 fullpath = None
118 118
119 119 if fullpath:
120 120 break
121 121
122 122 log.warning(
123 123 'Waiting {} sec for a valid file in {}: try {} ...'.format(
124 124 self.delay, self.path, nTries + 1),
125 125 self.name)
126 126 time.sleep(self.delay)
127 127
128 128 if not(fullpath):
129 129 raise schainpy.admin.SchainError(
130 130 'There isn\'t any valid file in {}'.format(self.path))
131 131
132 132 pathname, filename = os.path.split(fullpath)
133 133 self.year = int(filename[1:5])
134 134 self.doy = int(filename[5:8])
135 135 self.set = int(filename[8:11]) - 1
136 136 else:
137 137 log.log("Searching files in {}".format(self.path), self.name)
138 138 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
139 139 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt,self.filter)
140 140
141 141 self.setNextFile()
142 142
143 143 return
144 144
145 145 def readFirstHeader(self):
146 146 '''Read metadata and data'''
147 147
148 148 self.__readMetadata()
149 149 self.__readData()
150 150 self.__setBlockList()
151 151
152 152 if 'type' in self.meta:
153 153 self.dataOut = eval(self.meta['type'])()
154 154
155 155 if self.dparam:
156 156 setattr(self.dataOut, "dparam", 1)
157 157
158 158 for attr in self.meta:
159 159 setattr(self.dataOut, attr, self.meta[attr])
160 160
161 161 self.blockIndex = 0
162 162
163 163 return
164 164
165 165 def __setBlockList(self):
166 166 '''
167 167 Selects the data within the times defined
168 168
169 169 self.fp
170 170 self.startTime
171 171 self.endTime
172 172 self.blockList
173 173 self.blocksPerFile
174 174
175 175 '''
176 176
177 177 startTime = self.startTime
178 178 endTime = self.endTime
179 179 thisUtcTime = self.data['utctime'] + self.utcoffset
180 180 try:
181 181 self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
182 182 except:
183 183 self.interval = 0
184 184 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
185 185
186 186 thisDate = thisDatetime.date()
187 187 thisTime = thisDatetime.time()
188 188
189 189 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
190 190 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
191 191
192 192 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
193 193
194 194 self.blockList = ind
195 195 self.blocksPerFile = len(ind)
196 196 return
197 197
198 198 def __readMetadata(self):
199 199 '''
200 200 Reads Metadata
201 201 '''
202 202
203 203 meta = {}
204 204
205 205 if self.description:
206 206 for key, value in self.description['Metadata'].items():
207 207 meta[key] = self.fp[value][()]
208 208 else:
209 209 grp = self.fp['Metadata']
210 210 for name in grp:
211 211 meta[name] = grp[name][()]
212 212
213 213 if self.extras:
214 214 for key, value in self.extras.items():
215 215 meta[key] = value
216 216 self.meta = meta
217 217
218 218 return
219 219
220 220 def __readData(self):
221 221
222 222 data = {}
223 223
224 224 if self.description:
225 225 for key, value in self.description['Data'].items():
226 226 if isinstance(value, str):
227 227 if isinstance(self.fp[value], h5py.Dataset):
228 228 data[key] = self.fp[value][()]
229 229 elif isinstance(self.fp[value], h5py.Group):
230 230 array = []
231 231 for ch in self.fp[value]:
232 232 array.append(self.fp[value][ch][()])
233 233 data[key] = numpy.array(array)
234 234 elif isinstance(value, list):
235 235 array = []
236 236 for ch in value:
237 237 array.append(self.fp[ch][()])
238 238 data[key] = numpy.array(array)
239 239 else:
240 240 grp = self.fp['Data']
241 241 for name in grp:
242 242 if isinstance(grp[name], h5py.Dataset):
243 243 array = grp[name][()]
244 244 elif isinstance(grp[name], h5py.Group):
245 245 array = []
246 246 for ch in grp[name]:
247 247 array.append(grp[name][ch][()])
248 248 array = numpy.array(array)
249 249 else:
250 250 log.warning('Unknown type: {}'.format(name))
251 251
252 252 if name in self.description:
253 253 key = self.description[name]
254 254 else:
255 255 key = name
256 256 data[key] = array
257 257
258 258 self.data = data
259 259 return
260 260
261 261 def getData(self):
262 262
263 263 for attr in self.data:
264 264 if self.data[attr].ndim == 1:
265 265 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
266 266 else:
267 267 if self.dparam:
268 268 setattr(self.dataOut, attr, self.data[attr])
269 269 else:
270 270 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
271 271
272 272 self.dataOut.flagNoData = False
273 273 self.blockIndex += 1
274 274
275 275 log.log("Block No. {}/{} -> {}".format(
276 276 self.blockIndex,
277 277 self.blocksPerFile,
278 278 self.dataOut.datatime.ctime()), self.name)
279 279
280 280 return
281 281
282 282 def run(self, **kwargs):
283 283
284 284 if not(self.isConfig):
285 285 self.setup(**kwargs)
286 286 self.isConfig = True
287 287
288 288 if self.blockIndex == self.blocksPerFile:
289 289 self.setNextFile()
290 290
291 291 self.getData()
292 292
293 293 return
294 294
295 295 @MPDecorator
296 296 class HDFWriter(Operation):
297 297 """Operation to write HDF5 files.
298 298
299 299 The HDF5 file contains by default two groups Data and Metadata where
300 300 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
301 301 parameters, data attributes are normaly time dependent where the metadata
302 302 are not.
303 303 It is possible to customize the structure of the HDF5 file with the
304 304 optional description parameter see the examples.
305 305
306 306 Parameters:
307 307 -----------
308 308 path : str
309 309 Path where files will be saved.
310 310 blocksPerFile : int
311 311 Number of blocks per file
312 312 metadataList : list
313 313 List of the dataOut attributes that will be saved as metadata
314 314 dataList : int
315 315 List of the dataOut attributes that will be saved as data
316 316 setType : bool
317 317 If True the name of the files corresponds to the timestamp of the data
318 318 description : dict, optional
319 319 Dictionary with the desired description of the HDF5 file
320 320
321 321 Examples
322 322 --------
323 323
324 324 desc = {
325 325 'data_output': {'winds': ['z', 'w', 'v']},
326 326 'utctime': 'timestamps',
327 327 'heightList': 'heights'
328 328 }
329 329 desc = {
330 330 'data_output': ['z', 'w', 'v'],
331 331 'utctime': 'timestamps',
332 332 'heightList': 'heights'
333 333 }
334 334 desc = {
335 335 'Data': {
336 336 'data_output': 'winds',
337 337 'utctime': 'timestamps'
338 338 },
339 339 'Metadata': {
340 340 'heightList': 'heights'
341 341 }
342 342 }
343 343
344 344 writer = proc_unit.addOperation(name='HDFWriter')
345 345 writer.addParameter(name='path', value='/path/to/file')
346 346 writer.addParameter(name='blocksPerFile', value='32')
347 347 writer.addParameter(name='metadataList', value='heightList,timeZone')
348 348 writer.addParameter(name='dataList',value='data_output,utctime')
349 349 # writer.addParameter(name='description',value=json.dumps(desc))
350 350
351 351 """
352 352
353 353 ext = ".hdf5"
354 354 optchar = "D"
355 355 filename = None
356 356 path = None
357 357 setFile = None
358 358 fp = None
359 359 firsttime = True
360 360 #Configurations
361 361 blocksPerFile = None
362 362 blockIndex = None
363 363 dataOut = None
364 364 #Data Arrays
365 365 dataList = None
366 366 metadataList = None
367 367 currentDay = None
368 368 lastTime = None
369 369 last_Azipos = None
370 370 last_Elepos = None
371 371 mode = None
372 372 #-----------------------
373 373 Typename = None
374 374 mask = False
375 375
376 376 def __init__(self):
377 377
378 378 Operation.__init__(self)
379 379 return
380 380
381 381 def set_kwargs(self, **kwargs):
382 382
383 383 for key, value in kwargs.items():
384 384 setattr(self, key, value)
385 385
386 386 def set_kwargs_obj(self,obj, **kwargs):
387 387
388 388 for key, value in kwargs.items():
389 389 setattr(obj, key, value)
390 390
391 391 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None,type_data=None, localtime=True, **kwargs):
392 392 self.path = path
393 393 self.blocksPerFile = blocksPerFile
394 394 self.metadataList = metadataList
395 395 self.dataList = [s.strip() for s in dataList]
396 396 self.setType = setType
397 397 if self.setType == "weather":
398 398 self.set_kwargs(**kwargs)
399 399 self.set_kwargs_obj(self.dataOut,**kwargs)
400 400 self.weather_vars = {
401 401 'S' : 0,
402 402 'V' : 1,
403 403 'W' : 2,
404 404 'SNR' : 3,
405 405 'Z' : 4,
406 406 'D' : 5,
407 407 'P' : 6,
408 408 'R' : 7,
409 409 }
410 410
411 411 if localtime:
412 412 self.getDateTime = datetime.datetime.fromtimestamp
413 413 else:
414 414 self.getDateTime = datetime.datetime.utcfromtimestamp
415 415
416 416 self.description = description
417 417 self.type_data=type_data
418 418
419 419 if self.metadataList is None:
420 420 self.metadataList = self.dataOut.metadata_list
421 421
422 422 dsList = []
423 423
424 424 for i in range(len(self.dataList)):
425 425 dsDict = {}
426 426 if hasattr(self.dataOut, self.dataList[i]):
427 427 dataAux = getattr(self.dataOut, self.dataList[i])
428 428 if self.setType == 'weather' and self.dataList[i] == 'data_param':
429 429 dataAux = dataAux[:,self.weather_vars[self.weather_var],:]
430 430 dsDict['variable'] = self.dataList[i]
431 431 else:
432 432 log.warning('Attribute {} not found in dataOut'.format(self.dataList[i]), self.name)
433 433 continue
434 434
435 435 if dataAux is None:
436 436 continue
437 437 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
438 438 dsDict['nDim'] = 0
439 439 else:
440 440 dsDict['nDim'] = len(dataAux.shape)
441 441 dsDict['shape'] = dataAux.shape
442 442 dsDict['dsNumber'] = dataAux.shape[0]
443 443 dsDict['dtype'] = dataAux.dtype
444 444 dsList.append(dsDict)
445 445
446 446 self.dsList = dsList
447 447 self.currentDay = self.dataOut.datatime.date()
448 448
449 449 def timeFlag(self):
450 450 currentTime = self.dataOut.utctime
451 451 dt = self.getDateTime(currentTime)
452 452
453 453 dataDay = int(dt.strftime('%j'))
454 454
455 455 if self.lastTime is None:
456 456 self.lastTime = currentTime
457 457 self.currentDay = dataDay
458 458 return False
459 459
460 460 timeDiff = currentTime - self.lastTime
461 461
462 462 #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
463 463 if dataDay != self.currentDay:
464 464 self.currentDay = dataDay
465 465 return True
466 466 elif timeDiff > 3*60*60:
467 467 self.lastTime = currentTime
468 468 return True
469 469 else:
470 470 self.lastTime = currentTime
471 471 return False
472 472
473 473 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
474 474 dataList=[], setType=None, description={}, mode= None,
475 475 type_data=None, Reset = False, localtime=True, **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
486 486 if not(self.isConfig):
487 487 self.setup(path=path, blocksPerFile=blocksPerFile,
488 488 metadataList=metadataList, dataList=dataList,
489 489 setType=setType, description=description,type_data=type_data,
490 490 localtime=localtime, **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
500 500 ext = self.ext
501 501 path = self.path
502 502 setFile = self.setFile
503 503
504 504 dt = self.getDateTime(self.dataOut.utctime)
505 505
506 506 if self.setType == 'weather':
507 507 subfolder = dt.strftime('%Y-%m-%dT%H-00-00')
508 subfolder = ''
508 509 else:
509 510 subfolder = dt.strftime('d%Y%j')
510 511
511 512 fullpath = os.path.join(path, subfolder)
512 513
513 514 if os.path.exists(fullpath):
514 515 filesList = os.listdir(fullpath)
515 516 filesList = [k for k in filesList if k.startswith(self.optchar)]
516 517 if len( filesList ) > 0:
517 518 filesList = sorted(filesList, key=str.lower)
518 519 filen = filesList[-1]
519 520 # el filename debera tener el siguiente formato
520 521 # 0 1234 567 89A BCDE (hex)
521 522 # x YYYY DDD SSS .ext
522 523 if isNumber(filen[8:11]):
523 524 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
524 525 else:
525 526 setFile = -1
526 527 else:
527 528 setFile = -1 #inicializo mi contador de seteo
528 529 else:
529 530 os.makedirs(fullpath)
530 531 setFile = -1 #inicializo mi contador de seteo
531 532
532 533 if self.setType is None:
533 534 setFile += 1
534 535 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
535 536 dt.year,
536 537 int(dt.strftime('%j')),
537 538 setFile,
538 539 ext )
539 540 elif self.setType == "weather":
540 541
541 542 #SOPHY_20200505_140215_E10.0_Z.h5
542 543 #SOPHY_20200505_140215_A40.0_Z.h5
543 544 if self.dataOut.flagMode == 1: #'AZI' #PPI
544 ang_type = 'E'
545 ang_type = 'EL'
546 mode_type = 'PPI'
545 547 len_aux = int(self.dataOut.data_ele.shape[0]/4)
546 548 mean = numpy.mean(self.dataOut.data_ele[len_aux:-len_aux])
547 549 ang_ = round(mean,1)
548 550 elif self.dataOut.flagMode == 0: #'ELE' #RHI
549 ang_type = 'A'
551 ang_type = 'AZ'
552 mode_type = 'RHI'
550 553 len_aux = int(self.dataOut.data_azi.shape[0]/4)
551 554 mean = numpy.mean(self.dataOut.data_azi[len_aux:-len_aux])
552 555 ang_ = round(mean,1)
553 556
554 557 file = '%s_%2.2d%2.2d%2.2d_%2.2d%2.2d%2.2d_%s%2.1f_%s%s' % (
555 558 'SOPHY',
556 559 dt.year,
557 560 dt.month,
558 561 dt.day,
559 562 dt.hour,
560 563 dt.minute,
561 564 dt.second,
562 ang_type,
565 ang_type[0],
563 566 ang_,
564 567 self.weather_var,
565 568 ext )
566
569 subfolder = '{}_{}_{}_{:2.1f}'.format(self.weather_var, mode_type, ang_type, ang_)
570 fullpath = os.path.join(path, subfolder)
571 if not os.path.exists(fullpath):
572 os.makedirs(fullpath)
567 573 else:
568 574 setFile = dt.hour*60+dt.minute
569 575 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
570 576 dt.year,
571 577 int(dt.strftime('%j')),
572 578 setFile,
573 579 ext )
574 580
575 581 self.filename = os.path.join( path, subfolder, file )
576 582
577 583 self.fp = h5py.File(self.filename, 'w')
578 584 #write metadata
579 585 self.writeMetadata(self.fp)
580 586 #Write data
581 587 self.writeData(self.fp)
582 588
583 589 def getLabel(self, name, x=None):
584 590
585 591 if x is None:
586 592 if 'Data' in self.description:
587 593 data = self.description['Data']
588 594 if 'Metadata' in self.description:
589 595 data.update(self.description['Metadata'])
590 596 else:
591 597 data = self.description
592 598 if name in data:
593 599 if isinstance(data[name], str):
594 600 return data[name]
595 601 elif isinstance(data[name], list):
596 602 return None
597 603 elif isinstance(data[name], dict):
598 604 for key, value in data[name].items():
599 605 return key
600 606 return name
601 607 else:
602 608 if 'Data' in self.description:
603 609 data = self.description['Data']
604 610 if 'Metadata' in self.description:
605 611 data.update(self.description['Metadata'])
606 612 else:
607 613 data = self.description
608 614 if name in data:
609 615 if isinstance(data[name], list):
610 616 return data[name][x]
611 617 elif isinstance(data[name], dict):
612 618 for key, value in data[name].items():
613 619 return value[x]
614 620 if 'cspc' in name:
615 621 return 'pair{:02d}'.format(x)
616 622 else:
617 623 return 'channel{:02d}'.format(x)
618 624
619 625 def writeMetadata(self, fp):
620 626
621 627 if self.description:
622 628 if 'Metadata' in self.description:
623 629 grp = fp.create_group('Metadata')
624 630 else:
625 631 grp = fp
626 632 else:
627 633 grp = fp.create_group('Metadata')
628 634
629 635 for i in range(len(self.metadataList)):
630 636 if not hasattr(self.dataOut, self.metadataList[i]):
631 637 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
632 638 continue
633 639 value = getattr(self.dataOut, self.metadataList[i])
634 640 if isinstance(value, bool):
635 641 if value is True:
636 642 value = 1
637 643 else:
638 644 value = 0
639 645 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
640 646 return
641 647
642 648 def writeData(self, fp):
643 649
644 650 if self.description:
645 651 if 'Data' in self.description:
646 652 grp = fp.create_group('Data')
647 653 else:
648 654 grp = fp
649 655 else:
650 656 grp = fp.create_group('Data')
651 657
652 658 dtsets = []
653 659 data = []
654 660
655 661 for dsInfo in self.dsList:
656 662
657 663 if dsInfo['nDim'] == 0:
658 664 ds = grp.create_dataset(
659 665 self.getLabel(dsInfo['variable']),
660 666 (self.blocksPerFile, ),
661 667 chunks=True,
662 668 dtype=numpy.float64)
663 669 dtsets.append(ds)
664 670 data.append((dsInfo['variable'], -1))
665 671 else:
666 672 label = self.getLabel(dsInfo['variable'])
667 673 if label is not None:
668 674 sgrp = grp.create_group(label)
669 675 else:
670 676 sgrp = grp
671 677 if self.blocksPerFile == 1:
672 678 shape = dsInfo['shape'][1:]
673 679 else:
674 680 shape = (self.blocksPerFile, ) + dsInfo['shape'][1:]
675 681 for i in range(dsInfo['dsNumber']):
676 682 ds = sgrp.create_dataset(
677 683 self.getLabel(dsInfo['variable'], i),
678 684 shape,
679 685 chunks=True,
680 686 dtype=dsInfo['dtype'],
681 687 compression='gzip',
682 688 )
683 689 dtsets.append(ds)
684 690 data.append((dsInfo['variable'], i))
685 691 fp.flush()
686 692
687 693 log.log('Creating file: {}'.format(fp.filename), self.name)
688 694
689 695 self.ds = dtsets
690 696 self.data = data
691 697 self.firsttime = True
692 698 self.blockIndex = 0
693 699 return
694 700
695 701 def putData(self):
696 702
697 703 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
698 704 self.closeFile()
699 705 self.setNextFile()
700 706
701 707 for i, ds in enumerate(self.ds):
702 708 attr, ch = self.data[i]
703 709 if ch == -1:
704 710 ds[self.blockIndex] = getattr(self.dataOut, attr)
705 711 else:
706 712 if self.blocksPerFile == 1:
707 713 mask = self.dataOut.data_param[:,3,:][ch] < self.mask
708 714 tmp = getattr(self.dataOut, attr)[:,self.weather_vars[self.weather_var],:][ch]
709 715 if self.mask:
710 716 tmp[mask] = numpy.nan
711 717 ds[:] = tmp
712 718 else:
713 719 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
714 720
715 721 self.fp.flush()
716 722 self.blockIndex += 1
717 723 log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
718 724
719 725 return
720 726
721 727 def closeFile(self):
722 728
723 729 if self.blockIndex != self.blocksPerFile:
724 730 for ds in self.ds:
725 731 ds.resize(self.blockIndex, axis=0)
726 732
727 733 if self.fp:
728 734 self.fp.flush()
729 735 self.fp.close()
730 736
731 737 def close(self):
732 738
733 739 self.closeFile()
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now