##// END OF EJS Templates
schain-jc
RSS
Merge branch 'v3.0-devel' of http://jro-dev.igp.gob.pe/rhodecode/schain into v3.0-devel
George Yong -
r1190:29b77fb75833 merge
parent child
Show More
Side by Side Unified
Context file:
r1190:29b77fb75833 Loading diff...:
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (803 lines changed) Show them Hide them
@@ -0,0 +1,803
1
2 import os
3 import sys
4 import zmq
5 import time
6 import datetime
7 from functools import wraps
8 import numpy
9 import matplotlib
10
11 if 'BACKEND' in os.environ:
12 matplotlib.use(os.environ['BACKEND'])
13 elif 'linux' in sys.platform:
14 matplotlib.use("TkAgg")
15 elif 'darwin' in sys.platform:
16 matplotlib.use('TkAgg')
17 else:
18 from schainpy.utils import log
19 log.warning('Using default Backend="Agg"', 'INFO')
20 matplotlib.use('Agg')
21
22 import matplotlib.pyplot as plt
23 from matplotlib.patches import Polygon
24 from mpl_toolkits.axes_grid1 import make_axes_locatable
25 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
26
27 from schainpy.model.data.jrodata import PlotterData
28 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
29 from schainpy.utils import log
30
31 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
32 blu_values = matplotlib.pyplot.get_cmap(
33 'seismic_r', 20)(numpy.arange(20))[10:15]
34 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
35 'jro', numpy.vstack((blu_values, jet_values)))
36 matplotlib.pyplot.register_cmap(cmap=ncmap)
37
38 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
39 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
40
41 EARTH_RADIUS = 6.3710e3
42
43
44 def ll2xy(lat1, lon1, lat2, lon2):
45
46 p = 0.017453292519943295
47 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
48 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
49 r = 12742 * numpy.arcsin(numpy.sqrt(a))
50 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
51 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
52 theta = -theta + numpy.pi/2
53 return r*numpy.cos(theta), r*numpy.sin(theta)
54
55
56 def km2deg(km):
57 '''
58 Convert distance in km to degrees
59 '''
60
61 return numpy.rad2deg(km/EARTH_RADIUS)
62
63
64 def figpause(interval):
65 backend = plt.rcParams['backend']
66 if backend in matplotlib.rcsetup.interactive_bk:
67 figManager = matplotlib._pylab_helpers.Gcf.get_active()
68 if figManager is not None:
69 canvas = figManager.canvas
70 if canvas.figure.stale:
71 canvas.draw()
72 try:
73 canvas.start_event_loop(interval)
74 except:
75 pass
76 return
77
78
79 def popup(message):
80 '''
81 '''
82
83 fig = plt.figure(figsize=(12, 8), facecolor='r')
84 text = '\n'.join([s.strip() for s in message.split(':')])
85 fig.text(0.01, 0.5, text, ha='left', va='center',
86 size='20', weight='heavy', color='w')
87 fig.show()
88 figpause(1000)
89
90
91 class Throttle(object):
92 '''
93 Decorator that prevents a function from being called more than once every
94 time period.
95 To create a function that cannot be called more than once a minute, but
96 will sleep until it can be called:
97 @Throttle(minutes=1)
98 def foo():
99 pass
100
101 for i in range(10):
102 foo()
103 print "This function has run %s times." % i
104 '''
105
106 def __init__(self, seconds=0, minutes=0, hours=0):
107 self.throttle_period = datetime.timedelta(
108 seconds=seconds, minutes=minutes, hours=hours
109 )
110
111 self.time_of_last_call = datetime.datetime.min
112
113 def __call__(self, fn):
114 @wraps(fn)
115 def wrapper(*args, **kwargs):
116 coerce = kwargs.pop('coerce', None)
117 if coerce:
118 self.time_of_last_call = datetime.datetime.now()
119 return fn(*args, **kwargs)
120 else:
121 now = datetime.datetime.now()
122 time_since_last_call = now - self.time_of_last_call
123 time_left = self.throttle_period - time_since_last_call
124
125 if time_left > datetime.timedelta(seconds=0):
126 return
127
128 self.time_of_last_call = datetime.datetime.now()
129 return fn(*args, **kwargs)
130
131 return wrapper
132
133 def apply_throttle(value):
134
135 @Throttle(seconds=value)
136 def fnThrottled(fn):
137 fn()
138
139 return fnThrottled
140
141 @MPDecorator
142 class Plotter(ProcessingUnit):
143 '''
144 Proccessing unit to handle plot operations
145 '''
146
147 def __init__(self):
148
149 ProcessingUnit.__init__(self)
150
151 def setup(self, **kwargs):
152
153 self.connections = 0
154 self.web_address = kwargs.get('web_server', False)
155 self.realtime = kwargs.get('realtime', False)
156 self.localtime = kwargs.get('localtime', True)
157 self.buffering = kwargs.get('buffering', True)
158 self.throttle = kwargs.get('throttle', 2)
159 self.exp_code = kwargs.get('exp_code', None)
160 self.set_ready = apply_throttle(self.throttle)
161 self.dates = []
162 self.data = PlotterData(
163 self.plots, self.throttle, self.exp_code, self.buffering)
164 self.isConfig = True
165
166 def ready(self):
167 '''
168 Set dataOut ready
169 '''
170
171 self.data.ready = True
172 self.dataOut.data_plt = self.data
173
174 def run(self, realtime=True, localtime=True, buffering=True,
175 throttle=2, exp_code=None, web_server=None):
176
177 if not self.isConfig:
178 self.setup(realtime=realtime, localtime=localtime,
179 buffering=buffering, throttle=throttle, exp_code=exp_code,
180 web_server=web_server)
181
182 if self.web_address:
183 log.success(
184 'Sending to web: {}'.format(self.web_address),
185 self.name
186 )
187 self.context = zmq.Context()
188 self.sender_web = self.context.socket(zmq.REQ)
189 self.sender_web.connect(self.web_address)
190 self.poll = zmq.Poller()
191 self.poll.register(self.sender_web, zmq.POLLIN)
192 time.sleep(1)
193
194 # t = Thread(target=self.event_monitor, args=(monitor,))
195 # t.start()
196
197 self.dataOut = self.dataIn
198 self.data.ready = False
199
200 if self.dataOut.flagNoData:
201 coerce = True
202 else:
203 coerce = False
204
205 if self.dataOut.type == 'Parameters':
206 tm = self.dataOut.utctimeInit
207 else:
208 tm = self.dataOut.utctime
209 if self.dataOut.useLocalTime:
210 if not self.localtime:
211 tm += time.timezone
212 dt = datetime.datetime.fromtimestamp(tm).date()
213 else:
214 if self.localtime:
215 tm -= time.timezone
216 dt = datetime.datetime.utcfromtimestamp(tm).date()
217 if dt not in self.dates:
218 if self.data:
219 self.ready()
220 self.data.setup()
221 self.dates.append(dt)
222
223 self.data.update(self.dataOut, tm)
224
225 if False: # TODO check when publishers ends
226 self.connections -= 1
227 if self.connections == 0 and dt in self.dates:
228 self.data.ended = True
229 self.ready()
230 time.sleep(1)
231 else:
232 if self.realtime:
233 self.ready()
234 if self.web_address:
235 retries = 5
236 while True:
237 self.sender_web.send(self.data.jsonify())
238 socks = dict(self.poll.poll(5000))
239 if socks.get(self.sender_web) == zmq.POLLIN:
240 reply = self.sender_web.recv_string()
241 if reply == 'ok':
242 log.log("Response from server ok", self.name)
243 break
244 else:
245 log.warning(
246 "Malformed reply from server: {}".format(reply), self.name)
247
248 else:
249 log.warning(
250 "No response from server, retrying...", self.name)
251 self.sender_web.setsockopt(zmq.LINGER, 0)
252 self.sender_web.close()
253 self.poll.unregister(self.sender_web)
254 retries -= 1
255 if retries == 0:
256 log.error(
257 "Server seems to be offline, abandoning", self.name)
258 self.sender_web = self.context.socket(zmq.REQ)
259 self.sender_web.connect(self.web_address)
260 self.poll.register(self.sender_web, zmq.POLLIN)
261 time.sleep(1)
262 break
263 self.sender_web = self.context.socket(zmq.REQ)
264 self.sender_web.connect(self.web_address)
265 self.poll.register(self.sender_web, zmq.POLLIN)
266 time.sleep(1)
267 else:
268 self.set_ready(self.ready, coerce=coerce)
269
270 return
271
272 def close(self):
273 pass
274
275
276 @MPDecorator
277 class Plot(Operation):
278 '''
279 Base class for Schain plotting operations
280 '''
281
282 CODE = 'Figure'
283 colormap = 'jro'
284 bgcolor = 'white'
285 __missing = 1E30
286
287 __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax',
288 'zlimits', 'xlabel', 'ylabel', 'xaxis', 'cb_label', 'title',
289 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure',
290 'showprofile', 'decimation', 'pause']
291
292 def __init__(self):
293
294 Operation.__init__(self)
295 self.isConfig = False
296 self.isPlotConfig = False
297
298 def __fmtTime(self, x, pos):
299 '''
300 '''
301
302 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
303
304 def __setup(self, **kwargs):
305 '''
306 Initialize variables
307 '''
308
309 self.figures = []
310 self.axes = []
311 self.cb_axes = []
312 self.localtime = kwargs.pop('localtime', True)
313 self.show = kwargs.get('show', True)
314 self.save = kwargs.get('save', False)
315 self.ftp = kwargs.get('ftp', False)
316 self.colormap = kwargs.get('colormap', self.colormap)
317 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
318 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
319 self.colormaps = kwargs.get('colormaps', None)
320 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
321 self.showprofile = kwargs.get('showprofile', False)
322 self.title = kwargs.get('wintitle', self.CODE.upper())
323 self.cb_label = kwargs.get('cb_label', None)
324 self.cb_labels = kwargs.get('cb_labels', None)
325 self.labels = kwargs.get('labels', None)
326 self.xaxis = kwargs.get('xaxis', 'frequency')
327 self.zmin = kwargs.get('zmin', None)
328 self.zmax = kwargs.get('zmax', None)
329 self.zlimits = kwargs.get('zlimits', None)
330 self.xmin = kwargs.get('xmin', None)
331 self.xmax = kwargs.get('xmax', None)
332 self.xrange = kwargs.get('xrange', 12)
333 self.xscale = kwargs.get('xscale', None)
334 self.ymin = kwargs.get('ymin', None)
335 self.ymax = kwargs.get('ymax', None)
336 self.yscale = kwargs.get('yscale', None)
337 self.xlabel = kwargs.get('xlabel', None)
338 self.decimation = kwargs.get('decimation', None)
339 self.showSNR = kwargs.get('showSNR', False)
340 self.oneFigure = kwargs.get('oneFigure', True)
341 self.width = kwargs.get('width', None)
342 self.height = kwargs.get('height', None)
343 self.colorbar = kwargs.get('colorbar', True)
344 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
345 self.channels = kwargs.get('channels', None)
346 self.titles = kwargs.get('titles', [])
347 self.polar = False
348 self.grid = kwargs.get('grid', False)
349 self.pause = kwargs.get('pause', False)
350 self.save_labels = kwargs.get('save_labels', None)
351 self.realtime = kwargs.get('realtime', True)
352 self.buffering = kwargs.get('buffering', True)
353 self.throttle = kwargs.get('throttle', 2)
354 self.exp_code = kwargs.get('exp_code', None)
355 self.__throttle_plot = apply_throttle(self.throttle)
356 self.data = PlotterData(
357 self.CODE, self.throttle, self.exp_code, self.buffering)
358
359 def __setup_plot(self):
360 '''
361 Common setup for all figures, here figures and axes are created
362 '''
363
364 self.setup()
365
366 self.time_label = 'LT' if self.localtime else 'UTC'
367 if self.data.localtime:
368 self.getDateTime = datetime.datetime.fromtimestamp
369 else:
370 self.getDateTime = datetime.datetime.utcfromtimestamp
371
372 if self.width is None:
373 self.width = 8
374
375 self.figures = []
376 self.axes = []
377 self.cb_axes = []
378 self.pf_axes = []
379 self.cmaps = []
380
381 size = '15%' if self.ncols == 1 else '30%'
382 pad = '4%' if self.ncols == 1 else '8%'
383
384 if self.oneFigure:
385 if self.height is None:
386 self.height = 1.4 * self.nrows + 1
387 fig = plt.figure(figsize=(self.width, self.height),
388 edgecolor='k',
389 facecolor='w')
390 self.figures.append(fig)
391 for n in range(self.nplots):
392 ax = fig.add_subplot(self.nrows, self.ncols,
393 n + 1, polar=self.polar)
394 ax.tick_params(labelsize=8)
395 ax.firsttime = True
396 ax.index = 0
397 ax.press = None
398 self.axes.append(ax)
399 if self.showprofile:
400 cax = self.__add_axes(ax, size=size, pad=pad)
401 cax.tick_params(labelsize=8)
402 self.pf_axes.append(cax)
403 else:
404 if self.height is None:
405 self.height = 3
406 for n in range(self.nplots):
407 fig = plt.figure(figsize=(self.width, self.height),
408 edgecolor='k',
409 facecolor='w')
410 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
411 ax.tick_params(labelsize=8)
412 ax.firsttime = True
413 ax.index = 0
414 ax.press = None
415 self.figures.append(fig)
416 self.axes.append(ax)
417 if self.showprofile:
418 cax = self.__add_axes(ax, size=size, pad=pad)
419 cax.tick_params(labelsize=8)
420 self.pf_axes.append(cax)
421
422 for n in range(self.nrows):
423 if self.colormaps is not None:
424 cmap = plt.get_cmap(self.colormaps[n])
425 else:
426 cmap = plt.get_cmap(self.colormap)
427 cmap.set_bad(self.bgcolor, 1.)
428 self.cmaps.append(cmap)
429
430 for fig in self.figures:
431 fig.canvas.mpl_connect('key_press_event', self.OnKeyPress)
432 fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll)
433 fig.canvas.mpl_connect('button_press_event', self.onBtnPress)
434 fig.canvas.mpl_connect('motion_notify_event', self.onMotion)
435 fig.canvas.mpl_connect('button_release_event', self.onBtnRelease)
436 if self.show:
437 fig.show()
438
439 def OnKeyPress(self, event):
440 '''
441 Event for pressing keys (up, down) change colormap
442 '''
443 ax = event.inaxes
444 if ax in self.axes:
445 if event.key == 'down':
446 ax.index += 1
447 elif event.key == 'up':
448 ax.index -= 1
449 if ax.index < 0:
450 ax.index = len(CMAPS) - 1
451 elif ax.index == len(CMAPS):
452 ax.index = 0
453 cmap = CMAPS[ax.index]
454 ax.cbar.set_cmap(cmap)
455 ax.cbar.draw_all()
456 ax.plt.set_cmap(cmap)
457 ax.cbar.patch.figure.canvas.draw()
458 self.colormap = cmap.name
459
460 def OnBtnScroll(self, event):
461 '''
462 Event for scrolling, scale figure
463 '''
464 cb_ax = event.inaxes
465 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
466 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
467 pt = ax.cbar.ax.bbox.get_points()[:, 1]
468 nrm = ax.cbar.norm
469 vmin, vmax, p0, p1, pS = (
470 nrm.vmin, nrm.vmax, pt[0], pt[1], event.y)
471 scale = 2 if event.step == 1 else 0.5
472 point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0)
473 ax.cbar.norm.vmin = point - scale * (point - vmin)
474 ax.cbar.norm.vmax = point - scale * (point - vmax)
475 ax.plt.set_norm(ax.cbar.norm)
476 ax.cbar.draw_all()
477 ax.cbar.patch.figure.canvas.draw()
478
479 def onBtnPress(self, event):
480 '''
481 Event for mouse button press
482 '''
483 cb_ax = event.inaxes
484 if cb_ax is None:
485 return
486
487 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
488 cb_ax.press = event.x, event.y
489 else:
490 cb_ax.press = None
491
492 def onMotion(self, event):
493 '''
494 Event for move inside colorbar
495 '''
496 cb_ax = event.inaxes
497 if cb_ax is None:
498 return
499 if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]:
500 return
501 if cb_ax.press is None:
502 return
503
504 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
505 xprev, yprev = cb_ax.press
506 dx = event.x - xprev
507 dy = event.y - yprev
508 cb_ax.press = event.x, event.y
509 scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin
510 perc = 0.03
511
512 if event.button == 1:
513 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
514 ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy)
515 elif event.button == 3:
516 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
517 ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy)
518
519 ax.cbar.draw_all()
520 ax.plt.set_norm(ax.cbar.norm)
521 ax.cbar.patch.figure.canvas.draw()
522
523 def onBtnRelease(self, event):
524 '''
525 Event for mouse button release
526 '''
527 cb_ax = event.inaxes
528 if cb_ax is not None:
529 cb_ax.press = None
530
531 def __add_axes(self, ax, size='30%', pad='8%'):
532 '''
533 Add new axes to the given figure
534 '''
535 divider = make_axes_locatable(ax)
536 nax = divider.new_horizontal(size=size, pad=pad)
537 ax.figure.add_axes(nax)
538 return nax
539
540 def setup(self):
541 '''
542 This method should be implemented in the child class, the following
543 attributes should be set:
544
545 self.nrows: number of rows
546 self.ncols: number of cols
547 self.nplots: number of plots (channels or pairs)
548 self.ylabel: label for Y axes
549 self.titles: list of axes title
550
551 '''
552 raise NotImplementedError
553
554 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
555 '''
556 Create a masked array for missing data
557 '''
558 if x_buffer.shape[0] < 2:
559 return x_buffer, y_buffer, z_buffer
560
561 deltas = x_buffer[1:] - x_buffer[0:-1]
562 x_median = numpy.median(deltas)
563
564 index = numpy.where(deltas > 5 * x_median)
565
566 if len(index[0]) != 0:
567 z_buffer[::, index[0], ::] = self.__missing
568 z_buffer = numpy.ma.masked_inside(z_buffer,
569 0.99 * self.__missing,
570 1.01 * self.__missing)
571
572 return x_buffer, y_buffer, z_buffer
573
574 def decimate(self):
575
576 # dx = int(len(self.x)/self.__MAXNUMX) + 1
577 dy = int(len(self.y) / self.decimation) + 1
578
579 # x = self.x[::dx]
580 x = self.x
581 y = self.y[::dy]
582 z = self.z[::, ::, ::dy]
583
584 return x, y, z
585
586 def format(self):
587 '''
588 Set min and max values, labels, ticks and titles
589 '''
590
591 if self.xmin is None:
592 xmin = self.data.min_time
593 else:
594 if self.xaxis is 'time':
595 dt = self.getDateTime(self.data.min_time)
596 xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) -
597 datetime.datetime(1970, 1, 1)).total_seconds()
598 if self.data.localtime:
599 xmin += time.timezone
600 else:
601 xmin = self.xmin
602
603 if self.xmax is None:
604 xmax = xmin + self.xrange * 60 * 60
605 else:
606 if self.xaxis is 'time':
607 dt = self.getDateTime(self.data.max_time)
608 xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) -
609 datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds()
610 if self.data.localtime:
611 xmax += time.timezone
612 else:
613 xmax = self.xmax
614
615 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
616 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
617
618 Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])
619 i = 1 if numpy.where(
620 abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0]
621 ystep = Y[i] / 10.
622
623 if self.xaxis is not 'time':
624 X = numpy.array([1, 2, 5, 10, 20, 50, 100,
625 200, 500, 1000, 2000, 5000])/2.
626 i = 1 if numpy.where(
627 abs(xmax-xmin) <= X)[0][0] < 0 else numpy.where(abs(xmax-xmin) <= X)[0][0]
628 xstep = X[i] / 10.
629
630 for n, ax in enumerate(self.axes):
631 if ax.firsttime:
632 ax.set_facecolor(self.bgcolor)
633 ax.yaxis.set_major_locator(MultipleLocator(ystep))
634 if self.xscale:
635 ax.xaxis.set_major_formatter(FuncFormatter(
636 lambda x, pos: '{0:g}'.format(x*self.xscale)))
637 if self.xscale:
638 ax.yaxis.set_major_formatter(FuncFormatter(
639 lambda x, pos: '{0:g}'.format(x*self.yscale)))
640 if self.xaxis is 'time':
641 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
642 ax.xaxis.set_major_locator(LinearLocator(9))
643 else:
644 ax.xaxis.set_major_locator(MultipleLocator(xstep))
645 if self.xlabel is not None:
646 ax.set_xlabel(self.xlabel)
647 ax.set_ylabel(self.ylabel)
648 ax.firsttime = False
649 if self.showprofile:
650 self.pf_axes[n].set_ylim(ymin, ymax)
651 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
652 self.pf_axes[n].set_xlabel('dB')
653 self.pf_axes[n].grid(b=True, axis='x')
654 [tick.set_visible(False)
655 for tick in self.pf_axes[n].get_yticklabels()]
656 if self.colorbar:
657 ax.cbar = plt.colorbar(
658 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
659 ax.cbar.ax.tick_params(labelsize=8)
660 ax.cbar.ax.press = None
661 if self.cb_label:
662 ax.cbar.set_label(self.cb_label, size=8)
663 elif self.cb_labels:
664 ax.cbar.set_label(self.cb_labels[n], size=8)
665 else:
666 ax.cbar = None
667 if self.grid:
668 ax.grid(True)
669
670 if not self.polar:
671 ax.set_xlim(xmin, xmax)
672 ax.set_ylim(ymin, ymax)
673 ax.set_title('{} {} {}'.format(
674 self.titles[n],
675 self.getDateTime(self.data.max_time).strftime(
676 '%Y-%m-%dT%H:%M:%S'),
677 self.time_label),
678 size=8)
679 else:
680 ax.set_title('{}'.format(self.titles[n]), size=8)
681 ax.set_ylim(0, 90)
682 ax.set_yticks(numpy.arange(0, 90, 20))
683 ax.yaxis.labelpad = 40
684
685 def clear_figures(self):
686 '''
687 Reset axes for redraw plots
688 '''
689
690 for ax in self.axes:
691 ax.clear()
692 ax.firsttime = True
693 if ax.cbar:
694 ax.cbar.remove()
695
696 def __plot(self):
697 '''
698 Main function to plot, format and save figures
699 '''
700
701 #try:
702 self.plot()
703 self.format()
704 #except Exception as e:
705 # log.warning('{} Plot could not be updated... check data'.format(
706 # self.CODE), self.name)
707 # log.error(str(e), '')
708 # return
709
710 for n, fig in enumerate(self.figures):
711 if self.nrows == 0 or self.nplots == 0:
712 log.warning('No data', self.name)
713 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
714 fig.canvas.manager.set_window_title(self.CODE)
715 continue
716
717 fig.tight_layout()
718 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
719 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
720 fig.canvas.draw()
721
722 if self.save:
723
724 if self.save_labels:
725 labels = self.save_labels
726 else:
727 labels = list(range(self.nrows))
728
729 if self.oneFigure:
730 label = ''
731 else:
732 label = '-{}'.format(labels[n])
733 figname = os.path.join(
734 self.save,
735 self.CODE,
736 '{}{}_{}.png'.format(
737 self.CODE,
738 label,
739 self.getDateTime(self.data.max_time).strftime(
740 '%Y%m%d_%H%M%S'),
741 )
742 )
743 log.log('Saving figure: {}'.format(figname), self.name)
744 if not os.path.isdir(os.path.dirname(figname)):
745 os.makedirs(os.path.dirname(figname))
746 fig.savefig(figname)
747
748 def plot(self):
749 '''
750 Must be defined in the child class
751 '''
752 raise NotImplementedError
753
754 def run(self, dataOut, **kwargs):
755
756 if dataOut.flagNoData and not dataOut.error:
757 return dataOut
758
759 if dataOut.error:
760 coerce = True
761 else:
762 coerce = False
763
764 if self.isConfig is False:
765 self.__setup(**kwargs)
766 self.data.setup()
767 self.isConfig = True
768
769 if dataOut.type == 'Parameters':
770 tm = dataOut.utctimeInit
771 else:
772 tm = dataOut.utctime
773
774 if dataOut.useLocalTime:
775 if not self.localtime:
776 tm += time.timezone
777 else:
778 if self.localtime:
779 tm -= time.timezone
780
781 if self.data and (tm - self.data.min_time) >= self.xrange*60*60:
782 self.__plot()
783 self.data.setup()
784 self.clear_figures()
785
786 self.data.update(dataOut, tm)
787
788 if self.isPlotConfig is False:
789 self.__setup_plot()
790 self.isPlotConfig = True
791
792 if self.realtime:
793 self.__plot()
794 else:
795 self.__throttle_plot(self.__plot, coerce=coerce)
796
797 figpause(0.001)
798
799 def close(self):
800
801 if self.data and self.pause:
802 figpause(10)
803
Show file before | Show file after | Hide comments
@@ -112,5 +112,3 schainpy/scripts/
112 112 .vscode
113 113 trash
114 114 *.log
115 schainpy/scripts/testDigitalRF.py
116 schainpy/scripts/testDigitalRFWriter.py
Show file before | Show file after | Hide comments
@@ -1,5 +1,12
1 1 ## CHANGELOG:
2 2
3 ### 3.0
4 * Python 3.x compatible
5 * New architecture with multiprocessing and IPC communication
6 * Add @MPDecorator for multiprocessing Units and Operations
7 * Added new type of operation `external` for non-locking operations
8 * New plotting architecture with buffering/throttle capabilities to speed up plots
9
3 10 ### 2.3
4 11 * Added support for Madrigal formats (reading/writing).
5 12 * Added support for reading BLTR parameters (*.sswma).
Show file before | Show file after | Hide comments
@@ -24,7 +24,7 from email.mime.multipart import MIMEMultipart
24 24
25 25 import schainpy
26 26 from schainpy.utils import log
27 from schainpy.model.graphics.jroplot_data import popup
27 from schainpy.model.graphics.jroplot_base import popup
28 28
29 29 def get_path():
30 30 '''
Show file before | Show file after | Hide comments
@@ -97,9 +97,7 def wait(context):
97 97 receiver = c.socket(zmq.SUB)
98 98 receiver.connect('ipc:///tmp/schain_{}_pub'.format(self.id))
99 99 receiver.setsockopt(zmq.SUBSCRIBE, self.id.encode())
100 log.error('startinggg')
101 100 msg = receiver.recv_multipart()[1]
102 #log.error(msg)
103 101 context.terminate()
104 102
105 103 class ParameterConf():
@@ -1245,7 +1243,7 class Project(Process):
1245 1243
1246 1244 try:
1247 1245 zmq.proxy(xpub, xsub)
1248 except zmq.ContextTerminated:
1246 except: # zmq.ContextTerminated:
1249 1247 xpub.close()
1250 1248 xsub.close()
1251 1249
@@ -1260,6 +1258,6 class Project(Process):
1260 1258
1261 1259 # Iniciar todos los procesos .start(), monitoreo de procesos. ELiminar lo de abajo
1262 1260
1263 log.success('{} finished (time: {}s)'.format(
1261 log.success('{} Done (time: {}s)'.format(
1264 1262 self.name,
1265 1263 time.time()-self.start_time))
Show file before | Show file after | Hide comments
@@ -7,7 +7,9 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
7 7 import copy
8 8 import numpy
9 9 import datetime
10 import json
10 11
12 from schainpy.utils import log
11 13 from .jroheaderIO import SystemHeader, RadarControllerHeader
12 14
13 15
@@ -147,85 +149,48 class JROData(GenericData):
147 149 # m_ProcessingHeader = ProcessingHeader()
148 150
149 151 systemHeaderObj = SystemHeader()
150
151 152 radarControllerHeaderObj = RadarControllerHeader()
152
153 153 # data = None
154
155 154 type = None
156
157 155 datatype = None # dtype but in string
158
159 156 # dtype = None
160
161 157 # nChannels = None
162
163 158 # nHeights = None
164
165 159 nProfiles = None
166
167 160 heightList = None
168
169 161 channelList = None
170
171 162 flagDiscontinuousBlock = False
172
173 163 useLocalTime = False
174
175 164 utctime = None
176
177 165 timeZone = None
178
179 166 dstFlag = None
180
181 167 errorCount = None
182
183 168 blocksize = None
184
185 169 # nCode = None
186 #
187 170 # nBaud = None
188 #
189 171 # code = None
190
191 172 flagDecodeData = False # asumo q la data no esta decodificada
192
193 173 flagDeflipData = False # asumo q la data no esta sin flip
194
195 174 flagShiftFFT = False
196
197 175 # ippSeconds = None
198
199 176 # timeInterval = None
200
201 177 nCohInt = None
202
203 178 # noise = None
204
205 179 windowOfFilter = 1
206
207 180 # Speed of ligth
208 181 C = 3e8
209
210 182 frequency = 49.92e6
211
212 183 realtime = False
213
214 184 beacon_heiIndexList = None
215
216 185 last_block = None
217
218 186 blocknow = None
219
220 187 azimuth = None
221
222 188 zenith = None
223
224 189 beam = Beam()
225
226 190 profileIndex = None
227
228 error = (0, '')
191 error = None
192 data = None
193 nmodes = None
229 194
230 195 def __str__(self):
231 196
@@ -395,53 +360,29 class Voltage(JROData):
395 360 '''
396 361
397 362 self.useLocalTime = True
398
399 363 self.radarControllerHeaderObj = RadarControllerHeader()
400
401 364 self.systemHeaderObj = SystemHeader()
402
403 365 self.type = "Voltage"
404
405 366 self.data = None
406
407 367 # self.dtype = None
408
409 368 # self.nChannels = 0
410
411 369 # self.nHeights = 0
412
413 370 self.nProfiles = None
414
415 self.heightList = None
416
371 self.heightList = Non
417 372 self.channelList = None
418
419 373 # self.channelIndexList = None
420
421 374 self.flagNoData = True
422
423 375 self.flagDiscontinuousBlock = False
424
425 376 self.utctime = None
426
427 377 self.timeZone = None
428
429 378 self.dstFlag = None
430
431 379 self.errorCount = None
432
433 380 self.nCohInt = None
434
435 381 self.blocksize = None
436
437 382 self.flagDecodeData = False # asumo q la data no esta decodificada
438
439 383 self.flagDeflipData = False # asumo q la data no esta sin flip
440
441 384 self.flagShiftFFT = False
442
443 385 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
444
445 386 self.profileIndex = 0
446 387
447 388 def getNoisebyHildebrand(self, channel=None):
@@ -505,32 +446,20 class Spectra(JROData):
505 446
506 447 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
507 448 data_spc = None
508
509 449 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
510 450 data_cspc = None
511
512 451 # data dc es un numpy array de 2 dmensiones (canales, alturas)
513 452 data_dc = None
514
515 453 # data power
516 454 data_pwr = None
517
518 455 nFFTPoints = None
519
520 456 # nPairs = None
521
522 457 pairsList = None
523
524 458 nIncohInt = None
525
526 459 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
527
528 460 nCohInt = None # se requiere para determinar el valor de timeInterval
529
530 461 ippFactor = None
531
532 462 profileIndex = 0
533
534 463 plotting = "spectra"
535 464
536 465 def __init__(self):
@@ -539,59 +468,32 class Spectra(JROData):
539 468 '''
540 469
541 470 self.useLocalTime = True
542
543 471 self.radarControllerHeaderObj = RadarControllerHeader()
544
545 472 self.systemHeaderObj = SystemHeader()
546
547 473 self.type = "Spectra"
548
549 474 # self.data = None
550
551 475 # self.dtype = None
552
553 476 # self.nChannels = 0
554
555 477 # self.nHeights = 0
556
557 478 self.nProfiles = None
558
559 479 self.heightList = None
560
561 480 self.channelList = None
562
563 481 # self.channelIndexList = None
564
565 482 self.pairsList = None
566
567 483 self.flagNoData = True
568
569 484 self.flagDiscontinuousBlock = False
570
571 485 self.utctime = None
572
573 486 self.nCohInt = None
574
575 487 self.nIncohInt = None
576
577 488 self.blocksize = None
578
579 489 self.nFFTPoints = None
580
581 490 self.wavelength = None
582
583 491 self.flagDecodeData = False # asumo q la data no esta decodificada
584
585 492 self.flagDeflipData = False # asumo q la data no esta sin flip
586
587 493 self.flagShiftFFT = False
588
589 494 self.ippFactor = 1
590
591 495 #self.noise = None
592
593 496 self.beacon_heiIndexList = []
594
595 497 self.noise_estimation = None
596 498
597 499 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
@@ -652,8 +554,11 class Spectra(JROData):
652 554
653 555 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
654 556 velrange = deltav * (numpy.arange(self.nFFTPoints +
655 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
557 extrapoints) - self.nFFTPoints / 2.)
656 558
559 if self.nmodes:
560 return velrange/self.nmodes
561 else:
657 562 return velrange
658 563
659 564 def getNPairs(self):
@@ -692,7 +597,8 class Spectra(JROData):
692 597
693 598 def getTimeInterval(self):
694 599
695 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
600 timeInterval = self.ippSeconds * self.nCohInt * \
601 self.nIncohInt * self.nProfiles * self.ippFactor
696 602
697 603 return timeInterval
698 604
@@ -755,19 +661,12 class Spectra(JROData):
755 661 class SpectraHeis(Spectra):
756 662
757 663 data_spc = None
758
759 664 data_cspc = None
760
761 665 data_dc = None
762
763 666 nFFTPoints = None
764
765 667 # nPairs = None
766
767 668 pairsList = None
768
769 669 nCohInt = None
770
771 670 nIncohInt = None
772 671
773 672 def __init__(self):
@@ -830,36 +729,21 class SpectraHeis(Spectra):
830 729 class Fits(JROData):
831 730
832 731 heightList = None
833
834 732 channelList = None
835
836 733 flagNoData = True
837
838 734 flagDiscontinuousBlock = False
839
840 735 useLocalTime = False
841
842 736 utctime = None
843
844 737 timeZone = None
845
846 738 # ippSeconds = None
847
848 739 # timeInterval = None
849
850 740 nCohInt = None
851
852 741 nIncohInt = None
853
854 742 noise = None
855
856 743 windowOfFilter = 1
857
858 744 # Speed of ligth
859 745 C = 3e8
860
861 746 frequency = 49.92e6
862
863 747 realtime = False
864 748
865 749 def __init__(self):
@@ -978,33 +862,19 class Fits(JROData):
978 862 class Correlation(JROData):
979 863
980 864 noise = None
981
982 865 SNR = None
983
984 866 #--------------------------------------------------
985
986 867 mode = None
987
988 868 split = False
989
990 869 data_cf = None
991
992 870 lags = None
993
994 871 lagRange = None
995
996 872 pairsList = None
997
998 873 normFactor = None
999
1000 874 #--------------------------------------------------
1001
1002 875 # calculateVelocity = None
1003
1004 876 nLags = None
1005
1006 877 nPairs = None
1007
1008 878 nAvg = None
1009 879
1010 880 def __init__(self):
@@ -1068,7 +938,8 class Correlation(JROData):
1068 938 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
1069 939
1070 940 if ind_vel[0] < 0:
1071 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
941 ind_vel[list(range(0, 1))] = ind_vel[list(
942 range(0, 1))] + self.num_prof
1072 943
1073 944 if mode == 1:
1074 945 jspectra[:, freq_dc, :] = (
@@ -1154,55 +1025,30 class Correlation(JROData):
1154 1025 class Parameters(Spectra):
1155 1026
1156 1027 experimentInfo = None # Information about the experiment
1157
1158 1028 # Information from previous data
1159
1160 1029 inputUnit = None # Type of data to be processed
1161
1162 1030 operation = None # Type of operation to parametrize
1163
1164 1031 # normFactor = None #Normalization Factor
1165
1166 1032 groupList = None # List of Pairs, Groups, etc
1167
1168 1033 # Parameters
1169
1170 1034 data_param = None # Parameters obtained
1171
1172 1035 data_pre = None # Data Pre Parametrization
1173
1174 1036 data_SNR = None # Signal to Noise Ratio
1175
1176 1037 # heightRange = None #Heights
1177
1178 1038 abscissaList = None # Abscissa, can be velocities, lags or time
1179
1180 1039 # noise = None #Noise Potency
1181
1182 1040 utctimeInit = None # Initial UTC time
1183
1184 1041 paramInterval = None # Time interval to calculate Parameters in seconds
1185
1186 1042 useLocalTime = True
1187
1188 1043 # Fitting
1189
1190 1044 data_error = None # Error of the estimation
1191
1192 1045 constants = None
1193
1194 1046 library = None
1195
1196 1047 # Output signal
1197
1198 1048 outputInterval = None # Time interval to calculate output signal in seconds
1199
1200 1049 data_output = None # Out signal
1201
1202 1050 nAvg = None
1203
1204 1051 noise_estimation = None
1205
1206 1052 GauSPC = None # Fit gaussian SPC
1207 1053
1208 1054 def __init__(self):
@@ -1249,3 +1095,251 class Parameters(Spectra):
1249 1095
1250 1096 timeInterval = property(getTimeInterval)
1251 1097 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1098
1099
1100 class PlotterData(object):
1101 '''
1102 Object to hold data to be plotted
1103 '''
1104
1105 MAXNUMX = 100
1106 MAXNUMY = 100
1107
1108 def __init__(self, code, throttle_value, exp_code, buffering=True):
1109
1110 self.throttle = throttle_value
1111 self.exp_code = exp_code
1112 self.buffering = buffering
1113 self.ready = False
1114 self.localtime = False
1115 self.data = {}
1116 self.meta = {}
1117 self.__times = []
1118 self.__heights = []
1119
1120 if 'snr' in code:
1121 self.plottypes = ['snr']
1122 elif code == 'spc':
1123 self.plottypes = ['spc', 'noise', 'rti']
1124 elif code == 'rti':
1125 self.plottypes = ['noise', 'rti']
1126 else:
1127 self.plottypes = [code]
1128
1129 for plot in self.plottypes:
1130 self.data[plot] = {}
1131
1132 def __str__(self):
1133 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1134 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1135
1136 def __len__(self):
1137 return len(self.__times)
1138
1139 def __getitem__(self, key):
1140
1141 if key not in self.data:
1142 raise KeyError(log.error('Missing key: {}'.format(key)))
1143 if 'spc' in key or not self.buffering:
1144 ret = self.data[key]
1145 else:
1146 ret = numpy.array([self.data[key][x] for x in self.times])
1147 if ret.ndim > 1:
1148 ret = numpy.swapaxes(ret, 0, 1)
1149 return ret
1150
1151 def __contains__(self, key):
1152 return key in self.data
1153
1154 def setup(self):
1155 '''
1156 Configure object
1157 '''
1158
1159 self.type = ''
1160 self.ready = False
1161 self.data = {}
1162 self.__times = []
1163 self.__heights = []
1164 self.__all_heights = set()
1165 for plot in self.plottypes:
1166 if 'snr' in plot:
1167 plot = 'snr'
1168 self.data[plot] = {}
1169
1170 if 'spc' in self.data or 'rti' in self.data:
1171 self.data['noise'] = {}
1172 if 'noise' not in self.plottypes:
1173 self.plottypes.append('noise')
1174
1175 def shape(self, key):
1176 '''
1177 Get the shape of the one-element data for the given key
1178 '''
1179
1180 if len(self.data[key]):
1181 if 'spc' in key or not self.buffering:
1182 return self.data[key].shape
1183 return self.data[key][self.__times[0]].shape
1184 return (0,)
1185
1186 def update(self, dataOut, tm):
1187 '''
1188 Update data object with new dataOut
1189 '''
1190
1191 if tm in self.__times:
1192 return
1193
1194 self.type = dataOut.type
1195 self.parameters = getattr(dataOut, 'parameters', [])
1196 if hasattr(dataOut, 'pairsList'):
1197 self.pairs = dataOut.pairsList
1198 if hasattr(dataOut, 'meta'):
1199 self.meta = dataOut.meta
1200 self.channels = dataOut.channelList
1201 self.interval = dataOut.getTimeInterval()
1202 self.localtime = dataOut.useLocalTime
1203 if 'spc' in self.plottypes or 'cspc' in self.plottypes:
1204 self.xrange = (dataOut.getFreqRange(1)/1000.,
1205 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1206 self.__heights.append(dataOut.heightList)
1207 self.__all_heights.update(dataOut.heightList)
1208 self.__times.append(tm)
1209
1210 for plot in self.plottypes:
1211 if plot == 'spc':
1212 z = dataOut.data_spc/dataOut.normFactor
1213 buffer = 10*numpy.log10(z)
1214 if plot == 'cspc':
1215 buffer = dataOut.data_cspc
1216 if plot == 'noise':
1217 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1218 if plot == 'rti':
1219 buffer = dataOut.getPower()
1220 if plot == 'snr_db':
1221 buffer = dataOut.data_SNR
1222 if plot == 'snr':
1223 buffer = 10*numpy.log10(dataOut.data_SNR)
1224 if plot == 'dop':
1225 buffer = 10*numpy.log10(dataOut.data_DOP)
1226 if plot == 'mean':
1227 buffer = dataOut.data_MEAN
1228 if plot == 'std':
1229 buffer = dataOut.data_STD
1230 if plot == 'coh':
1231 buffer = dataOut.getCoherence()
1232 if plot == 'phase':
1233 buffer = dataOut.getCoherence(phase=True)
1234 if plot == 'output':
1235 buffer = dataOut.data_output
1236 if plot == 'param':
1237 buffer = dataOut.data_param
1238
1239 if 'spc' in plot:
1240 self.data[plot] = buffer
1241 else:
1242 if self.buffering:
1243 self.data[plot][tm] = buffer
1244 else:
1245 self.data[plot] = buffer
1246
1247 def normalize_heights(self):
1248 '''
1249 Ensure same-dimension of the data for different heighList
1250 '''
1251
1252 H = numpy.array(list(self.__all_heights))
1253 H.sort()
1254 for key in self.data:
1255 shape = self.shape(key)[:-1] + H.shape
1256 for tm, obj in list(self.data[key].items()):
1257 h = self.__heights[self.__times.index(tm)]
1258 if H.size == h.size:
1259 continue
1260 index = numpy.where(numpy.in1d(H, h))[0]
1261 dummy = numpy.zeros(shape) + numpy.nan
1262 if len(shape) == 2:
1263 dummy[:, index] = obj
1264 else:
1265 dummy[index] = obj
1266 self.data[key][tm] = dummy
1267
1268 self.__heights = [H for tm in self.__times]
1269
1270 def jsonify(self, decimate=False):
1271 '''
1272 Convert data to json
1273 '''
1274
1275 data = {}
1276 tm = self.times[-1]
1277 dy = int(self.heights.size/self.MAXNUMY) + 1
1278 for key in self.data:
1279 if key in ('spc', 'cspc') or not self.buffering:
1280 dx = int(self.data[key].shape[1]/self.MAXNUMX) + 1
1281 data[key] = self.roundFloats(
1282 self.data[key][::, ::dx, ::dy].tolist())
1283 else:
1284 data[key] = self.roundFloats(self.data[key][tm].tolist())
1285
1286 ret = {'data': data}
1287 ret['exp_code'] = self.exp_code
1288 ret['time'] = float(tm)
1289 ret['interval'] = float(self.interval)
1290 ret['localtime'] = self.localtime
1291 ret['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1292 if 'spc' in self.data or 'cspc' in self.data:
1293 ret['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1294 else:
1295 ret['xrange'] = []
1296 if hasattr(self, 'pairs'):
1297 ret['pairs'] = [(int(p[0]), int(p[1])) for p in self.pairs]
1298 else:
1299 ret['pairs'] = []
1300
1301 for key, value in list(self.meta.items()):
1302 ret[key] = value
1303
1304 return json.dumps(ret)
1305
1306 @property
1307 def times(self):
1308 '''
1309 Return the list of times of the current data
1310 '''
1311
1312 ret = numpy.array(self.__times)
1313 ret.sort()
1314 return ret
1315
1316 @property
1317 def min_time(self):
1318 '''
1319 Return the minimun time value
1320 '''
1321
1322 return self.times[0]
1323
1324 @property
1325 def max_time(self):
1326 '''
1327 Return the maximun time value
1328 '''
1329
1330 return self.times[-1]
1331
1332 @property
1333 def heights(self):
1334 '''
1335 Return the list of heights of the current data
1336 '''
1337
1338 return numpy.array(self.__heights[-1])
1339
1340 @staticmethod
1341 def roundFloats(obj):
1342 if isinstance(obj, list):
1343 return list(map(PlotterData.roundFloats, obj))
1344 elif isinstance(obj, float):
1345 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -4,4 +4,3 from .jroplot_heispectra import *
4 4 from .jroplot_correlation import *
5 5 from .jroplot_parameters import *
6 6 from .jroplot_data import *
7 from .jroplotter import *
Show file before | Show file after | Hide comments
@@ -5,7 +5,7 import copy
5 5 from schainpy.model import *
6 6 from .figure import Figure, isRealtime
7 7
8 class CorrelationPlot(Figure):
8 class CorrelationPlot_(Figure):
9 9 isConfig = None
10 10 __nsubplots = None
11 11
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (639 lines changed) Show them Hide them
@@ -1,41 +1,32
1 '''
2 New Plots Operations
3
4 @author: juan.espinoza@jro.igp.gob.pe
5 '''
6
1 7
2 import os
3 8 import time
4 import glob
5 9 import datetime
6 from multiprocessing import Process
7
8 import zmq
9 10 import numpy
10 import matplotlib
11 import matplotlib.pyplot as plt
12 from matplotlib.patches import Polygon
13 from mpl_toolkits.axes_grid1 import make_axes_locatable
14 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
15 11
16 from schainpy.model.proc.jroproc_base import Operation
12 from schainpy.model.graphics.jroplot_base import Plot, plt
17 13 from schainpy.utils import log
18 14
19 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
20 blu_values = matplotlib.pyplot.get_cmap(
21 'seismic_r', 20)(numpy.arange(20))[10:15]
22 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
23 'jro', numpy.vstack((blu_values, jet_values)))
24 matplotlib.pyplot.register_cmap(cmap=ncmap)
25
26 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis', 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
27
28 15 EARTH_RADIUS = 6.3710e3
29 16
17
30 18 def ll2xy(lat1, lon1, lat2, lon2):
31 19
32 20 p = 0.017453292519943295
33 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
21 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
22 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
34 23 r = 12742 * numpy.arcsin(numpy.sqrt(a))
35 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)*numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
24 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
25 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
36 26 theta = -theta + numpy.pi/2
37 27 return r*numpy.cos(theta), r*numpy.sin(theta)
38 28
29
39 30 def km2deg(km):
40 31 '''
41 32 Convert distance in km to degrees
@@ -43,536 +34,8 def km2deg(km):
43 34
44 35 return numpy.rad2deg(km/EARTH_RADIUS)
45 36
46 def figpause(interval):
47 backend = plt.rcParams['backend']
48 if backend in matplotlib.rcsetup.interactive_bk:
49 figManager = matplotlib._pylab_helpers.Gcf.get_active()
50 if figManager is not None:
51 canvas = figManager.canvas
52 if canvas.figure.stale:
53 canvas.draw()
54 try:
55 canvas.start_event_loop(interval)
56 except:
57 pass
58 return
59
60 def popup(message):
61 '''
62 '''
63
64 fig = plt.figure(figsize=(12, 8), facecolor='r')
65 text = '\n'.join([s.strip() for s in message.split(':')])
66 fig.text(0.01, 0.5, text, ha='left', va='center', size='20', weight='heavy', color='w')
67 fig.show()
68 figpause(1000)
69
70
71 class PlotData(Operation, Process):
72 '''
73 Base class for Schain plotting operations
74 '''
75
76 CODE = 'Figure'
77 colormap = 'jro'
78 bgcolor = 'white'
79 CONFLATE = False
80 __missing = 1E30
81
82 __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax',
83 'zlimits', 'xlabel', 'ylabel', 'xaxis','cb_label', 'title',
84 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure',
85 'showprofile', 'decimation', 'ftp']
86
87 def __init__(self, **kwargs):
88
89 Operation.__init__(self, plot=True, **kwargs)
90 Process.__init__(self)
91
92 self.kwargs['code'] = self.CODE
93 self.mp = False
94 self.data = None
95 self.isConfig = False
96 self.figures = []
97 self.axes = []
98 self.cb_axes = []
99 self.localtime = kwargs.pop('localtime', True)
100 self.show = kwargs.get('show', True)
101 self.save = kwargs.get('save', False)
102 self.ftp = kwargs.get('ftp', False)
103 self.colormap = kwargs.get('colormap', self.colormap)
104 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
105 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
106 self.colormaps = kwargs.get('colormaps', None)
107 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
108 self.showprofile = kwargs.get('showprofile', False)
109 self.title = kwargs.get('wintitle', self.CODE.upper())
110 self.cb_label = kwargs.get('cb_label', None)
111 self.cb_labels = kwargs.get('cb_labels', None)
112 self.labels = kwargs.get('labels', None)
113 self.xaxis = kwargs.get('xaxis', 'frequency')
114 self.zmin = kwargs.get('zmin', None)
115 self.zmax = kwargs.get('zmax', None)
116 self.zlimits = kwargs.get('zlimits', None)
117 self.xmin = kwargs.get('xmin', None)
118 self.xmax = kwargs.get('xmax', None)
119 self.xrange = kwargs.get('xrange', 24)
120 self.xscale = kwargs.get('xscale', None)
121 self.ymin = kwargs.get('ymin', None)
122 self.ymax = kwargs.get('ymax', None)
123 self.yscale = kwargs.get('yscale', None)
124 self.xlabel = kwargs.get('xlabel', None)
125 self.decimation = kwargs.get('decimation', None)
126 self.showSNR = kwargs.get('showSNR', False)
127 self.oneFigure = kwargs.get('oneFigure', True)
128 self.width = kwargs.get('width', None)
129 self.height = kwargs.get('height', None)
130 self.colorbar = kwargs.get('colorbar', True)
131 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
132 self.channels = kwargs.get('channels', None)
133 self.titles = kwargs.get('titles', [])
134 self.polar = False
135 self.grid = kwargs.get('grid', False)
136
137 def __fmtTime(self, x, pos):
138 '''
139 '''
140
141 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
142
143 def __setup(self):
144 '''
145 Common setup for all figures, here figures and axes are created
146 '''
147
148 if self.CODE not in self.data:
149 raise ValueError(log.error('Missing data for {}'.format(self.CODE),
150 self.name))
151
152 self.setup()
153
154 self.time_label = 'LT' if self.localtime else 'UTC'
155 if self.data.localtime:
156 self.getDateTime = datetime.datetime.fromtimestamp
157 else:
158 self.getDateTime = datetime.datetime.utcfromtimestamp
159
160 if self.width is None:
161 self.width = 8
162
163 self.figures = []
164 self.axes = []
165 self.cb_axes = []
166 self.pf_axes = []
167 self.cmaps = []
168
169 size = '15%' if self.ncols == 1 else '30%'
170 pad = '4%' if self.ncols == 1 else '8%'
171
172 if self.oneFigure:
173 if self.height is None:
174 self.height = 1.4 * self.nrows + 1
175 fig = plt.figure(figsize=(self.width, self.height),
176 edgecolor='k',
177 facecolor='w')
178 self.figures.append(fig)
179 for n in range(self.nplots):
180 ax = fig.add_subplot(self.nrows, self.ncols,
181 n + 1, polar=self.polar)
182 ax.tick_params(labelsize=8)
183 ax.firsttime = True
184 ax.index = 0
185 ax.press = None
186 self.axes.append(ax)
187 if self.showprofile:
188 cax = self.__add_axes(ax, size=size, pad=pad)
189 cax.tick_params(labelsize=8)
190 self.pf_axes.append(cax)
191 else:
192 if self.height is None:
193 self.height = 3
194 for n in range(self.nplots):
195 fig = plt.figure(figsize=(self.width, self.height),
196 edgecolor='k',
197 facecolor='w')
198 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
199 ax.tick_params(labelsize=8)
200 ax.firsttime = True
201 ax.index = 0
202 ax.press = None
203 self.figures.append(fig)
204 self.axes.append(ax)
205 if self.showprofile:
206 cax = self.__add_axes(ax, size=size, pad=pad)
207 cax.tick_params(labelsize=8)
208 self.pf_axes.append(cax)
209
210 for n in range(self.nrows):
211 if self.colormaps is not None:
212 cmap = plt.get_cmap(self.colormaps[n])
213 else:
214 cmap = plt.get_cmap(self.colormap)
215 cmap.set_bad(self.bgcolor, 1.)
216 self.cmaps.append(cmap)
217
218 for fig in self.figures:
219 fig.canvas.mpl_connect('key_press_event', self.OnKeyPress)
220 fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll)
221 fig.canvas.mpl_connect('button_press_event', self.onBtnPress)
222 fig.canvas.mpl_connect('motion_notify_event', self.onMotion)
223 fig.canvas.mpl_connect('button_release_event', self.onBtnRelease)
224 if self.show:
225 fig.show()
226
227 def OnKeyPress(self, event):
228 '''
229 Event for pressing keys (up, down) change colormap
230 '''
231 ax = event.inaxes
232 if ax in self.axes:
233 if event.key == 'down':
234 ax.index += 1
235 elif event.key == 'up':
236 ax.index -= 1
237 if ax.index < 0:
238 ax.index = len(CMAPS) - 1
239 elif ax.index == len(CMAPS):
240 ax.index = 0
241 cmap = CMAPS[ax.index]
242 ax.cbar.set_cmap(cmap)
243 ax.cbar.draw_all()
244 ax.plt.set_cmap(cmap)
245 ax.cbar.patch.figure.canvas.draw()
246 self.colormap = cmap.name
247
248 def OnBtnScroll(self, event):
249 '''
250 Event for scrolling, scale figure
251 '''
252 cb_ax = event.inaxes
253 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
254 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
255 pt = ax.cbar.ax.bbox.get_points()[:, 1]
256 nrm = ax.cbar.norm
257 vmin, vmax, p0, p1, pS = (
258 nrm.vmin, nrm.vmax, pt[0], pt[1], event.y)
259 scale = 2 if event.step == 1 else 0.5
260 point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0)
261 ax.cbar.norm.vmin = point - scale * (point - vmin)
262 ax.cbar.norm.vmax = point - scale * (point - vmax)
263 ax.plt.set_norm(ax.cbar.norm)
264 ax.cbar.draw_all()
265 ax.cbar.patch.figure.canvas.draw()
266
267 def onBtnPress(self, event):
268 '''
269 Event for mouse button press
270 '''
271 cb_ax = event.inaxes
272 if cb_ax is None:
273 return
274
275 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
276 cb_ax.press = event.x, event.y
277 else:
278 cb_ax.press = None
279
280 def onMotion(self, event):
281 '''
282 Event for move inside colorbar
283 '''
284 cb_ax = event.inaxes
285 if cb_ax is None:
286 return
287 if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]:
288 return
289 if cb_ax.press is None:
290 return
291
292 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
293 xprev, yprev = cb_ax.press
294 dx = event.x - xprev
295 dy = event.y - yprev
296 cb_ax.press = event.x, event.y
297 scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin
298 perc = 0.03
299
300 if event.button == 1:
301 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
302 ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy)
303 elif event.button == 3:
304 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
305 ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy)
306
307 ax.cbar.draw_all()
308 ax.plt.set_norm(ax.cbar.norm)
309 ax.cbar.patch.figure.canvas.draw()
310
311 def onBtnRelease(self, event):
312 '''
313 Event for mouse button release
314 '''
315 cb_ax = event.inaxes
316 if cb_ax is not None:
317 cb_ax.press = None
318
319 def __add_axes(self, ax, size='30%', pad='8%'):
320 '''
321 Add new axes to the given figure
322 '''
323 divider = make_axes_locatable(ax)
324 nax = divider.new_horizontal(size=size, pad=pad)
325 ax.figure.add_axes(nax)
326 return nax
327
328 self.setup()
329
330 def setup(self):
331 '''
332 This method should be implemented in the child class, the following
333 attributes should be set:
334
335 self.nrows: number of rows
336 self.ncols: number of cols
337 self.nplots: number of plots (channels or pairs)
338 self.ylabel: label for Y axes
339 self.titles: list of axes title
340
341 '''
342 raise NotImplementedError
343
344 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
345 '''
346 Create a masked array for missing data
347 '''
348 if x_buffer.shape[0] < 2:
349 return x_buffer, y_buffer, z_buffer
350
351 deltas = x_buffer[1:] - x_buffer[0:-1]
352 x_median = numpy.median(deltas)
353 37
354 index = numpy.where(deltas > 5 * x_median)
355
356 if len(index[0]) != 0:
357 z_buffer[::, index[0], ::] = self.__missing
358 z_buffer = numpy.ma.masked_inside(z_buffer,
359 0.99 * self.__missing,
360 1.01 * self.__missing)
361
362 return x_buffer, y_buffer, z_buffer
363
364 def decimate(self):
365
366 # dx = int(len(self.x)/self.__MAXNUMX) + 1
367 dy = int(len(self.y) / self.decimation) + 1
368
369 # x = self.x[::dx]
370 x = self.x
371 y = self.y[::dy]
372 z = self.z[::, ::, ::dy]
373
374 return x, y, z
375
376 def format(self):
377 '''
378 Set min and max values, labels, ticks and titles
379 '''
380
381 if self.xmin is None:
382 xmin = self.min_time
383 else:
384 if self.xaxis is 'time':
385 dt = self.getDateTime(self.min_time)
386 xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) -
387 datetime.datetime(1970, 1, 1)).total_seconds()
388 if self.data.localtime:
389 xmin += time.timezone
390 else:
391 xmin = self.xmin
392
393 if self.xmax is None:
394 xmax = xmin + self.xrange * 60 * 60
395 else:
396 if self.xaxis is 'time':
397 dt = self.getDateTime(self.max_time)
398 xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) -
399 datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds()
400 if self.data.localtime:
401 xmax += time.timezone
402 else:
403 xmax = self.xmax
404
405 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
406 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
407
408 Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])
409 i = 1 if numpy.where(abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0]
410 ystep = Y[i] / 10.
411
412 if self.xaxis is not 'time':
413 X = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])/2.
414 i = 1 if numpy.where(abs(xmax-xmin) <= X)[0][0] < 0 else numpy.where(abs(xmax-xmin) <= X)[0][0]
415 xstep = X[i] / 10.
416
417 for n, ax in enumerate(self.axes):
418 if ax.firsttime:
419 ax.set_facecolor(self.bgcolor)
420 ax.yaxis.set_major_locator(MultipleLocator(ystep))
421 if self.xscale:
422 ax.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{0:g}'.format(x*self.xscale)))
423 if self.xscale:
424 ax.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{0:g}'.format(x*self.yscale)))
425 if self.xaxis is 'time':
426 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
427 ax.xaxis.set_major_locator(LinearLocator(9))
428 else:
429 ax.xaxis.set_major_locator(MultipleLocator(xstep))
430 if self.xlabel is not None:
431 ax.set_xlabel(self.xlabel)
432 ax.set_ylabel(self.ylabel)
433 ax.firsttime = False
434 if self.showprofile:
435 self.pf_axes[n].set_ylim(ymin, ymax)
436 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
437 self.pf_axes[n].set_xlabel('dB')
438 self.pf_axes[n].grid(b=True, axis='x')
439 [tick.set_visible(False)
440 for tick in self.pf_axes[n].get_yticklabels()]
441 if self.colorbar:
442 ax.cbar = plt.colorbar(
443 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
444 ax.cbar.ax.tick_params(labelsize=8)
445 ax.cbar.ax.press = None
446 if self.cb_label:
447 ax.cbar.set_label(self.cb_label, size=8)
448 elif self.cb_labels:
449 ax.cbar.set_label(self.cb_labels[n], size=8)
450 else:
451 ax.cbar = None
452 if self.grid:
453 ax.grid(True)
454
455 if not self.polar:
456 ax.set_xlim(xmin, xmax)
457 ax.set_ylim(ymin, ymax)
458 ax.set_title('{} {} {}'.format(
459 self.titles[n],
460 self.getDateTime(self.max_time).strftime('%Y-%m-%dT%H:%M:%S'),
461 self.time_label),
462 size=8)
463 else:
464 ax.set_title('{}'.format(self.titles[n]), size=8)
465 ax.set_ylim(0, 90)
466 ax.set_yticks(numpy.arange(0, 90, 20))
467 ax.yaxis.labelpad = 40
468
469 def __plot(self):
470 '''
471 '''
472 log.log('Plotting', self.name)
473
474 try:
475 self.plot()
476 self.format()
477 except Exception as e:
478 log.warning('{} Plot could not be updated... check data'.format(self.CODE), self.name)
479 log.error(str(e), '')
480 return
481
482 for n, fig in enumerate(self.figures):
483 if self.nrows == 0 or self.nplots == 0:
484 log.warning('No data', self.name)
485 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
486 fig.canvas.manager.set_window_title(self.CODE)
487 continue
488
489 fig.tight_layout()
490 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
491 self.getDateTime(self.max_time).strftime('%Y/%m/%d')))
492 fig.canvas.draw()
493
494 if self.save and (self.data.ended or not self.data.buffering):
495
496 if self.save_labels:
497 labels = self.save_labels
498 else:
499 labels = list(range(self.nrows))
500
501 if self.oneFigure:
502 label = ''
503 else:
504 label = '-{}'.format(labels[n])
505 figname = os.path.join(
506 self.save,
507 self.CODE,
508 '{}{}_{}.png'.format(
509 self.CODE,
510 label,
511 self.getDateTime(self.saveTime).strftime(
512 '%Y%m%d_%H%M%S'),
513 )
514 )
515 log.log('Saving figure: {}'.format(figname), self.name)
516 if not os.path.isdir(os.path.dirname(figname)):
517 os.makedirs(os.path.dirname(figname))
518 fig.savefig(figname)
519
520 def plot(self):
521 '''
522 '''
523 raise NotImplementedError
524
525 def run(self):
526
527 log.log('Starting', self.name)
528
529 context = zmq.Context()
530 receiver = context.socket(zmq.SUB)
531 receiver.setsockopt(zmq.SUBSCRIBE, '')
532 receiver.setsockopt(zmq.CONFLATE, self.CONFLATE)
533
534 if 'server' in self.kwargs['parent']:
535 receiver.connect(
536 'ipc:///tmp/{}.plots'.format(self.kwargs['parent']['server']))
537 else:
538 receiver.connect("ipc:///tmp/zmq.plots")
539
540 while True:
541 try:
542 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)
543 if self.data.localtime and self.localtime:
544 self.times = self.data.times
545 elif self.data.localtime and not self.localtime:
546 self.times = self.data.times + time.timezone
547 elif not self.data.localtime and self.localtime:
548 self.times = self.data.times - time.timezone
549 else:
550 self.times = self.data.times
551
552 self.min_time = self.times[0]
553 self.max_time = self.times[-1]
554
555 if self.isConfig is False:
556 self.__setup()
557 self.isConfig = True
558
559 self.__plot()
560
561 except zmq.Again as e:
562 if self.data and self.data.ended:
563 break
564 log.log('Waiting for data...')
565 if self.data:
566 figpause(self.data.throttle)
567 else:
568 time.sleep(2)
569
570 def close(self):
571 if self.data:
572 self.__plot()
573
574
575 class PlotSpectraData(PlotData):
38 class SpectraPlot(Plot):
576 39 '''
577 40 Plot for Spectra data
578 41 '''
@@ -644,10 +107,9 class PlotSpectraData(PlotData):
644 107 ax.plt_mean.set_data(mean, y)
645 108
646 109 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
647 self.saveTime = self.max_time
648 110
649 111
650 class PlotCrossSpectraData(PlotData):
112 class CrossSpectraPlot(Plot):
651 113
652 114 CODE = 'cspc'
653 115 zmin_coh = None
@@ -741,10 +203,8 class PlotCrossSpectraData(PlotData):
741 203 ax.plt.set_array(phase.T.ravel())
742 204 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
743 205
744 self.saveTime = self.max_time
745
746 206
747 class PlotSpectraMeanData(PlotSpectraData):
207 class SpectraMeanPlot(SpectraPlot):
748 208 '''
749 209 Plot for Spectra and Mean
750 210 '''
@@ -752,7 +212,7 class PlotSpectraMeanData(PlotSpectraData):
752 212 colormap = 'jro'
753 213
754 214
755 class PlotRTIData(PlotData):
215 class RTIPlot(Plot):
756 216 '''
757 217 Plot for RTI data
758 218 '''
@@ -771,7 +231,7 class PlotRTIData(PlotData):
771 231 self.CODE.upper(), x) for x in range(self.nrows)]
772 232
773 233 def plot(self):
774 self.x = self.times
234 self.x = self.data.times
775 235 self.y = self.data.heights
776 236 self.z = self.data[self.CODE]
777 237 self.z = numpy.ma.masked_invalid(self.z)
@@ -807,10 +267,8 class PlotRTIData(PlotData):
807 267 ax.plot_noise.set_data(numpy.repeat(
808 268 self.data['noise'][n][-1], len(self.y)), self.y)
809 269
810 self.saveTime = self.min_time
811 270
812
813 class PlotCOHData(PlotRTIData):
271 class CoherencePlot(RTIPlot):
814 272 '''
815 273 Plot for Coherence data
816 274 '''
@@ -833,7 +291,7 class PlotCOHData(PlotRTIData):
833 291 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
834 292
835 293
836 class PlotPHASEData(PlotCOHData):
294 class PhasePlot(CoherencePlot):
837 295 '''
838 296 Plot for Phase map data
839 297 '''
@@ -842,7 +300,7 class PlotPHASEData(PlotCOHData):
842 300 colormap = 'seismic'
843 301
844 302
845 class PlotNoiseData(PlotData):
303 class NoisePlot(Plot):
846 304 '''
847 305 Plot for noise
848 306 '''
@@ -860,8 +318,8 class PlotNoiseData(PlotData):
860 318
861 319 def plot(self):
862 320
863 x = self.times
864 xmin = self.min_time
321 x = self.data.times
322 xmin = self.data.min_time
865 323 xmax = xmin + self.xrange * 60 * 60
866 324 Y = self.data[self.CODE]
867 325
@@ -877,10 +335,9 class PlotNoiseData(PlotData):
877 335
878 336 self.ymin = numpy.nanmin(Y) - 5
879 337 self.ymax = numpy.nanmax(Y) + 5
880 self.saveTime = self.min_time
881 338
882 339
883 class PlotSNRData(PlotRTIData):
340 class SnrPlot(RTIPlot):
884 341 '''
885 342 Plot for SNR Data
886 343 '''
@@ -889,7 +346,7 class PlotSNRData(PlotRTIData):
889 346 colormap = 'jet'
890 347
891 348
892 class PlotDOPData(PlotRTIData):
349 class DopplerPlot(RTIPlot):
893 350 '''
894 351 Plot for DOPPLER Data
895 352 '''
@@ -898,7 +355,7 class PlotDOPData(PlotRTIData):
898 355 colormap = 'jet'
899 356
900 357
901 class PlotSkyMapData(PlotData):
358 class SkyMapPlot(Plot):
902 359 '''
903 360 Plot for meteors detection data
904 361 '''
@@ -938,16 +395,15 class PlotSkyMapData(PlotData):
938 395 else:
939 396 ax.plot.set_data(x, y)
940 397
941 dt1 = self.getDateTime(self.min_time).strftime('%y/%m/%d %H:%M:%S')
942 dt2 = self.getDateTime(self.max_time).strftime('%y/%m/%d %H:%M:%S')
398 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
399 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
943 400 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
944 401 dt2,
945 402 len(x))
946 403 self.titles[0] = title
947 self.saveTime = self.max_time
948 404
949 405
950 class PlotParamData(PlotRTIData):
406 class ParametersPlot(RTIPlot):
951 407 '''
952 408 Plot for data_param object
953 409 '''
@@ -973,7 +429,7 class PlotParamData(PlotRTIData):
973 429
974 430 def plot(self):
975 431 self.data.normalize_heights()
976 self.x = self.times
432 self.x = self.data.times
977 433 self.y = self.data.heights
978 434 if self.showSNR:
979 435 self.z = numpy.concatenate(
@@ -1015,10 +471,8 class PlotParamData(PlotRTIData):
1015 471 cmap=self.cmaps[n]
1016 472 )
1017 473
1018 self.saveTime = self.min_time
1019
1020 474
1021 class PlotOutputData(PlotParamData):
475 class OutputPlot(ParametersPlot):
1022 476 '''
1023 477 Plot data_output object
1024 478 '''
@@ -1027,9 +481,9 class PlotOutputData(PlotParamData):
1027 481 colormap = 'seismic'
1028 482
1029 483
1030 class PlotPolarMapData(PlotData):
484 class PolarMapPlot(Plot):
1031 485 '''
1032 Plot for meteors detection data
486 Plot for weather radar
1033 487 '''
1034 488
1035 489 CODE = 'param'
@@ -1058,8 +512,10 class PlotPolarMapData(PlotData):
1058 512 self.cb_labels = self.data.meta['units']
1059 513 self.lat = self.data.meta['latitude']
1060 514 self.lon = self.data.meta['longitude']
1061 self.xmin, self.xmax = float(km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
1062 self.ymin, self.ymax = float(km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
515 self.xmin, self.xmax = float(
516 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
517 self.ymin, self.ymax = float(
518 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
1063 519 # self.polar = True
1064 520
1065 521 def plot(self):
@@ -1067,11 +523,13 class PlotPolarMapData(PlotData):
1067 523 for n, ax in enumerate(self.axes):
1068 524 data = self.data['param'][self.channels[n]]
1069 525
1070 zeniths = numpy.linspace(0, self.data.meta['max_range'], data.shape[1])
526 zeniths = numpy.linspace(
527 0, self.data.meta['max_range'], data.shape[1])
1071 528 if self.mode == 'E':
1072 529 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
1073 530 r, theta = numpy.meshgrid(zeniths, azimuths)
1074 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
531 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
532 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
1075 533 x = km2deg(x) + self.lon
1076 534 y = km2deg(y) + self.lat
1077 535 else:
@@ -1108,7 +566,8 class PlotPolarMapData(PlotData):
1108 566 lat = float(lat)
1109 567 lon = float(lon)
1110 568 # ax.plot(lon, lat, '.b', ms=2)
1111 ax.text(lon, lat, label.decode('utf8'), ha='center', va='bottom', size='8', color='black')
569 ax.text(lon, lat, label.decode('utf8'), ha='center',
570 va='bottom', size='8', color='black')
1112 571
1113 572 # plot limites
1114 573 limites =[]
@@ -1122,7 +581,8 class PlotPolarMapData(PlotData):
1122 581 values = line.strip().split(',')
1123 582 tmp.append((float(values[0]), float(values[1])))
1124 583 for points in limites:
1125 ax.add_patch(Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
584 ax.add_patch(
585 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
1126 586
1127 587 # plot Cuencas
1128 588 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
@@ -1133,7 +593,8 class PlotPolarMapData(PlotData):
1133 593
1134 594 # plot grid
1135 595 for r in (15, 30, 45, 60):
1136 ax.add_artist(plt.Circle((self.lon, self.lat), km2deg(r), color='0.6', fill=False, lw=0.2))
596 ax.add_artist(plt.Circle((self.lon, self.lat),
597 km2deg(r), color='0.6', fill=False, lw=0.2))
1137 598 ax.text(
1138 599 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
1139 600 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
@@ -1148,7 +609,5 class PlotPolarMapData(PlotData):
1148 609 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
1149 610
1150 611 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
1151 self.titles = ['{} {}'.format(self.data.parameters[x], title) for x in self.channels]
1152 self.saveTime = self.max_time
1153
1154 No newline at end of file
612 self.titles = ['{} {}'.format(
613 self.data.parameters[x], title) for x in self.channels]
Show file before | Show file after | Hide comments
@@ -10,7 +10,7 import numpy
10 10 from .figure import Figure, isRealtime
11 11 from .plotting_codes import *
12 12
13 class SpectraHeisScope(Figure):
13 class SpectraHeisScope_(Figure):
14 14
15 15
16 16 isConfig = None
@@ -173,7 +173,7 class SpectraHeisScope(Figure):
173 173 wr_period=wr_period,
174 174 thisDatetime=thisDatetime)
175 175
176 class RTIfromSpectraHeis(Figure):
176 class RTIfromSpectraHeis_(Figure):
177 177
178 178 isConfig = None
179 179 __nsubplots = None
Show file before | Show file after | Hide comments
@@ -7,14 +7,190 from .plotting_codes import *
7 7 from schainpy.model.proc.jroproc_base import MPDecorator
8 8 from schainpy.utils import log
9 9
10 class FitGauPlot(Figure):
10 class ParamLine_(Figure):
11
12 isConfig = None
13
14 def __init__(self):
15
16 self.isConfig = False
17 self.WIDTH = 300
18 self.HEIGHT = 200
19 self.counter_imagwr = 0
20
21 def getSubplots(self):
22
23 nrow = self.nplots
24 ncol = 3
25 return nrow, ncol
26
27 def setup(self, id, nplots, wintitle, show):
28
29 self.nplots = nplots
30
31 self.createFigure(id=id,
32 wintitle=wintitle,
33 show=show)
34
35 nrow,ncol = self.getSubplots()
36 colspan = 3
37 rowspan = 1
38
39 for i in range(nplots):
40 self.addAxes(nrow, ncol, i, 0, colspan, rowspan)
41
42 def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
43 yreal = y[channelIndexList,:].real
44 yimag = y[channelIndexList,:].imag
45
46 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
47 xlabel = "Range (Km)"
48 ylabel = "Intensity - IQ"
49
50 if not self.isConfig:
51 nplots = len(channelIndexList)
52
53 self.setup(id=id,
54 nplots=nplots,
55 wintitle='',
56 show=show)
57
58 if xmin == None: xmin = numpy.nanmin(x)
59 if xmax == None: xmax = numpy.nanmax(x)
60 if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag))
61 if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag))
62
63 self.isConfig = True
64
65 self.setWinTitle(title)
66
67 for i in range(len(self.axesList)):
68 title = "Channel %d" %(i)
69 axes = self.axesList[i]
70
71 axes.pline(x, yreal[i,:],
72 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
73 xlabel=xlabel, ylabel=ylabel, title=title)
74
75 axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2)
76
77 def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
78 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
79 yreal = y.real
80
81 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
82 xlabel = "Range (Km)"
83 ylabel = "Intensity"
84
85 if not self.isConfig:
86 nplots = len(channelIndexList)
87
88 self.setup(id=id,
89 nplots=nplots,
90 wintitle='',
91 show=show)
92
93 if xmin == None: xmin = numpy.nanmin(x)
94 if xmax == None: xmax = numpy.nanmax(x)
95 if ymin == None: ymin = numpy.nanmin(yreal)
96 if ymax == None: ymax = numpy.nanmax(yreal)
97
98 self.isConfig = True
99
100 self.setWinTitle(title)
101
102 for i in range(len(self.axesList)):
103 title = "Channel %d" %(i)
104 axes = self.axesList[i]
105 ychannel = yreal[i,:]
106 axes.pline(x, ychannel,
107 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
108 xlabel=xlabel, ylabel=ylabel, title=title)
109
110
111 def run(self, dataOut, id, wintitle="", channelList=None,
112 xmin=None, xmax=None, ymin=None, ymax=None, save=False,
113 figpath='./', figfile=None, show=True, wr_period=1,
114 ftp=False, server=None, folder=None, username=None, password=None):
115
116 """
117
118 Input:
119 dataOut :
120 id :
121 wintitle :
122 channelList :
123 xmin : None,
124 xmax : None,
125 ymin : None,
126 ymax : None,
127 """
128
129 if channelList == None:
130 channelIndexList = dataOut.channelIndexList
131 else:
132 channelIndexList = []
133 for channel in channelList:
134 if channel not in dataOut.channelList:
135 raise ValueError("Channel %d is not in dataOut.channelList" % channel)
136 channelIndexList.append(dataOut.channelList.index(channel))
137
138 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
139
140 y = dataOut.RR
141
142 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
143 xlabel = "Range (Km)"
144 ylabel = "Intensity"
145
146 if not self.isConfig:
147 nplots = len(channelIndexList)
148
149 self.setup(id=id,
150 nplots=nplots,
151 wintitle='',
152 show=show)
153
154 if xmin == None: xmin = numpy.nanmin(x)
155 if xmax == None: xmax = numpy.nanmax(x)
156 if ymin == None: ymin = numpy.nanmin(y)
157 if ymax == None: ymax = numpy.nanmax(y)
158
159 self.isConfig = True
160
161 self.setWinTitle(title)
162
163 for i in range(len(self.axesList)):
164 title = "Channel %d" %(i)
165 axes = self.axesList[i]
166 ychannel = y[i,:]
167 axes.pline(x, ychannel,
168 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
169 xlabel=xlabel, ylabel=ylabel, title=title)
170
171
172 self.draw()
173
174 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + "_" + str(dataOut.profileIndex)
175 figfile = self.getFilename(name = str_datetime)
176
177 self.save(figpath=figpath,
178 figfile=figfile,
179 save=save,
180 ftp=ftp,
181 wr_period=wr_period,
182 thisDatetime=thisDatetime)
183
184
185
186 class SpcParamPlot_(Figure):
11 187
12 188 isConfig = None
13 189 __nsubplots = None
14 190
15 191 WIDTHPROF = None
16 192 HEIGHTPROF = None
17 PREFIX = 'fitgau'
193 PREFIX = 'SpcParam'
18 194
19 195 def __init__(self, **kwargs):
20 196 Figure.__init__(self, **kwargs)
@@ -83,7 +259,7 class FitGauPlot(Figure):
83 259 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
84 260 server=None, folder=None, username=None, password=None,
85 261 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
86 xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 1):
262 xaxis="frequency", colormap='jet', normFactor=None , Selector = 0):
87 263
88 264 """
89 265
@@ -119,23 +295,22 class FitGauPlot(Figure):
119 295 # else:
120 296 # factor = normFactor
121 297 if xaxis == "frequency":
122 x = dataOut.spc_range[0]
298 x = dataOut.spcparam_range[0]
123 299 xlabel = "Frequency (kHz)"
124 300
125 301 elif xaxis == "time":
126 x = dataOut.spc_range[1]
302 x = dataOut.spcparam_range[1]
127 303 xlabel = "Time (ms)"
128 304
129 305 else:
130 x = dataOut.spc_range[2]
306 x = dataOut.spcparam_range[2]
131 307 xlabel = "Velocity (m/s)"
132 308
133 ylabel = "Range (Km)"
309 ylabel = "Range (km)"
134 310
135 311 y = dataOut.getHeiRange()
136 312
137 z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor
138 print('GausSPC', z[0,32,10:40])
313 z = dataOut.SPCparam[Selector] /1966080.0#/ dataOut.normFactor#GauSelector] #dataOut.data_spc/factor
139 314 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
140 315 zdB = 10*numpy.log10(z)
141 316
@@ -218,7 +393,7 class FitGauPlot(Figure):
218 393
219 394
220 395
221 class MomentsPlot(Figure):
396 class MomentsPlot_(Figure):
222 397
223 398 isConfig = None
224 399 __nsubplots = None
@@ -405,7 +580,7 class MomentsPlot(Figure):
405 580 thisDatetime=thisDatetime)
406 581
407 582
408 class SkyMapPlot(Figure):
583 class SkyMapPlot_(Figure):
409 584
410 585 __isConfig = None
411 586 __nsubplots = None
@@ -561,7 +736,7 class SkyMapPlot(Figure):
561 736
562 737
563 738
564 class WindProfilerPlot(Figure):
739 class WindProfilerPlot_(Figure):
565 740
566 741 __isConfig = None
567 742 __nsubplots = None
@@ -666,13 +841,14 class WindProfilerPlot(Figure):
666 841 #If there is a SNR function defined
667 842 if dataOut.data_SNR is not None:
668 843 nplots += 1
669 SNR = dataOut.data_SNR
670 SNRavg = numpy.average(SNR, axis=0)
844 SNR = dataOut.data_SNR[0]
845 SNRavg = SNR#numpy.average(SNR, axis=0)
671 846
672 847 SNRdB = 10*numpy.log10(SNR)
673 848 SNRavgdB = 10*numpy.log10(SNRavg)
674 849
675 if SNRthresh == None: SNRthresh = -5.0
850 if SNRthresh == None:
851 SNRthresh = -5.0
676 852 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
677 853
678 854 for i in range(nplotsw):
@@ -741,7 +917,6 class WindProfilerPlot(Figure):
741 917 axes = self.axesList[i*self.__nsubplots]
742 918
743 919 z1 = z[i,:].reshape((1,-1))*windFactor[i]
744 #z1=numpy.ma.masked_where(z1==0.,z1)
745 920
746 921 axes.pcolorbuffer(x, y, z1,
747 922 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
@@ -774,7 +949,7 class WindProfilerPlot(Figure):
774 949 update_figfile = True
775 950
776 951 @MPDecorator
777 class ParametersPlot(Figure):
952 class ParametersPlot_(Figure):
778 953
779 954 __isConfig = None
780 955 __nsubplots = None
@@ -792,8 +967,8 class ParametersPlot(Figure):
792 967 self.isConfig = False
793 968 self.__nsubplots = 1
794 969
795 self.WIDTH = 800
796 self.HEIGHT = 180
970 self.WIDTH = 300
971 self.HEIGHT = 550
797 972 self.WIDTHPROF = 120
798 973 self.HEIGHTPROF = 0
799 974 self.counter_imagwr = 0
@@ -905,7 +1080,7 class ParametersPlot(Figure):
905 1080 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
906 1081 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
907 1082 xlabel = ""
908 ylabel = "Range (Km)"
1083 ylabel = "Range (km)"
909 1084
910 1085 update_figfile = False
911 1086
@@ -949,9 +1124,28 class ParametersPlot(Figure):
949 1124
950 1125 self.setWinTitle(title)
951 1126
952 for i in range(self.nchan):
1127 # for i in range(self.nchan):
1128 # index = channelIndexList[i]
1129 # title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1130 # axes = self.axesList[i*self.plotFact]
1131 # z1 = z[i,:].reshape((1,-1))
1132 # axes.pcolorbuffer(x, y, z1,
1133 # xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1134 # xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1135 # ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
1136 #
1137 # if showSNR:
1138 # title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1139 # axes = self.axesList[i*self.plotFact + 1]
1140 # SNRdB1 = SNRdB[i,:].reshape((1,-1))
1141 # axes.pcolorbuffer(x, y, SNRdB1,
1142 # xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1143 # xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1144 # ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1145
1146 i=0
953 1147 index = channelIndexList[i]
954 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1148 title = "Factor de reflectividad Z [dBZ]"
955 1149 axes = self.axesList[i*self.plotFact]
956 1150 z1 = z[i,:].reshape((1,-1))
957 1151 axes.pcolorbuffer(x, y, z1,
@@ -968,6 +1162,44 class ParametersPlot(Figure):
968 1162 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
969 1163 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
970 1164
1165 i=1
1166 index = channelIndexList[i]
1167 title = "Velocidad vertical Doppler [m/s]"
1168 axes = self.axesList[i*self.plotFact]
1169 z1 = z[i,:].reshape((1,-1))
1170 axes.pcolorbuffer(x, y, z1,
1171 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-10, zmax=10,
1172 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1173 ticksize=9, cblabel='', cbsize="1%",colormap='seismic_r')
1174
1175 if showSNR:
1176 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1177 axes = self.axesList[i*self.plotFact + 1]
1178 SNRdB1 = SNRdB[i,:].reshape((1,-1))
1179 axes.pcolorbuffer(x, y, SNRdB1,
1180 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1181 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1182 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1183
1184 i=2
1185 index = channelIndexList[i]
1186 title = "Intensidad de lluvia [mm/h]"
1187 axes = self.axesList[i*self.plotFact]
1188 z1 = z[i,:].reshape((1,-1))
1189 axes.pcolorbuffer(x, y, z1,
1190 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=40,
1191 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1192 ticksize=9, cblabel='', cbsize="1%",colormap='ocean_r')
1193
1194 if showSNR:
1195 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1196 axes = self.axesList[i*self.plotFact + 1]
1197 SNRdB1 = SNRdB[i,:].reshape((1,-1))
1198 axes.pcolorbuffer(x, y, SNRdB1,
1199 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1200 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1201 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1202
971 1203
972 1204 self.draw()
973 1205
@@ -986,7 +1218,7 class ParametersPlot(Figure):
986 1218
987 1219 return dataOut
988 1220 @MPDecorator
989 class Parameters1Plot(Figure):
1221 class Parameters1Plot_(Figure):
990 1222
991 1223 __isConfig = None
992 1224 __nsubplots = None
@@ -1067,9 +1299,8 class Parameters1Plot(Figure):
1067 1299 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1068 1300 server=None, folder=None, username=None, password=None,
1069 1301 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1070 #print inspect.getargspec(self.run).args
1071 """
1072 1302
1303 """
1073 1304 Input:
1074 1305 dataOut :
1075 1306 id :
@@ -1237,7 +1468,7 class Parameters1Plot(Figure):
1237 1468 update_figfile=False)
1238 1469 return dataOut
1239 1470
1240 class SpectralFittingPlot(Figure):
1471 class SpectralFittingPlot_(Figure):
1241 1472
1242 1473 __isConfig = None
1243 1474 __nsubplots = None
@@ -1415,7 +1646,7 class SpectralFittingPlot(Figure):
1415 1646 thisDatetime=thisDatetime)
1416 1647
1417 1648
1418 class EWDriftsPlot(Figure):
1649 class EWDriftsPlot_(Figure):
1419 1650
1420 1651 __isConfig = None
1421 1652 __nsubplots = None
@@ -1621,7 +1852,7 class EWDriftsPlot(Figure):
1621 1852
1622 1853
1623 1854
1624 class PhasePlot(Figure):
1855 class PhasePlot_(Figure):
1625 1856
1626 1857 __isConfig = None
1627 1858 __nsubplots = None
@@ -1785,7 +2016,7 class PhasePlot(Figure):
1785 2016
1786 2017
1787 2018
1788 class NSMeteorDetection1Plot(Figure):
2019 class NSMeteorDetection1Plot_(Figure):
1789 2020
1790 2021 isConfig = None
1791 2022 __nsubplots = None
@@ -1969,7 +2200,7 class NSMeteorDetection1Plot(Figure):
1969 2200 thisDatetime=thisDatetime)
1970 2201
1971 2202
1972 class NSMeteorDetection2Plot(Figure):
2203 class NSMeteorDetection2Plot_(Figure):
1973 2204
1974 2205 isConfig = None
1975 2206 __nsubplots = None
Show file before | Show file after | Hide comments
@@ -14,7 +14,7 from schainpy.model.proc.jroproc_base import MPDecorator
14 14 from schainpy.utils import log
15 15
16 16 @MPDecorator
17 class SpectraPlot(Figure):
17 class SpectraPlot_(Figure):
18 18
19 19 isConfig = None
20 20 __nsubplots = None
@@ -43,6 +43,8 class SpectraPlot(Figure):
43 43 self.__xfilter_ena = False
44 44 self.__yfilter_ena = False
45 45
46 self.indice=1
47
46 48 def getSubplots(self):
47 49
48 50 ncol = int(numpy.sqrt(self.nplots)+0.9)
@@ -139,10 +141,9 class SpectraPlot(Figure):
139 141 x = dataOut.getVelRange(1)
140 142 xlabel = "Velocity (m/s)"
141 143
142 ylabel = "Range (Km)"
144 ylabel = "Range (km)"
143 145
144 146 y = dataOut.getHeiRange()
145
146 147 z = dataOut.data_spc/factor
147 148 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
148 149 zdB = 10*numpy.log10(z)
@@ -155,6 +156,7 class SpectraPlot(Figure):
155 156
156 157 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
157 158 title = wintitle + " Spectra"
159
158 160 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
159 161 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
160 162
@@ -224,9 +226,10 class SpectraPlot(Figure):
224 226 wr_period=wr_period,
225 227 thisDatetime=thisDatetime)
226 228
229
227 230 return dataOut
228 231 @MPDecorator
229 class CrossSpectraPlot(Figure):
232 class CrossSpectraPlot_(Figure):
230 233
231 234 isConfig = None
232 235 __nsubplots = None
@@ -253,6 +256,8 class CrossSpectraPlot(Figure):
253 256 self.SUB_EXP_CODE = None
254 257 self.PLOT_POS = None
255 258
259 self.indice=0
260
256 261 def getSubplots(self):
257 262
258 263 ncol = 4
@@ -397,6 +402,7 class CrossSpectraPlot(Figure):
397 402
398 403 self.setWinTitle(title)
399 404
405
400 406 for i in range(self.nplots):
401 407 pair = dataOut.pairsList[pairsIndexList[i]]
402 408
@@ -439,8 +445,6 class CrossSpectraPlot(Figure):
439 445 xlabel=xlabel, ylabel=ylabel, title=title,
440 446 ticksize=9, colormap=phase_cmap, cblabel='')
441 447
442
443
444 448 self.draw()
445 449
446 450 self.save(figpath=figpath,
@@ -453,7 +457,7 class CrossSpectraPlot(Figure):
453 457 return dataOut
454 458
455 459 @MPDecorator
456 class RTIPlot(Figure):
460 class RTIPlot_(Figure):
457 461
458 462 __isConfig = None
459 463 __nsubplots = None
@@ -470,7 +474,7 class RTIPlot(Figure):
470 474 self.__nsubplots = 1
471 475
472 476 self.WIDTH = 800
473 self.HEIGHT = 180
477 self.HEIGHT = 250
474 478 self.WIDTHPROF = 120
475 479 self.HEIGHTPROF = 0
476 480 self.counter_imagwr = 0
@@ -667,7 +671,7 class RTIPlot(Figure):
667 671 return dataOut
668 672
669 673 @MPDecorator
670 class CoherenceMap(Figure):
674 class CoherenceMap_(Figure):
671 675 isConfig = None
672 676 __nsubplots = None
673 677
@@ -878,7 +882,7 class CoherenceMap(Figure):
878 882 return dataOut
879 883
880 884 @MPDecorator
881 class PowerProfilePlot(Figure):
885 class PowerProfilePlot_(Figure):
882 886
883 887 isConfig = None
884 888 __nsubplots = None
@@ -1008,7 +1012,7 class PowerProfilePlot(Figure):
1008 1012 return dataOut
1009 1013
1010 1014 @MPDecorator
1011 class SpectraCutPlot(Figure):
1015 class SpectraCutPlot_(Figure):
1012 1016
1013 1017 isConfig = None
1014 1018 __nsubplots = None
@@ -1145,7 +1149,7 class SpectraCutPlot(Figure):
1145 1149 return dataOut
1146 1150
1147 1151 @MPDecorator
1148 class Noise(Figure):
1152 class Noise_(Figure):
1149 1153
1150 1154 isConfig = None
1151 1155 __nsubplots = None
@@ -1352,7 +1356,7 class Noise(Figure):
1352 1356 return dataOut
1353 1357
1354 1358 @MPDecorator
1355 class BeaconPhase(Figure):
1359 class BeaconPhase_(Figure):
1356 1360
1357 1361 __isConfig = None
1358 1362 __nsubplots = None
@@ -1497,9 +1501,6 class BeaconPhase(Figure):
1497 1501 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1498 1502 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1499 1503
1500 #print "Phase %d%d" %(pair[0], pair[1])
1501 #print phase[dataOut.beacon_heiIndexList]
1502
1503 1504 if dataOut.beacon_heiIndexList:
1504 1505 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1505 1506 else:
Show file before | Show file after | Hide comments
@@ -12,7 +12,7 from .figure import Figure
12 12
13 13
14 14 @MPDecorator
15 class Scope(Figure):
15 class Scope_(Figure):
16 16
17 17 isConfig = None
18 18
Show file before | Show file after | Hide comments
@@ -434,7 +434,6 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel=''
434 434 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
435 435
436 436 ax = iplot.axes
437
438 437 printLabels(ax, xlabel, ylabel, title)
439 438
440 439 for i in range(len(ax.lines)):
Show file before | Show file after | Hide comments
@@ -13,7 +13,7 import datetime
13 13
14 14 import numpy
15 15
16 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator
17 17 from schainpy.model.data.jrodata import Parameters
18 18 from schainpy.model.io.jroIO_base import JRODataReader, isNumber
19 19 from schainpy.utils import log
@@ -111,7 +111,6 class BLTRParamReader(JRODataReader, ProcessingUnit):
111 111 timezone=0,
112 112 status_value=0,
113 113 **kwargs):
114
115 114 self.path = path
116 115 self.startDate = startDate
117 116 self.endDate = endDate
@@ -185,7 +184,6 class BLTRParamReader(JRODataReader, ProcessingUnit):
185 184 file_id = self.fileIndex
186 185
187 186 if file_id == len(self.fileList):
188 log.success('No more files in the folder', 'BLTRParamReader')
189 187 self.flagNoMoreFiles = 1
190 188 return 0
191 189
@@ -240,7 +238,7 class BLTRParamReader(JRODataReader, ProcessingUnit):
240 238
241 239 pointer = self.fp.tell()
242 240 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
243 self.nchannels = header_rec['nchan'][0] / 2
241 self.nchannels = int(header_rec['nchan'][0] / 2)
244 242 self.kchan = header_rec['nrxs'][0]
245 243 self.nmodes = header_rec['nmodes'][0]
246 244 self.nranges = header_rec['nranges'][0]
@@ -358,8 +356,7 class BLTRParamReader(JRODataReader, ProcessingUnit):
358 356 '''
359 357 if self.flagNoMoreFiles:
360 358 self.dataOut.flagNoData = True
361 log.success('No file left to process', 'BLTRParamReader')
362 return 0
359 self.dataOut.error = (1, 'No More files to read')
363 360
364 361 if not self.readNextBlock():
365 362 self.dataOut.flagNoData = True
Show file before | Show file after | Hide comments
@@ -1815,7 +1815,7 class JRODataWriter(JRODataIO):
1815 1815
1816 1816 return 1
1817 1817
1818 def run(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1818 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1819 1819
1820 1820 if not(self.isConfig):
1821 1821
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (507 lines changed) Show them Hide them
@@ -63,9 +63,6 class Header(object):
63 63 if attr:
64 64 message += "%s = %s" % ("size", attr) + "\n"
65 65
66 # print message
67
68
69 66 FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes
70 67 ('FileMgcNumber', '<u4'), # 0x23020100
71 68 # No Of FDT data records in this file (0 or more)
@@ -94,29 +91,6 class FileHeaderBLTR(Header):
94 91
95 92 header = numpy.fromfile(startFp, FILE_STRUCTURE, 1)
96 93
97 print(' ')
98 print('puntero file header', startFp.tell())
99 print(' ')
100
101 ''' numpy.fromfile(file, dtype, count, sep='')
102 file : file or str
103 Open file object or filename.
104
105 dtype : data-type
106 Data type of the returned array. For binary files, it is used to determine
107 the size and byte-order of the items in the file.
108
109 count : int
110 Number of items to read. -1 means all items (i.e., the complete file).
111
112 sep : str
113 Separator between items if file is a text file. Empty ("") separator means
114 the file should be treated as binary. Spaces (" ") in the separator match zero
115 or more whitespace characters. A separator consisting only of spaces must match
116 at least one whitespace.
117
118 '''
119
120 94 self.FileMgcNumber = hex(header['FileMgcNumber'][0])
121 95 # No Of FDT data records in this file (0 or more)
122 96 self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
@@ -124,8 +98,6 class FileHeaderBLTR(Header):
124 98 self.OffsetStartHeader = int(header['OffsetStartHeader'][0])
125 99 self.SiteName = str(header['SiteName'][0])
126 100
127 # print 'Numero de bloques', self.nFDTdataRecors
128
129 101 if self.size < 48:
130 102 return 0
131 103
@@ -316,36 +288,10 class RecordHeaderBLTR(Header):
316 288 self.OffsetStartHeader = 48
317 289
318 290 def RHread(self, fp):
319 # print fp
320 # startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
321 # The method tell() returns the current position of the file read/write pointer within the file.
322 291 startFp = open(fp, "rb")
323 # RecCounter=0
324 # Off2StartNxtRec=811248
325 292 OffRHeader = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
326 print(' ')
327 print('puntero Record Header', startFp.tell())
328 print(' ')
329
330 293 startFp.seek(OffRHeader, os.SEEK_SET)
331
332 print(' ')
333 print('puntero Record Header con seek', startFp.tell())
334 print(' ')
335
336 # print 'Posicion del bloque: ',OffRHeader
337
338 294 header = numpy.fromfile(startFp, RECORD_STRUCTURE, 1)
339
340 print(' ')
341 print('puntero Record Header con seek', startFp.tell())
342 print(' ')
343
344 print(' ')
345 #
346 # print 'puntero Record Header despues de seek', header.tell()
347 print(' ')
348
349 295 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) # 0x23030001
350 296 self.RecCounter = int(header['RecCounter'][0])
351 297 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
@@ -397,51 +343,8 class RecordHeaderBLTR(Header):
397 343
398 344 self.RHsize = 180 + 20 * self.nChannels
399 345 self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
400 # print 'Datasize',self.Datasize
401 346 endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
402 347
403 print('==============================================')
404 print('RecMgcNumber ', self.RecMgcNumber)
405 print('RecCounter ', self.RecCounter)
406 print('Off2StartNxtRec ', self.Off2StartNxtRec)
407 print('Off2StartData ', self.Off2StartData)
408 print('Range Resolution ', self.SampResolution)
409 print('First Height ', self.StartRangeSamp)
410 print('PRF (Hz) ', self.PRFhz)
411 print('Heights (K) ', self.nHeights)
412 print('Channels (N) ', self.nChannels)
413 print('Profiles (J) ', self.nProfiles)
414 print('iCoh ', self.nCohInt)
415 print('iInCoh ', self.nIncohInt)
416 print('BeamAngleAzim ', self.BeamAngleAzim)
417 print('BeamAngleZen ', self.BeamAngleZen)
418
419 # print 'ModoEnUso ',self.DualModeIndex
420 # print 'UtcTime ',self.nUtime
421 # print 'MiliSec ',self.nMilisec
422 # print 'Exp TagName ',self.ExpTagName
423 # print 'Exp Comment ',self.ExpComment
424 # print 'FFT Window Index ',self.FFTwindowingInd
425 # print 'N Dig. Channels ',self.nDigChannels
426 print('Size de bloque ', self.RHsize)
427 print('DataSize ', self.Datasize)
428 print('BeamAngleAzim ', self.BeamAngleAzim)
429 # print 'AntennaCoord0 ',self.AntennaCoord0
430 # print 'AntennaAngl0 ',self.AntennaAngl0
431 # print 'AntennaCoord1 ',self.AntennaCoord1
432 # print 'AntennaAngl1 ',self.AntennaAngl1
433 # print 'AntennaCoord2 ',self.AntennaCoord2
434 # print 'AntennaAngl2 ',self.AntennaAngl2
435 print('RecPhaseCalibr0 ', self.RecPhaseCalibr0)
436 print('RecPhaseCalibr1 ', self.RecPhaseCalibr1)
437 print('RecPhaseCalibr2 ', self.RecPhaseCalibr2)
438 print('RecAmpCalibr0 ', self.RecAmpCalibr0)
439 print('RecAmpCalibr1 ', self.RecAmpCalibr1)
440 print('RecAmpCalibr2 ', self.RecAmpCalibr2)
441 print('ReceiverGaindB0 ', self.ReceiverGaindB0)
442 print('ReceiverGaindB1 ', self.ReceiverGaindB1)
443 print('ReceiverGaindB2 ', self.ReceiverGaindB2)
444 print('==============================================')
445 348
446 349 if OffRHeader > endFp:
447 350 sys.stderr.write(
@@ -537,9 +440,6 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
537 440 FileList.append(IndexFile)
538 441 nFiles += 1
539 442
540 # print 'Files2Read'
541 # print 'Existen '+str(nFiles)+' archivos .fdt'
542
543 443 self.filenameList = FileList # List of files from least to largest by names
544 444
545 445 def run(self, **kwargs):
@@ -553,7 +453,6 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
553 453 self.isConfig = True
554 454
555 455 self.getData()
556 # print 'running'
557 456
558 457 def setup(self, path=None,
559 458 startDate=None,
@@ -590,7 +489,6 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
590 489
591 490 if self.flagNoMoreFiles:
592 491 self.dataOut.flagNoData = True
593 print('NoData se vuelve true')
594 492 return 0
595 493
596 494 self.fp = self.path
@@ -600,11 +498,9 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
600 498 self.dataOut.data_cspc = self.data_cspc
601 499 self.dataOut.data_output = self.data_output
602 500
603 print('self.dataOut.data_output', shape(self.dataOut.data_output))
604
605 # self.removeDC()
606 501 return self.dataOut.data_spc
607 502
503
608 504 def readFile(self, fp):
609 505 '''
610 506 You must indicate if you are reading in Online or Offline mode and load the
@@ -616,21 +512,16 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
616 512 2. Start reading the first block.
617 513 '''
618 514
619 # The address of the folder is generated the name of the .fdt file that will be read
620 print("File: ", self.fileSelector + 1)
621
622 515 if self.fileSelector < len(self.filenameList):
623 516
624 517 self.fpFile = str(fp) + '/' + \
625 518 str(self.filenameList[self.fileSelector])
626 # print self.fpFile
627 519 fheader = FileHeaderBLTR()
628 520 fheader.FHread(self.fpFile) # Bltr FileHeader Reading
629 521 self.nFDTdataRecors = fheader.nFDTdataRecors
630 522
631 523 self.readBlock() # Block reading
632 524 else:
633 print('readFile FlagNoData becomes true')
634 525 self.flagNoMoreFiles = True
635 526 self.dataOut.flagNoData = True
636 527 return 0
@@ -659,8 +550,7 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
659 550
660 551 '''
661 552
662 if self.BlockCounter < self.nFDTdataRecors - 2:
663 print(self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
553 if self.BlockCounter < self.nFDTdataRecors-1:
664 554 if self.ReadMode == 1:
665 555 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter + 1)
666 556 elif self.ReadMode == 0:
@@ -683,12 +573,10 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
683 573
684 574 self.nRdPairs = len(self.dataOut.pairsList)
685 575 self.dataOut.nRdPairs = self.nRdPairs
686
687 576 self.__firstHeigth = rheader.StartRangeSamp
688 577 self.__deltaHeigth = rheader.SampResolution
689 self.dataOut.heightList = self.__firstHeigth + \
690 numpy.array(list(range(self.nHeights))) * self.__deltaHeigth
691 self.dataOut.channelList = list(range(self.nChannels))
578 self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth
579 self.dataOut.channelList = range(self.nChannels)
692 580 self.dataOut.nProfiles = rheader.nProfiles
693 581 self.dataOut.nIncohInt = rheader.nIncohInt
694 582 self.dataOut.nCohInt = rheader.nCohInt
@@ -697,13 +585,10 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
697 585 self.dataOut.nFFTPoints = rheader.nProfiles
698 586 self.dataOut.utctime = rheader.nUtime
699 587 self.dataOut.timeZone = 0
700 self.dataOut.normFactor = self.dataOut.nProfiles * \
701 self.dataOut.nIncohInt * self.dataOut.nCohInt
702 self.dataOut.outputInterval = self.dataOut.ippSeconds * \
703 self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
588 self.dataOut.normFactor= self.dataOut.nProfiles*self.dataOut.nIncohInt*self.dataOut.nCohInt
589 self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
704 590
705 591 self.data_output = numpy.ones([3, rheader.nHeights]) * numpy.NaN
706 print('self.data_output', shape(self.data_output))
707 592 self.dataOut.velocityX = []
708 593 self.dataOut.velocityY = []
709 594 self.dataOut.velocityV = []
@@ -736,13 +621,12 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
736 621
737 622 if self.DualModeIndex == self.ReadMode:
738 623
739 self.data_fft = numpy.fromfile(
740 startDATA, [('complex', '<c8')], self.nProfiles * self.nChannels * self.nHeights)
624 self.data_fft = numpy.fromfile( startDATA, [('complex','<c8')],self.nProfiles*self.nChannels*self.nHeights )
625 self.data_fft = numpy.empty(101376)
741 626
742 627 self.data_fft = self.data_fft.astype(numpy.dtype('complex'))
743 628
744 self.data_block = numpy.reshape(
745 self.data_fft, (self.nHeights, self.nChannels, self.nProfiles))
629 self.data_block=numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles ))
746 630
747 631 self.data_block = numpy.transpose(self.data_block, (1, 2, 0))
748 632
@@ -756,18 +640,8 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
756 640
757 641 z = self.data_spc.copy() # /factor
758 642 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
759 #zdB = 10*numpy.log10(z)
760 print(' ')
761 print('Z: ')
762 print(shape(z))
763 print(' ')
764 print(' ')
765
766 643 self.dataOut.data_spc = self.data_spc
767
768 self.noise = self.dataOut.getNoise(
769 ymin_index=80, ymax_index=132) # /factor
770 #noisedB = 10*numpy.log10(self.noise)
644 self.noise = self.dataOut.getNoise(ymin_index=80, ymax_index=132)#/factor
771 645
772 646 ySamples = numpy.ones([3, self.nProfiles])
773 647 phase = numpy.ones([3, self.nProfiles])
@@ -784,12 +658,8 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
784 658 xFrec = self.getVelRange(1)
785 659 VelRange = self.getVelRange(1)
786 660 self.dataOut.VelRange = VelRange
787 # print ' '
788 # print ' '
789 # print 'xFrec',xFrec
790 # print ' '
791 # print ' '
792 # Height=35
661
662
793 663 for i in range(self.nRdPairs):
794 664
795 665 chan_index0 = self.dataOut.pairsList[i][0]
@@ -820,361 +690,8 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
820 690
821 691 self.dataOut.ChanDist = self.ChanDist
822 692
823
824 # for Height in range(self.nHeights):
825 #
826 # for i in range(self.nRdPairs):
827 #
828 # '''****** Line of Data SPC ******'''
829 # zline=z[i,:,Height]
830 #
831 # '''****** DC is removed ******'''
832 # DC=Find(zline,numpy.amax(zline))
833 # zline[DC]=(zline[DC-1]+zline[DC+1])/2
834 #
835 #
836 # '''****** SPC is normalized ******'''
837 # FactNorm= zline.copy() / numpy.sum(zline.copy())
838 # FactNorm= FactNorm/numpy.sum(FactNorm)
839 #
840 # SmoothSPC=moving_average(FactNorm,N=3)
841 #
842 # xSamples = ar(range(len(SmoothSPC)))
843 # ySamples[i] = SmoothSPC-self.noise[i]
844 #
845 # for i in range(self.nRdPairs):
846 #
847 # '''****** Line of Data CSPC ******'''
848 # cspcLine=self.data_cspc[i,:,Height].copy()
849 #
850 #
851 #
852 # '''****** CSPC is normalized ******'''
853 # chan_index0 = self.dataOut.pairsList[i][0]
854 # chan_index1 = self.dataOut.pairsList[i][1]
855 # CSPCFactor= numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])
856 #
857 #
858 # CSPCNorm= cspcLine.copy() / numpy.sqrt(CSPCFactor)
859 #
860 #
861 # CSPCSamples[i] = CSPCNorm-self.noise[i]
862 # coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
863 #
864 # '''****** DC is removed ******'''
865 # DC=Find(coherence[i],numpy.amax(coherence[i]))
866 # coherence[i][DC]=(coherence[i][DC-1]+coherence[i][DC+1])/2
867 # coherence[i]= moving_average(coherence[i],N=2)
868 #
869 # phase[i] = moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
870 #
871 #
872 # '''****** Getting fij width ******'''
873 #
874 # yMean=[]
875 # yMean2=[]
876 #
877 # for j in range(len(ySamples[1])):
878 # yMean=numpy.append(yMean,numpy.average([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
879 #
880 # '''******* Getting fitting Gaussian ******'''
881 # meanGauss=sum(xSamples*yMean) / len(xSamples)
882 # sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
883 # #print 'Height',Height,'SNR', meanGauss/sigma**2
884 #
885 # if (abs(meanGauss/sigma**2) > 0.0001) :
886 #
887 # try:
888 # popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
889 #
890 # if numpy.amax(popt)>numpy.amax(yMean)*0.3:
891 # FitGauss=gaus(xSamples,*popt)
892 #
893 # else:
894 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
895 # print 'Verificador: Dentro', Height
896 # except RuntimeError:
897 #
898 # try:
899 # for j in range(len(ySamples[1])):
900 # yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]]))
901 # popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma])
902 # FitGauss=gaus(xSamples,*popt)
903 # print 'Verificador: Exepcion1', Height
904 # except RuntimeError:
905 #
906 # try:
907 # popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma])
908 # FitGauss=gaus(xSamples,*popt)
909 # print 'Verificador: Exepcion2', Height
910 # except RuntimeError:
911 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
912 # print 'Verificador: Exepcion3', Height
913 # else:
914 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
915 # #print 'Verificador: Fuera', Height
916 #
917 #
918 #
919 # Maximun=numpy.amax(yMean)
920 # eMinus1=Maximun*numpy.exp(-1)
921 #
922 # HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
923 # HalfWidth= xFrec[HWpos]
924 # GCpos=Find(FitGauss, numpy.amax(FitGauss))
925 # Vpos=Find(FactNorm, numpy.amax(FactNorm))
926 # #Vpos=numpy.sum(FactNorm)/len(FactNorm)
927 # #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) )))
928 # #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos
929 # '''****** Getting Fij ******'''
930 #
931 # GaussCenter=xFrec[GCpos]
932 # if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
933 # Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
934 # else:
935 # Fij=abs(GaussCenter-HalfWidth)+0.0000001
936 #
937 # '''****** Getting Frecuency range of significant data ******'''
938 #
939 # Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
940 #
941 # if Rangpos<GCpos:
942 # Range=numpy.array([Rangpos,2*GCpos-Rangpos])
943 # else:
944 # Range=numpy.array([2*GCpos-Rangpos,Rangpos])
945 #
946 # FrecRange=xFrec[Range[0]:Range[1]]
947 #
948 # #print 'FrecRange', FrecRange
949 # '''****** Getting SCPC Slope ******'''
950 #
951 # for i in range(self.nRdPairs):
952 #
953 # if len(FrecRange)>5 and len(FrecRange)<self.nProfiles*0.5:
954 # PhaseRange=moving_average(phase[i,Range[0]:Range[1]],N=3)
955 #
956 # slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
957 # PhaseSlope[i]=slope
958 # PhaseInter[i]=intercept
959 # else:
960 # PhaseSlope[i]=0
961 # PhaseInter[i]=0
962 #
963 # # plt.figure(i+15)
964 # # plt.title('FASE ( CH%s*CH%s )' %(self.dataOut.pairsList[i][0],self.dataOut.pairsList[i][1]))
965 # # plt.xlabel('Frecuencia (KHz)')
966 # # plt.ylabel('Magnitud')
967 # # #plt.subplot(311+i)
968 # # plt.plot(FrecRange,PhaseRange,'b')
969 # # plt.plot(FrecRange,FrecRange*PhaseSlope[i]+PhaseInter[i],'r')
970 #
971 # #plt.axis([-0.6, 0.2, -3.2, 3.2])
972 #
973 #
974 # '''Getting constant C'''
975 # cC=(Fij*numpy.pi)**2
976 #
977 # # '''Getting Eij and Nij'''
978 # # (AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180)
979 # # (AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180)
980 # # (AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180)
981 # #
982 # # E01=AntennaX0-AntennaX1
983 # # N01=AntennaY0-AntennaY1
984 # #
985 # # E02=AntennaX0-AntennaX2
986 # # N02=AntennaY0-AntennaY2
987 # #
988 # # E12=AntennaX1-AntennaX2
989 # # N12=AntennaY1-AntennaY2
990 #
991 # '''****** Getting constants F and G ******'''
992 # MijEijNij=numpy.array([[E02,N02], [E12,N12]])
993 # MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
994 # MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
995 # MijResults=numpy.array([MijResult0,MijResult1])
996 # (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
997 #
998 # '''****** Getting constants A, B and H ******'''
999 # W01=numpy.amax(coherence[0])
1000 # W02=numpy.amax(coherence[1])
1001 # W12=numpy.amax(coherence[2])
1002 #
1003 # WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1004 # WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1005 # WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1006 #
1007 # WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1008 #
1009 # WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
1010 # (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1011 #
1012 # VxVy=numpy.array([[cA,cH],[cH,cB]])
1013 #
1014 # VxVyResults=numpy.array([-cF,-cG])
1015 # (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1016 # Vzon = Vy
1017 # Vmer = Vx
1018 # Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1019 # Vang=numpy.arctan2(Vmer,Vzon)
1020 #
1021 # if abs(Vy)<100 and abs(Vy)> 0.:
1022 # self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag
1023 # #print 'Vmag',Vmag
1024 # else:
1025 # self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN)
1026 #
1027 # if abs(Vx)<100 and abs(Vx) > 0.:
1028 # self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang
1029 # #print 'Vang',Vang
1030 # else:
1031 # self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN)
1032 #
1033 # if abs(GaussCenter)<2:
1034 # self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos])
1035 #
1036 # else:
1037 # self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN)
1038 #
1039 #
1040 # # print '********************************************'
1041 # # print 'HalfWidth ', HalfWidth
1042 # # print 'Maximun ', Maximun
1043 # # print 'eMinus1 ', eMinus1
1044 # # print 'Rangpos ', Rangpos
1045 # # print 'GaussCenter ',GaussCenter
1046 # # print 'E01 ',E01
1047 # # print 'N01 ',N01
1048 # # print 'E02 ',E02
1049 # # print 'N02 ',N02
1050 # # print 'E12 ',E12
1051 # # print 'N12 ',N12
1052 # #print 'self.dataOut.velocityX ', self.dataOut.velocityX
1053 # # print 'Fij ', Fij
1054 # # print 'cC ', cC
1055 # # print 'cF ', cF
1056 # # print 'cG ', cG
1057 # # print 'cA ', cA
1058 # # print 'cB ', cB
1059 # # print 'cH ', cH
1060 # # print 'Vx ', Vx
1061 # # print 'Vy ', Vy
1062 # # print 'Vmag ', Vmag
1063 # # print 'Vang ', Vang*180/numpy.pi
1064 # # print 'PhaseSlope ',PhaseSlope[0]
1065 # # print 'PhaseSlope ',PhaseSlope[1]
1066 # # print 'PhaseSlope ',PhaseSlope[2]
1067 # # print '********************************************'
1068 # #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY)
1069 #
1070 # #print 'self.dataOut.velocityX', len(self.dataOut.velocityX)
1071 # #print 'self.dataOut.velocityY', len(self.dataOut.velocityY)
1072 # #print 'self.dataOut.velocityV', self.dataOut.velocityV
1073 #
1074 # self.data_output[0]=numpy.array(self.dataOut.velocityX)
1075 # self.data_output[1]=numpy.array(self.dataOut.velocityY)
1076 # self.data_output[2]=numpy.array(self.dataOut.velocityV)
1077 #
1078 # prin= self.data_output[0][~numpy.isnan(self.data_output[0])]
1079 # print ' '
1080 # print 'VmagAverage',numpy.mean(prin)
1081 # print ' '
1082 # # plt.figure(5)
1083 # # plt.subplot(211)
1084 # # plt.plot(self.dataOut.velocityX,'yo:')
1085 # # plt.subplot(212)
1086 # # plt.plot(self.dataOut.velocityY,'yo:')
1087 #
1088 # # plt.figure(1)
1089 # # # plt.subplot(121)
1090 # # # plt.plot(xFrec,ySamples[0],'k',label='Ch0')
1091 # # # plt.plot(xFrec,ySamples[1],'g',label='Ch1')
1092 # # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
1093 # # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
1094 # # # plt.legend()
1095 # # plt.title('DATOS A ALTURA DE 2850 METROS')
1096 # #
1097 # # plt.xlabel('Frecuencia (KHz)')
1098 # # plt.ylabel('Magnitud')
1099 # # # plt.subplot(122)
1100 # # # plt.title('Fit for Time Constant')
1101 # # #plt.plot(xFrec,zline)
1102 # # #plt.plot(xFrec,SmoothSPC,'g')
1103 # # plt.plot(xFrec,FactNorm)
1104 # # plt.axis([-4, 4, 0, 0.15])
1105 # # # plt.xlabel('SelfSpectra KHz')
1106 # #
1107 # # plt.figure(10)
1108 # # # plt.subplot(121)
1109 # # plt.plot(xFrec,ySamples[0],'b',label='Ch0')
1110 # # plt.plot(xFrec,ySamples[1],'y',label='Ch1')
1111 # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
1112 # # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
1113 # # plt.legend()
1114 # # plt.title('SELFSPECTRA EN CANALES')
1115 # #
1116 # # plt.xlabel('Frecuencia (KHz)')
1117 # # plt.ylabel('Magnitud')
1118 # # # plt.subplot(122)
1119 # # # plt.title('Fit for Time Constant')
1120 # # #plt.plot(xFrec,zline)
1121 # # #plt.plot(xFrec,SmoothSPC,'g')
1122 # # # plt.plot(xFrec,FactNorm)
1123 # # # plt.axis([-4, 4, 0, 0.15])
1124 # # # plt.xlabel('SelfSpectra KHz')
1125 # #
1126 # # plt.figure(9)
1127 # #
1128 # #
1129 # # plt.title('DATOS SUAVIZADOS')
1130 # # plt.xlabel('Frecuencia (KHz)')
1131 # # plt.ylabel('Magnitud')
1132 # # plt.plot(xFrec,SmoothSPC,'g')
1133 # #
1134 # # #plt.plot(xFrec,FactNorm)
1135 # # plt.axis([-4, 4, 0, 0.15])
1136 # # # plt.xlabel('SelfSpectra KHz')
1137 # # #
1138 # # plt.figure(2)
1139 # # # #plt.subplot(121)
1140 # # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra')
1141 # # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano')
1142 # # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo')
1143 # # # #plt.plot(xFrec,phase)
1144 # # # plt.xlabel('Suavizado, promediado KHz')
1145 # # plt.title('SELFSPECTRA PROMEDIADO')
1146 # # # #plt.subplot(122)
1147 # # # #plt.plot(xSamples,zline)
1148 # # plt.xlabel('Frecuencia (KHz)')
1149 # # plt.ylabel('Magnitud')
1150 # # plt.legend()
1151 # # #
1152 # # # plt.figure(3)
1153 # # # plt.subplot(311)
1154 # # # #plt.plot(xFrec,phase[0])
1155 # # # plt.plot(xFrec,phase[0],'g')
1156 # # # plt.subplot(312)
1157 # # # plt.plot(xFrec,phase[1],'g')
1158 # # # plt.subplot(313)
1159 # # # plt.plot(xFrec,phase[2],'g')
1160 # # # #plt.plot(xFrec,phase[2])
1161 # # #
1162 # # # plt.figure(4)
1163 # # #
1164 # # # plt.plot(xSamples,coherence[0],'b')
1165 # # # plt.plot(xSamples,coherence[1],'r')
1166 # # # plt.plot(xSamples,coherence[2],'g')
1167 # # plt.show()
1168 # # #
1169 # # # plt.clf()
1170 # # # plt.cla()
1171 # # # plt.close()
1172 #
1173 # print ' '
1174
1175 693 self.BlockCounter += 2
1176 694
1177 695 else:
1178 696 self.fileSelector += 1
1179 697 self.BlockCounter = 0
1180 print("Next File") No newline at end of file
Show file before | Show file after | Hide comments
@@ -179,9 +179,6 class ParamReader(JRODataReader,ProcessingUnit):
179 179 print("[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime))
180 180 print()
181 181
182 # for i in range(len(filenameList)):
183 # print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
184
185 182 self.filenameList = filenameList
186 183 self.datetimeList = datetimeList
187 184
@@ -504,20 +501,11 class ParamReader(JRODataReader,ProcessingUnit):
504 501
505 502 def getData(self):
506 503
507 # if self.flagNoMoreFiles:
508 # self.dataOut.flagNoData = True
509 # print 'Process finished'
510 # return 0
511 #
512 504 if self.blockIndex==self.blocksPerFile:
513 505 if not( self.__setNextFileOffline() ):
514 506 self.dataOut.flagNoData = True
515 507 return 0
516 508
517 # if self.datablock == None: # setear esta condicion cuando no hayan datos por leers
518 # self.dataOut.flagNoData = True
519 # return 0
520 # self.__readData()
521 509 self.__setDataOut()
522 510 self.dataOut.flagNoData = False
523 511
@@ -637,7 +625,10 class ParamWriter(Operation):
637 625 dsDict['variable'] = self.dataList[i]
638 626 #--------------------- Conditionals ------------------------
639 627 #There is no data
628
629
640 630 if dataAux is None:
631
641 632 return 0
642 633
643 634 #Not array, just a number
@@ -1095,7 +1086,6 class ParamWriter(Operation):
1095 1086 return
1096 1087
1097 1088 self.isConfig = True
1098 # self.putMetadata()
1099 1089 self.setNextFile()
1100 1090
1101 1091 self.putData()
Show file before | Show file after | Hide comments
@@ -413,9 +413,7 class SpectraWriter(JRODataWriter, Operation):
413 413
414 414 data_dc = None
415 415
416 # dataOut = None
417
418 def __init__(self):#, **kwargs):
416 def __init__(self):
419 417 """
420 418 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
421 419
@@ -429,9 +427,7 class SpectraWriter(JRODataWriter, Operation):
429 427 Return: None
430 428 """
431 429
432 Operation.__init__(self)#, **kwargs)
433
434 #self.isConfig = False
430 Operation.__init__(self)
435 431
436 432 self.nTotalBlocks = 0
437 433
@@ -496,7 +492,7 class SpectraWriter(JRODataWriter, Operation):
496 492
497 493
498 494 def writeBlock(self):
499 """
495 """processingHeaderObj
500 496 Escribe el buffer en el file designado
501 497
502 498 Affected:
@@ -519,8 +515,10 class SpectraWriter(JRODataWriter, Operation):
519 515 data.tofile(self.fp)
520 516
521 517 if self.data_cspc is not None:
522 data = numpy.zeros( self.shape_cspc_Buffer, self.dtype )
518
523 519 cspc = numpy.transpose( self.data_cspc, (0,2,1) )
520 #data = numpy.zeros( numpy.shape(cspc), self.dtype )
521 #print 'data.shape', self.shape_cspc_Buffer
524 522 if not self.processingHeaderObj.shif_fft:
525 523 cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
526 524 data['real'] = cspc.real
@@ -529,8 +527,9 class SpectraWriter(JRODataWriter, Operation):
529 527 data.tofile(self.fp)
530 528
531 529 if self.data_dc is not None:
532 data = numpy.zeros( self.shape_dc_Buffer, self.dtype )
530
533 531 dc = self.data_dc
532 data = numpy.zeros( numpy.shape(dc), self.dtype )
534 533 data['real'] = dc.real
535 534 data['imag'] = dc.imag
536 535 data = data.reshape((-1))
Show file before | Show file after | Hide comments
@@ -12,13 +12,11 from time import gmtime
12 12
13 13 from numpy import transpose
14 14
15 from .jroproc_base import ProcessingUnit, MPDecorator, Operation
15 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
16 16 from schainpy.model.data.jrodata import Parameters
17 17
18 18 @MPDecorator
19 19 class BLTRParametersProc(ProcessingUnit):
20
21 METHODS = {}
22 20 '''
23 21 Processing unit for BLTR parameters data (winds)
24 22
@@ -46,9 +44,7 class BLTRParametersProc(ProcessingUnit):
46 44 Inputs: None
47 45 '''
48 46 ProcessingUnit.__init__(self)
49 self.setupReq = False
50 47 self.dataOut = Parameters()
51 self.isConfig = False
52 48
53 49 def setup(self, mode):
54 50 '''
Show file before | Show file after | Hide comments
@@ -11,13 +11,15 Based on:
11 11 $Author: murco $
12 12 $Id: jroproc_base.py 1 2012-11-12 18:56:07Z murco $
13 13 '''
14 from platform import python_version
14
15 15 import inspect
16 16 import zmq
17 17 import time
18 18 import pickle
19 19 import os
20 20 from multiprocessing import Process
21 from zmq.utils.monitor import recv_monitor_message
22
21 23 from schainpy.utils import log
22 24
23 25
@@ -36,21 +38,13 class ProcessingUnit(object):
36 38
37 39 """
38 40
39 METHODS = {}
40 dataIn = None
41 dataInList = []
42 id = None
43 inputId = None
44 dataOut = None
45 dictProcs = None
46 isConfig = False
47
48 41 def __init__(self):
49 42
50 43 self.dataIn = None
51 44 self.dataOut = None
52 45 self.isConfig = False
53 46 self.operations = []
47 self.plots = []
54 48
55 49 def getAllowedArgs(self):
56 50 if hasattr(self, '__attrs__'):
@@ -59,7 +53,6 class ProcessingUnit(object):
59 53 return inspect.getargspec(self.run).args
60 54
61 55 def addOperation(self, conf, operation):
62
63 56 """
64 57 This method is used in the controller, and update the dictionary containing the operations to execute. The dict
65 58 posses the id of the operation process (IPC purposes)
@@ -76,7 +69,11 class ProcessingUnit(object):
76 69 objId : identificador del objeto, necesario para comunicar con master(procUnit)
77 70 """
78 71
79 self.operations.append((operation, conf.type, conf.id, conf.getKwargs()))
72 self.operations.append(
73 (operation, conf.type, conf.id, conf.getKwargs()))
74
75 if 'plot' in self.name.lower():
76 self.plots.append(operation.CODE)
80 77
81 78 def getOperationObj(self, objId):
82 79
@@ -86,7 +83,6 class ProcessingUnit(object):
86 83 return self.operations[objId]
87 84
88 85 def operation(self, **kwargs):
89
90 86 """
91 87 Operacion directa sobre la data (dataOut.data). Es necesario actualizar los valores de los
92 88 atributos del objeto dataOut
@@ -107,9 +103,10 class ProcessingUnit(object):
107 103 raise NotImplementedError
108 104
109 105 def close(self):
110 #Close every thread, queue or any other object here is it is neccesary.
106
111 107 return
112 108
109
113 110 class Operation(object):
114 111
115 112 """
@@ -126,18 +123,11 class Operation(object):
126 123 Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer)
127 124
128 125 """
129 id = None
130 __buffer = None
131 dest = None
132 isConfig = False
133 readyFlag = None
134 126
135 127 def __init__(self):
136 128
137 self.buffer = None
138 self.dest = None
129 self.id = None
139 130 self.isConfig = False
140 self.readyFlag = False
141 131
142 132 if not hasattr(self, 'name'):
143 133 self.name = self.__class__.__name__
@@ -154,9 +144,7 class Operation(object):
154 144
155 145 raise NotImplementedError
156 146
157
158 147 def run(self, dataIn, **kwargs):
159
160 148 """
161 149 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
162 150 atributos del objeto dataIn.
@@ -180,11 +168,10 class Operation(object):
180 168
181 169 def close(self):
182 170
183 pass
171 return
184 172
185 173
186 174 def MPDecorator(BaseClass):
187
188 175 """
189 176 Multiprocessing class decorator
190 177
@@ -203,25 +190,19 def MPDecorator(BaseClass):
203 190 self.sender = None
204 191 self.receiver = None
205 192 self.name = BaseClass.__name__
193 self.start_time = time.time()
206 194
207 195 if len(self.args) is 3:
208 196 self.typeProc = "ProcUnit"
209 197 self.id = args[0]
210 198 self.inputId = args[1]
211 199 self.project_id = args[2]
212 else:
200 elif len(self.args) is 2:
213 201 self.id = args[0]
214 202 self.inputId = args[0]
215 203 self.project_id = args[1]
216 204 self.typeProc = "Operation"
217 205
218 def getAllowedArgs(self):
219
220 if hasattr(self, '__attrs__'):
221 return self.__attrs__
222 else:
223 return inspect.getargspec(BaseClass.run).args
224
225 206 def subscribe(self):
226 207 '''
227 208 This function create a socket to receive objects from the
@@ -230,7 +211,8 def MPDecorator(BaseClass):
230 211
231 212 c = zmq.Context()
232 213 self.receiver = c.socket(zmq.SUB)
233 self.receiver.connect('ipc:///tmp/schain/{}_pub'.format(self.project_id))
214 self.receiver.connect(
215 'ipc:///tmp/schain/{}_pub'.format(self.project_id))
234 216 self.receiver.setsockopt(zmq.SUBSCRIBE, self.inputId.encode())
235 217
236 218 def listen(self):
@@ -239,6 +221,7 def MPDecorator(BaseClass):
239 221 '''
240 222
241 223 data = pickle.loads(self.receiver.recv_multipart()[1])
224
242 225 return data
243 226
244 227 def set_publisher(self):
@@ -249,13 +232,13 def MPDecorator(BaseClass):
249 232 time.sleep(1)
250 233 c = zmq.Context()
251 234 self.sender = c.socket(zmq.PUB)
252 self.sender.connect('ipc:///tmp/schain/{}_sub'.format(self.project_id))
235 self.sender.connect(
236 'ipc:///tmp/schain/{}_sub'.format(self.project_id))
253 237
254 238 def publish(self, data, id):
255 239 '''
256 240 This function publish an object, to a specific topic.
257 241 '''
258
259 242 self.sender.send_multipart([str(id).encode(), pickle.dumps(data)])
260 243
261 244 def runReader(self):
@@ -266,13 +249,7 def MPDecorator(BaseClass):
266 249
267 250 BaseClass.run(self, **self.kwargs)
268 251
269 if self.dataOut.error[0] == -1:
270 log.error(self.dataOut.error[1])
271 self.publish('end', self.id)
272 #self.sender.send_multipart([str(self.project_id).encode(), 'end'.encode()])
273 break
274
275 for op, optype, id, kwargs in self.operations:
252 for op, optype, opId, kwargs in self.operations:
276 253 if optype=='self':
277 254 op(**kwargs)
278 255 elif optype=='other':
@@ -280,11 +257,21 def MPDecorator(BaseClass):
280 257 elif optype=='external':
281 258 self.publish(self.dataOut, opId)
282 259
283 if self.dataOut.flagNoData:
260 if self.dataOut.flagNoData and self.dataOut.error is None:
284 261 continue
285 262
286 263 self.publish(self.dataOut, self.id)
287 264
265 if self.dataOut.error:
266 if self.dataOut.error[0] == -1:
267 log.error(self.dataOut.error[1], self.name)
268 if self.dataOut.error[0] == 1:
269 log.success(self.dataOut.error[1], self.name)
270 # self.sender.send_multipart([str(self.project_id).encode(), 'end'.encode()])
271 break
272
273 time.sleep(1)
274
288 275 def runProc(self):
289 276 '''
290 277 Run function for proccessing units
@@ -293,14 +280,7 def MPDecorator(BaseClass):
293 280 while True:
294 281 self.dataIn = self.listen()
295 282
296 if self.dataIn == 'end':
297 self.publish('end', self.id)
298 for op, optype, opId, kwargs in self.operations:
299 if optype == 'external':
300 self.publish('end', opId)
301 break
302
303 if self.dataIn.flagNoData:
283 if self.dataIn.flagNoData and self.dataIn.error is None:
304 284 continue
305 285
306 286 BaseClass.run(self, **self.kwargs)
@@ -313,25 +293,28 def MPDecorator(BaseClass):
313 293 elif optype=='external':
314 294 self.publish(self.dataOut, opId)
315 295
316 if self.dataOut.flagNoData:
317 continue
318
319 296 self.publish(self.dataOut, self.id)
297 if self.dataIn.error:
298 break
299
300 time.sleep(1)
320 301
321 302 def runOp(self):
322 303 '''
323 Run function for operations
304 Run function for external operations (this operations just receive data
305 ex: plots, writers, publishers)
324 306 '''
325 307
326 308 while True:
327 309
328 310 dataOut = self.listen()
329 311
330 if dataOut == 'end':
331 break
332
333 312 BaseClass.run(self, dataOut, **self.kwargs)
334 313
314 if dataOut.error:
315 break
316 time.sleep(1)
317
335 318 def run(self):
336 319
337 320 if self.typeProc is "ProcUnit":
@@ -353,15 +336,41 def MPDecorator(BaseClass):
353 336 else:
354 337 raise ValueError("Unknown type")
355 338
356 print("%s done" % BaseClass.__name__)
357 339 self.close()
358 340
341 def event_monitor(self, monitor):
342
343 events = {}
344
345 for name in dir(zmq):
346 if name.startswith('EVENT_'):
347 value = getattr(zmq, name)
348 events[value] = name
349
350 while monitor.poll():
351 evt = recv_monitor_message(monitor)
352 if evt['event'] == 32:
353 self.connections += 1
354 if evt['event'] == 512:
355 pass
356
357 evt.update({'description': events[evt['event']]})
358
359 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
360 break
361 monitor.close()
362 print('event monitor thread done!')
363
359 364 def close(self):
360 365
366 BaseClass.close(self)
367
361 368 if self.sender:
362 369 self.sender.close()
363 370
364 371 if self.receiver:
365 372 self.receiver.close()
366 373
374 log.success('Done...(Time:{:4.2f} secs)'.format(time.time()-self.start_time), self.name)
375
367 376 return MPClass
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1151 lines changed) Show them Hide them
@@ -10,11 +10,7 import importlib
10 10 import itertools
11 11 from multiprocessing import Pool, TimeoutError
12 12 from multiprocessing.pool import ThreadPool
13 import types
14 from functools import partial
15 13 import time
16 #from sklearn.cluster import KMeans
17
18 14
19 15 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
20 16 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
@@ -128,6 +124,7 class ParametersProc(ProcessingUnit):
128 124 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
129 125 self.dataOut.spc_noise = self.dataIn.getNoise()
130 126 self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
127 # self.dataOut.normFactor = self.dataIn.normFactor
131 128 self.dataOut.pairsList = self.dataIn.pairsList
132 129 self.dataOut.groupList = self.dataIn.pairsList
133 130 self.dataOut.flagNoData = False
@@ -136,9 +133,9 class ParametersProc(ProcessingUnit):
136 133 self.dataOut.ChanDist = self.dataIn.ChanDist
137 134 else: self.dataOut.ChanDist = None
138 135
139 if hasattr(self.dataIn, 'VelRange'): #Velocities range
140 self.dataOut.VelRange = self.dataIn.VelRange
141 else: self.dataOut.VelRange = None
136 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
137 # self.dataOut.VelRange = self.dataIn.VelRange
138 #else: self.dataOut.VelRange = None
142 139
143 140 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
144 141 self.dataOut.RadarConst = self.dataIn.RadarConst
@@ -184,9 +181,112 class ParametersProc(ProcessingUnit):
184 181 def target(tups):
185 182
186 183 obj, args = tups
187 #print 'TARGETTT', obj, args
184
188 185 return obj.FitGau(args)
189 186
187
188 class SpectralFilters(Operation):
189
190 '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR
191
192 LimitR : It is the limit in m/s of Rainfall
193 LimitW : It is the limit in m/s for Winds
194
195 Input:
196
197 self.dataOut.data_pre : SPC and CSPC
198 self.dataOut.spc_range : To select wind and rainfall velocities
199
200 Affected:
201
202 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
203 self.dataOut.spcparam_range : Used in SpcParamPlot
204 self.dataOut.SPCparam : Used in PrecipitationProc
205
206
207 '''
208
209 def __init__(self):
210 Operation.__init__(self)
211 self.i=0
212
213 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
214
215
216 #Limite de vientos
217 LimitR = PositiveLimit
218 LimitN = NegativeLimit
219
220 self.spc = dataOut.data_pre[0].copy()
221 self.cspc = dataOut.data_pre[1].copy()
222
223 self.Num_Hei = self.spc.shape[2]
224 self.Num_Bin = self.spc.shape[1]
225 self.Num_Chn = self.spc.shape[0]
226
227 VelRange = dataOut.spc_range[2]
228 TimeRange = dataOut.spc_range[1]
229 FrecRange = dataOut.spc_range[0]
230
231 Vmax= 2*numpy.max(dataOut.spc_range[2])
232 Tmax= 2*numpy.max(dataOut.spc_range[1])
233 Fmax= 2*numpy.max(dataOut.spc_range[0])
234
235 Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
236 Breaker1R=numpy.where(VelRange == Breaker1R)
237
238 Delta = self.Num_Bin/2 - Breaker1R[0]
239
240
241 '''Reacomodando SPCrange'''
242
243 VelRange=numpy.roll(VelRange,-(self.Num_Bin/2) ,axis=0)
244
245 VelRange[-(self.Num_Bin/2):]+= Vmax
246
247 FrecRange=numpy.roll(FrecRange,-(self.Num_Bin/2),axis=0)
248
249 FrecRange[-(self.Num_Bin/2):]+= Fmax
250
251 TimeRange=numpy.roll(TimeRange,-(self.Num_Bin/2),axis=0)
252
253 TimeRange[-(self.Num_Bin/2):]+= Tmax
254
255 ''' ------------------ '''
256
257 Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
258 Breaker2R=numpy.where(VelRange == Breaker2R)
259
260
261 SPCroll = numpy.roll(self.spc,-(self.Num_Bin/2) ,axis=1)
262
263 SPCcut = SPCroll.copy()
264 for i in range(self.Num_Chn):
265
266 SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
267 SPCcut[i,-int(Delta):,:] = dataOut.noise[i]
268
269 SPCcut[i]=SPCcut[i]- dataOut.noise[i]
270 SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20
271
272 SPCroll[i]=SPCroll[i]-dataOut.noise[i]
273 SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20
274
275 SPC_ch1 = SPCroll
276
277 SPC_ch2 = SPCcut
278
279 SPCparam = (SPC_ch1, SPC_ch2, self.spc)
280 dataOut.SPCparam = numpy.asarray(SPCparam)
281
282
283 dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1])
284
285 dataOut.spcparam_range[2]=VelRange
286 dataOut.spcparam_range[1]=TimeRange
287 dataOut.spcparam_range[0]=FrecRange
288
289
190 290 class GaussianFit(Operation):
191 291
192 292 '''
@@ -198,15 +298,15 class GaussianFit(Operation):
198 298 self.dataOut.data_pre : SelfSpectra
199 299
200 300 Output:
201 self.dataOut.GauSPC : SPC_ch1, SPC_ch2
301 self.dataOut.SPCparam : SPC_ch1, SPC_ch2
202 302
203 303 '''
204 def __init__(self, **kwargs):
205 Operation.__init__(self, **kwargs)
304 def __init__(self):
305 Operation.__init__(self)
206 306 self.i=0
207 307
208 308
209 def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
309 def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
210 310 """This routine will find a couple of generalized Gaussians to a power spectrum
211 311 input: spc
212 312 output:
@@ -214,31 +314,12 class GaussianFit(Operation):
214 314 """
215 315
216 316 self.spc = dataOut.data_pre[0].copy()
217
218
219 print('SelfSpectra Shape', numpy.asarray(self.spc).shape)
220
221
222 #plt.figure(50)
223 #plt.subplot(121)
224 #plt.plot(self.spc,'k',label='spc(66)')
225 #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
226 #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
227 #plt.plot(xFrec,FitGauss,'yo:',label='fit')
228 #plt.legend()
229 #plt.title('DATOS A ALTURA DE 7500 METROS')
230 #plt.show()
231
232 317 self.Num_Hei = self.spc.shape[2]
233 #self.Num_Bin = len(self.spc)
234 318 self.Num_Bin = self.spc.shape[1]
235 319 self.Num_Chn = self.spc.shape[0]
236
237 320 Vrange = dataOut.abscissaList
238 321
239 #print 'self.spc2', numpy.asarray(self.spc).shape
240
241 GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei])
322 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
242 323 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
243 324 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
244 325 SPC_ch1[:] = numpy.NaN
@@ -250,272 +331,12 class GaussianFit(Operation):
250 331 noise_ = dataOut.spc_noise[0].copy()
251 332
252 333
253
254 334 pool = Pool(processes=self.Num_Chn)
255 335 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
256 336 objs = [self for __ in range(self.Num_Chn)]
257 337 attrs = list(zip(objs, args))
258 338 gauSPC = pool.map(target, attrs)
259 dataOut.GauSPC = numpy.asarray(gauSPC)
260 # ret = []
261 # for n in range(self.Num_Chn):
262 # self.FitGau(args[n])
263 # dataOut.GauSPC = ret
264
265
266
267 # for ch in range(self.Num_Chn):
268 #
269 # for ht in range(self.Num_Hei):
270 # #print (numpy.asarray(self.spc).shape)
271 # spc = numpy.asarray(self.spc)[ch,:,ht]
272 #
273 # #############################################
274 # # normalizing spc and noise
275 # # This part differs from gg1
276 # spc_norm_max = max(spc)
277 # spc = spc / spc_norm_max
278 # pnoise = pnoise / spc_norm_max
279 # #############################################
280 #
281 # if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's
282 # fatspectra=1.0
283 # else:
284 # fatspectra=0.5
285 #
286 # wnoise = noise_ / spc_norm_max
287 # #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise
288 # #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
289 # #if wnoise>1.1*pnoise: # to be tested later
290 # # wnoise=pnoise
291 # noisebl=wnoise*0.9; noisebh=wnoise*1.1
292 # spc=spc-wnoise
293 #
294 # minx=numpy.argmin(spc)
295 # spcs=numpy.roll(spc,-minx)
296 # cum=numpy.cumsum(spcs)
297 # tot_noise=wnoise * self.Num_Bin #64;
298 # #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
299 # #snr=tot_signal/tot_noise
300 # #snr=cum[-1]/tot_noise
301 #
302 # #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
303 #
304 # snr = sum(spcs)/tot_noise
305 # snrdB=10.*numpy.log10(snr)
306 #
307 # #if snrdB < -9 :
308 # # snrdB = numpy.NaN
309 # # continue
310 #
311 # #print 'snr',snrdB # , sum(spcs) , tot_noise
312 #
313 #
314 # #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
315 # # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
316 #
317 # cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
318 # cumlo=cummax*epsi;
319 # cumhi=cummax*(1-epsi)
320 # powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
321 #
322 # #if len(powerindex)==1:
323 # ##return [numpy.mod(powerindex[0]+minx,64),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
324 # #return [numpy.mod(powerindex[0]+minx, self.Num_Bin ),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
325 # #elif len(powerindex)<4*fatspectra:
326 # #return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
327 #
328 # if len(powerindex) < 1:# case for powerindex 0
329 # continue
330 # powerlo=powerindex[0]
331 # powerhi=powerindex[-1]
332 # powerwidth=powerhi-powerlo
333 #
334 # firstpeak=powerlo+powerwidth/10.# first gaussian energy location
335 # secondpeak=powerhi-powerwidth/10.#second gaussian energy location
336 # midpeak=(firstpeak+secondpeak)/2.
337 # firstamp=spcs[int(firstpeak)]
338 # secondamp=spcs[int(secondpeak)]
339 # midamp=spcs[int(midpeak)]
340 # #x=numpy.spc.shape[1]
341 #
342 # #x=numpy.arange(64)
343 # x=numpy.arange( self.Num_Bin )
344 # y_data=spc+wnoise
345 #
346 # # single gaussian
347 # #shift0=numpy.mod(midpeak+minx,64)
348 # shift0=numpy.mod(midpeak+minx, self.Num_Bin )
349 # width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
350 # power0=2.
351 # amplitude0=midamp
352 # state0=[shift0,width0,amplitude0,power0,wnoise]
353 # #bnds=((0,63),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
354 # bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
355 # #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(0.1,0.5))
356 # # bnds = range of fft, power width, amplitude, power, noise
357 # lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
358 #
359 # chiSq1=lsq1[1];
360 # jack1= self.y_jacobian1(x,lsq1[0])
361 #
362 #
363 # try:
364 # sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
365 # except:
366 # std1=32.; sigmas1=numpy.ones(5)
367 # else:
368 # std1=sigmas1[0]
369 #
370 #
371 # if fatspectra<1.0 and powerwidth<4:
372 # choice=0
373 # Amplitude0=lsq1[0][2]
374 # shift0=lsq1[0][0]
375 # width0=lsq1[0][1]
376 # p0=lsq1[0][3]
377 # Amplitude1=0.
378 # shift1=0.
379 # width1=0.
380 # p1=0.
381 # noise=lsq1[0][4]
382 # #return (numpy.array([shift0,width0,Amplitude0,p0]),
383 # # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
384 #
385 # # two gaussians
386 # #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
387 # shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
388 # shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
389 # width0=powerwidth/6.;
390 # width1=width0
391 # power0=2.;
392 # power1=power0
393 # amplitude0=firstamp;
394 # amplitude1=secondamp
395 # state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
396 # #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
397 # bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
398 # #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
399 #
400 # lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
401 #
402 #
403 # chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
404 #
405 #
406 # try:
407 # sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
408 # except:
409 # std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
410 # else:
411 # std2a=sigmas2[0]; std2b=sigmas2[4]
412 #
413 #
414 #
415 # oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
416 #
417 # if snrdB>-9: # when SNR is strong pick the peak with least shift (LOS velocity) error
418 # if oneG:
419 # choice=0
420 # else:
421 # w1=lsq2[0][1]; w2=lsq2[0][5]
422 # a1=lsq2[0][2]; a2=lsq2[0][6]
423 # p1=lsq2[0][3]; p2=lsq2[0][7]
424 # s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
425 # gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
426 #
427 # if gp1>gp2:
428 # if a1>0.7*a2:
429 # choice=1
430 # else:
431 # choice=2
432 # elif gp2>gp1:
433 # if a2>0.7*a1:
434 # choice=2
435 # else:
436 # choice=1
437 # else:
438 # choice=numpy.argmax([a1,a2])+1
439 # #else:
440 # #choice=argmin([std2a,std2b])+1
441 #
442 # else: # with low SNR go to the most energetic peak
443 # choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
444 #
445 # #print 'choice',choice
446 #
447 # if choice==0: # pick the single gaussian fit
448 # Amplitude0=lsq1[0][2]
449 # shift0=lsq1[0][0]
450 # width0=lsq1[0][1]
451 # p0=lsq1[0][3]
452 # Amplitude1=0.
453 # shift1=0.
454 # width1=0.
455 # p1=0.
456 # noise=lsq1[0][4]
457 # elif choice==1: # take the first one of the 2 gaussians fitted
458 # Amplitude0 = lsq2[0][2]
459 # shift0 = lsq2[0][0]
460 # width0 = lsq2[0][1]
461 # p0 = lsq2[0][3]
462 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
463 # shift1 = lsq2[0][4] # This is 0 in gg1
464 # width1 = lsq2[0][5] # This is 0 in gg1
465 # p1 = lsq2[0][7] # This is 0 in gg1
466 # noise = lsq2[0][8]
467 # else: # the second one
468 # Amplitude0 = lsq2[0][6]
469 # shift0 = lsq2[0][4]
470 # width0 = lsq2[0][5]
471 # p0 = lsq2[0][7]
472 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
473 # shift1 = lsq2[0][0] # This is 0 in gg1
474 # width1 = lsq2[0][1] # This is 0 in gg1
475 # p1 = lsq2[0][3] # This is 0 in gg1
476 # noise = lsq2[0][8]
477 #
478 # #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
479 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
480 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
481 # #print 'SPC_ch1.shape',SPC_ch1.shape
482 # #print 'SPC_ch2.shape',SPC_ch2.shape
483 # #dataOut.data_param = SPC_ch1
484 # GauSPC[0] = SPC_ch1
485 # GauSPC[1] = SPC_ch2
486
487 # #plt.gcf().clear()
488 # plt.figure(50+self.i)
489 # self.i=self.i+1
490 # #plt.subplot(121)
491 # plt.plot(self.spc,'k')#,label='spc(66)')
492 # plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1')
493 # #plt.plot(SPC_ch2,'r')#,label='gg2')
494 # #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
495 # #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
496 # #plt.plot(xFrec,FitGauss,'yo:',label='fit')
497 # plt.legend()
498 # plt.title('DATOS A ALTURA DE 7500 METROS')
499 # plt.show()
500 # print 'shift0', shift0
501 # print 'Amplitude0', Amplitude0
502 # print 'width0', width0
503 # print 'p0', p0
504 # print '========================'
505 # print 'shift1', shift1
506 # print 'Amplitude1', Amplitude1
507 # print 'width1', width1
508 # print 'p1', p1
509 # print 'noise', noise
510 # print 's_noise', wnoise
511
512 print('========================================================')
513 print('total_time: ', time.time()-start_time)
514
515 # re-normalizing spc and noise
516 # This part differs from gg1
517
518
339 dataOut.SPCparam = numpy.asarray(SPCparam)
519 340
520 341 ''' Parameters:
521 342 1. Amplitude
@@ -524,16 +345,11 class GaussianFit(Operation):
524 345 4. Power
525 346 '''
526 347
527
528 ###############################################################################
529 348 def FitGau(self, X):
530 349
531 350 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
532 #print 'VARSSSS', ch, pnoise, noise, num_intg
533
534 #print 'HEIGHTS', self.Num_Hei
535 351
536 GauSPC = []
352 SPCparam = []
537 353 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
538 354 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
539 355 SPC_ch1[:] = 0#numpy.NaN
@@ -542,10 +358,6 class GaussianFit(Operation):
542 358
543 359
544 360 for ht in range(self.Num_Hei):
545 #print (numpy.asarray(self.spc).shape)
546
547 #print 'TTTTT', ch , ht
548 #print self.spc.shape
549 361
550 362
551 363 spc = numpy.asarray(self.spc)[ch,:,ht]
@@ -554,27 +366,26 class GaussianFit(Operation):
554 366 # normalizing spc and noise
555 367 # This part differs from gg1
556 368 spc_norm_max = max(spc)
557 spc = spc / spc_norm_max
558 pnoise = pnoise / spc_norm_max
369 #spc = spc / spc_norm_max
370 pnoise = pnoise #/ spc_norm_max
559 371 #############################################
560 372
561 373 fatspectra=1.0
562 374
563 wnoise = noise_ / spc_norm_max
375 wnoise = noise_ #/ spc_norm_max
564 376 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
565 377 #if wnoise>1.1*pnoise: # to be tested later
566 378 # wnoise=pnoise
567 noisebl=wnoise*0.9; noisebh=wnoise*1.1
379 noisebl=wnoise*0.9;
380 noisebh=wnoise*1.1
568 381 spc=spc-wnoise
569 # print 'wnoise', noise_[0], spc_norm_max, wnoise
382
570 383 minx=numpy.argmin(spc)
384 #spcs=spc.copy()
571 385 spcs=numpy.roll(spc,-minx)
572 386 cum=numpy.cumsum(spcs)
573 387 tot_noise=wnoise * self.Num_Bin #64;
574 #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
575 #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
576 #snr=tot_signal/tot_noise
577 #snr=cum[-1]/tot_noise
388
578 389 snr = sum(spcs)/tot_noise
579 390 snrdB=10.*numpy.log10(snr)
580 391
@@ -582,16 +393,15 class GaussianFit(Operation):
582 393 snr = numpy.NaN
583 394 SPC_ch1[:,ht] = 0#numpy.NaN
584 395 SPC_ch1[:,ht] = 0#numpy.NaN
585 GauSPC = (SPC_ch1,SPC_ch2)
396 SPCparam = (SPC_ch1,SPC_ch2)
586 397 continue
587 #print 'snr',snrdB #, sum(spcs) , tot_noise
588
589 398
590 399
591 400 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
592 401 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
593 402
594 cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
403 cummax=max(cum);
404 epsi=0.08*fatspectra # cumsum to narrow down the energy region
595 405 cumlo=cummax*epsi;
596 406 cumhi=cummax*(1-epsi)
597 407 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
@@ -613,7 +423,7 class GaussianFit(Operation):
613 423 x=numpy.arange( self.Num_Bin )
614 424 y_data=spc+wnoise
615 425
616 # single gaussian
426 ''' single Gaussian '''
617 427 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
618 428 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
619 429 power0=2.
@@ -623,15 +433,6 class GaussianFit(Operation):
623 433 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
624 434
625 435 chiSq1=lsq1[1];
626 jack1= self.y_jacobian1(x,lsq1[0])
627
628
629 try:
630 sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
631 except:
632 std1=32.; sigmas1=numpy.ones(5)
633 else:
634 std1=sigmas1[0]
635 436
636 437
637 438 if fatspectra<1.0 and powerwidth<4:
@@ -648,7 +449,7 class GaussianFit(Operation):
648 449 #return (numpy.array([shift0,width0,Amplitude0,p0]),
649 450 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
650 451
651 # two gaussians
452 ''' two gaussians '''
652 453 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
653 454 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
654 455 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
@@ -666,21 +467,13 class GaussianFit(Operation):
666 467 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
667 468
668 469
669 chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
670
671
672 try:
673 sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
674 except:
675 std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
676 else:
677 std2a=sigmas2[0]; std2b=sigmas2[4]
470 chiSq2=lsq2[1];
678 471
679 472
680 473
681 474 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
682 475
683 if snrdB>-6: # when SNR is strong pick the peak with least shift (LOS velocity) error
476 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
684 477 if oneG:
685 478 choice=0
686 479 else:
@@ -710,13 +503,15 class GaussianFit(Operation):
710 503 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
711 504
712 505
713 shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
714 shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
506 shift0=lsq2[0][0];
507 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
508 shift1=lsq2[0][4];
509 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
715 510
716 max_vel = 20
511 max_vel = 1.0
717 512
718 513 #first peak will be 0, second peak will be 1
719 if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range
514 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
720 515 shift0=lsq2[0][0]
721 516 width0=lsq2[0][1]
722 517 Amplitude0=lsq2[0][2]
@@ -739,120 +534,18 class GaussianFit(Operation):
739 534 p0=lsq2[0][7]
740 535 noise=lsq2[0][8]
741 536
742 if Amplitude0<0.1: # in case the peak is noise
743 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
744 if Amplitude1<0.1:
745 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
746
747
748 # if choice==0: # pick the single gaussian fit
749 # Amplitude0=lsq1[0][2]
750 # shift0=lsq1[0][0]
751 # width0=lsq1[0][1]
752 # p0=lsq1[0][3]
753 # Amplitude1=0.
754 # shift1=0.
755 # width1=0.
756 # p1=0.
757 # noise=lsq1[0][4]
758 # elif choice==1: # take the first one of the 2 gaussians fitted
759 # Amplitude0 = lsq2[0][2]
760 # shift0 = lsq2[0][0]
761 # width0 = lsq2[0][1]
762 # p0 = lsq2[0][3]
763 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
764 # shift1 = lsq2[0][4] # This is 0 in gg1
765 # width1 = lsq2[0][5] # This is 0 in gg1
766 # p1 = lsq2[0][7] # This is 0 in gg1
767 # noise = lsq2[0][8]
768 # else: # the second one
769 # Amplitude0 = lsq2[0][6]
770 # shift0 = lsq2[0][4]
771 # width0 = lsq2[0][5]
772 # p0 = lsq2[0][7]
773 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
774 # shift1 = lsq2[0][0] # This is 0 in gg1
775 # width1 = lsq2[0][1] # This is 0 in gg1
776 # p1 = lsq2[0][3] # This is 0 in gg1
777 # noise = lsq2[0][8]
778
779 #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
780 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
781 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
782 #print 'SPC_ch1.shape',SPC_ch1.shape
783 #print 'SPC_ch2.shape',SPC_ch2.shape
784 #dataOut.data_param = SPC_ch1
785 GauSPC = (SPC_ch1,SPC_ch2)
786 #GauSPC[1] = SPC_ch2
787
788 # print 'shift0', shift0
789 # print 'Amplitude0', Amplitude0
790 # print 'width0', width0
791 # print 'p0', p0
792 # print '========================'
793 # print 'shift1', shift1
794 # print 'Amplitude1', Amplitude1
795 # print 'width1', width1
796 # print 'p1', p1
797 # print 'noise', noise
798 # print 's_noise', wnoise
799
800 return GauSPC
801
802
803 def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination.
804 y_model=self.y_model1(x,state)
805 s0,w0,a0,p0,n=state
806 e0=((x-s0)/w0)**2;
807
808 e0u=((x-s0-self.Num_Bin)/w0)**2;
809
810 e0d=((x-s0+self.Num_Bin)/w0)**2
811 m0=numpy.exp(-0.5*e0**(p0/2.));
812 m0u=numpy.exp(-0.5*e0u**(p0/2.));
813 m0d=numpy.exp(-0.5*e0d**(p0/2.))
814 JA=m0+m0u+m0d
815 JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
537 if Amplitude0<0.05: # in case the peak is noise
538 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
539 if Amplitude1<0.05:
540 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
816 541
817 JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
818 542
819 JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
820 jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model])
821 return jack1.T
822
823 def y_jacobian2(self,x,state):
824 y_model=self.y_model2(x,state)
825 s0,w0,a0,p0,s1,w1,a1,p1,n=state
826 e0=((x-s0)/w0)**2;
827
828 e0u=((x-s0- self.Num_Bin )/w0)**2;
829
830 e0d=((x-s0+ self.Num_Bin )/w0)**2
831 e1=((x-s1)/w1)**2;
832
833 e1u=((x-s1- self.Num_Bin )/w1)**2;
834
835 e1d=((x-s1+ self.Num_Bin )/w1)**2
836 m0=numpy.exp(-0.5*e0**(p0/2.));
837 m0u=numpy.exp(-0.5*e0u**(p0/2.));
838 m0d=numpy.exp(-0.5*e0d**(p0/2.))
839 m1=numpy.exp(-0.5*e1**(p1/2.));
840 m1u=numpy.exp(-0.5*e1u**(p1/2.));
841 m1d=numpy.exp(-0.5*e1d**(p1/2.))
842 JA=m0+m0u+m0d
843 JA1=m1+m1u+m1d
844 JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
845 JP1=(-1/4.)*a1*m1*e1**(p1/2.)*numpy.log(e1)+(-1/4.)*a1*m1u*e1u**(p1/2.)*numpy.log(e1u)+(-1/4.)*a1*m1d*e1d**(p1/2.)*numpy.log(e1d)
846
847 JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
848
849 JS1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)
543 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
544 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
545 SPCparam = (SPC_ch1,SPC_ch2)
850 546
851 JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
852 547
853 JW1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)**2+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)**2+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)**2
854 jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model])
855 return jack2.T
548 return GauSPC
856 549
857 550 def y_model1(self,x,state):
858 551 shift0,width0,amplitude0,power0,noise=state
@@ -885,6 +578,7 class GaussianFit(Operation):
885 578 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
886 579
887 580
581
888 582 class PrecipitationProc(Operation):
889 583
890 584 '''
@@ -901,22 +595,59 class PrecipitationProc(Operation):
901 595 Parameters affected:
902 596 '''
903 597
598 def __init__(self):
599 Operation.__init__(self)
600 self.i=0
601
602
603 def gaus(self,xSamples,Amp,Mu,Sigma):
604 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
904 605
905 def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None,
906 tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None):
606
607
608 def Moments(self, ySamples, xSamples):
609 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
610 yNorm = ySamples / Pot
611
612 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
613 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
614 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
615
616 return numpy.array([Pot, Vr, Desv])
617
618 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
619 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
620
621
622 Velrange = dataOut.spcparam_range[2]
623 FrecRange = dataOut.spcparam_range[0]
624
625 dV= Velrange[1]-Velrange[0]
626 dF= FrecRange[1]-FrecRange[0]
627
628 if radar == "MIRA35C" :
907 629
908 630 self.spc = dataOut.data_pre[0].copy()
909 631 self.Num_Hei = self.spc.shape[2]
910 632 self.Num_Bin = self.spc.shape[1]
911 633 self.Num_Chn = self.spc.shape[0]
634 Ze = self.dBZeMODE2(dataOut)
635
636 else:
912 637
913 Velrange = dataOut.abscissaList
638 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
914 639
915 if radar == "MIRA35C" :
640 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
916 641
917 Ze = self.dBZeMODE2(dataOut)
642 self.spc[:,:,0:7]= numpy.NaN
918 643
919 else:
644 """##########################################"""
645
646 self.Num_Hei = self.spc.shape[2]
647 self.Num_Bin = self.spc.shape[1]
648 self.Num_Chn = self.spc.shape[0]
649
650 ''' Se obtiene la constante del RADAR '''
920 651
921 652 self.Pt = Pt
922 653 self.Gt = Gt
@@ -927,47 +658,100 class PrecipitationProc(Operation):
927 658 self.ThetaT = ThetaT
928 659 self.ThetaR = ThetaR
929 660
930 RadarConstant = GetRadarConstant()
931 SPCmean = numpy.mean(self.spc,0)
932 ETA = numpy.zeros(self.Num_Hei)
933 Pr = numpy.sum(SPCmean,0)
661 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
662 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
663 RadarConstant = 5e-26 * Numerator / Denominator #
664
665 ''' ============================= '''
666
667 self.spc[0] = (self.spc[0]-dataOut.noise[0])
668 self.spc[1] = (self.spc[1]-dataOut.noise[1])
669 self.spc[2] = (self.spc[2]-dataOut.noise[2])
934 670
935 #for R in range(self.Num_Hei):
936 # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
671 self.spc[ numpy.where(self.spc < 0)] = 0
937 672
938 D_range = numpy.zeros(self.Num_Hei)
939 EqSec = numpy.zeros(self.Num_Hei)
673 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
674 SPCmean[ numpy.where(SPCmean < 0)] = 0
675
676 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
677 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
678 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
679
680 Pr = SPCmean[:,:]
681
682 VelMeteoro = numpy.mean(SPCmean,axis=0)
683
684 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
685 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
686 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
687 V_mean = numpy.zeros(self.Num_Hei)
940 688 del_V = numpy.zeros(self.Num_Hei)
689 Z = numpy.zeros(self.Num_Hei)
690 Ze = numpy.zeros(self.Num_Hei)
691 RR = numpy.zeros(self.Num_Hei)
692
693 Range = dataOut.heightList*1000.
941 694
942 695 for R in range(self.Num_Hei):
943 ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
944 696
945 h = R + Altitude #Range from ground to radar pulse altitude
697 h = Range[R] + Altitude #Range from ground to radar pulse altitude
946 698 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
947 699
948 D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity
949 SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma)
700 D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3
701
702 '''NOTA: ETA(n) dn = ETA(f) df
703
704 dn = 1 Diferencial de muestreo
705 df = ETA(n) / ETA(f)
950 706
951 N_dist[R] = ETA[R] / SIGMA[R]
707 '''
708
709 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
710
711 ETAv[:,R]=ETAn[:,R]/dV
712
713 ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
714
715 SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
716
717 N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
718
719 DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin])
720
721 try:
722 popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments)
723 except:
724 popt01=numpy.zeros(3)
725 popt01[1]= DMoments[1]
726
727 if popt01[1]<0 or popt01[1]>20:
728 popt01[1]=numpy.NaN
729
730
731 V_mean[R]=popt01[1]
732
733 Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18
734
735 RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
736
737 Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km)
952 738
953 Ze = (ETA * Lambda**4) / (numpy.pi * Km)
954 Z = numpy.sum( N_dist * D_range**6 )
955 RR = 6*10**-4*numpy.pi * numpy.sum( D_range**3 * N_dist * Velrange ) #Rainfall rate
956 739
957 740
958 RR = (Ze/200)**(1/1.6)
741 RR2 = (Z/200)**(1/1.6)
959 742 dBRR = 10*numpy.log10(RR)
743 dBRR2 = 10*numpy.log10(RR2)
960 744
961 745 dBZe = 10*numpy.log10(Ze)
962 dataOut.data_output = Ze
963 dataOut.data_param = numpy.ones([2,self.Num_Hei])
964 dataOut.channelList = [0,1]
965 print('channelList', dataOut.channelList)
966 dataOut.data_param[0]=dBZe
967 dataOut.data_param[1]=dBRR
968 print('RR SHAPE', dBRR.shape)
969 print('Ze SHAPE', dBZe.shape)
970 print('dataOut.data_param SHAPE', dataOut.data_param.shape)
746 dBZ = 10*numpy.log10(Z)
747
748 dataOut.data_output = RR[8]
749 dataOut.data_param = numpy.ones([3,self.Num_Hei])
750 dataOut.channelList = [0,1,2]
751
752 dataOut.data_param[0]=dBZ
753 dataOut.data_param[1]=V_mean
754 dataOut.data_param[2]=RR
971 755
972 756
973 757 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
@@ -983,7 +767,7 class PrecipitationProc(Operation):
983 767 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
984 768
985 769 ETA = numpy.sum(SNR,1)
986 print('ETA' , ETA)
770
987 771 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
988 772
989 773 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
@@ -995,26 +779,27 class PrecipitationProc(Operation):
995 779
996 780 return Ze
997 781
998 def GetRadarConstant(self):
999
1000 """
1001 Constants:
1002
1003 Pt: Transmission Power dB
1004 Gt: Transmission Gain dB
1005 Gr: Reception Gain dB
1006 Lambda: Wavelenght m
1007 aL: Attenuation loses dB
1008 tauW: Width of transmission pulse s
1009 ThetaT: Transmission antenna bean angle rad
1010 ThetaR: Reception antenna beam angle rad
1011
1012 """
1013 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
1014 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
1015 RadarConstant = Numerator / Denominator
1016
1017 return RadarConstant
782 # def GetRadarConstant(self):
783 #
784 # """
785 # Constants:
786 #
787 # Pt: Transmission Power dB 5kW 5000
788 # Gt: Transmission Gain dB 24.7 dB 295.1209
789 # Gr: Reception Gain dB 18.5 dB 70.7945
790 # Lambda: Wavelenght m 0.6741 m 0.6741
791 # aL: Attenuation loses dB 4dB 2.5118
792 # tauW: Width of transmission pulse s 4us 4e-6
793 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317
794 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087
795 #
796 # """
797 #
798 # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
799 # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
800 # RadarConstant = Numerator / Denominator
801 #
802 # return RadarConstant
1018 803
1019 804
1020 805
@@ -1037,10 +822,20 class FullSpectralAnalysis(Operation):
1037 822 Parameters affected: Winds, height range, SNR
1038 823
1039 824 """
1040 def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7):
825 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7):
826
827 self.indice=int(numpy.random.rand()*1000)
1041 828
1042 829 spc = dataOut.data_pre[0].copy()
1043 cspc = dataOut.data_pre[1].copy()
830 cspc = dataOut.data_pre[1]
831
832 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
833
834 SNRspc = spc.copy()
835 SNRspc[:,:,0:7]= numpy.NaN
836
837 """##########################################"""
838
1044 839
1045 840 nChannel = spc.shape[0]
1046 841 nProfiles = spc.shape[1]
@@ -1050,14 +845,9 class FullSpectralAnalysis(Operation):
1050 845 if dataOut.ChanDist is not None :
1051 846 ChanDist = dataOut.ChanDist
1052 847 else:
1053 ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]])
1054
1055 #print 'ChanDist', ChanDist
848 ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
1056 849
1057 if dataOut.VelRange is not None:
1058 VelRange= dataOut.VelRange
1059 else:
1060 VelRange= dataOut.abscissaList
850 FrecRange = dataOut.spc_range[0]
1061 851
1062 852 ySamples=numpy.ones([nChannel,nProfiles])
1063 853 phase=numpy.ones([nChannel,nProfiles])
@@ -1065,113 +855,108 class FullSpectralAnalysis(Operation):
1065 855 coherence=numpy.ones([nChannel,nProfiles])
1066 856 PhaseSlope=numpy.ones(nChannel)
1067 857 PhaseInter=numpy.ones(nChannel)
1068 dataSNR = dataOut.data_SNR
1069
1070
858 data_SNR=numpy.zeros([nProfiles])
1071 859
1072 860 data = dataOut.data_pre
1073 861 noise = dataOut.noise
1074 print('noise',noise)
1075 #SNRdB = 10*numpy.log10(dataOut.data_SNR)
1076 862
1077 FirstMoment = numpy.average(dataOut.data_param[:,1,:],0)
1078 #SNRdBMean = []
863 dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
1079 864
865 dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20
1080 866
1081 #for j in range(nHeights):
1082 # FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]]))
1083 # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
1084 867
1085 data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN
868 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
1086 869
1087 870 velocityX=[]
1088 871 velocityY=[]
1089 872 velocityV=[]
873 PhaseLine=[]
1090 874
1091 dbSNR = 10*numpy.log10(dataSNR)
875 dbSNR = 10*numpy.log10(dataOut.data_SNR)
1092 876 dbSNR = numpy.average(dbSNR,0)
877
1093 878 for Height in range(nHeights):
1094 879
1095 [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit)
880 [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range.copy(), dbSNR[Height], SNRlimit)
881 PhaseLine = numpy.append(PhaseLine, PhaseSlope)
1096 882
1097 883 if abs(Vzon)<100. and abs(Vzon)> 0.:
1098 884 velocityX=numpy.append(velocityX, Vzon)#Vmag
1099 885
1100 886 else:
1101 print('Vzon',Vzon)
1102 887 velocityX=numpy.append(velocityX, numpy.NaN)
1103 888
1104 889 if abs(Vmer)<100. and abs(Vmer) > 0.:
1105 velocityY=numpy.append(velocityY, Vmer)#Vang
890 velocityY=numpy.append(velocityY, -Vmer)#Vang
1106 891
1107 892 else:
1108 print('Vmer',Vmer)
1109 893 velocityY=numpy.append(velocityY, numpy.NaN)
1110 894
1111 895 if dbSNR[Height] > SNRlimit:
1112 velocityV=numpy.append(velocityV, FirstMoment[Height])
896 velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height])
1113 897 else:
1114 898 velocityV=numpy.append(velocityV, numpy.NaN)
1115 #FirstMoment[Height]= numpy.NaN
1116 # if SNRdBMean[Height] <12:
1117 # FirstMoment[Height] = numpy.NaN
1118 # velocityX[Height] = numpy.NaN
1119 # velocityY[Height] = numpy.NaN
1120
1121
1122 data_output[0]=numpy.array(velocityX)
1123 data_output[1]=numpy.array(velocityY)
1124 data_output[2]=-velocityV#FirstMoment
1125
1126 print(' ')
1127 #print 'FirstMoment'
1128 #print FirstMoment
1129 print('velocityX',data_output[0])
1130 print(' ')
1131 print('velocityY',data_output[1])
1132 #print numpy.array(velocityY)
1133 print(' ')
1134 #print 'SNR'
1135 #print 10*numpy.log10(dataOut.data_SNR)
1136 #print numpy.shape(10*numpy.log10(dataOut.data_SNR))
1137 print(' ')
1138 899
1139 900
901
902 '''Nota: Cambiar el signo de numpy.array(velocityX) cuando se intente procesar datos de BLTR'''
903 data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1)
904 data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1)
905 data_output[2] = velocityV#FirstMoment
906
907 xFrec=FrecRange[0:spc.shape[1]]
908
1140 909 dataOut.data_output=data_output
910
1141 911 return
1142 912
1143 913
1144 914 def moving_average(self,x, N=2):
1145 915 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
1146 916
1147 def gaus(self,xSamples,a,x0,sigma):
1148 return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2))
917 def gaus(self,xSamples,Amp,Mu,Sigma):
918 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
919
920
921
922 def Moments(self, ySamples, xSamples):
923 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
924 yNorm = ySamples / Pot
925 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
926 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
927 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
928
929 return numpy.array([Pot, Vr, Desv])
930
931 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
1149 932
1150 def Find(self,x,value):
1151 for index in range(len(x)):
1152 if x[index]==value:
1153 return index
1154 933
1155 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit):
1156 934
1157 935 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
1158 936 phase=numpy.ones([spc.shape[0],spc.shape[1]])
1159 937 CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
1160 938 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
1161 PhaseSlope=numpy.ones(spc.shape[0])
939 PhaseSlope=numpy.zeros(spc.shape[0])
1162 940 PhaseInter=numpy.ones(spc.shape[0])
1163 xFrec=VelRange
941 xFrec=AbbsisaRange[0][0:spc.shape[1]]
942 xVel =AbbsisaRange[2][0:spc.shape[1]]
943 Vv=numpy.empty(spc.shape[2])*0
944 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]#
945
946 SPCmoments = self.Moments(SPCav[:,Height], xVel )
947 CSPCmoments = []
948 cspcNoise = numpy.empty(3)
1164 949
1165 950 '''Getting Eij and Nij'''
1166 951
1167 E01=ChanDist[0][0]
1168 N01=ChanDist[0][1]
952 Xi01=ChanDist[0][0]
953 Eta01=ChanDist[0][1]
1169 954
1170 E02=ChanDist[1][0]
1171 N02=ChanDist[1][1]
955 Xi02=ChanDist[1][0]
956 Eta02=ChanDist[1][1]
1172 957
1173 E12=ChanDist[2][0]
1174 N12=ChanDist[2][1]
958 Xi12=ChanDist[2][0]
959 Eta12=ChanDist[2][1]
1175 960
1176 961 z = spc.copy()
1177 962 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
@@ -1179,87 +964,94 class FullSpectralAnalysis(Operation):
1179 964 for i in range(spc.shape[0]):
1180 965
1181 966 '''****** Line of Data SPC ******'''
1182 zline=z[i,:,Height]
967 zline=z[i,:,Height].copy() - noise[i] # Se resta ruido
1183 968
1184 969 '''****** SPC is normalized ******'''
1185 FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy())
1186 FactNorm= FactNorm/numpy.sum(FactNorm)
970 SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido
971 FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado
1187 972
1188 SmoothSPC=self.moving_average(FactNorm,N=3)
1189
1190 xSamples = ar(list(range(len(SmoothSPC))))
1191 ySamples[i] = SmoothSPC
1192
1193 #dbSNR=10*numpy.log10(dataSNR)
1194 print(' ')
1195 print(' ')
1196 print(' ')
1197
1198 #print 'dataSNR', dbSNR.shape, dbSNR[0,40:120]
1199 print('SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20])
1200 print('noise',noise)
1201 print('zline',zline.shape, zline[0:20])
1202 print('FactNorm',FactNorm.shape, FactNorm[0:20])
1203 print('FactNorm suma', numpy.sum(FactNorm))
973 xSamples = xFrec # Se toma el rango de frecuncias
974 ySamples[i] = FactNorm # Se toman los valores de SPC normalizado
1204 975
1205 976 for i in range(spc.shape[0]):
1206 977
1207 978 '''****** Line of Data CSPC ******'''
1208 cspcLine=cspc[i,:,Height].copy()
979 cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido
980 SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido
981 cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado
1209 982
1210 '''****** CSPC is normalized ******'''
983 '''****** CSPC is normalized with respect to Briggs and Vincent ******'''
1211 984 chan_index0 = pairsList[i][0]
1212 985 chan_index1 = pairsList[i][1]
1213 CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) #
1214 986
1215 CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor)
987 CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2
988 CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor)
1216 989
1217 990 CSPCSamples[i] = CSPCNorm
991
1218 992 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
1219 993
1220 coherence[i]= self.moving_average(coherence[i],N=2)
994 #coherence[i]= self.moving_average(coherence[i],N=1)
1221 995
1222 996 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
1223 997
1224 print('cspcLine', cspcLine.shape, cspcLine[0:20])
1225 print('CSPCFactor', CSPCFactor)#, CSPCFactor[0:20]
1226 print(numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i])
1227 print('CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20])
1228 print('CSPCNorm suma', numpy.sum(CSPCNorm))
1229 print('CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20])
998 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples),
999 self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
1000 self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
1230 1001
1231 '''****** Getting fij width ******'''
1232 1002
1233 yMean=[]
1234 yMean2=[]
1003 popt=[1e-10,0,1e-10]
1004 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1005 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1006
1007 CSPCMask01 = numpy.abs(CSPCSamples[0])
1008 CSPCMask02 = numpy.abs(CSPCSamples[1])
1009 CSPCMask12 = numpy.abs(CSPCSamples[2])
1235 1010
1236 for j in range(len(ySamples[1])):
1237 yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
1011 mask01 = ~numpy.isnan(CSPCMask01)
1012 mask02 = ~numpy.isnan(CSPCMask02)
1013 mask12 = ~numpy.isnan(CSPCMask12)
1238 1014
1239 '''******* Getting fitting Gaussian ******'''
1240 meanGauss=sum(xSamples*yMean) / len(xSamples)
1241 sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
1015 #mask = ~numpy.isnan(CSPCMask01)
1016 CSPCMask01 = CSPCMask01[mask01]
1017 CSPCMask02 = CSPCMask02[mask02]
1018 CSPCMask12 = CSPCMask12[mask12]
1019 #CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01)
1242 1020
1243 print('****************************')
1244 print('len(xSamples): ',len(xSamples))
1245 print('yMean: ', yMean.shape, yMean[0:20])
1246 print('ySamples', ySamples.shape, ySamples[0,0:20])
1247 print('xSamples: ',xSamples.shape, xSamples[0:20])
1248 1021
1249 print('meanGauss',meanGauss)
1250 print('sigma',sigma)
1251 1022
1252 #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001):
1253 if dbSNR > SNRlimit :
1023 '''***Fit Gauss CSPC01***'''
1024 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 :
1254 1025 try:
1255 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
1026 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1027 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1028 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1029 FitGauss01 = self.gaus(xSamples,*popt01)
1030 FitGauss02 = self.gaus(xSamples,*popt02)
1031 FitGauss12 = self.gaus(xSamples,*popt12)
1032 except:
1033 FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0]))
1034 FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1]))
1035 FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2]))
1036
1037
1038 CSPCopt = numpy.vstack([popt01,popt02,popt12])
1039
1040 '''****** Getting fij width ******'''
1256 1041
1257 if numpy.amax(popt)>numpy.amax(yMean)*0.3:
1042 yMean = numpy.average(ySamples, axis=0) # ySamples[0]
1043
1044 '''******* Getting fitting Gaussian *******'''
1045 meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia)
1046 sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia)
1047
1048 yMoments = self.Moments(yMean, xSamples)
1049
1050 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001:
1051 try:
1052 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
1258 1053 FitGauss=self.gaus(xSamples,*popt)
1259 1054
1260 else:
1261 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1262 print('Verificador: Dentro', Height)
1263 1055 except :#RuntimeError:
1264 1056 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1265 1057
@@ -1267,55 +1059,69 class FullSpectralAnalysis(Operation):
1267 1059 else:
1268 1060 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1269 1061
1270 Maximun=numpy.amax(yMean)
1271 eMinus1=Maximun*numpy.exp(-1)#*0.8
1272 1062
1273 HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
1274 HalfWidth= xFrec[HWpos]
1275 GCpos=self.Find(FitGauss, numpy.amax(FitGauss))
1276 Vpos=self.Find(FactNorm, numpy.amax(FactNorm))
1277
1278 #Vpos=FirstMoment[]
1279 1063
1280 1064 '''****** Getting Fij ******'''
1065 Fijcspc = CSPCopt[:,2]/2*3
1281 1066
1282 GaussCenter=xFrec[GCpos]
1283 if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
1284 Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
1285 else:
1286 Fij=abs(GaussCenter-HalfWidth)+0.0000001
1287 1067
1288 '''****** Getting Frecuency range of significant data ******'''
1068 GaussCenter = popt[1] #xFrec[GCpos]
1069 #Punto en Eje X de la Gaussiana donde se encuentra el centro
1070 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1071 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1289 1072
1290 Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
1073 #Punto e^-1 hubicado en la Gaussiana
1074 PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1)
1075 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1076 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1077
1078 if xSamples[PointFij] > xSamples[PointGauCenter]:
1079 Fij = xSamples[PointFij] - xSamples[PointGauCenter]
1291 1080
1292 if Rangpos<GCpos:
1293 Range=numpy.array([Rangpos,2*GCpos-Rangpos])
1294 elif Rangpos< ( len(xFrec)- len(xFrec)*0.1):
1295 Range=numpy.array([2*GCpos-Rangpos,Rangpos])
1296 1081 else:
1297 Range = numpy.array([0,0])
1082 Fij = xSamples[PointGauCenter] - xSamples[PointFij]
1083
1084
1085 '''****** Taking frequency ranges from SPCs ******'''
1086
1087
1088 #GaussCenter = popt[1] #Primer momento 01
1089 GauWidth = popt[2] *3/2 #Ancho de banda de Gau01
1090 Range = numpy.empty(2)
1091 Range[0] = GaussCenter - GauWidth
1092 Range[1] = GaussCenter + GauWidth
1093 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max)
1094 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1095 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1096
1097 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1098 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1099
1100 Range=numpy.array([ PointRangeMin, PointRangeMax ])
1298 1101
1299 print(' ')
1300 print('GCpos',GCpos, ( len(xFrec)- len(xFrec)*0.1))
1301 print('Rangpos',Rangpos)
1302 print('RANGE: ', Range)
1303 1102 FrecRange = xFrec[ Range[0] : Range[1] ]
1103 VelRange = xVel[ Range[0] : Range[1] ]
1104
1304 1105
1305 1106 '''****** Getting SCPC Slope ******'''
1306 1107
1307 1108 for i in range(spc.shape[0]):
1308 1109
1309 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.5:
1110 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.3:
1310 1111 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1311 1112
1312 print('FrecRange', len(FrecRange) , FrecRange)
1313 print('PhaseRange', len(PhaseRange), PhaseRange)
1314 print(' ')
1113 '''***********************VelRange******************'''
1114
1115 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1116
1315 1117 if len(FrecRange) == len(PhaseRange):
1316 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
1118 try:
1119 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask])
1317 1120 PhaseSlope[i]=slope
1318 1121 PhaseInter[i]=intercept
1122 except:
1123 PhaseSlope[i]=0
1124 PhaseInter[i]=0
1319 1125 else:
1320 1126 PhaseSlope[i]=0
1321 1127 PhaseInter[i]=0
@@ -1323,32 +1129,32 class FullSpectralAnalysis(Operation):
1323 1129 PhaseSlope[i]=0
1324 1130 PhaseInter[i]=0
1325 1131
1132
1326 1133 '''Getting constant C'''
1327 1134 cC=(Fij*numpy.pi)**2
1328 1135
1329 1136 '''****** Getting constants F and G ******'''
1330 MijEijNij=numpy.array([[E02,N02], [E12,N12]])
1137 MijEijNij=numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1331 1138 MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
1332 1139 MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
1333 1140 MijResults=numpy.array([MijResult0,MijResult1])
1334 1141 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1335 1142
1336 1143 '''****** Getting constants A, B and H ******'''
1337 W01=numpy.amax(coherence[0])
1338 W02=numpy.amax(coherence[1])
1339 W12=numpy.amax(coherence[2])
1144 W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0]))
1145 W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1]))
1146 W12=numpy.nanmax( FitGauss12 ) #numpy.abs(CSPCSamples[2]))
1340 1147
1341 WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1342 WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1343 WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1148 WijResult0=((cF*Xi01+cG*Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1149 WijResult1=((cF*Xi02+cG*Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1150 WijResult2=((cF*Xi12+cG*Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1344 1151
1345 1152 WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1346 1153
1347 WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
1154 WijEijNij=numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1348 1155 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1349 1156
1350 1157 VxVy=numpy.array([[cA,cH],[cH,cB]])
1351
1352 1158 VxVyResults=numpy.array([-cF,-cG])
1353 1159 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1354 1160
@@ -1356,9 +1162,14 class FullSpectralAnalysis(Operation):
1356 1162 Vmer = Vx
1357 1163 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1358 1164 Vang=numpy.arctan2(Vmer,Vzon)
1359 Vver=xFrec[Vpos]
1360 print('vzon y vmer', Vzon, Vmer)
1361 return Vzon, Vmer, Vver, GaussCenter
1165 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange)>4:
1166 Vver=popt[1]
1167 else:
1168 Vver=numpy.NaN
1169 FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
1170
1171
1172 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1362 1173
1363 1174 class SpectralMoments(Operation):
1364 1175
@@ -1384,7 +1195,7 class SpectralMoments(Operation):
1384 1195 self.dataOut.noise : Noise level per channel
1385 1196
1386 1197 Affected:
1387 self.dataOut.data_param : Parameters per channel
1198 self.dataOut.moments : Parameters per channel
1388 1199 self.dataOut.data_SNR : SNR per channel
1389 1200
1390 1201 '''
@@ -1401,7 +1212,7 class SpectralMoments(Operation):
1401 1212 for ind in range(nChannel):
1402 1213 data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
1403 1214
1404 dataOut.data_param = data_param[:,1:,:]
1215 dataOut.moments = data_param[:,1:,:]
1405 1216 dataOut.data_SNR = data_param[:,0]
1406 1217 dataOut.data_DOP = data_param[:,1]
1407 1218 dataOut.data_MEAN = data_param[:,2]
@@ -1432,6 +1243,8 class SpectralMoments(Operation):
1432 1243 vec_w = numpy.zeros(oldspec.shape[1])
1433 1244 vec_snr = numpy.zeros(oldspec.shape[1])
1434 1245
1246 oldspec = numpy.ma.masked_invalid(oldspec)
1247
1435 1248 for ind in range(oldspec.shape[1]):
1436 1249
1437 1250 spec = oldspec[:,ind]
@@ -1675,7 +1488,6 class SpectralFitting(Operation):
1675 1488 dataCross = dataCross**2/K
1676 1489
1677 1490 for h in range(nHeights):
1678 # print self.dataOut.heightList[h]
1679 1491
1680 1492 #Input
1681 1493 d = data[:,h]
@@ -1768,8 +1580,8 class WindProfiler(Operation):
1768 1580
1769 1581 n = None
1770 1582
1771 def __init__(self, **kwargs):
1772 Operation.__init__(self, **kwargs)
1583 def __init__(self):
1584 Operation.__init__(self)
1773 1585
1774 1586 def __calculateCosDir(self, elev, azim):
1775 1587 zen = (90 - elev)*numpy.pi/180
@@ -2071,12 +1883,9 class WindProfiler(Operation):
2071 1883
2072 1884 Parameters affected: Winds
2073 1885 '''
2074 # print arrayMeteor.shape
2075 1886 #Settings
2076 1887 nInt = (heightMax - heightMin)/2
2077 # print nInt
2078 1888 nInt = int(nInt)
2079 # print nInt
2080 1889 winds = numpy.zeros((2,nInt))*numpy.nan
2081 1890
2082 1891 #Filter errors
@@ -2475,8 +2284,8 class WindProfiler(Operation):
2475 2284
2476 2285 class EWDriftsEstimation(Operation):
2477 2286
2478 def __init__(self, **kwargs):
2479 Operation.__init__(self, **kwargs)
2287 def __init__(self):
2288 Operation.__init__(self)
2480 2289
2481 2290 def __correctValues(self, heiRang, phi, velRadial, SNR):
2482 2291 listPhi = phi.tolist()
Show file before | Show file after | Hide comments
@@ -159,9 +159,7 class SpectraProc(ProcessingUnit):
159 159 dtype='complex')
160 160
161 161 if self.dataIn.flagDataAsBlock:
162 # data dimension: [nChannels, nProfiles, nSamples]
163 162 nVoltProfiles = self.dataIn.data.shape[1]
164 # nVoltProfiles = self.dataIn.nProfiles
165 163
166 164 if nVoltProfiles == nProfiles:
167 165 self.buffer = self.dataIn.data.copy()
@@ -300,6 +298,56 class SpectraProc(ProcessingUnit):
300 298
301 299 return 1
302 300
301
302 def selectFFTs(self, minFFT, maxFFT ):
303 """
304 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
305 minFFT<= FFT <= maxFFT
306 """
307
308 if (minFFT > maxFFT):
309 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
310
311 if (minFFT < self.dataOut.getFreqRange()[0]):
312 minFFT = self.dataOut.getFreqRange()[0]
313
314 if (maxFFT > self.dataOut.getFreqRange()[-1]):
315 maxFFT = self.dataOut.getFreqRange()[-1]
316
317 minIndex = 0
318 maxIndex = 0
319 FFTs = self.dataOut.getFreqRange()
320
321 inda = numpy.where(FFTs >= minFFT)
322 indb = numpy.where(FFTs <= maxFFT)
323
324 try:
325 minIndex = inda[0][0]
326 except:
327 minIndex = 0
328
329 try:
330 maxIndex = indb[0][-1]
331 except:
332 maxIndex = len(FFTs)
333
334 self.selectFFTsByIndex(minIndex, maxIndex)
335
336 return 1
337
338
339 def setH0(self, h0, deltaHeight = None):
340
341 if not deltaHeight:
342 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
343
344 nHeights = self.dataOut.nHeights
345
346 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
347
348 self.dataOut.heightList = newHeiRange
349
350
303 351 def selectHeights(self, minHei, maxHei):
304 352 """
305 353 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
@@ -316,9 +364,9 class SpectraProc(ProcessingUnit):
316 364 1 si el metodo se ejecuto con exito caso contrario devuelve 0
317 365 """
318 366
367
319 368 if (minHei > maxHei):
320 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (
321 minHei, maxHei))
369 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei))
322 370
323 371 if (minHei < self.dataOut.heightList[0]):
324 372 minHei = self.dataOut.heightList[0]
@@ -345,6 +393,7 class SpectraProc(ProcessingUnit):
345 393
346 394 self.selectHeightsByIndex(minIndex, maxIndex)
347 395
396
348 397 return 1
349 398
350 399 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
@@ -389,6 +438,40 class SpectraProc(ProcessingUnit):
389 438
390 439 return 1
391 440
441 def selectFFTsByIndex(self, minIndex, maxIndex):
442 """
443
444 """
445
446 if (minIndex < 0) or (minIndex > maxIndex):
447 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
448
449 if (maxIndex >= self.dataOut.nProfiles):
450 maxIndex = self.dataOut.nProfiles-1
451
452 #Spectra
453 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
454
455 data_cspc = None
456 if self.dataOut.data_cspc is not None:
457 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
458
459 data_dc = None
460 if self.dataOut.data_dc is not None:
461 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
462
463 self.dataOut.data_spc = data_spc
464 self.dataOut.data_cspc = data_cspc
465 self.dataOut.data_dc = data_dc
466
467 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
468 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
469 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
470
471 return 1
472
473
474
392 475 def selectHeightsByIndex(self, minIndex, maxIndex):
393 476 """
394 477 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
@@ -494,6 +577,31 class SpectraProc(ProcessingUnit):
494 577
495 578 return 1
496 579
580 def removeInterference2(self):
581
582 cspc = self.dataOut.data_cspc
583 spc = self.dataOut.data_spc
584 Heights = numpy.arange(cspc.shape[2])
585 realCspc = numpy.abs(cspc)
586
587 for i in range(cspc.shape[0]):
588 LinePower= numpy.sum(realCspc[i], axis=0)
589 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
590 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
591 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
592 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
593 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
594
595
596 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
597 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
598 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
599 cspc[i,InterferenceRange,:] = numpy.NaN
600
601
602
603 self.dataOut.data_cspc = cspc
604
497 605 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
498 606
499 607 jspectra = self.dataOut.data_spc
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (566 lines changed) Show them Hide them
@@ -7,11 +7,9 import glob
7 7 import time
8 8 import json
9 9 import numpy
10 import paho.mqtt.client as mqtt
11 10 import zmq
12 11 import datetime
13 12 import ftplib
14 from zmq.utils.monitor import recv_monitor_message
15 13 from functools import wraps
16 14 from threading import Thread
17 15 from multiprocessing import Process
@@ -54,330 +52,30 def get_plot_code(s):
54 52 else:
55 53 return 24
56 54
57 def roundFloats(obj):
58 if isinstance(obj, list):
59 return list(map(roundFloats, obj))
60 elif isinstance(obj, float):
61 return round(obj, 2)
62
63 55 def decimate(z, MAXNUMY):
64 56 dy = int(len(z[0])/MAXNUMY) + 1
65 57
66 58 return z[::, ::dy]
67 59
68 class throttle(object):
69 '''
70 Decorator that prevents a function from being called more than once every
71 time period.
72 To create a function that cannot be called more than once a minute, but
73 will sleep until it can be called:
74 @throttle(minutes=1)
75 def foo():
76 pass
77
78 for i in range(10):
79 foo()
80 print "This function has run %s times." % i
81 '''
82
83 def __init__(self, seconds=0, minutes=0, hours=0):
84 self.throttle_period = datetime.timedelta(
85 seconds=seconds, minutes=minutes, hours=hours
86 )
87
88 self.time_of_last_call = datetime.datetime.min
89
90 def __call__(self, fn):
91 @wraps(fn)
92 def wrapper(*args, **kwargs):
93 coerce = kwargs.pop('coerce', None)
94 if coerce:
95 self.time_of_last_call = datetime.datetime.now()
96 return fn(*args, **kwargs)
97 else:
98 now = datetime.datetime.now()
99 time_since_last_call = now - self.time_of_last_call
100 time_left = self.throttle_period - time_since_last_call
101
102 if time_left > datetime.timedelta(seconds=0):
103 return
104
105 self.time_of_last_call = datetime.datetime.now()
106 return fn(*args, **kwargs)
107
108 return wrapper
109
110 class Data(object):
111 '''
112 Object to hold data to be plotted
113 '''
114
115 def __init__(self, plottypes, throttle_value, exp_code, buffering=True):
116 self.plottypes = plottypes
117 self.throttle = throttle_value
118 self.exp_code = exp_code
119 self.buffering = buffering
120 self.ended = False
121 self.localtime = False
122 self.meta = {}
123 self.__times = []
124 self.__heights = []
125
126 def __str__(self):
127 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
128 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
129
130 def __len__(self):
131 return len(self.__times)
132
133 def __getitem__(self, key):
134 if key not in self.data:
135 raise KeyError(log.error('Missing key: {}'.format(key)))
136
137 if 'spc' in key or not self.buffering:
138 ret = self.data[key]
139 else:
140 ret = numpy.array([self.data[key][x] for x in self.times])
141 if ret.ndim > 1:
142 ret = numpy.swapaxes(ret, 0, 1)
143 return ret
144
145 def __contains__(self, key):
146 return key in self.data
147
148 def setup(self):
149 '''
150 Configure object
151 '''
152
153 self.type = ''
154 self.ended = False
155 self.data = {}
156 self.__times = []
157 self.__heights = []
158 self.__all_heights = set()
159 for plot in self.plottypes:
160 if 'snr' in plot:
161 plot = 'snr'
162 self.data[plot] = {}
163
164 def shape(self, key):
165 '''
166 Get the shape of the one-element data for the given key
167 '''
168
169 if len(self.data[key]):
170 if 'spc' in key or not self.buffering:
171 return self.data[key].shape
172 return self.data[key][self.__times[0]].shape
173 return (0,)
174
175 def update(self, dataOut, tm):
176 '''
177 Update data object with new dataOut
178 '''
179
180 if tm in self.__times:
181 return
182
183 self.type = dataOut.type
184 self.parameters = getattr(dataOut, 'parameters', [])
185 if hasattr(dataOut, 'pairsList'):
186 self.pairs = dataOut.pairsList
187 if hasattr(dataOut, 'meta'):
188 self.meta = dataOut.meta
189 self.channels = dataOut.channelList
190 self.interval = dataOut.getTimeInterval()
191 self.localtime = dataOut.useLocalTime
192 if 'spc' in self.plottypes or 'cspc' in self.plottypes:
193 self.xrange = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
194 self.__heights.append(dataOut.heightList)
195 self.__all_heights.update(dataOut.heightList)
196 self.__times.append(tm)
197
198 for plot in self.plottypes:
199 if plot == 'spc':
200 z = dataOut.data_spc/dataOut.normFactor
201 buffer = 10*numpy.log10(z)
202 if plot == 'cspc':
203 buffer = dataOut.data_cspc
204 if plot == 'noise':
205 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
206 if plot == 'rti':
207 buffer = dataOut.getPower()
208 if plot == 'snr_db':
209 buffer = dataOut.data_SNR
210 if plot == 'snr':
211 buffer = 10*numpy.log10(dataOut.data_SNR)
212 if plot == 'dop':
213 buffer = 10*numpy.log10(dataOut.data_DOP)
214 if plot == 'mean':
215 buffer = dataOut.data_MEAN
216 if plot == 'std':
217 buffer = dataOut.data_STD
218 if plot == 'coh':
219 buffer = dataOut.getCoherence()
220 if plot == 'phase':
221 buffer = dataOut.getCoherence(phase=True)
222 if plot == 'output':
223 buffer = dataOut.data_output
224 if plot == 'param':
225 buffer = dataOut.data_param
226
227 if 'spc' in plot:
228 self.data[plot] = buffer
229 else:
230 if self.buffering:
231 self.data[plot][tm] = buffer
232 else:
233 self.data[plot] = buffer
234
235 def normalize_heights(self):
236 '''
237 Ensure same-dimension of the data for different heighList
238 '''
239
240 H = numpy.array(list(self.__all_heights))
241 H.sort()
242 for key in self.data:
243 shape = self.shape(key)[:-1] + H.shape
244 for tm, obj in list(self.data[key].items()):
245 h = self.__heights[self.__times.index(tm)]
246 if H.size == h.size:
247 continue
248 index = numpy.where(numpy.in1d(H, h))[0]
249 dummy = numpy.zeros(shape) + numpy.nan
250 if len(shape) == 2:
251 dummy[:, index] = obj
252 else:
253 dummy[index] = obj
254 self.data[key][tm] = dummy
255
256 self.__heights = [H for tm in self.__times]
257
258 def jsonify(self, decimate=False):
259 '''
260 Convert data to json
261 '''
262
263 data = {}
264 tm = self.times[-1]
265 dy = int(self.heights.size/MAXNUMY) + 1
266 for key in self.data:
267 if key in ('spc', 'cspc') or not self.buffering:
268 dx = int(self.data[key].shape[1]/MAXNUMX) + 1
269 data[key] = roundFloats(self.data[key][::, ::dx, ::dy].tolist())
270 else:
271 data[key] = roundFloats(self.data[key][tm].tolist())
272
273 ret = {'data': data}
274 ret['exp_code'] = self.exp_code
275 ret['time'] = tm
276 ret['interval'] = self.interval
277 ret['localtime'] = self.localtime
278 ret['yrange'] = roundFloats(self.heights[::dy].tolist())
279 if 'spc' in self.data or 'cspc' in self.data:
280 ret['xrange'] = roundFloats(self.xrange[2][::dx].tolist())
281 else:
282 ret['xrange'] = []
283 if hasattr(self, 'pairs'):
284 ret['pairs'] = self.pairs
285 else:
286 ret['pairs'] = []
287
288 for key, value in list(self.meta.items()):
289 ret[key] = value
290
291 return json.dumps(ret)
292
293 @property
294 def times(self):
295 '''
296 Return the list of times of the current data
297 '''
298
299 ret = numpy.array(self.__times)
300 ret.sort()
301 return ret
302
303 @property
304 def heights(self):
305 '''
306 Return the list of heights of the current data
307 '''
308
309 return numpy.array(self.__heights[-1])
310 60
311 61 class PublishData(Operation):
312 62 '''
313 63 Operation to send data over zmq.
314 64 '''
315 65
316 __attrs__ = ['host', 'port', 'delay', 'zeromq', 'mqtt', 'verbose']
66 __attrs__ = ['host', 'port', 'delay', 'verbose']
317 67
318 68 def __init__(self, **kwargs):
319 69 """Inicio."""
320 70 Operation.__init__(self, **kwargs)
321 71 self.isConfig = False
322 self.client = None
323 self.zeromq = None
324 self.mqtt = None
325
326 def on_disconnect(self, client, userdata, rc):
327 if rc != 0:
328 log.warning('Unexpected disconnection.')
329 self.connect()
330 72
331 def connect(self):
332 log.warning('trying to connect')
333 try:
334 self.client.connect(
335 host=self.host,
336 port=self.port,
337 keepalive=60*10,
338 bind_address='')
339 self.client.loop_start()
340 # self.client.publish(
341 # self.topic + 'SETUP',
342 # json.dumps(setup),
343 # retain=True
344 # )
345 except:
346 log.error('MQTT Conection error.')
347 self.client = False
348
349 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, verbose=True, **kwargs):
73 def setup(self, server='zmq.pipe', delay=0, verbose=True, **kwargs):
350 74 self.counter = 0
351 self.topic = kwargs.get('topic', 'schain')
352 75 self.delay = kwargs.get('delay', 0)
353 self.plottype = kwargs.get('plottype', 'spectra')
354 self.host = kwargs.get('host', "10.10.10.82")
355 self.port = kwargs.get('port', 3000)
356 self.clientId = clientId
357 76 self.cnt = 0
358 self.zeromq = zeromq
359 self.mqtt = kwargs.get('plottype', 0)
360 self.client = None
361 77 self.verbose = verbose
362 78 setup = []
363 if mqtt is 1:
364 self.client = mqtt.Client(
365 client_id=self.clientId + self.topic + 'SCHAIN',
366 clean_session=True)
367 self.client.on_disconnect = self.on_disconnect
368 self.connect()
369 for plot in self.plottype:
370 setup.append({
371 'plot': plot,
372 'topic': self.topic + plot,
373 'title': getattr(self, plot + '_' + 'title', False),
374 'xlabel': getattr(self, plot + '_' + 'xlabel', False),
375 'ylabel': getattr(self, plot + '_' + 'ylabel', False),
376 'xrange': getattr(self, plot + '_' + 'xrange', False),
377 'yrange': getattr(self, plot + '_' + 'yrange', False),
378 'zrange': getattr(self, plot + '_' + 'zrange', False),
379 })
380 if zeromq is 1:
381 79 context = zmq.Context()
382 80 self.zmq_socket = context.socket(zmq.PUSH)
383 81 server = kwargs.get('server', 'zmq.pipe')
@@ -393,83 +91,7 class PublishData(Operation):
393 91
394 92 def publish_data(self):
395 93 self.dataOut.finished = False
396 if self.mqtt is 1:
397 yData = self.dataOut.heightList[:2].tolist()
398 if self.plottype == 'spectra':
399 data = getattr(self.dataOut, 'data_spc')
400 z = data/self.dataOut.normFactor
401 zdB = 10*numpy.log10(z)
402 xlen, ylen = zdB[0].shape
403 dx = int(xlen/MAXNUMX) + 1
404 dy = int(ylen/MAXNUMY) + 1
405 Z = [0 for i in self.dataOut.channelList]
406 for i in self.dataOut.channelList:
407 Z[i] = zdB[i][::dx, ::dy].tolist()
408 payload = {
409 'timestamp': self.dataOut.utctime,
410 'data': roundFloats(Z),
411 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
412 'interval': self.dataOut.getTimeInterval(),
413 'type': self.plottype,
414 'yData': yData
415 }
416
417 elif self.plottype in ('rti', 'power'):
418 data = getattr(self.dataOut, 'data_spc')
419 z = data/self.dataOut.normFactor
420 avg = numpy.average(z, axis=1)
421 avgdB = 10*numpy.log10(avg)
422 xlen, ylen = z[0].shape
423 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
424 AVG = [0 for i in self.dataOut.channelList]
425 for i in self.dataOut.channelList:
426 AVG[i] = avgdB[i][::dy].tolist()
427 payload = {
428 'timestamp': self.dataOut.utctime,
429 'data': roundFloats(AVG),
430 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
431 'interval': self.dataOut.getTimeInterval(),
432 'type': self.plottype,
433 'yData': yData
434 }
435 elif self.plottype == 'noise':
436 noise = self.dataOut.getNoise()/self.dataOut.normFactor
437 noisedB = 10*numpy.log10(noise)
438 payload = {
439 'timestamp': self.dataOut.utctime,
440 'data': roundFloats(noisedB.reshape(-1, 1).tolist()),
441 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
442 'interval': self.dataOut.getTimeInterval(),
443 'type': self.plottype,
444 'yData': yData
445 }
446 elif self.plottype == 'snr':
447 data = getattr(self.dataOut, 'data_SNR')
448 avgdB = 10*numpy.log10(data)
449
450 ylen = data[0].size
451 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
452 AVG = [0 for i in self.dataOut.channelList]
453 for i in self.dataOut.channelList:
454 AVG[i] = avgdB[i][::dy].tolist()
455 payload = {
456 'timestamp': self.dataOut.utctime,
457 'data': roundFloats(AVG),
458 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
459 'type': self.plottype,
460 'yData': yData
461 }
462 else:
463 print("Tipo de grafico invalido")
464 payload = {
465 'data': 'None',
466 'timestamp': 'None',
467 'type': None
468 }
469 94
470 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
471
472 if self.zeromq is 1:
473 95 if self.verbose:
474 96 log.log(
475 97 'Sending {} - {}'.format(self.dataOut.type, self.dataOut.datatime),
@@ -487,14 +109,11 class PublishData(Operation):
487 109 time.sleep(self.delay)
488 110
489 111 def close(self):
490 if self.zeromq is 1:
112
491 113 self.dataOut.finished = True
492 114 self.zmq_socket.send_pyobj(self.dataOut)
493 115 time.sleep(0.1)
494 116 self.zmq_socket.close()
495 if self.client:
496 self.client.loop_stop()
497 self.client.disconnect()
498 117
499 118
500 119 class ReceiverData(ProcessingUnit):
@@ -536,185 +155,6 class ReceiverData(ProcessingUnit):
536 155 'Receiving')
537 156
538 157
539 class PlotterReceiver(ProcessingUnit, Process):
540
541 throttle_value = 5
542 __attrs__ = ['server', 'plottypes', 'realtime', 'localtime', 'throttle',
543 'exp_code', 'web_server', 'buffering']
544
545 def __init__(self, **kwargs):
546
547 ProcessingUnit.__init__(self, **kwargs)
548 Process.__init__(self)
549 self.mp = False
550 self.isConfig = False
551 self.isWebConfig = False
552 self.connections = 0
553 server = kwargs.get('server', 'zmq.pipe')
554 web_server = kwargs.get('web_server', None)
555 if 'tcp://' in server:
556 address = server
557 else:
558 address = 'ipc:///tmp/%s' % server
559 self.address = address
560 self.web_address = web_server
561 self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')]
562 self.realtime = kwargs.get('realtime', False)
563 self.localtime = kwargs.get('localtime', True)
564 self.buffering = kwargs.get('buffering', True)
565 self.throttle_value = kwargs.get('throttle', 5)
566 self.exp_code = kwargs.get('exp_code', None)
567 self.sendData = self.initThrottle(self.throttle_value)
568 self.dates = []
569 self.setup()
570
571 def setup(self):
572
573 self.data = Data(self.plottypes, self.throttle_value, self.exp_code, self.buffering)
574 self.isConfig = True
575
576 def event_monitor(self, monitor):
577
578 events = {}
579
580 for name in dir(zmq):
581 if name.startswith('EVENT_'):
582 value = getattr(zmq, name)
583 events[value] = name
584
585 while monitor.poll():
586 evt = recv_monitor_message(monitor)
587 if evt['event'] == 32:
588 self.connections += 1
589 if evt['event'] == 512:
590 pass
591
592 evt.update({'description': events[evt['event']]})
593
594 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
595 break
596 monitor.close()
597 print('event monitor thread done!')
598
599 def initThrottle(self, throttle_value):
600
601 @throttle(seconds=throttle_value)
602 def sendDataThrottled(fn_sender, data):
603 fn_sender(data)
604
605 return sendDataThrottled
606
607 def send(self, data):
608 log.log('Sending {}'.format(data), self.name)
609 self.sender.send_pyobj(data)
610
611 def run(self):
612
613 log.log(
614 'Starting from {}'.format(self.address),
615 self.name
616 )
617
618 self.context = zmq.Context()
619 self.receiver = self.context.socket(zmq.PULL)
620 self.receiver.bind(self.address)
621 monitor = self.receiver.get_monitor_socket()
622 self.sender = self.context.socket(zmq.PUB)
623 if self.web_address:
624 log.success(
625 'Sending to web: {}'.format(self.web_address),
626 self.name
627 )
628 self.sender_web = self.context.socket(zmq.REQ)
629 self.sender_web.connect(self.web_address)
630 self.poll = zmq.Poller()
631 self.poll.register(self.sender_web, zmq.POLLIN)
632 time.sleep(1)
633
634 if 'server' in self.kwargs:
635 self.sender.bind("ipc:///tmp/{}.plots".format(self.kwargs['server']))
636 else:
637 self.sender.bind("ipc:///tmp/zmq.plots")
638
639 time.sleep(2)
640
641 t = Thread(target=self.event_monitor, args=(monitor,))
642 t.start()
643
644 while True:
645 dataOut = self.receiver.recv_pyobj()
646 if not dataOut.flagNoData:
647 if dataOut.type == 'Parameters':
648 tm = dataOut.utctimeInit
649 else:
650 tm = dataOut.utctime
651 if dataOut.useLocalTime:
652 if not self.localtime:
653 tm += time.timezone
654 dt = datetime.datetime.fromtimestamp(tm).date()
655 else:
656 if self.localtime:
657 tm -= time.timezone
658 dt = datetime.datetime.utcfromtimestamp(tm).date()
659 coerce = False
660 if dt not in self.dates:
661 if self.data:
662 self.data.ended = True
663 self.send(self.data)
664 coerce = True
665 self.data.setup()
666 self.dates.append(dt)
667
668 self.data.update(dataOut, tm)
669
670 if dataOut.finished is True:
671 self.connections -= 1
672 if self.connections == 0 and dt in self.dates:
673 self.data.ended = True
674 self.send(self.data)
675 # self.data.setup()
676 time.sleep(1)
677 break
678 else:
679 if self.realtime:
680 self.send(self.data)
681 if self.web_address:
682 retries = 5
683 while True:
684 self.sender_web.send(self.data.jsonify())
685 socks = dict(self.poll.poll(5000))
686 if socks.get(self.sender_web) == zmq.POLLIN:
687 reply = self.sender_web.recv_string()
688 if reply == 'ok':
689 log.log("Response from server ok", self.name)
690 break
691 else:
692 log.warning("Malformed reply from server: {}".format(reply), self.name)
693
694 else:
695 log.warning("No response from server, retrying...", self.name)
696 self.sender_web.setsockopt(zmq.LINGER, 0)
697 self.sender_web.close()
698 self.poll.unregister(self.sender_web)
699 retries -= 1
700 if retries == 0:
701 log.error("Server seems to be offline, abandoning", self.name)
702 self.sender_web = self.context.socket(zmq.REQ)
703 self.sender_web.connect(self.web_address)
704 self.poll.register(self.sender_web, zmq.POLLIN)
705 time.sleep(1)
706 break
707 self.sender_web = self.context.socket(zmq.REQ)
708 self.sender_web.connect(self.web_address)
709 self.poll.register(self.sender_web, zmq.POLLIN)
710 time.sleep(1)
711 else:
712 self.sendData(self.send, self.data, coerce=coerce)
713 coerce = False
714
715 return
716
717
718 158 class SendToFTP(Operation, Process):
719 159
720 160 '''
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
General Comments 0
You need to be logged in to leave comments. Login now