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