##// END OF EJS Templates
schain
RSS
clean Rayleigh funcionando, tiempo de ejecución elevado usando todos los pares
joabAM -
r1391:fba03565a781
parent child
Show More
Side by Side Unified
Context file:
r1391:fba03565a781 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,693 +1,695
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Base class to create plot operations
6 6
7 7 """
8 8
9 9 import os
10 10 import sys
11 11 import zmq
12 12 import time
13 13 import numpy
14 14 import datetime
15 15 from collections import deque
16 16 from functools import wraps
17 17 from threading import Thread
18 18 import matplotlib
19 19
20 20 if 'BACKEND' in os.environ:
21 21 matplotlib.use(os.environ['BACKEND'])
22 22 elif 'linux' in sys.platform:
23 23 matplotlib.use("TkAgg")
24 24 elif 'darwin' in sys.platform:
25 25 matplotlib.use('MacOSX')
26 26 else:
27 27 from schainpy.utils import log
28 28 log.warning('Using default Backend="Agg"', 'INFO')
29 29 matplotlib.use('Agg')
30 30
31 31 import matplotlib.pyplot as plt
32 32 from matplotlib.patches import Polygon
33 33 from mpl_toolkits.axes_grid1 import make_axes_locatable
34 34 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
35 35
36 36 from schainpy.model.data.jrodata import PlotterData
37 37 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
38 38 from schainpy.utils import log
39 39
40 40 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
41 41 blu_values = matplotlib.pyplot.get_cmap(
42 42 'seismic_r', 20)(numpy.arange(20))[10:15]
43 43 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
44 44 'jro', numpy.vstack((blu_values, jet_values)))
45 45 matplotlib.pyplot.register_cmap(cmap=ncmap)
46 46
47 47 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
48 48 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
49 49
50 50 EARTH_RADIUS = 6.3710e3
51 51
52 52 def ll2xy(lat1, lon1, lat2, lon2):
53 53
54 54 p = 0.017453292519943295
55 55 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
56 56 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
57 57 r = 12742 * numpy.arcsin(numpy.sqrt(a))
58 58 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
59 59 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
60 60 theta = -theta + numpy.pi/2
61 61 return r*numpy.cos(theta), r*numpy.sin(theta)
62 62
63 63
64 64 def km2deg(km):
65 65 '''
66 66 Convert distance in km to degrees
67 67 '''
68 68
69 69 return numpy.rad2deg(km/EARTH_RADIUS)
70 70
71 71
72 72 def figpause(interval):
73 73 backend = plt.rcParams['backend']
74 74 if backend in matplotlib.rcsetup.interactive_bk:
75 75 figManager = matplotlib._pylab_helpers.Gcf.get_active()
76 76 if figManager is not None:
77 77 canvas = figManager.canvas
78 78 if canvas.figure.stale:
79 79 canvas.draw()
80 80 try:
81 81 canvas.start_event_loop(interval)
82 82 except:
83 83 pass
84 84 return
85 85
86 86 def popup(message):
87 87 '''
88 88 '''
89 89
90 90 fig = plt.figure(figsize=(12, 8), facecolor='r')
91 91 text = '\n'.join([s.strip() for s in message.split(':')])
92 92 fig.text(0.01, 0.5, text, ha='left', va='center',
93 93 size='20', weight='heavy', color='w')
94 94 fig.show()
95 95 figpause(1000)
96 96
97 97
98 98 class Throttle(object):
99 99 '''
100 100 Decorator that prevents a function from being called more than once every
101 101 time period.
102 102 To create a function that cannot be called more than once a minute, but
103 103 will sleep until it can be called:
104 104 @Throttle(minutes=1)
105 105 def foo():
106 106 pass
107 107
108 108 for i in range(10):
109 109 foo()
110 110 print "This function has run %s times." % i
111 111 '''
112 112
113 113 def __init__(self, seconds=0, minutes=0, hours=0):
114 114 self.throttle_period = datetime.timedelta(
115 115 seconds=seconds, minutes=minutes, hours=hours
116 116 )
117 117
118 118 self.time_of_last_call = datetime.datetime.min
119 119
120 120 def __call__(self, fn):
121 121 @wraps(fn)
122 122 def wrapper(*args, **kwargs):
123 123 coerce = kwargs.pop('coerce', None)
124 124 if coerce:
125 125 self.time_of_last_call = datetime.datetime.now()
126 126 return fn(*args, **kwargs)
127 127 else:
128 128 now = datetime.datetime.now()
129 129 time_since_last_call = now - self.time_of_last_call
130 130 time_left = self.throttle_period - time_since_last_call
131 131
132 132 if time_left > datetime.timedelta(seconds=0):
133 133 return
134 134
135 135 self.time_of_last_call = datetime.datetime.now()
136 136 return fn(*args, **kwargs)
137 137
138 138 return wrapper
139 139
140 140 def apply_throttle(value):
141 141
142 142 @Throttle(seconds=value)
143 143 def fnThrottled(fn):
144 144 fn()
145 145
146 146 return fnThrottled
147 147
148 148
149 149 @MPDecorator
150 150 class Plot(Operation):
151 151 """Base class for Schain plotting operations
152 152
153 153 This class should never be use directtly you must subclass a new operation,
154 154 children classes must be defined as follow:
155 155
156 156 ExamplePlot(Plot):
157 157
158 158 CODE = 'code'
159 159 colormap = 'jet'
160 160 plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
161 161
162 162 def setup(self):
163 163 pass
164 164
165 165 def plot(self):
166 166 pass
167 167
168 168 """
169 169
170 170 CODE = 'Figure'
171 171 colormap = 'jet'
172 172 bgcolor = 'white'
173 173 buffering = True
174 174 __missing = 1E30
175 175
176 176 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
177 177 'showprofile']
178 178
179 179 def __init__(self):
180 180
181 181 Operation.__init__(self)
182 182 self.isConfig = False
183 183 self.isPlotConfig = False
184 184 self.save_time = 0
185 185 self.sender_time = 0
186 186 self.data = None
187 187 self.firsttime = True
188 188 self.sender_queue = deque(maxlen=10)
189 189 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
190 190
191 191 def __fmtTime(self, x, pos):
192 192 '''
193 193 '''
194 194
195 195 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
196 196
197 197 def __setup(self, **kwargs):
198 198 '''
199 199 Initialize variables
200 200 '''
201 201
202 202 self.figures = []
203 203 self.axes = []
204 204 self.cb_axes = []
205 self.pf_axes = []
205 206 self.localtime = kwargs.pop('localtime', True)
206 207 self.show = kwargs.get('show', True)
207 208 self.save = kwargs.get('save', False)
208 209 self.save_period = kwargs.get('save_period', 0)
209 210 self.colormap = kwargs.get('colormap', self.colormap)
210 211 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
211 212 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
212 213 self.colormaps = kwargs.get('colormaps', None)
213 214 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
214 215 self.showprofile = kwargs.get('showprofile', False)
215 216 self.title = kwargs.get('wintitle', self.CODE.upper())
216 217 self.cb_label = kwargs.get('cb_label', None)
217 218 self.cb_labels = kwargs.get('cb_labels', None)
218 219 self.labels = kwargs.get('labels', None)
219 220 self.xaxis = kwargs.get('xaxis', 'frequency')
220 221 self.zmin = kwargs.get('zmin', None)
221 222 self.zmax = kwargs.get('zmax', None)
222 223 self.zlimits = kwargs.get('zlimits', None)
223 224 self.xmin = kwargs.get('xmin', None)
224 225 self.xmax = kwargs.get('xmax', None)
225 226 self.xrange = kwargs.get('xrange', 12)
226 227 self.xscale = kwargs.get('xscale', None)
227 228 self.ymin = kwargs.get('ymin', None)
228 229 self.ymax = kwargs.get('ymax', None)
229 230 self.yscale = kwargs.get('yscale', None)
230 231 self.xlabel = kwargs.get('xlabel', None)
231 232 self.attr_time = kwargs.get('attr_time', 'utctime')
232 233 self.attr_data = kwargs.get('attr_data', 'data_param')
233 234 self.decimation = kwargs.get('decimation', None)
234 235 self.oneFigure = kwargs.get('oneFigure', True)
235 236 self.width = kwargs.get('width', None)
236 237 self.height = kwargs.get('height', None)
237 238 self.colorbar = kwargs.get('colorbar', True)
238 239 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
239 240 self.channels = kwargs.get('channels', None)
240 241 self.titles = kwargs.get('titles', [])
241 242 self.polar = False
242 243 self.type = kwargs.get('type', 'iq')
243 244 self.grid = kwargs.get('grid', False)
244 245 self.pause = kwargs.get('pause', False)
245 246 self.save_code = kwargs.get('save_code', self.CODE)
246 247 self.throttle = kwargs.get('throttle', 0)
247 248 self.exp_code = kwargs.get('exp_code', None)
248 249 self.server = kwargs.get('server', False)
249 250 self.sender_period = kwargs.get('sender_period', 60)
250 251 self.tag = kwargs.get('tag', '')
251 252 self.height_index = kwargs.get('height_index', None)
252 253 self.__throttle_plot = apply_throttle(self.throttle)
253 254 code = self.attr_data if self.attr_data else self.CODE
254 255 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
256 self.tmin = kwargs.get('tmin', None)
255 257
256 258 if self.server:
257 259 if not self.server.startswith('tcp://'):
258 260 self.server = 'tcp://{}'.format(self.server)
259 261 log.success(
260 262 'Sending to server: {}'.format(self.server),
261 263 self.name
262 264 )
263 265
264 266 if isinstance(self.attr_data, str):
265 267 self.attr_data = [self.attr_data]
266 268
267 269 def __setup_plot(self):
268 270 '''
269 271 Common setup for all figures, here figures and axes are created
270 272 '''
271 273
272 274 self.setup()
273 275
274 276 self.time_label = 'LT' if self.localtime else 'UTC'
275 277
276 278 if self.width is None:
277 279 self.width = 8
278 280
279 281 self.figures = []
280 282 self.axes = []
281 283 self.cb_axes = []
282 284 self.pf_axes = []
283 285 self.cmaps = []
284 286
285 287 size = '15%' if self.ncols == 1 else '30%'
286 288 pad = '4%' if self.ncols == 1 else '8%'
287 289
288 290 if self.oneFigure:
289 291 if self.height is None:
290 292 self.height = 1.4 * self.nrows + 1
291 293 fig = plt.figure(figsize=(self.width, self.height),
292 294 edgecolor='k',
293 295 facecolor='w')
294 296 self.figures.append(fig)
295 297 for n in range(self.nplots):
296 298 ax = fig.add_subplot(self.nrows, self.ncols,
297 299 n + 1, polar=self.polar)
298 300 ax.tick_params(labelsize=8)
299 301 ax.firsttime = True
300 302 ax.index = 0
301 303 ax.press = None
302 304 self.axes.append(ax)
303 305 if self.showprofile:
304 306 cax = self.__add_axes(ax, size=size, pad=pad)
305 307 cax.tick_params(labelsize=8)
306 308 self.pf_axes.append(cax)
307 309 else:
308 310 if self.height is None:
309 311 self.height = 3
310 312 for n in range(self.nplots):
311 313 fig = plt.figure(figsize=(self.width, self.height),
312 314 edgecolor='k',
313 315 facecolor='w')
314 316 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
315 317 ax.tick_params(labelsize=8)
316 318 ax.firsttime = True
317 319 ax.index = 0
318 320 ax.press = None
319 321 self.figures.append(fig)
320 322 self.axes.append(ax)
321 323 if self.showprofile:
322 324 cax = self.__add_axes(ax, size=size, pad=pad)
323 325 cax.tick_params(labelsize=8)
324 326 self.pf_axes.append(cax)
325 327
326 328 for n in range(self.nrows):
327 329 if self.colormaps is not None:
328 330 cmap = plt.get_cmap(self.colormaps[n])
329 331 else:
330 332 cmap = plt.get_cmap(self.colormap)
331 333 cmap.set_bad(self.bgcolor, 1.)
332 334 self.cmaps.append(cmap)
333 335
334 336 def __add_axes(self, ax, size='30%', pad='8%'):
335 337 '''
336 338 Add new axes to the given figure
337 339 '''
338 340 divider = make_axes_locatable(ax)
339 341 nax = divider.new_horizontal(size=size, pad=pad)
340 342 ax.figure.add_axes(nax)
341 343 return nax
342 344
343 345 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
344 346 '''
345 347 Create a masked array for missing data
346 348 '''
347 349 if x_buffer.shape[0] < 2:
348 350 return x_buffer, y_buffer, z_buffer
349 351
350 352 deltas = x_buffer[1:] - x_buffer[0:-1]
351 353 x_median = numpy.median(deltas)
352 354
353 355 index = numpy.where(deltas > 5 * x_median)
354 356
355 357 if len(index[0]) != 0:
356 358 z_buffer[::, index[0], ::] = self.__missing
357 359 z_buffer = numpy.ma.masked_inside(z_buffer,
358 360 0.99 * self.__missing,
359 361 1.01 * self.__missing)
360 362
361 363 return x_buffer, y_buffer, z_buffer
362 364
363 365 def decimate(self):
364 366
365 367 # dx = int(len(self.x)/self.__MAXNUMX) + 1
366 368 dy = int(len(self.y) / self.decimation) + 1
367 369
368 370 # x = self.x[::dx]
369 371 x = self.x
370 372 y = self.y[::dy]
371 373 z = self.z[::, ::, ::dy]
372 374
373 375 return x, y, z
374 376
375 377 def format(self):
376 378 '''
377 379 Set min and max values, labels, ticks and titles
378 380 '''
379 381
380 382 for n, ax in enumerate(self.axes):
381 383 if ax.firsttime:
382 384 if self.xaxis != 'time':
383 385 xmin = self.xmin
384 386 xmax = self.xmax
385 387 else:
386 388 xmin = self.tmin
387 389 xmax = self.tmin + self.xrange*60*60
388 390 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
389 391 ax.xaxis.set_major_locator(LinearLocator(9))
390 392 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
391 393 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
392 394 ax.set_facecolor(self.bgcolor)
393 395 if self.xscale:
394 396 ax.xaxis.set_major_formatter(FuncFormatter(
395 397 lambda x, pos: '{0:g}'.format(x*self.xscale)))
396 398 if self.yscale:
397 399 ax.yaxis.set_major_formatter(FuncFormatter(
398 400 lambda x, pos: '{0:g}'.format(x*self.yscale)))
399 401 if self.xlabel is not None:
400 402 ax.set_xlabel(self.xlabel)
401 403 if self.ylabel is not None:
402 404 ax.set_ylabel(self.ylabel)
403 405 if self.showprofile:
404 406 self.pf_axes[n].set_ylim(ymin, ymax)
405 407 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
406 408 self.pf_axes[n].set_xlabel('dB')
407 409 self.pf_axes[n].grid(b=True, axis='x')
408 410 [tick.set_visible(False)
409 411 for tick in self.pf_axes[n].get_yticklabels()]
410 412 if self.colorbar:
411 413 ax.cbar = plt.colorbar(
412 414 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
413 415 ax.cbar.ax.tick_params(labelsize=8)
414 416 ax.cbar.ax.press = None
415 417 if self.cb_label:
416 418 ax.cbar.set_label(self.cb_label, size=8)
417 419 elif self.cb_labels:
418 420 ax.cbar.set_label(self.cb_labels[n], size=8)
419 421 else:
420 422 ax.cbar = None
421 423 ax.set_xlim(xmin, xmax)
422 424 ax.set_ylim(ymin, ymax)
423 425 ax.firsttime = False
424 426 if self.grid:
425 427 ax.grid(True)
426 428 if not self.polar:
427 429 ax.set_title('{} {} {}'.format(
428 430 self.titles[n],
429 431 self.getDateTime(self.data.max_time).strftime(
430 432 '%Y-%m-%d %H:%M:%S'),
431 433 self.time_label),
432 434 size=8)
433 435 else:
434 436 ax.set_title('{}'.format(self.titles[n]), size=8)
435 437 ax.set_ylim(0, 90)
436 438 ax.set_yticks(numpy.arange(0, 90, 20))
437 439 ax.yaxis.labelpad = 40
438 440
439 441 if self.firsttime:
440 442 for n, fig in enumerate(self.figures):
441 443 fig.subplots_adjust(**self.plots_adjust)
442 444 self.firsttime = False
443 445
444 446 def clear_figures(self):
445 447 '''
446 448 Reset axes for redraw plots
447 449 '''
448 450
449 451 for ax in self.axes+self.pf_axes+self.cb_axes:
450 452 ax.clear()
451 453 ax.firsttime = True
452 454 if hasattr(ax, 'cbar') and ax.cbar:
453 455 ax.cbar.remove()
454 456
455 457 def __plot(self):
456 458 '''
457 459 Main function to plot, format and save figures
458 460 '''
459 461
460 462 self.plot()
461 463 self.format()
462 464
463 465 for n, fig in enumerate(self.figures):
464 466 if self.nrows == 0 or self.nplots == 0:
465 467 log.warning('No data', self.name)
466 468 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
467 469 fig.canvas.manager.set_window_title(self.CODE)
468 470 continue
469 471
470 472 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
471 473 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
472 474 fig.canvas.draw()
473 475 if self.show:
474 476 fig.show()
475 477 figpause(0.01)
476 478
477 479 if self.save:
478 480 self.save_figure(n)
479 481
480 482 if self.server:
481 483 self.send_to_server()
482 484
483 485 def __update(self, dataOut, timestamp):
484 486 '''
485 487 '''
486 488
487 489 metadata = {
488 490 'yrange': dataOut.heightList,
489 491 'interval': dataOut.timeInterval,
490 492 'channels': dataOut.channelList
491 493 }
492 494
493 495 data, meta = self.update(dataOut)
494 496 metadata.update(meta)
495 497 self.data.update(data, timestamp, metadata)
496 498
497 499 def save_figure(self, n):
498 500 '''
499 501 '''
500 502
501 503 if (self.data.max_time - self.save_time) <= self.save_period:
502 504 return
503 505
504 506 self.save_time = self.data.max_time
505 507
506 508 fig = self.figures[n]
507 509
508 510 if self.throttle == 0:
509 511 figname = os.path.join(
510 512 self.save,
511 513 self.save_code,
512 514 '{}_{}.png'.format(
513 515 self.save_code,
514 516 self.getDateTime(self.data.max_time).strftime(
515 517 '%Y%m%d_%H%M%S'
516 518 ),
517 519 )
518 520 )
519 521 log.log('Saving figure: {}'.format(figname), self.name)
520 522 if not os.path.isdir(os.path.dirname(figname)):
521 523 os.makedirs(os.path.dirname(figname))
522 524 fig.savefig(figname)
523 525
524 526 figname = os.path.join(
525 527 self.save,
526 528 '{}_{}.png'.format(
527 529 self.save_code,
528 530 self.getDateTime(self.data.min_time).strftime(
529 531 '%Y%m%d'
530 532 ),
531 533 )
532 534 )
533 535
534 536 log.log('Saving figure: {}'.format(figname), self.name)
535 537 if not os.path.isdir(os.path.dirname(figname)):
536 538 os.makedirs(os.path.dirname(figname))
537 539 fig.savefig(figname)
538 540
539 541 def send_to_server(self):
540 542 '''
541 543 '''
542 544
543 545 if self.exp_code == None:
544 546 log.warning('Missing `exp_code` skipping sending to server...')
545 547
546 548 last_time = self.data.max_time
547 549 interval = last_time - self.sender_time
548 550 if interval < self.sender_period:
549 551 return
550 552
551 553 self.sender_time = last_time
552 554
553 555 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
554 556 for attr in attrs:
555 557 value = getattr(self, attr)
556 558 if value:
557 559 if isinstance(value, (numpy.float32, numpy.float64)):
558 560 value = round(float(value), 2)
559 561 self.data.meta[attr] = value
560 562 if self.colormap == 'jet':
561 563 self.data.meta['colormap'] = 'Jet'
562 564 elif 'RdBu' in self.colormap:
563 565 self.data.meta['colormap'] = 'RdBu'
564 566 else:
565 567 self.data.meta['colormap'] = 'Viridis'
566 568 self.data.meta['interval'] = int(interval)
567 569
568 570 self.sender_queue.append(last_time)
569 571
570 572 while True:
571 573 try:
572 574 tm = self.sender_queue.popleft()
573 575 except IndexError:
574 576 break
575 577 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
576 578 self.socket.send_string(msg)
577 579 socks = dict(self.poll.poll(2000))
578 580 if socks.get(self.socket) == zmq.POLLIN:
579 581 reply = self.socket.recv_string()
580 582 if reply == 'ok':
581 583 log.log("Response from server ok", self.name)
582 584 time.sleep(0.1)
583 585 continue
584 586 else:
585 587 log.warning(
586 588 "Malformed reply from server: {}".format(reply), self.name)
587 589 else:
588 590 log.warning(
589 591 "No response from server, retrying...", self.name)
590 592 self.sender_queue.appendleft(tm)
591 593 self.socket.setsockopt(zmq.LINGER, 0)
592 594 self.socket.close()
593 595 self.poll.unregister(self.socket)
594 596 self.socket = self.context.socket(zmq.REQ)
595 597 self.socket.connect(self.server)
596 598 self.poll.register(self.socket, zmq.POLLIN)
597 599 break
598 600
599 601 def setup(self):
600 602 '''
601 603 This method should be implemented in the child class, the following
602 604 attributes should be set:
603 605
604 606 self.nrows: number of rows
605 607 self.ncols: number of cols
606 608 self.nplots: number of plots (channels or pairs)
607 609 self.ylabel: label for Y axes
608 610 self.titles: list of axes title
609 611
610 612 '''
611 613 raise NotImplementedError
612 614
613 615 def plot(self):
614 616 '''
615 617 Must be defined in the child class, the actual plotting method
616 618 '''
617 619 raise NotImplementedError
618 620
619 621 def update(self, dataOut):
620 622 '''
621 623 Must be defined in the child class, update self.data with new data
622 624 '''
623 625
624 626 data = {
625 627 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
626 628 }
627 629 meta = {}
628 630
629 631 return data, meta
630 632
631 633 def run(self, dataOut, **kwargs):
632 634 '''
633 635 Main plotting routine
634 636 '''
635 637
636 638 if self.isConfig is False:
637 639 self.__setup(**kwargs)
638 640
639 641 if self.localtime:
640 642 self.getDateTime = datetime.datetime.fromtimestamp
641 643 else:
642 644 self.getDateTime = datetime.datetime.utcfromtimestamp
643 645
644 646 self.data.setup()
645 647 self.isConfig = True
646 648 if self.server:
647 649 self.context = zmq.Context()
648 650 self.socket = self.context.socket(zmq.REQ)
649 651 self.socket.connect(self.server)
650 652 self.poll = zmq.Poller()
651 653 self.poll.register(self.socket, zmq.POLLIN)
652 654
653 655 tm = getattr(dataOut, self.attr_time)
654 656
655 657 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
656 658 self.save_time = tm
657 659 self.__plot()
658 660 self.tmin += self.xrange*60*60
659 661 self.data.setup()
660 662 self.clear_figures()
661 663
662 664 self.__update(dataOut, tm)
663 665
664 666 if self.isPlotConfig is False:
665 667 self.__setup_plot()
666 668 self.isPlotConfig = True
667 669 if self.xaxis == 'time':
668 670 dt = self.getDateTime(tm)
669 671 if self.xmin is None:
670 672 self.tmin = tm
671 673 self.xmin = dt.hour
672 674 minutes = (self.xmin-int(self.xmin)) * 60
673 675 seconds = (minutes - int(minutes)) * 60
674 676 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
675 677 datetime.datetime(1970, 1, 1)).total_seconds()
676 678 if self.localtime:
677 679 self.tmin += time.timezone
678 680
679 681 if self.xmin is not None and self.xmax is not None:
680 682 self.xrange = self.xmax - self.xmin
681 683
682 684 if self.throttle == 0:
683 685 self.__plot()
684 686 else:
685 687 self.__throttle_plot(self.__plot)#, coerce=coerce)
686 688
687 689 def close(self):
688 690
689 691 if self.data and not self.data.flagNoData:
690 692 self.save_time = 0
691 693 self.__plot()
692 694 if self.data and not self.data.flagNoData and self.pause:
693 695 figpause(10)
Show file before | Show file after | Hide comments
@@ -1,357 +1,356
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot
7 7 from schainpy.utils import log
8 8
9 9 EARTH_RADIUS = 6.3710e3
10 10
11 11
12 12 def ll2xy(lat1, lon1, lat2, lon2):
13 13
14 14 p = 0.017453292519943295
15 15 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
16 16 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
17 17 r = 12742 * numpy.arcsin(numpy.sqrt(a))
18 18 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
19 19 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
20 20 theta = -theta + numpy.pi/2
21 21 return r*numpy.cos(theta), r*numpy.sin(theta)
22 22
23 23
24 24 def km2deg(km):
25 25 '''
26 26 Convert distance in km to degrees
27 27 '''
28 28
29 29 return numpy.rad2deg(km/EARTH_RADIUS)
30 30
31 31
32 32
33 33 class SpectralMomentsPlot(SpectraPlot):
34 34 '''
35 35 Plot for Spectral Moments
36 36 '''
37 37 CODE = 'spc_moments'
38 38 colormap = 'jet'
39 39 plot_type = 'pcolor'
40 40
41 41
42 42 class SnrPlot(RTIPlot):
43 43 '''
44 44 Plot for SNR Data
45 45 '''
46 46
47 47 CODE = 'snr'
48 48 colormap = 'jet'
49 49
50 50 def update(self, dataOut):
51 51
52 52 data = {
53 53 'snr': 10*numpy.log10(dataOut.data_snr)
54 54 }
55 55
56 56 return data, {}
57 57
58 58 class DopplerPlot(RTIPlot):
59 59 '''
60 60 Plot for DOPPLER Data (1st moment)
61 61 '''
62 62
63 63 CODE = 'dop'
64 64 colormap = 'jet'
65 65
66 66 def update(self, dataOut):
67 67
68 68 data = {
69 69 'dop': 10*numpy.log10(dataOut.data_dop)
70 70 }
71 71
72 72 return data, {}
73 73
74 74 class PowerPlot(RTIPlot):
75 75 '''
76 76 Plot for Power Data (0 moment)
77 77 '''
78 78
79 79 CODE = 'pow'
80 80 colormap = 'jet'
81 81
82 82 def update(self, dataOut):
83 83
84 84 data = {
85 85 'pow': 10*numpy.log10(dataOut.data_pow)
86 86 }
87
87 print("data",data)
88 88 return data, {}
89 89
90 90 class SpectralWidthPlot(RTIPlot):
91 91 '''
92 92 Plot for Spectral Width Data (2nd moment)
93 93 '''
94 94
95 95 CODE = 'width'
96 96 colormap = 'jet'
97 97
98 98 def update(self, dataOut):
99 99
100 100 data = {
101 101 'width': dataOut.data_width
102 102 }
103 103
104 104 return data, {}
105 105
106 106 class SkyMapPlot(Plot):
107 107 '''
108 108 Plot for meteors detection data
109 109 '''
110 110
111 111 CODE = 'param'
112 112
113 113 def setup(self):
114 114
115 115 self.ncols = 1
116 116 self.nrows = 1
117 117 self.width = 7.2
118 118 self.height = 7.2
119 119 self.nplots = 1
120 120 self.xlabel = 'Zonal Zenith Angle (deg)'
121 121 self.ylabel = 'Meridional Zenith Angle (deg)'
122 122 self.polar = True
123 123 self.ymin = -180
124 124 self.ymax = 180
125 125 self.colorbar = False
126 126
127 127 def plot(self):
128 128
129 129 arrayParameters = numpy.concatenate(self.data['param'])
130 130 error = arrayParameters[:, -1]
131 131 indValid = numpy.where(error == 0)[0]
132 132 finalMeteor = arrayParameters[indValid, :]
133 133 finalAzimuth = finalMeteor[:, 3]
134 134 finalZenith = finalMeteor[:, 4]
135 135
136 136 x = finalAzimuth * numpy.pi / 180
137 137 y = finalZenith
138 138
139 139 ax = self.axes[0]
140 140
141 141 if ax.firsttime:
142 142 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
143 143 else:
144 144 ax.plot.set_data(x, y)
145 145
146 146 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
147 147 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
148 148 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
149 149 dt2,
150 150 len(x))
151 151 self.titles[0] = title
152 152
153 153
154 154 class GenericRTIPlot(Plot):
155 155 '''
156 156 Plot for data_xxxx object
157 157 '''
158 158
159 159 CODE = 'param'
160 160 colormap = 'viridis'
161 161 plot_type = 'pcolorbuffer'
162 162
163 163 def setup(self):
164 164 self.xaxis = 'time'
165 165 self.ncols = 1
166 166 self.nrows = self.data.shape('param')[0]
167 167 self.nplots = self.nrows
168 168 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
169 169
170 170 if not self.xlabel:
171 171 self.xlabel = 'Time'
172 172
173 173 self.ylabel = 'Height [km]'
174 174 if not self.titles:
175 175 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
176 176
177 177 def update(self, dataOut):
178 178
179 179 data = {
180 180 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
181 181 }
182 182
183 183 meta = {}
184 184
185 185 return data, meta
186 186
187 187 def plot(self):
188 188 # self.data.normalize_heights()
189 189 self.x = self.data.times
190 190 self.y = self.data.yrange
191 191 self.z = self.data['param']
192 192
193 193 self.z = numpy.ma.masked_invalid(self.z)
194 194
195 195 if self.decimation is None:
196 196 x, y, z = self.fill_gaps(self.x, self.y, self.z)
197 197 else:
198 198 x, y, z = self.fill_gaps(*self.decimate())
199 199
200 200 for n, ax in enumerate(self.axes):
201 201
202 202 self.zmax = self.zmax if self.zmax is not None else numpy.max(
203 203 self.z[n])
204 204 self.zmin = self.zmin if self.zmin is not None else numpy.min(
205 205 self.z[n])
206 206
207 207 if ax.firsttime:
208 208 if self.zlimits is not None:
209 209 self.zmin, self.zmax = self.zlimits[n]
210 210
211 211 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
212 212 vmin=self.zmin,
213 213 vmax=self.zmax,
214 214 cmap=self.cmaps[n]
215 215 )
216 216 else:
217 217 if self.zlimits is not None:
218 218 self.zmin, self.zmax = self.zlimits[n]
219 219 ax.collections.remove(ax.collections[0])
220 220 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
221 221 vmin=self.zmin,
222 222 vmax=self.zmax,
223 223 cmap=self.cmaps[n]
224 224 )
225 225
226 226
227 227 class PolarMapPlot(Plot):
228 228 '''
229 229 Plot for weather radar
230 230 '''
231 231
232 232 CODE = 'param'
233 233 colormap = 'seismic'
234 234
235 235 def setup(self):
236 236 self.ncols = 1
237 237 self.nrows = 1
238 238 self.width = 9
239 239 self.height = 8
240 240 self.mode = self.data.meta['mode']
241 241 if self.channels is not None:
242 242 self.nplots = len(self.channels)
243 243 self.nrows = len(self.channels)
244 244 else:
245 245 self.nplots = self.data.shape(self.CODE)[0]
246 246 self.nrows = self.nplots
247 247 self.channels = list(range(self.nplots))
248 248 if self.mode == 'E':
249 249 self.xlabel = 'Longitude'
250 250 self.ylabel = 'Latitude'
251 251 else:
252 252 self.xlabel = 'Range (km)'
253 253 self.ylabel = 'Height (km)'
254 254 self.bgcolor = 'white'
255 255 self.cb_labels = self.data.meta['units']
256 256 self.lat = self.data.meta['latitude']
257 257 self.lon = self.data.meta['longitude']
258 258 self.xmin, self.xmax = float(
259 259 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
260 260 self.ymin, self.ymax = float(
261 261 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
262 262 # self.polar = True
263 263
264 264 def plot(self):
265 265
266 266 for n, ax in enumerate(self.axes):
267 267 data = self.data['param'][self.channels[n]]
268 268
269 269 zeniths = numpy.linspace(
270 270 0, self.data.meta['max_range'], data.shape[1])
271 271 if self.mode == 'E':
272 272 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
273 273 r, theta = numpy.meshgrid(zeniths, azimuths)
274 274 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
275 275 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
276 276 x = km2deg(x) + self.lon
277 277 y = km2deg(y) + self.lat
278 278 else:
279 279 azimuths = numpy.radians(self.data.yrange)
280 280 r, theta = numpy.meshgrid(zeniths, azimuths)
281 281 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
282 282 self.y = zeniths
283 283
284 284 if ax.firsttime:
285 285 if self.zlimits is not None:
286 286 self.zmin, self.zmax = self.zlimits[n]
287 287 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
288 288 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
289 289 vmin=self.zmin,
290 290 vmax=self.zmax,
291 291 cmap=self.cmaps[n])
292 292 else:
293 293 if self.zlimits is not None:
294 294 self.zmin, self.zmax = self.zlimits[n]
295 295 ax.collections.remove(ax.collections[0])
296 296 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
297 297 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
298 298 vmin=self.zmin,
299 299 vmax=self.zmax,
300 300 cmap=self.cmaps[n])
301 301
302 302 if self.mode == 'A':
303 303 continue
304 304
305 305 # plot district names
306 306 f = open('/data/workspace/schain_scripts/distrito.csv')
307 307 for line in f:
308 308 label, lon, lat = [s.strip() for s in line.split(',') if s]
309 309 lat = float(lat)
310 310 lon = float(lon)
311 311 # ax.plot(lon, lat, '.b', ms=2)
312 312 ax.text(lon, lat, label.decode('utf8'), ha='center',
313 313 va='bottom', size='8', color='black')
314 314
315 315 # plot limites
316 316 limites = []
317 317 tmp = []
318 318 for line in open('/data/workspace/schain_scripts/lima.csv'):
319 319 if '#' in line:
320 320 if tmp:
321 321 limites.append(tmp)
322 322 tmp = []
323 323 continue
324 324 values = line.strip().split(',')
325 325 tmp.append((float(values[0]), float(values[1])))
326 326 for points in limites:
327 327 ax.add_patch(
328 328 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
329 329
330 330 # plot Cuencas
331 331 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
332 332 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
333 333 values = [line.strip().split(',') for line in f]
334 334 points = [(float(s[0]), float(s[1])) for s in values]
335 335 ax.add_patch(Polygon(points, ec='b', fc='none'))
336 336
337 337 # plot grid
338 338 for r in (15, 30, 45, 60):
339 339 ax.add_artist(plt.Circle((self.lon, self.lat),
340 340 km2deg(r), color='0.6', fill=False, lw=0.2))
341 341 ax.text(
342 342 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
343 343 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
344 344 '{}km'.format(r),
345 345 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
346 346
347 347 if self.mode == 'E':
348 348 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
349 349 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
350 350 else:
351 351 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
352 352 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
353 353
354 354 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
355 355 self.titles = ['{} {}'.format(
356 356 self.data.parameters[x], title) for x in self.channels]
357
Show file before | Show file after | Hide comments
@@ -1,711 +1,712
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13
14 14
15 15 class SpectraPlot(Plot):
16 16 '''
17 17 Plot for Spectra data
18 18 '''
19 19
20 20 CODE = 'spc'
21 21 colormap = 'jet'
22 22 plot_type = 'pcolor'
23 23 buffering = False
24 channelList = None
24 channelList = []
25 25
26 26 def setup(self):
27 27 self.nplots = len(self.data.channels)
28 28 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
29 29 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
30 30 self.height = 2.6 * self.nrows
31 31
32 32 self.cb_label = 'dB'
33 33 if self.showprofile:
34 34 self.width = 4 * self.ncols
35 35 else:
36 36 self.width = 3.5 * self.ncols
37 37 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
38 38 self.ylabel = 'Range [km]'
39 39
40 40 def update(self, dataOut):
41 41 if self.channelList == None:
42 42 self.channelList = dataOut.channelList
43 43 data = {}
44 44 meta = {}
45 45 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
46 46 data['spc'] = spc
47 47 data['rti'] = dataOut.getPower()
48 48 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
49 49 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
50 50 if self.CODE == 'spc_moments':
51 51 data['moments'] = dataOut.moments
52 52
53 53 return data, meta
54 54
55 55 def plot(self):
56 56 if self.xaxis == "frequency":
57 57 x = self.data.xrange[0]
58 58 self.xlabel = "Frequency (kHz)"
59 59 elif self.xaxis == "time":
60 60 x = self.data.xrange[1]
61 61 self.xlabel = "Time (ms)"
62 62 else:
63 63 x = self.data.xrange[2]
64 64 self.xlabel = "Velocity (m/s)"
65 65
66 66 if self.CODE == 'spc_moments':
67 67 x = self.data.xrange[2]
68 68 self.xlabel = "Velocity (m/s)"
69 69
70 70 self.titles = []
71 71
72 72 y = self.data.yrange
73 73 self.y = y
74 74
75 75 data = self.data[-1]
76 76 z = data['spc']
77 77
78 78 for n, ax in enumerate(self.axes):
79 79 noise = data['noise'][n]
80 80 if self.CODE == 'spc_moments':
81 81 mean = data['moments'][n, 1]
82 82 if ax.firsttime:
83 83 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
84 84 self.xmin = self.xmin if self.xmin else -self.xmax
85 85 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
86 86 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
87 87 ax.plt = ax.pcolormesh(x, y, z[n].T,
88 88 vmin=self.zmin,
89 89 vmax=self.zmax,
90 90 cmap=plt.get_cmap(self.colormap)
91 91 )
92 92
93 93 if self.showprofile:
94 94 ax.plt_profile = self.pf_axes[n].plot(
95 95 data['rti'][n], y)[0]
96 96 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
97 97 color="k", linestyle="dashed", lw=1)[0]
98 98 if self.CODE == 'spc_moments':
99 99 ax.plt_mean = ax.plot(mean, y, color='k')[0]
100 100 else:
101 101 ax.plt.set_array(z[n].T.ravel())
102 102 if self.showprofile:
103 103 ax.plt_profile.set_data(data['rti'][n], y)
104 104 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
105 105 if self.CODE == 'spc_moments':
106 106 ax.plt_mean.set_data(mean, y)
107 107 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
108 108
109 109
110 110 class CrossSpectraPlot(Plot):
111 111
112 112 CODE = 'cspc'
113 113 colormap = 'jet'
114 114 plot_type = 'pcolor'
115 115 zmin_coh = None
116 116 zmax_coh = None
117 117 zmin_phase = None
118 118 zmax_phase = None
119 119
120 120 def setup(self):
121 121
122 122 self.ncols = 4
123 123 self.nplots = len(self.data.pairs) * 2
124 124 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
125 125 self.width = 3.1 * self.ncols
126 126 self.height = 2.6 * self.nrows
127 127 self.ylabel = 'Range [km]'
128 128 self.showprofile = False
129 129 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
130 130
131 131 def update(self, dataOut):
132 132
133 133 data = {}
134 134 meta = {}
135 135
136 136 spc = dataOut.data_spc
137 137 cspc = dataOut.data_cspc
138 138 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
139 139 meta['pairs'] = dataOut.pairsList
140 140
141 141 tmp = []
142 142
143 143 for n, pair in enumerate(meta['pairs']):
144 144 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
145 145 coh = numpy.abs(out)
146 146 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
147 147 tmp.append(coh)
148 148 tmp.append(phase)
149 149
150 150 data['cspc'] = numpy.array(tmp)
151 151
152 152 return data, meta
153 153
154 154 def plot(self):
155 155
156 156 if self.xaxis == "frequency":
157 157 x = self.data.xrange[0]
158 158 self.xlabel = "Frequency (kHz)"
159 159 elif self.xaxis == "time":
160 160 x = self.data.xrange[1]
161 161 self.xlabel = "Time (ms)"
162 162 else:
163 163 x = self.data.xrange[2]
164 164 self.xlabel = "Velocity (m/s)"
165 165
166 166 self.titles = []
167 167
168 168 y = self.data.yrange
169 169 self.y = y
170 170
171 171 data = self.data[-1]
172 172 cspc = data['cspc']
173 173
174 174 for n in range(len(self.data.pairs)):
175 175 pair = self.data.pairs[n]
176 176 coh = cspc[n*2]
177 177 phase = cspc[n*2+1]
178 178 ax = self.axes[2 * n]
179 179 if ax.firsttime:
180 180 ax.plt = ax.pcolormesh(x, y, coh.T,
181 181 vmin=0,
182 182 vmax=1,
183 183 cmap=plt.get_cmap(self.colormap_coh)
184 184 )
185 185 else:
186 186 ax.plt.set_array(coh.T.ravel())
187 187 self.titles.append(
188 188 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
189 189
190 190 ax = self.axes[2 * n + 1]
191 191 if ax.firsttime:
192 192 ax.plt = ax.pcolormesh(x, y, phase.T,
193 193 vmin=-180,
194 194 vmax=180,
195 195 cmap=plt.get_cmap(self.colormap_phase)
196 196 )
197 197 else:
198 198 ax.plt.set_array(phase.T.ravel())
199 199 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
200 200
201 201
202 202 class RTIPlot(Plot):
203 203 '''
204 204 Plot for RTI data
205 205 '''
206 206
207 207 CODE = 'rti'
208 208 colormap = 'jet'
209 209 plot_type = 'pcolorbuffer'
210 210 titles = None
211 channelList = None
211 channelList = []
212 212
213 213 def setup(self):
214 214 self.xaxis = 'time'
215 215 self.ncols = 1
216 print("dataChannels ",self.data.channels)
216 217 self.nrows = len(self.data.channels)
217 218 self.nplots = len(self.data.channels)
218 219 self.ylabel = 'Range [km]'
219 220 self.xlabel = 'Time'
220 221 self.cb_label = 'dB'
221 222 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
222 223 self.titles = ['{} Channel {}'.format(
223 224 self.CODE.upper(), x) for x in range(self.nplots)]
224
225 print("SETUP")
225 226 def update(self, dataOut):
226 if self.channelList == None:
227 if len(self.channelList) == 0:
227 228 self.channelList = dataOut.channelList
228 229 data = {}
229 230 meta = {}
230 231 data['rti'] = dataOut.getPower()
231 232 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
232 233
233 234 return data, meta
234 235
235 236 def plot(self):
236 237 self.x = self.data.times
237 238 self.y = self.data.yrange
238 239 self.z = self.data[self.CODE]
239 240 self.z = numpy.ma.masked_invalid(self.z)
240 241 if self.channelList != None:
241 242 self.titles = ['{} Channel {}'.format(
242 243 self.CODE.upper(), x) for x in self.channelList]
243 244
244 245 if self.decimation is None:
245 246 x, y, z = self.fill_gaps(self.x, self.y, self.z)
246 247 else:
247 248 x, y, z = self.fill_gaps(*self.decimate())
248 249
249 250 for n, ax in enumerate(self.axes):
250 251 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
251 252 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
252 253 data = self.data[-1]
253 254 if ax.firsttime:
254 255 ax.plt = ax.pcolormesh(x, y, z[n].T,
255 256 vmin=self.zmin,
256 257 vmax=self.zmax,
257 258 cmap=plt.get_cmap(self.colormap)
258 259 )
259 260 if self.showprofile:
260 261 ax.plot_profile = self.pf_axes[n].plot(
261 262 data['rti'][n], self.y)[0]
262 263 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
263 264 color="k", linestyle="dashed", lw=1)[0]
264 265 else:
265 266 ax.collections.remove(ax.collections[0])
266 267 ax.plt = ax.pcolormesh(x, y, z[n].T,
267 268 vmin=self.zmin,
268 269 vmax=self.zmax,
269 270 cmap=plt.get_cmap(self.colormap)
270 271 )
271 272 if self.showprofile:
272 273 ax.plot_profile.set_data(data['rti'][n], self.y)
273 274 ax.plot_noise.set_data(numpy.repeat(
274 275 data['noise'][n], len(self.y)), self.y)
275 276
276 277
277 278 class CoherencePlot(RTIPlot):
278 279 '''
279 280 Plot for Coherence data
280 281 '''
281 282
282 283 CODE = 'coh'
283 284
284 285 def setup(self):
285 286 self.xaxis = 'time'
286 287 self.ncols = 1
287 288 self.nrows = len(self.data.pairs)
288 289 self.nplots = len(self.data.pairs)
289 290 self.ylabel = 'Range [km]'
290 291 self.xlabel = 'Time'
291 292 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
292 293 if self.CODE == 'coh':
293 294 self.cb_label = ''
294 295 self.titles = [
295 296 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
296 297 else:
297 298 self.cb_label = 'Degrees'
298 299 self.titles = [
299 300 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
300 301
301 302 def update(self, dataOut):
302 303
303 304 data = {}
304 305 meta = {}
305 306 data['coh'] = dataOut.getCoherence()
306 307 meta['pairs'] = dataOut.pairsList
307 308
308 309 return data, meta
309 310
310 311 class PhasePlot(CoherencePlot):
311 312 '''
312 313 Plot for Phase map data
313 314 '''
314 315
315 316 CODE = 'phase'
316 317 colormap = 'seismic'
317 318
318 319 def update(self, dataOut):
319 320
320 321 data = {}
321 322 meta = {}
322 323 data['phase'] = dataOut.getCoherence(phase=True)
323 324 meta['pairs'] = dataOut.pairsList
324 325
325 326 return data, meta
326 327
327 328 class NoisePlot(Plot):
328 329 '''
329 330 Plot for noise
330 331 '''
331 332
332 333 CODE = 'noise'
333 334 plot_type = 'scatterbuffer'
334 335
335 336 def setup(self):
336 337 self.xaxis = 'time'
337 338 self.ncols = 1
338 339 self.nrows = 1
339 340 self.nplots = 1
340 341 self.ylabel = 'Intensity [dB]'
341 342 self.xlabel = 'Time'
342 343 self.titles = ['Noise']
343 344 self.colorbar = False
344 345 self.plots_adjust.update({'right': 0.85 })
345 346
346 347 def update(self, dataOut):
347 348
348 349 data = {}
349 350 meta = {}
350 351 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
351 352 meta['yrange'] = numpy.array([])
352 353
353 354 return data, meta
354 355
355 356 def plot(self):
356 357
357 358 x = self.data.times
358 359 xmin = self.data.min_time
359 360 xmax = xmin + self.xrange * 60 * 60
360 361 Y = self.data['noise']
361 362
362 363 if self.axes[0].firsttime:
363 364 self.ymin = numpy.nanmin(Y) - 5
364 365 self.ymax = numpy.nanmax(Y) + 5
365 366 for ch in self.data.channels:
366 367 y = Y[ch]
367 368 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
368 369 plt.legend(bbox_to_anchor=(1.18, 1.0))
369 370 else:
370 371 for ch in self.data.channels:
371 372 y = Y[ch]
372 373 self.axes[0].lines[ch].set_data(x, y)
373 374
374 375
375 376 class PowerProfilePlot(Plot):
376 377
377 378 CODE = 'pow_profile'
378 379 plot_type = 'scatter'
379 380
380 381 def setup(self):
381 382
382 383 self.ncols = 1
383 384 self.nrows = 1
384 385 self.nplots = 1
385 386 self.height = 4
386 387 self.width = 3
387 388 self.ylabel = 'Range [km]'
388 389 self.xlabel = 'Intensity [dB]'
389 390 self.titles = ['Power Profile']
390 391 self.colorbar = False
391 392
392 393 def update(self, dataOut):
393 394
394 395 data = {}
395 396 meta = {}
396 397 data[self.CODE] = dataOut.getPower()
397 398
398 399 return data, meta
399 400
400 401 def plot(self):
401 402
402 403 y = self.data.yrange
403 404 self.y = y
404 405
405 406 x = self.data[-1][self.CODE]
406 407
407 408 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
408 409 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
409 410
410 411 if self.axes[0].firsttime:
411 412 for ch in self.data.channels:
412 413 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
413 414 plt.legend()
414 415 else:
415 416 for ch in self.data.channels:
416 417 self.axes[0].lines[ch].set_data(x[ch], y)
417 418
418 419
419 420 class SpectraCutPlot(Plot):
420 421
421 422 CODE = 'spc_cut'
422 423 plot_type = 'scatter'
423 424 buffering = False
424 425
425 426 def setup(self):
426 427
427 428 self.nplots = len(self.data.channels)
428 429 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
429 430 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
430 431 self.width = 3.4 * self.ncols + 1.5
431 432 self.height = 3 * self.nrows
432 433 self.ylabel = 'Power [dB]'
433 434 self.colorbar = False
434 435 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
435 436
436 437 def update(self, dataOut):
437 438
438 439 data = {}
439 440 meta = {}
440 441 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
441 442 data['spc'] = spc
442 443 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
443 444
444 445 return data, meta
445 446
446 447 def plot(self):
447 448 if self.xaxis == "frequency":
448 449 x = self.data.xrange[0][1:]
449 450 self.xlabel = "Frequency (kHz)"
450 451 elif self.xaxis == "time":
451 452 x = self.data.xrange[1]
452 453 self.xlabel = "Time (ms)"
453 454 else:
454 455 x = self.data.xrange[2]
455 456 self.xlabel = "Velocity (m/s)"
456 457
457 458 self.titles = []
458 459
459 460 y = self.data.yrange
460 461 z = self.data[-1]['spc']
461 462
462 463 if self.height_index:
463 464 index = numpy.array(self.height_index)
464 465 else:
465 466 index = numpy.arange(0, len(y), int((len(y))/9))
466 467
467 468 for n, ax in enumerate(self.axes):
468 469 if ax.firsttime:
469 470 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
470 471 self.xmin = self.xmin if self.xmin else -self.xmax
471 472 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
472 473 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
473 474 ax.plt = ax.plot(x, z[n, :, index].T)
474 475 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
475 476 self.figures[0].legend(ax.plt, labels, loc='center right')
476 477 else:
477 478 for i, line in enumerate(ax.plt):
478 479 line.set_data(x, z[n, :, index[i]])
479 480 self.titles.append('CH {}'.format(n))
480 481
481 482
482 483 class BeaconPhase(Plot):
483 484
484 485 __isConfig = None
485 486 __nsubplots = None
486 487
487 488 PREFIX = 'beacon_phase'
488 489
489 490 def __init__(self):
490 491 Plot.__init__(self)
491 492 self.timerange = 24*60*60
492 493 self.isConfig = False
493 494 self.__nsubplots = 1
494 495 self.counter_imagwr = 0
495 496 self.WIDTH = 800
496 497 self.HEIGHT = 400
497 498 self.WIDTHPROF = 120
498 499 self.HEIGHTPROF = 0
499 500 self.xdata = None
500 501 self.ydata = None
501 502
502 503 self.PLOT_CODE = BEACON_CODE
503 504
504 505 self.FTP_WEI = None
505 506 self.EXP_CODE = None
506 507 self.SUB_EXP_CODE = None
507 508 self.PLOT_POS = None
508 509
509 510 self.filename_phase = None
510 511
511 512 self.figfile = None
512 513
513 514 self.xmin = None
514 515 self.xmax = None
515 516
516 517 def getSubplots(self):
517 518
518 519 ncol = 1
519 520 nrow = 1
520 521
521 522 return nrow, ncol
522 523
523 524 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
524 525
525 526 self.__showprofile = showprofile
526 527 self.nplots = nplots
527 528
528 529 ncolspan = 7
529 530 colspan = 6
530 531 self.__nsubplots = 2
531 532
532 533 self.createFigure(id = id,
533 534 wintitle = wintitle,
534 535 widthplot = self.WIDTH+self.WIDTHPROF,
535 536 heightplot = self.HEIGHT+self.HEIGHTPROF,
536 537 show=show)
537 538
538 539 nrow, ncol = self.getSubplots()
539 540
540 541 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
541 542
542 543 def save_phase(self, filename_phase):
543 544 f = open(filename_phase,'w+')
544 545 f.write('\n\n')
545 546 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
546 547 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
547 548 f.close()
548 549
549 550 def save_data(self, filename_phase, data, data_datetime):
550 551 f=open(filename_phase,'a')
551 552 timetuple_data = data_datetime.timetuple()
552 553 day = str(timetuple_data.tm_mday)
553 554 month = str(timetuple_data.tm_mon)
554 555 year = str(timetuple_data.tm_year)
555 556 hour = str(timetuple_data.tm_hour)
556 557 minute = str(timetuple_data.tm_min)
557 558 second = str(timetuple_data.tm_sec)
558 559 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
559 560 f.close()
560 561
561 562 def plot(self):
562 563 log.warning('TODO: Not yet implemented...')
563 564
564 565 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
565 566 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
566 567 timerange=None,
567 568 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
568 569 server=None, folder=None, username=None, password=None,
569 570 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
570 571
571 572 if dataOut.flagNoData:
572 573 return dataOut
573 574
574 575 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
575 576 return
576 577
577 578 if pairsList == None:
578 579 pairsIndexList = dataOut.pairsIndexList[:10]
579 580 else:
580 581 pairsIndexList = []
581 582 for pair in pairsList:
582 583 if pair not in dataOut.pairsList:
583 584 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
584 585 pairsIndexList.append(dataOut.pairsList.index(pair))
585 586
586 587 if pairsIndexList == []:
587 588 return
588 589
589 590 # if len(pairsIndexList) > 4:
590 591 # pairsIndexList = pairsIndexList[0:4]
591 592
592 593 hmin_index = None
593 594 hmax_index = None
594 595
595 596 if hmin != None and hmax != None:
596 597 indexes = numpy.arange(dataOut.nHeights)
597 598 hmin_list = indexes[dataOut.heightList >= hmin]
598 599 hmax_list = indexes[dataOut.heightList <= hmax]
599 600
600 601 if hmin_list.any():
601 602 hmin_index = hmin_list[0]
602 603
603 604 if hmax_list.any():
604 605 hmax_index = hmax_list[-1]+1
605 606
606 607 x = dataOut.getTimeRange()
607 608
608 609 thisDatetime = dataOut.datatime
609 610
610 611 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
611 612 xlabel = "Local Time"
612 613 ylabel = "Phase (degrees)"
613 614
614 615 update_figfile = False
615 616
616 617 nplots = len(pairsIndexList)
617 618 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
618 619 phase_beacon = numpy.zeros(len(pairsIndexList))
619 620 for i in range(nplots):
620 621 pair = dataOut.pairsList[pairsIndexList[i]]
621 622 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
622 623 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
623 624 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
624 625 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
625 626 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
626 627
627 628 if dataOut.beacon_heiIndexList:
628 629 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
629 630 else:
630 631 phase_beacon[i] = numpy.average(phase)
631 632
632 633 if not self.isConfig:
633 634
634 635 nplots = len(pairsIndexList)
635 636
636 637 self.setup(id=id,
637 638 nplots=nplots,
638 639 wintitle=wintitle,
639 640 showprofile=showprofile,
640 641 show=show)
641 642
642 643 if timerange != None:
643 644 self.timerange = timerange
644 645
645 646 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
646 647
647 648 if ymin == None: ymin = 0
648 649 if ymax == None: ymax = 360
649 650
650 651 self.FTP_WEI = ftp_wei
651 652 self.EXP_CODE = exp_code
652 653 self.SUB_EXP_CODE = sub_exp_code
653 654 self.PLOT_POS = plot_pos
654 655
655 656 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
656 657 self.isConfig = True
657 658 self.figfile = figfile
658 659 self.xdata = numpy.array([])
659 660 self.ydata = numpy.array([])
660 661
661 662 update_figfile = True
662 663
663 664 #open file beacon phase
664 665 path = '%s%03d' %(self.PREFIX, self.id)
665 666 beacon_file = os.path.join(path,'%s.txt'%self.name)
666 667 self.filename_phase = os.path.join(figpath,beacon_file)
667 668 #self.save_phase(self.filename_phase)
668 669
669 670
670 671 #store data beacon phase
671 672 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
672 673
673 674 self.setWinTitle(title)
674 675
675 676
676 677 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
677 678
678 679 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
679 680
680 681 axes = self.axesList[0]
681 682
682 683 self.xdata = numpy.hstack((self.xdata, x[0:1]))
683 684
684 685 if len(self.ydata)==0:
685 686 self.ydata = phase_beacon.reshape(-1,1)
686 687 else:
687 688 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
688 689
689 690
690 691 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
691 692 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
692 693 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
693 694 XAxisAsTime=True, grid='both'
694 695 )
695 696
696 697 self.draw()
697 698
698 699 if dataOut.ltctime >= self.xmax:
699 700 self.counter_imagwr = wr_period
700 701 self.isConfig = False
701 702 update_figfile = True
702 703
703 704 self.save(figpath=figpath,
704 705 figfile=figfile,
705 706 save=save,
706 707 ftp=ftp,
707 708 wr_period=wr_period,
708 709 thisDatetime=thisDatetime,
709 710 update_figfile=update_figfile)
710 711
711 712 return dataOut
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
Show file before | Show file after | Hide comments
@@ -1,663 +1,661
1 1 '''
2 2 Created on Set 9, 2015
3 3
4 4 @author: roj-idl71 Karim Kuyeng
5 5
6 6 @update: 2021, Joab Apaza
7 7 '''
8 8
9 9 import os
10 10 import sys
11 11 import glob
12 12 import fnmatch
13 13 import datetime
14 14 import time
15 15 import re
16 16 import h5py
17 17 import numpy
18 18
19 19 try:
20 20 from gevent import sleep
21 21 except:
22 22 from time import sleep
23 23
24 24 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
25 25 from schainpy.model.data.jrodata import Voltage
26 26 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
27 27 from numpy import imag
28 28
29 29
30 30 class AMISRReader(ProcessingUnit):
31 31 '''
32 32 classdocs
33 33 '''
34 34
35 35 def __init__(self):
36 36 '''
37 37 Constructor
38 38 '''
39 39
40 40 ProcessingUnit.__init__(self)
41 41
42 42 self.set = None
43 43 self.subset = None
44 44 self.extension_file = '.h5'
45 45 self.dtc_str = 'dtc'
46 46 self.dtc_id = 0
47 47 self.status = True
48 48 self.isConfig = False
49 49 self.dirnameList = []
50 50 self.filenameList = []
51 51 self.fileIndex = None
52 52 self.flagNoMoreFiles = False
53 53 self.flagIsNewFile = 0
54 54 self.filename = ''
55 55 self.amisrFilePointer = None
56 56 self.realBeamCode = []
57 57 self.beamCodeMap = None
58 58 self.azimuthList = []
59 59 self.elevationList = []
60 60 self.dataShape = None
61 61
62 62
63 63
64 64 self.profileIndex = 0
65 65
66 66
67 67 self.beamCodeByFrame = None
68 68 self.radacTimeByFrame = None
69 69
70 70 self.dataset = None
71 71
72 72 self.__firstFile = True
73 73
74 74 self.buffer = None
75 75
76 76 self.timezone = 'ut'
77 77
78 78 self.__waitForNewFile = 20
79 79 self.__filename_online = None
80 80 #Is really necessary create the output object in the initializer
81 81 self.dataOut = Voltage()
82 82 self.dataOut.error=False
83 83
84 84
85 85 def setup(self,path=None,
86 86 startDate=None,
87 87 endDate=None,
88 88 startTime=None,
89 89 endTime=None,
90 90 walk=True,
91 91 timezone='ut',
92 92 all=0,
93 93 code = None,
94 94 nCode = 0,
95 95 nBaud = 0,
96 96 online=False):
97 97
98 98
99 99
100 100 self.timezone = timezone
101 101 self.all = all
102 102 self.online = online
103 103
104 104 self.code = code
105 105 self.nCode = int(nCode)
106 106 self.nBaud = int(nBaud)
107 107
108 108
109 109
110 110 #self.findFiles()
111 111 if not(online):
112 112 #Busqueda de archivos offline
113 113 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk)
114 114 else:
115 115 self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
116 116
117 117 if not(self.filenameList):
118 118 print("There is no files into the folder: %s"%(path))
119 119 sys.exit()
120 120
121 121 self.fileIndex = 0
122 122
123 123 self.readNextFile(online)
124 124
125 125 '''
126 126 Add code
127 127 '''
128 128 self.isConfig = True
129 129 # print("Setup Done")
130 130 pass
131 131
132 132
133 133 def readAMISRHeader(self,fp):
134 134
135 135 if self.isConfig and (not self.flagNoMoreFiles):
136 136 newShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
137 137 if self.dataShape != newShape and newShape != None:
138 138 print("\nNEW FILE HAS A DIFFERENT SHAPE")
139 139 print(self.dataShape,newShape,"\n")
140 140 return 0
141 141 else:
142 142 self.dataShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
143 143
144 144
145 145 header = 'Raw11/Data/RadacHeader'
146 146 self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
147 147 if (self.startDate> datetime.date(2021, 7, 15)): #Se cambió la forma de extracción de Apuntes el 17
148 148 self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode()
149 149 self.trueBeams = self.beamcodeFile.split("\n")
150 150 self.trueBeams.pop()#remove last
151 151 [self.realBeamCode.append(x) for x in self.trueBeams if x not in self.realBeamCode]
152 152 self.beamCode = [int(x, 16) for x in self.realBeamCode]
153 153 else:
154 154 _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes
155 155 self.beamCode = _beamCode[0,:]
156 156
157 157 if self.beamCodeMap == None:
158 158 self.beamCodeMap = fp['Setup/BeamcodeMap']
159 159 for beam in self.beamCode:
160 160 beamAziElev = numpy.where(self.beamCodeMap[:,0]==beam)
161 161 beamAziElev = beamAziElev[0].squeeze()
162 162 self.azimuthList.append(self.beamCodeMap[beamAziElev,1])
163 163 self.elevationList.append(self.beamCodeMap[beamAziElev,2])
164 164 #print("Beamssss: ",self.beamCodeMap[beamAziElev,1],self.beamCodeMap[beamAziElev,2])
165 165 #print(self.beamCode)
166 166 #self.code = fp.get(header+'/Code') # NOT USE FOR THIS
167 167 self.frameCount = fp.get(header+'/FrameCount')# NOT USE FOR THIS
168 168 self.modeGroup = fp.get(header+'/ModeGroup')# NOT USE FOR THIS
169 169 self.nsamplesPulse = fp.get(header+'/NSamplesPulse')# TO GET NSA OR USING DATA FOR THAT
170 170 self.pulseCount = fp.get(header+'/PulseCount')# NOT USE FOR THIS
171 171 self.radacTime = fp.get(header+'/RadacTime')# 1st TIME ON FILE ANDE CALCULATE THE REST WITH IPP*nindexprofile
172 172 self.timeCount = fp.get(header+'/TimeCount')# NOT USE FOR THIS
173 173 self.timeStatus = fp.get(header+'/TimeStatus')# NOT USE FOR THIS
174 174 self.rangeFromFile = fp.get('Raw11/Data/Samples/Range')
175 175 self.frequency = fp.get('Rx/Frequency')
176 176 txAus = fp.get('Raw11/Data/Pulsewidth')
177 177
178 178
179 179 self.nblocks = self.pulseCount.shape[0] #nblocks
180 180
181 181 self.nprofiles = self.pulseCount.shape[1] #nprofile
182 182 self.nsa = self.nsamplesPulse[0,0] #ngates
183 183 self.nchannels = len(self.beamCode)
184 184 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
185 185 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
186 186 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
187 187
188 188 #filling radar controller header parameters
189 189 self.__ippKm = self.ippSeconds *.15*1e6 # in km
190 190 self.__txA = (txAus.value)*.15 #(ipp[us]*.15km/1us) in km
191 191 self.__txB = 0
192 192 nWindows=1
193 193 self.__nSamples = self.nsa
194 194 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
195 195 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
196 196
197 197 #for now until understand why the code saved is different (code included even though code not in tuf file)
198 198 #self.__codeType = 0
199 199 # self.__nCode = None
200 200 # self.__nBaud = None
201 201 self.__code = self.code
202 202 self.__codeType = 0
203 203 if self.code != None:
204 204 self.__codeType = 1
205 205 self.__nCode = self.nCode
206 206 self.__nBaud = self.nBaud
207 207 #self.__code = 0
208 208
209 209 #filling system header parameters
210 210 self.__nSamples = self.nsa
211 211 self.newProfiles = self.nprofiles/self.nchannels
212 212 self.__channelList = list(range(self.nchannels))
213 213
214 214 self.__frequency = self.frequency[0][0]
215 215
216 216
217 217 return 1
218 218
219 219
220 220 def createBuffers(self):
221 221
222 222 pass
223 223
224 224 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
225 225 self.path = path
226 226 self.startDate = startDate
227 227 self.endDate = endDate
228 228 self.startTime = startTime
229 229 self.endTime = endTime
230 230 self.walk = walk
231 231
232 232 def __checkPath(self):
233 233 if os.path.exists(self.path):
234 234 self.status = 1
235 235 else:
236 236 self.status = 0
237 237 print('Path:%s does not exists'%self.path)
238 238
239 239 return
240 240
241 241
242 242 def __selDates(self, amisr_dirname_format):
243 243 try:
244 244 year = int(amisr_dirname_format[0:4])
245 245 month = int(amisr_dirname_format[4:6])
246 246 dom = int(amisr_dirname_format[6:8])
247 247 thisDate = datetime.date(year,month,dom)
248 248
249 249 if (thisDate>=self.startDate and thisDate <= self.endDate):
250 250 return amisr_dirname_format
251 251 except:
252 252 return None
253 253
254 254
255 255 def __findDataForDates(self,online=False):
256 256
257 257 if not(self.status):
258 258 return None
259 259
260 260 pat = '\d+.\d+'
261 261 dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
262 262 dirnameList = [x for x in dirnameList if x!=None]
263 263 dirnameList = [x.string for x in dirnameList]
264 264 if not(online):
265 265 dirnameList = [self.__selDates(x) for x in dirnameList]
266 266 dirnameList = [x for x in dirnameList if x!=None]
267 267 if len(dirnameList)>0:
268 268 self.status = 1
269 269 self.dirnameList = dirnameList
270 270 self.dirnameList.sort()
271 271 else:
272 272 self.status = 0
273 273 return None
274 274
275 275 def __getTimeFromData(self):
276 276 startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
277 277 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
278 278
279 279 print('Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader))
280 280 print('........................................')
281 281 filter_filenameList = []
282 282 self.filenameList.sort()
283 283 #for i in range(len(self.filenameList)-1):
284 284 for i in range(len(self.filenameList)):
285 285 filename = self.filenameList[i]
286 286 fp = h5py.File(filename,'r')
287 287 time_str = fp.get('Time/RadacTimeString')
288 288
289 289 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
290 290 #startDateTimeStr_File = "2019-12-16 09:21:11"
291 291 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
292 292 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
293 293
294 294 #endDateTimeStr_File = "2019-12-16 11:10:11"
295 295 endDateTimeStr_File = time_str[-1][-1].decode('UTF-8').split('.')[0]
296 296 junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
297 297 endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
298 298
299 299 fp.close()
300 300
301 301 #print("check time", startDateTime_File)
302 302 if self.timezone == 'lt':
303 303 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
304 304 endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
305 305 if (endDateTime_File>=startDateTime_Reader and endDateTime_File<=endDateTime_Reader):
306 306 filter_filenameList.append(filename)
307 307
308 308 if (endDateTime_File>endDateTime_Reader):
309 309 break
310 310
311 311
312 312 filter_filenameList.sort()
313 313 self.filenameList = filter_filenameList
314 314 return 1
315 315
316 316 def __filterByGlob1(self, dirName):
317 317 filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
318 318 filter_files.sort()
319 319 filterDict = {}
320 320 filterDict.setdefault(dirName)
321 321 filterDict[dirName] = filter_files
322 322 return filterDict
323 323
324 324 def __getFilenameList(self, fileListInKeys, dirList):
325 325 for value in fileListInKeys:
326 326 dirName = list(value.keys())[0]
327 327 for file in value[dirName]:
328 328 filename = os.path.join(dirName, file)
329 329 self.filenameList.append(filename)
330 330
331 331
332 332 def __selectDataForTimes(self, online=False):
333 333 #aun no esta implementado el filtro for tiempo
334 334 if not(self.status):
335 335 return None
336 336
337 337 dirList = [os.path.join(self.path,x) for x in self.dirnameList]
338 338
339 339 fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
340 340
341 341 self.__getFilenameList(fileListInKeys, dirList)
342 342 if not(online):
343 343 #filtro por tiempo
344 344 if not(self.all):
345 345 self.__getTimeFromData()
346 346
347 347 if len(self.filenameList)>0:
348 348 self.status = 1
349 349 self.filenameList.sort()
350 350 else:
351 351 self.status = 0
352 352 return None
353 353
354 354 else:
355 355 #get the last file - 1
356 356 self.filenameList = [self.filenameList[-2]]
357 357 new_dirnameList = []
358 358 for dirname in self.dirnameList:
359 359 junk = numpy.array([dirname in x for x in self.filenameList])
360 360 junk_sum = junk.sum()
361 361 if junk_sum > 0:
362 362 new_dirnameList.append(dirname)
363 363 self.dirnameList = new_dirnameList
364 364 return 1
365 365
366 366 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
367 367 endTime=datetime.time(23,59,59),walk=True):
368 368
369 369 if endDate ==None:
370 370 startDate = datetime.datetime.utcnow().date()
371 371 endDate = datetime.datetime.utcnow().date()
372 372
373 373 self.__setParameters(path=path, startDate=startDate, endDate=endDate,startTime = startTime,endTime=endTime, walk=walk)
374 374
375 375 self.__checkPath()
376 376
377 377 self.__findDataForDates(online=True)
378 378
379 379 self.dirnameList = [self.dirnameList[-1]]
380 380
381 381 self.__selectDataForTimes(online=True)
382 382
383 383 return
384 384
385 385
386 386 def searchFilesOffLine(self,
387 387 path,
388 388 startDate,
389 389 endDate,
390 390 startTime=datetime.time(0,0,0),
391 391 endTime=datetime.time(23,59,59),
392 392 walk=True):
393 393
394 394 self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
395 395
396 396 self.__checkPath()
397 397
398 398 self.__findDataForDates()
399 399
400 400 self.__selectDataForTimes()
401 401
402 402 for i in range(len(self.filenameList)):
403 403 print("%s" %(self.filenameList[i]))
404 404
405 405 return
406 406
407 407 def __setNextFileOffline(self):
408 408
409 409 try:
410 410 self.filename = self.filenameList[self.fileIndex]
411 411 self.amisrFilePointer = h5py.File(self.filename,'r')
412 412 self.fileIndex += 1
413 413 except:
414 414 self.flagNoMoreFiles = 1
415 415 print("No more Files")
416 416 return 0
417 417
418 418 self.flagIsNewFile = 1
419 419 print("Setting the file: %s"%self.filename)
420 420
421 421 return 1
422 422
423 423
424 424 def __setNextFileOnline(self):
425 425 filename = self.filenameList[0]
426 426 if self.__filename_online != None:
427 427 self.__selectDataForTimes(online=True)
428 428 filename = self.filenameList[0]
429 429 wait = 0
430 430 self.__waitForNewFile=300 ## DEBUG:
431 431 while self.__filename_online == filename:
432 432 print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
433 433 if wait == 5:
434 434 self.flagNoMoreFiles = 1
435 435 return 0
436 436 sleep(self.__waitForNewFile)
437 437 self.__selectDataForTimes(online=True)
438 438 filename = self.filenameList[0]
439 439 wait += 1
440 440
441 441 self.__filename_online = filename
442 442
443 443 self.amisrFilePointer = h5py.File(filename,'r')
444 444 self.flagIsNewFile = 1
445 445 self.filename = filename
446 446 print("Setting the file: %s"%self.filename)
447 447 return 1
448 448
449 449
450 450 def readData(self):
451 451 buffer = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
452 452 re = buffer[:,:,:,0]
453 453 im = buffer[:,:,:,1]
454 454 dataset = re + im*1j
455 455
456 456 self.radacTime = self.amisrFilePointer.get('Raw11/Data/RadacHeader/RadacTime')
457 457 timeset = self.radacTime[:,0]
458 458
459 459 return dataset,timeset
460 460
461 461 def reshapeData(self):
462 462 #self.beamCodeByPulse, self.beamCode, self.nblocks, self.nprofiles, self.nsa,
463 463 channels = self.beamCodeByPulse[0,:]
464 464 nchan = self.nchannels
465 465 #self.newProfiles = self.nprofiles/nchan #must be defined on filljroheader
466 466 nblocks = self.nblocks
467 467 nsamples = self.nsa
468 468
469 469 #Dimensions : nChannels, nProfiles, nSamples
470 470 new_block = numpy.empty((nblocks, nchan, numpy.int_(self.newProfiles), nsamples), dtype="complex64")
471 471 ############################################
472 472
473 473 for thisChannel in range(nchan):
474 474 new_block[:,thisChannel,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[thisChannel])[0],:]
475 475
476 476
477 477 new_block = numpy.transpose(new_block, (1,0,2,3))
478 478 new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
479 479
480 480 return new_block
481 481
482 482 def updateIndexes(self):
483 483
484 484 pass
485 485
486 486 def fillJROHeader(self):
487 487
488 488 #fill radar controller header
489 489 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(ipp=self.__ippKm,
490 490 txA=self.__txA,
491 491 txB=0,
492 492 nWindows=1,
493 493 nHeights=self.__nSamples,
494 494 firstHeight=self.__firstHeight,
495 495 deltaHeight=self.__deltaHeight,
496 496 codeType=self.__codeType,
497 497 nCode=self.__nCode, nBaud=self.__nBaud,
498 498 code = self.__code,
499 499 fClock=1)
500 500
501 501 #fill system header
502 502 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
503 503 nProfiles=self.newProfiles,
504 504 nChannels=len(self.__channelList),
505 505 adcResolution=14,
506 506 pciDioBusWidth=32)
507 507
508 508 self.dataOut.type = "Voltage"
509 509 self.dataOut.data = None
510 510 self.dataOut.dtype = numpy.dtype([('real','<i8'),('imag','<i8')])
511 511 # self.dataOut.nChannels = 0
512 512
513 513 # self.dataOut.nHeights = 0
514 514
515 515 self.dataOut.nProfiles = self.newProfiles*self.nblocks
516 516 #self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
517 517 ranges = numpy.reshape(self.rangeFromFile.value,(-1))
518 518 self.dataOut.heightList = ranges/1000.0 #km
519 519 self.dataOut.channelList = self.__channelList
520 520 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
521 521
522 522 # self.dataOut.channelIndexList = None
523 523
524 524
525 525 self.dataOut.azimuthList = numpy.array(self.azimuthList)
526 526 self.dataOut.elevationList = numpy.array(self.elevationList)
527 527 self.dataOut.codeList = numpy.array(self.beamCode)
528 528 #print(self.dataOut.elevationList)
529 529 self.dataOut.flagNoData = True
530 530
531 531 #Set to TRUE if the data is discontinuous
532 532 self.dataOut.flagDiscontinuousBlock = False
533 533
534 534 self.dataOut.utctime = None
535 535
536 536 #self.dataOut.timeZone = -5 #self.__timezone/60 #timezone like jroheader, difference in minutes between UTC and localtime
537 537 if self.timezone == 'lt':
538 538 self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
539 539 else:
540 540 self.dataOut.timeZone = 0 #by default time is UTC
541 541
542 542 self.dataOut.dstFlag = 0
543 543 self.dataOut.errorCount = 0
544 544 self.dataOut.nCohInt = 1
545 545 self.dataOut.flagDecodeData = False #asumo que la data esta decodificada
546 546 self.dataOut.flagDeflipData = False #asumo que la data esta sin flip
547 547 self.dataOut.flagShiftFFT = False
548 548 self.dataOut.ippSeconds = self.ippSeconds
549 549
550 550 #Time interval between profiles
551 551 #self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
552 552
553 553 self.dataOut.frequency = self.__frequency
554 554 self.dataOut.realtime = self.online
555 555 pass
556 556
557 557 def readNextFile(self,online=False):
558 558
559 559 if not(online):
560 560 newFile = self.__setNextFileOffline()
561 561 else:
562 562 newFile = self.__setNextFileOnline()
563 563
564 564 if not(newFile):
565 565 self.dataOut.error = True
566 566 return 0
567 567
568 568 if not self.readAMISRHeader(self.amisrFilePointer):
569 569 self.dataOut.error = True
570 570 return 0
571 571
572 572 self.createBuffers()
573 573 self.fillJROHeader()
574 574
575 575 #self.__firstFile = False
576 576
577 577
578 578
579 579 self.dataset,self.timeset = self.readData()
580 580
581 581 if self.endDate!=None:
582 582 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
583 583 time_str = self.amisrFilePointer.get('Time/RadacTimeString')
584 584 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
585 585 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
586 586 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
587 587 if self.timezone == 'lt':
588 588 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
589 589 if (startDateTime_File>endDateTime_Reader):
590 590 return 0
591 591
592 592 self.jrodataset = self.reshapeData()
593 593 #----self.updateIndexes()
594 594 self.profileIndex = 0
595 595
596 596 return 1
597 597
598 598
599 599 def __hasNotDataInBuffer(self):
600 600 if self.profileIndex >= (self.newProfiles*self.nblocks):
601 601 return 1
602 602 return 0
603 603
604 604
605 605 def getData(self):
606 606
607 607 if self.flagNoMoreFiles:
608 608 self.dataOut.flagNoData = True
609 609 return 0
610 610
611 611 if self.__hasNotDataInBuffer():
612 612 if not (self.readNextFile(self.online)):
613 613 return 0
614 614
615 615
616 616 if self.dataset is None: # setear esta condicion cuando no hayan datos por leer
617 617 self.dataOut.flagNoData = True
618 618 return 0
619 619
620 620 #self.dataOut.data = numpy.reshape(self.jrodataset[self.profileIndex,:],(1,-1))
621 621
622 622 self.dataOut.data = self.jrodataset[:,self.profileIndex,:]
623 623
624 624 #print("R_t",self.timeset)
625 625
626 626 #self.dataOut.utctime = self.jrotimeset[self.profileIndex]
627 627 #verificar basic header de jro data y ver si es compatible con este valor
628 628 #self.dataOut.utctime = self.timeset + (self.profileIndex * self.ippSeconds * self.nchannels)
629 629 indexprof = numpy.mod(self.profileIndex, self.newProfiles)
630 630 indexblock = self.profileIndex/self.newProfiles
631 631 #print (indexblock, indexprof)
632 632 diffUTC = 1.8e4 #UTC diference from peru in seconds --Joab
633 633 diffUTC = 0
634 634 t_comp = (indexprof * self.ippSeconds * self.nchannels) + diffUTC #
635 635
636 636 #print("utc :",indexblock," __ ",t_comp)
637 637 #print(numpy.shape(self.timeset))
638 638 self.dataOut.utctime = self.timeset[numpy.int_(indexblock)] + t_comp
639 639 #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
640 640 #print(self.dataOut.utctime)
641 641 self.dataOut.profileIndex = self.profileIndex
642 642 #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
643 643 self.dataOut.flagNoData = False
644 644 # if indexprof == 0:
645 645 # print self.dataOut.utctime
646 646
647 647 self.profileIndex += 1
648 648
649 #return self.dataOut.data
649 return self.dataOut.data
650 650
651 651
652 652 def run(self, **kwargs):
653 653 '''
654 654 This method will be called many times so here you should put all your code
655 655 '''
656 656 #print("running kamisr")
657 657 if not self.isConfig:
658 658 self.setup(**kwargs)
659 659 self.isConfig = True
660 660
661 661 self.getData()
662 #return(self.dataOut.data)
663 return(self.dataOut)
Show file before | Show file after | Hide comments
@@ -1,626 +1,651
1 1 import os
2 2 import time
3 3 import datetime
4 4
5 5 import numpy
6 6 import h5py
7 7
8 8 import schainpy.admin
9 9 from schainpy.model.data.jrodata import *
10 10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 11 from schainpy.model.io.jroIO_base import *
12 12 from schainpy.utils import log
13 13
14 14
15 15 class HDFReader(Reader, ProcessingUnit):
16 16 """Processing unit to read HDF5 format files
17 17
18 18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 20 attributes.
21 21 It is possible to read any HDF5 file by given the structure in the `description`
22 22 parameter, also you can add extra values to metadata with the parameter `extras`.
23 23
24 24 Parameters:
25 25 -----------
26 26 path : str
27 27 Path where files are located.
28 28 startDate : date
29 29 Start date of the files
30 30 endDate : list
31 31 End date of the files
32 32 startTime : time
33 33 Start time of the files
34 34 endTime : time
35 35 End time of the files
36 36 description : dict, optional
37 37 Dictionary with the description of the HDF5 file
38 38 extras : dict, optional
39 39 Dictionary with extra metadata to be be added to `dataOut`
40 40
41 41 Examples
42 42 --------
43 43
44 44 desc = {
45 45 'Data': {
46 46 'data_output': ['u', 'v', 'w'],
47 47 'utctime': 'timestamps',
48 48 } ,
49 49 'Metadata': {
50 50 'heightList': 'heights'
51 51 }
52 52 }
53 53
54 54 desc = {
55 55 'Data': {
56 56 'data_output': 'winds',
57 57 'utctime': 'timestamps'
58 58 },
59 59 'Metadata': {
60 60 'heightList': 'heights'
61 61 }
62 62 }
63 63
64 64 extras = {
65 65 'timeZone': 300
66 66 }
67 67
68 68 reader = project.addReadUnit(
69 69 name='HDFReader',
70 70 path='/path/to/files',
71 71 startDate='2019/01/01',
72 72 endDate='2019/01/31',
73 73 startTime='00:00:00',
74 74 endTime='23:59:59',
75 75 # description=json.dumps(desc),
76 76 # extras=json.dumps(extras),
77 77 )
78 78
79 79 """
80 80
81 81 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
82 82
83 83 def __init__(self):
84 84 ProcessingUnit.__init__(self)
85 85 self.dataOut = Parameters()
86 86 self.ext = ".hdf5"
87 87 self.optchar = "D"
88 88 self.meta = {}
89 89 self.data = {}
90 90 self.open_file = h5py.File
91 91 self.open_mode = 'r'
92 92 self.description = {}
93 93 self.extras = {}
94 94 self.filefmt = "*%Y%j***"
95 95 self.folderfmt = "*%Y%j"
96 96 self.utcoffset = 0
97 97
98 98 def setup(self, **kwargs):
99 99
100 100 self.set_kwargs(**kwargs)
101 101 if not self.ext.startswith('.'):
102 102 self.ext = '.{}'.format(self.ext)
103 103
104 104 if self.online:
105 105 log.log("Searching files in online mode...", self.name)
106 106
107 107 for nTries in range(self.nTries):
108 108 fullpath = self.searchFilesOnLine(self.path, self.startDate,
109 109 self.endDate, self.expLabel, self.ext, self.walk,
110 110 self.filefmt, self.folderfmt)
111 pathname, filename = os.path.split(fullpath)
112 print(pathname,filename)
111 113 try:
112 114 fullpath = next(fullpath)
115
113 116 except:
114 117 fullpath = None
115 118
116 119 if fullpath:
117 120 break
118 121
119 122 log.warning(
120 123 'Waiting {} sec for a valid file in {}: try {} ...'.format(
121 124 self.delay, self.path, nTries + 1),
122 125 self.name)
123 126 time.sleep(self.delay)
124 127
125 128 if not(fullpath):
126 129 raise schainpy.admin.SchainError(
127 130 'There isn\'t any valid file in {}'.format(self.path))
128 131
129 132 pathname, filename = os.path.split(fullpath)
130 133 self.year = int(filename[1:5])
131 134 self.doy = int(filename[5:8])
132 135 self.set = int(filename[8:11]) - 1
133 136 else:
134 137 log.log("Searching files in {}".format(self.path), self.name)
135 138 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
136 139 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
137 140
138 141 self.setNextFile()
139 142
140 143 return
141 144
145
142 146 def readFirstHeader(self):
143 147 '''Read metadata and data'''
144 148
145 149 self.__readMetadata()
146 150 self.__readData()
147 151 self.__setBlockList()
148 152
149 153 if 'type' in self.meta:
150 154 self.dataOut = eval(self.meta['type'])()
151 155
152 156 for attr in self.meta:
157 print("attr: ", attr)
153 158 setattr(self.dataOut, attr, self.meta[attr])
154 159
160
155 161 self.blockIndex = 0
156 162
157 163 return
158 164
159 165 def __setBlockList(self):
160 166 '''
161 167 Selects the data within the times defined
162 168
163 169 self.fp
164 170 self.startTime
165 171 self.endTime
166 172 self.blockList
167 173 self.blocksPerFile
168 174
169 175 '''
170 176
171 177 startTime = self.startTime
172 178 endTime = self.endTime
173 179 thisUtcTime = self.data['utctime'] + self.utcoffset
174 180 self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
175 181 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
176
182 self.startFileDatetime = thisDatetime
183 print("datee ",self.startFileDatetime)
177 184 thisDate = thisDatetime.date()
178 185 thisTime = thisDatetime.time()
179 186
180 187 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
181 188 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
182 189
183 190 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
184 191
185 192 self.blockList = ind
186 193 self.blocksPerFile = len(ind)
194 self.blocksPerFile = len(thisUtcTime)
187 195 return
188 196
189 197 def __readMetadata(self):
190 198 '''
191 199 Reads Metadata
192 200 '''
193 201
194 202 meta = {}
195 203
196 204 if self.description:
197 205 for key, value in self.description['Metadata'].items():
198 206 meta[key] = self.fp[value][()]
199 207 else:
200 208 grp = self.fp['Metadata']
201 209 for name in grp:
202 210 meta[name] = grp[name][()]
203 211
204 212 if self.extras:
205 213 for key, value in self.extras.items():
206 214 meta[key] = value
207 215 self.meta = meta
208 216
209 217 return
210 218
219
220
221 def checkForRealPath(self, nextFile, nextDay):
222
223 # print("check FRP")
224 # dt = self.startFileDatetime + datetime.timedelta(1)
225 # filename = '{}.{}{}'.format(self.path, dt.strftime('%Y%m%d'), self.ext)
226 # fullfilename = os.path.join(self.path, filename)
227 # print("check Path ",fullfilename,filename)
228 # if os.path.exists(fullfilename):
229 # return fullfilename, filename
230 # return None, filename
231 return None,None
232
211 233 def __readData(self):
212 234
213 235 data = {}
214 236
215 237 if self.description:
216 238 for key, value in self.description['Data'].items():
217 239 if isinstance(value, str):
218 240 if isinstance(self.fp[value], h5py.Dataset):
219 241 data[key] = self.fp[value][()]
220 242 elif isinstance(self.fp[value], h5py.Group):
221 243 array = []
222 244 for ch in self.fp[value]:
223 245 array.append(self.fp[value][ch][()])
224 246 data[key] = numpy.array(array)
225 247 elif isinstance(value, list):
226 248 array = []
227 249 for ch in value:
228 250 array.append(self.fp[ch][()])
229 251 data[key] = numpy.array(array)
230 252 else:
231 253 grp = self.fp['Data']
232 254 for name in grp:
233 255 if isinstance(grp[name], h5py.Dataset):
234 256 array = grp[name][()]
235 257 elif isinstance(grp[name], h5py.Group):
236 258 array = []
237 259 for ch in grp[name]:
238 260 array.append(grp[name][ch][()])
239 261 array = numpy.array(array)
240 262 else:
241 263 log.warning('Unknown type: {}'.format(name))
242 264
243 265 if name in self.description:
244 266 key = self.description[name]
245 267 else:
246 268 key = name
247 269 data[key] = array
248 270
249 271 self.data = data
250 272 return
251 273
252 274 def getData(self):
253
275 if not self.isDateTimeInRange(self.startFileDatetime, self.startDate, self.endDate, self.startTime, self.endTime):
276 self.dataOut.flagNoData = True
277 self.dataOut.error = True
278 return
254 279 for attr in self.data:
255 280 if self.data[attr].ndim == 1:
256 281 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
257 282 else:
258 283 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
259 284
260 285 self.dataOut.flagNoData = False
261 286 self.blockIndex += 1
262 287
263 288 log.log("Block No. {}/{} -> {}".format(
264 289 self.blockIndex,
265 290 self.blocksPerFile,
266 291 self.dataOut.datatime.ctime()), self.name)
267 292
268 293 return
269 294
270 295 def run(self, **kwargs):
271 296
272 297 if not(self.isConfig):
273 298 self.setup(**kwargs)
274 299 self.isConfig = True
275 300
276 301 if self.blockIndex == self.blocksPerFile:
277 302 self.setNextFile()
278 303
279 304 self.getData()
280 305
281 306 return
282 307
283 308 @MPDecorator
284 309 class HDFWriter(Operation):
285 310 """Operation to write HDF5 files.
286 311
287 312 The HDF5 file contains by default two groups Data and Metadata where
288 313 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
289 314 parameters, data attributes are normaly time dependent where the metadata
290 315 are not.
291 316 It is possible to customize the structure of the HDF5 file with the
292 317 optional description parameter see the examples.
293 318
294 319 Parameters:
295 320 -----------
296 321 path : str
297 322 Path where files will be saved.
298 323 blocksPerFile : int
299 324 Number of blocks per file
300 325 metadataList : list
301 326 List of the dataOut attributes that will be saved as metadata
302 327 dataList : int
303 328 List of the dataOut attributes that will be saved as data
304 329 setType : bool
305 330 If True the name of the files corresponds to the timestamp of the data
306 331 description : dict, optional
307 332 Dictionary with the desired description of the HDF5 file
308 333
309 334 Examples
310 335 --------
311 336
312 337 desc = {
313 338 'data_output': {'winds': ['z', 'w', 'v']},
314 339 'utctime': 'timestamps',
315 340 'heightList': 'heights'
316 341 }
317 342 desc = {
318 343 'data_output': ['z', 'w', 'v'],
319 344 'utctime': 'timestamps',
320 345 'heightList': 'heights'
321 346 }
322 347 desc = {
323 348 'Data': {
324 349 'data_output': 'winds',
325 350 'utctime': 'timestamps'
326 351 },
327 352 'Metadata': {
328 353 'heightList': 'heights'
329 354 }
330 355 }
331 356
332 357 writer = proc_unit.addOperation(name='HDFWriter')
333 358 writer.addParameter(name='path', value='/path/to/file')
334 359 writer.addParameter(name='blocksPerFile', value='32')
335 360 writer.addParameter(name='metadataList', value='heightList,timeZone')
336 361 writer.addParameter(name='dataList',value='data_output,utctime')
337 362 # writer.addParameter(name='description',value=json.dumps(desc))
338 363
339 364 """
340 365
341 366 ext = ".hdf5"
342 367 optchar = "D"
343 368 filename = None
344 369 path = None
345 370 setFile = None
346 371 fp = None
347 372 firsttime = True
348 373 #Configurations
349 374 blocksPerFile = None
350 375 blockIndex = None
351 376 dataOut = None
352 377 #Data Arrays
353 378 dataList = None
354 379 metadataList = None
355 380 currentDay = None
356 381 lastTime = None
357 382
358 383 def __init__(self):
359 384
360 385 Operation.__init__(self)
361 386 return
362 387
363 388 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None):
364 389 self.path = path
365 390 self.blocksPerFile = blocksPerFile
366 391 self.metadataList = metadataList
367 392 self.dataList = [s.strip() for s in dataList]
368 393 self.setType = setType
369 394 self.description = description
370 395
371 396 if self.metadataList is None:
372 397 self.metadataList = self.dataOut.metadata_list
373 398
374 399 tableList = []
375 400 dsList = []
376 401
377 402 for i in range(len(self.dataList)):
378 403 dsDict = {}
379 404 if hasattr(self.dataOut, self.dataList[i]):
380 405 dataAux = getattr(self.dataOut, self.dataList[i])
381 406 dsDict['variable'] = self.dataList[i]
382 407 else:
383 408 log.warning('Attribute {} not found in dataOut', self.name)
384 409 continue
385 410
386 411 if dataAux is None:
387 412 continue
388 413 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
389 414 dsDict['nDim'] = 0
390 415 else:
391 416 dsDict['nDim'] = len(dataAux.shape)
392 417 dsDict['shape'] = dataAux.shape
393 418 dsDict['dsNumber'] = dataAux.shape[0]
394 419 dsDict['dtype'] = dataAux.dtype
395 420
396 421 dsList.append(dsDict)
397 422
398 423 self.dsList = dsList
399 424 self.currentDay = self.dataOut.datatime.date()
400 425
401 426 def timeFlag(self):
402 427 currentTime = self.dataOut.utctime
403 428 timeTuple = time.localtime(currentTime)
404 429 dataDay = timeTuple.tm_yday
405 430
406 431 if self.lastTime is None:
407 432 self.lastTime = currentTime
408 433 self.currentDay = dataDay
409 434 return False
410 435
411 436 timeDiff = currentTime - self.lastTime
412 437
413 438 #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
414 439 if dataDay != self.currentDay:
415 440 self.currentDay = dataDay
416 441 return True
417 442 elif timeDiff > 3*60*60:
418 443 self.lastTime = currentTime
419 444 return True
420 445 else:
421 446 self.lastTime = currentTime
422 447 return False
423 448
424 449 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
425 450 dataList=[], setType=None, description={}):
426 451
427 452 self.dataOut = dataOut
428 453 if not(self.isConfig):
429 454 self.setup(path=path, blocksPerFile=blocksPerFile,
430 455 metadataList=metadataList, dataList=dataList,
431 456 setType=setType, description=description)
432 457
433 458 self.isConfig = True
434 459 self.setNextFile()
435 460
436 461 self.putData()
437 462 return
438 463
439 464 def setNextFile(self):
440 465
441 466 ext = self.ext
442 467 path = self.path
443 468 setFile = self.setFile
444 469
445 470 timeTuple = time.localtime(self.dataOut.utctime)
446 471 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
447 472 fullpath = os.path.join(path, subfolder)
448 473
449 474 if os.path.exists(fullpath):
450 475 filesList = os.listdir(fullpath)
451 476 filesList = [k for k in filesList if k.startswith(self.optchar)]
452 477 if len( filesList ) > 0:
453 478 filesList = sorted(filesList, key=str.lower)
454 479 filen = filesList[-1]
455 480 # el filename debera tener el siguiente formato
456 481 # 0 1234 567 89A BCDE (hex)
457 482 # x YYYY DDD SSS .ext
458 483 if isNumber(filen[8:11]):
459 484 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
460 485 else:
461 486 setFile = -1
462 487 else:
463 488 setFile = -1 #inicializo mi contador de seteo
464 489 else:
465 490 os.makedirs(fullpath)
466 491 setFile = -1 #inicializo mi contador de seteo
467 492
468 493 if self.setType is None:
469 494 setFile += 1
470 495 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
471 496 timeTuple.tm_year,
472 497 timeTuple.tm_yday,
473 498 setFile,
474 499 ext )
475 500 else:
476 501 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
477 502 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
478 503 timeTuple.tm_year,
479 504 timeTuple.tm_yday,
480 505 setFile,
481 506 ext )
482 507
483 508 self.filename = os.path.join( path, subfolder, file )
484 509
485 510 #Setting HDF5 File
486 511 self.fp = h5py.File(self.filename, 'w')
487 512 #write metadata
488 513 self.writeMetadata(self.fp)
489 514 #Write data
490 515 self.writeData(self.fp)
491 516
492 517 def getLabel(self, name, x=None):
493 518
494 519 if x is None:
495 520 if 'Data' in self.description:
496 521 data = self.description['Data']
497 522 if 'Metadata' in self.description:
498 523 data.update(self.description['Metadata'])
499 524 else:
500 525 data = self.description
501 526 if name in data:
502 527 if isinstance(data[name], str):
503 528 return data[name]
504 529 elif isinstance(data[name], list):
505 530 return None
506 531 elif isinstance(data[name], dict):
507 532 for key, value in data[name].items():
508 533 return key
509 534 return name
510 535 else:
511 536 if 'Metadata' in self.description:
512 537 meta = self.description['Metadata']
513 538 else:
514 539 meta = self.description
515 540 if name in meta:
516 541 if isinstance(meta[name], list):
517 542 return meta[name][x]
518 543 elif isinstance(meta[name], dict):
519 544 for key, value in meta[name].items():
520 545 return value[x]
521 546 if 'cspc' in name:
522 547 return 'pair{:02d}'.format(x)
523 548 else:
524 549 return 'channel{:02d}'.format(x)
525 550
526 551 def writeMetadata(self, fp):
527 552
528 553 if self.description:
529 554 if 'Metadata' in self.description:
530 555 grp = fp.create_group('Metadata')
531 556 else:
532 557 grp = fp
533 558 else:
534 559 grp = fp.create_group('Metadata')
535 560
536 561 for i in range(len(self.metadataList)):
537 562 if not hasattr(self.dataOut, self.metadataList[i]):
538 563 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
539 564 continue
540 565 value = getattr(self.dataOut, self.metadataList[i])
541 566 if isinstance(value, bool):
542 567 if value is True:
543 568 value = 1
544 569 else:
545 570 value = 0
546 571 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
547 572 return
548 573
549 574 def writeData(self, fp):
550 575
551 576 if self.description:
552 577 if 'Data' in self.description:
553 578 grp = fp.create_group('Data')
554 579 else:
555 580 grp = fp
556 581 else:
557 582 grp = fp.create_group('Data')
558 583
559 584 dtsets = []
560 585 data = []
561 586
562 587 for dsInfo in self.dsList:
563 588 if dsInfo['nDim'] == 0:
564 589 ds = grp.create_dataset(
565 590 self.getLabel(dsInfo['variable']),
566 591 (self.blocksPerFile, ),
567 592 chunks=True,
568 593 dtype=numpy.float64)
569 594 dtsets.append(ds)
570 595 data.append((dsInfo['variable'], -1))
571 596 else:
572 597 label = self.getLabel(dsInfo['variable'])
573 598 if label is not None:
574 599 sgrp = grp.create_group(label)
575 600 else:
576 601 sgrp = grp
577 602 for i in range(dsInfo['dsNumber']):
578 603 ds = sgrp.create_dataset(
579 604 self.getLabel(dsInfo['variable'], i),
580 605 (self.blocksPerFile, ) + dsInfo['shape'][1:],
581 606 chunks=True,
582 607 dtype=dsInfo['dtype'])
583 608 dtsets.append(ds)
584 609 data.append((dsInfo['variable'], i))
585 610 fp.flush()
586 611
587 612 log.log('Creating file: {}'.format(fp.filename), self.name)
588 613
589 614 self.ds = dtsets
590 615 self.data = data
591 616 self.firsttime = True
592 617 self.blockIndex = 0
593 618 return
594 619
595 620 def putData(self):
596 621
597 622 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
598 623 self.closeFile()
599 624 self.setNextFile()
600 625
601 626 for i, ds in enumerate(self.ds):
602 627 attr, ch = self.data[i]
603 628 if ch == -1:
604 629 ds[self.blockIndex] = getattr(self.dataOut, attr)
605 630 else:
606 631 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
607 632
608 633 self.fp.flush()
609 634 self.blockIndex += 1
610 635 log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
611 636
612 637 return
613 638
614 639 def closeFile(self):
615 640
616 641 if self.blockIndex != self.blocksPerFile:
617 642 for ds in self.ds:
618 643 ds.resize(self.blockIndex, axis=0)
619 644
620 645 if self.fp:
621 646 self.fp.flush()
622 647 self.fp.close()
623 648
624 649 def close(self):
625 650
626 651 self.closeFile()
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,1411 +1,1411
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15 import math
16 16
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 18 from schainpy.model.data.jrodata import Spectra
19 19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 20 from schainpy.utils import log
21 21
22 22 from scipy.optimize import curve_fit
23 23
24 24
25 25 class SpectraProc(ProcessingUnit):
26 26
27 27 def __init__(self):
28 28
29 29 ProcessingUnit.__init__(self)
30 30
31 31 self.buffer = None
32 32 self.firstdatatime = None
33 33 self.profIndex = 0
34 34 self.dataOut = Spectra()
35 35 self.id_min = None
36 36 self.id_max = None
37 37 self.setupReq = False #Agregar a todas las unidades de proc
38 38
39 39 def __updateSpecFromVoltage(self):
40 40
41 41 self.dataOut.timeZone = self.dataIn.timeZone
42 42 self.dataOut.dstFlag = self.dataIn.dstFlag
43 43 self.dataOut.errorCount = self.dataIn.errorCount
44 44 self.dataOut.useLocalTime = self.dataIn.useLocalTime
45 45 try:
46 46 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
47 47 except:
48 48 pass
49 49 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
50 50 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
51 51 self.dataOut.channelList = self.dataIn.channelList
52 52 self.dataOut.heightList = self.dataIn.heightList
53 53 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
54 54 self.dataOut.nProfiles = self.dataOut.nFFTPoints
55 55 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
56 56 self.dataOut.utctime = self.firstdatatime
57 57 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
58 58 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
59 59 self.dataOut.flagShiftFFT = False
60 60 self.dataOut.nCohInt = self.dataIn.nCohInt
61 61 self.dataOut.nIncohInt = 1
62 62 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
63 63 self.dataOut.frequency = self.dataIn.frequency
64 64 self.dataOut.realtime = self.dataIn.realtime
65 65 self.dataOut.azimuth = self.dataIn.azimuth
66 66 self.dataOut.zenith = self.dataIn.zenith
67 67 self.dataOut.codeList = self.dataIn.codeList
68 68 self.dataOut.azimuthList = self.dataIn.azimuthList
69 69 self.dataOut.elevationList = self.dataIn.elevationList
70 70
71 71 def __getFft(self):
72 72 """
73 73 Convierte valores de Voltaje a Spectra
74 74
75 75 Affected:
76 76 self.dataOut.data_spc
77 77 self.dataOut.data_cspc
78 78 self.dataOut.data_dc
79 79 self.dataOut.heightList
80 80 self.profIndex
81 81 self.buffer
82 82 self.dataOut.flagNoData
83 83 """
84 84 fft_volt = numpy.fft.fft(
85 85 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
86 86 fft_volt = fft_volt.astype(numpy.dtype('complex'))
87 87 dc = fft_volt[:, 0, :]
88 88
89 89 # calculo de self-spectra
90 90 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
91 91 spc = fft_volt * numpy.conjugate(fft_volt)
92 92 spc = spc.real
93 93
94 94 blocksize = 0
95 95 blocksize += dc.size
96 96 blocksize += spc.size
97 97
98 98 cspc = None
99 99 pairIndex = 0
100 100 if self.dataOut.pairsList != None:
101 101 # calculo de cross-spectra
102 102 cspc = numpy.zeros(
103 103 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
104 104 for pair in self.dataOut.pairsList:
105 105 if pair[0] not in self.dataOut.channelList:
106 106 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
107 107 str(pair), str(self.dataOut.channelList)))
108 108 if pair[1] not in self.dataOut.channelList:
109 109 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
110 110 str(pair), str(self.dataOut.channelList)))
111 111
112 112 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
113 113 numpy.conjugate(fft_volt[pair[1], :, :])
114 114 pairIndex += 1
115 115 blocksize += cspc.size
116 116
117 117 self.dataOut.data_spc = spc
118 118 self.dataOut.data_cspc = cspc
119 119 self.dataOut.data_dc = dc
120 120 self.dataOut.blockSize = blocksize
121 121 self.dataOut.flagShiftFFT = False
122 122
123 123 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
124 124
125 125 if self.dataIn.type == "Spectra":
126 126 self.dataOut.copy(self.dataIn)
127 127 if shift_fft:
128 128 #desplaza a la derecha en el eje 2 determinadas posiciones
129 129 shift = int(self.dataOut.nFFTPoints/2)
130 130 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
131 131
132 132 if self.dataOut.data_cspc is not None:
133 133 #desplaza a la derecha en el eje 2 determinadas posiciones
134 134 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
135 135 if pairsList:
136 136 self.__selectPairs(pairsList)
137 137
138 138 elif self.dataIn.type == "Voltage":
139 139
140 140 self.dataOut.flagNoData = True
141 141
142 142 if nFFTPoints == None:
143 143 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
144 144
145 145 if nProfiles == None:
146 146 nProfiles = nFFTPoints
147 147
148 148 if ippFactor == None:
149 149 self.dataOut.ippFactor = 1
150 150
151 151 self.dataOut.nFFTPoints = nFFTPoints
152 152
153 153 if self.buffer is None:
154 154 self.buffer = numpy.zeros((self.dataIn.nChannels,
155 155 nProfiles,
156 156 self.dataIn.nHeights),
157 157 dtype='complex')
158 158
159 159 if self.dataIn.flagDataAsBlock:
160 160 nVoltProfiles = self.dataIn.data.shape[1]
161 161
162 162 if nVoltProfiles == nProfiles:
163 163 self.buffer = self.dataIn.data.copy()
164 164 self.profIndex = nVoltProfiles
165 165
166 166 elif nVoltProfiles < nProfiles:
167 167
168 168 if self.profIndex == 0:
169 169 self.id_min = 0
170 170 self.id_max = nVoltProfiles
171 171
172 172 self.buffer[:, self.id_min:self.id_max,
173 173 :] = self.dataIn.data
174 174 self.profIndex += nVoltProfiles
175 175 self.id_min += nVoltProfiles
176 176 self.id_max += nVoltProfiles
177 177 else:
178 178 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
179 179 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
180 180 self.dataOut.flagNoData = True
181 181 else:
182 182 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
183 183 self.profIndex += 1
184 184
185 185 if self.firstdatatime == None:
186 186 self.firstdatatime = self.dataIn.utctime
187 187
188 188 if self.profIndex == nProfiles:
189 189 self.__updateSpecFromVoltage()
190 190 if pairsList == None:
191 191 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
192 192 else:
193 193 self.dataOut.pairsList = pairsList
194 194 self.__getFft()
195 195 self.dataOut.flagNoData = False
196 196 self.firstdatatime = None
197 197 self.profIndex = 0
198 198 else:
199 199 raise ValueError("The type of input object '%s' is not valid".format(
200 200 self.dataIn.type))
201 201
202 202 def __selectPairs(self, pairsList):
203 203
204 204 if not pairsList:
205 205 return
206 206
207 207 pairs = []
208 208 pairsIndex = []
209 209
210 210 for pair in pairsList:
211 211 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
212 212 continue
213 213 pairs.append(pair)
214 214 pairsIndex.append(pairs.index(pair))
215 215
216 216 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
217 217 self.dataOut.pairsList = pairs
218 218
219 219 return
220 220
221 221 def selectFFTs(self, minFFT, maxFFT ):
222 222 """
223 223 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
224 224 minFFT<= FFT <= maxFFT
225 225 """
226 226
227 227 if (minFFT > maxFFT):
228 228 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
229 229
230 230 if (minFFT < self.dataOut.getFreqRange()[0]):
231 231 minFFT = self.dataOut.getFreqRange()[0]
232 232
233 233 if (maxFFT > self.dataOut.getFreqRange()[-1]):
234 234 maxFFT = self.dataOut.getFreqRange()[-1]
235 235
236 236 minIndex = 0
237 237 maxIndex = 0
238 238 FFTs = self.dataOut.getFreqRange()
239 239
240 240 inda = numpy.where(FFTs >= minFFT)
241 241 indb = numpy.where(FFTs <= maxFFT)
242 242
243 243 try:
244 244 minIndex = inda[0][0]
245 245 except:
246 246 minIndex = 0
247 247
248 248 try:
249 249 maxIndex = indb[0][-1]
250 250 except:
251 251 maxIndex = len(FFTs)
252 252
253 253 self.selectFFTsByIndex(minIndex, maxIndex)
254 254
255 255 return 1
256 256
257 257 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
258 258 newheis = numpy.where(
259 259 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
260 260
261 261 if hei_ref != None:
262 262 newheis = numpy.where(self.dataOut.heightList > hei_ref)
263 263
264 264 minIndex = min(newheis[0])
265 265 maxIndex = max(newheis[0])
266 266 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
267 267 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
268 268
269 269 # determina indices
270 270 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
271 271 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
272 272 avg_dB = 10 * \
273 273 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
274 274 beacon_dB = numpy.sort(avg_dB)[-nheis:]
275 275 beacon_heiIndexList = []
276 276 for val in avg_dB.tolist():
277 277 if val >= beacon_dB[0]:
278 278 beacon_heiIndexList.append(avg_dB.tolist().index(val))
279 279
280 280 #data_spc = data_spc[:,:,beacon_heiIndexList]
281 281 data_cspc = None
282 282 if self.dataOut.data_cspc is not None:
283 283 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
284 284 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
285 285
286 286 data_dc = None
287 287 if self.dataOut.data_dc is not None:
288 288 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
289 289 #data_dc = data_dc[:,beacon_heiIndexList]
290 290
291 291 self.dataOut.data_spc = data_spc
292 292 self.dataOut.data_cspc = data_cspc
293 293 self.dataOut.data_dc = data_dc
294 294 self.dataOut.heightList = heightList
295 295 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
296 296
297 297 return 1
298 298
299 299 def selectFFTsByIndex(self, minIndex, maxIndex):
300 300 """
301 301
302 302 """
303 303
304 304 if (minIndex < 0) or (minIndex > maxIndex):
305 305 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
306 306
307 307 if (maxIndex >= self.dataOut.nProfiles):
308 308 maxIndex = self.dataOut.nProfiles-1
309 309
310 310 #Spectra
311 311 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
312 312
313 313 data_cspc = None
314 314 if self.dataOut.data_cspc is not None:
315 315 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
316 316
317 317 data_dc = None
318 318 if self.dataOut.data_dc is not None:
319 319 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
320 320
321 321 self.dataOut.data_spc = data_spc
322 322 self.dataOut.data_cspc = data_cspc
323 323 self.dataOut.data_dc = data_dc
324 324
325 325 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
326 326 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
327 327 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
328 328
329 329 return 1
330 330
331 331 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
332 332 # validacion de rango
333 333 if minHei == None:
334 334 minHei = self.dataOut.heightList[0]
335 335
336 336 if maxHei == None:
337 337 maxHei = self.dataOut.heightList[-1]
338 338
339 339 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
340 340 print('minHei: %.2f is out of the heights range' % (minHei))
341 341 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
342 342 minHei = self.dataOut.heightList[0]
343 343
344 344 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
345 345 print('maxHei: %.2f is out of the heights range' % (maxHei))
346 346 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
347 347 maxHei = self.dataOut.heightList[-1]
348 348
349 349 # validacion de velocidades
350 350 velrange = self.dataOut.getVelRange(1)
351 351
352 352 if minVel == None:
353 353 minVel = velrange[0]
354 354
355 355 if maxVel == None:
356 356 maxVel = velrange[-1]
357 357
358 358 if (minVel < velrange[0]) or (minVel > maxVel):
359 359 print('minVel: %.2f is out of the velocity range' % (minVel))
360 360 print('minVel is setting to %.2f' % (velrange[0]))
361 361 minVel = velrange[0]
362 362
363 363 if (maxVel > velrange[-1]) or (maxVel < minVel):
364 364 print('maxVel: %.2f is out of the velocity range' % (maxVel))
365 365 print('maxVel is setting to %.2f' % (velrange[-1]))
366 366 maxVel = velrange[-1]
367 367
368 368 # seleccion de indices para rango
369 369 minIndex = 0
370 370 maxIndex = 0
371 371 heights = self.dataOut.heightList
372 372
373 373 inda = numpy.where(heights >= minHei)
374 374 indb = numpy.where(heights <= maxHei)
375 375
376 376 try:
377 377 minIndex = inda[0][0]
378 378 except:
379 379 minIndex = 0
380 380
381 381 try:
382 382 maxIndex = indb[0][-1]
383 383 except:
384 384 maxIndex = len(heights)
385 385
386 386 if (minIndex < 0) or (minIndex > maxIndex):
387 387 raise ValueError("some value in (%d,%d) is not valid" % (
388 388 minIndex, maxIndex))
389 389
390 390 if (maxIndex >= self.dataOut.nHeights):
391 391 maxIndex = self.dataOut.nHeights - 1
392 392
393 393 # seleccion de indices para velocidades
394 394 indminvel = numpy.where(velrange >= minVel)
395 395 indmaxvel = numpy.where(velrange <= maxVel)
396 396 try:
397 397 minIndexVel = indminvel[0][0]
398 398 except:
399 399 minIndexVel = 0
400 400
401 401 try:
402 402 maxIndexVel = indmaxvel[0][-1]
403 403 except:
404 404 maxIndexVel = len(velrange)
405 405
406 406 # seleccion del espectro
407 407 data_spc = self.dataOut.data_spc[:,
408 408 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
409 409 # estimacion de ruido
410 410 noise = numpy.zeros(self.dataOut.nChannels)
411 411
412 412 for channel in range(self.dataOut.nChannels):
413 413 daux = data_spc[channel, :, :]
414 414 sortdata = numpy.sort(daux, axis=None)
415 415 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
416 416
417 417 self.dataOut.noise_estimation = noise.copy()
418 418
419 419 return 1
420 420
421 421 class removeDC(Operation):
422 422
423 423 def run(self, dataOut, mode=2):
424 424 self.dataOut = dataOut
425 425 jspectra = self.dataOut.data_spc
426 426 jcspectra = self.dataOut.data_cspc
427 427
428 428 num_chan = jspectra.shape[0]
429 429 num_hei = jspectra.shape[2]
430 430
431 431 if jcspectra is not None:
432 432 jcspectraExist = True
433 433 num_pairs = jcspectra.shape[0]
434 434 else:
435 435 jcspectraExist = False
436 436
437 437 freq_dc = int(jspectra.shape[1] / 2)
438 438 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
439 439 ind_vel = ind_vel.astype(int)
440 440
441 441 if ind_vel[0] < 0:
442 442 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
443 443
444 444 if mode == 1:
445 445 jspectra[:, freq_dc, :] = (
446 446 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
447 447
448 448 if jcspectraExist:
449 449 jcspectra[:, freq_dc, :] = (
450 450 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
451 451
452 452 if mode == 2:
453 453
454 454 vel = numpy.array([-2, -1, 1, 2])
455 455 xx = numpy.zeros([4, 4])
456 456
457 457 for fil in range(4):
458 458 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
459 459
460 460 xx_inv = numpy.linalg.inv(xx)
461 461 xx_aux = xx_inv[0, :]
462 462
463 463 for ich in range(num_chan):
464 464 yy = jspectra[ich, ind_vel, :]
465 465 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
466 466
467 467 junkid = jspectra[ich, freq_dc, :] <= 0
468 468 cjunkid = sum(junkid)
469 469
470 470 if cjunkid.any():
471 471 jspectra[ich, freq_dc, junkid.nonzero()] = (
472 472 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
473 473
474 474 if jcspectraExist:
475 475 for ip in range(num_pairs):
476 476 yy = jcspectra[ip, ind_vel, :]
477 477 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
478 478
479 479 self.dataOut.data_spc = jspectra
480 480 self.dataOut.data_cspc = jcspectra
481 481
482 482 return self.dataOut
483 483
484 484 # import matplotlib.pyplot as plt
485 485
486 486 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
487 487 z = (x - a1) / a2
488 488 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
489 489 return y
490 490 class CleanRayleigh(Operation):
491 491
492 492 def __init__(self):
493 493
494 494 Operation.__init__(self)
495 495 self.i=0
496 496 self.isConfig = False
497 497 self.__dataReady = False
498 498 self.__profIndex = 0
499 499 self.byTime = False
500 500 self.byProfiles = False
501 501
502 502 self.bloques = None
503 503 self.bloque0 = None
504 504
505 505 self.index = 0
506 506
507 507 self.buffer = 0
508 508 self.buffer2 = 0
509 509 self.buffer3 = 0
510 510
511 511
512 512 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
513 513
514 514 self.nChannels = dataOut.nChannels
515 515 self.nProf = dataOut.nProfiles
516 516 self.nPairs = dataOut.data_cspc.shape[0]
517 517 self.pairsArray = numpy.array(dataOut.pairsList)
518 518 self.spectra = dataOut.data_spc
519 519 self.cspectra = dataOut.data_cspc
520 520 self.heights = dataOut.heightList #alturas totales
521 521 self.nHeights = len(self.heights)
522 522 self.min_hei = min_hei
523 523 self.max_hei = max_hei
524 524 if (self.min_hei == None):
525 525 self.min_hei = 0
526 526 if (self.max_hei == None):
527 527 self.max_hei = dataOut.heightList[-1]
528 528 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
529 529 self.heightsClean = self.heights[self.hval] #alturas filtradas
530 530 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
531 531 self.nHeightsClean = len(self.heightsClean)
532 532 self.channels = dataOut.channelList
533 533 self.nChan = len(self.channels)
534 534 self.nIncohInt = dataOut.nIncohInt
535 535 self.__initime = dataOut.utctime
536 536 self.maxAltInd = self.hval[-1]+1
537 537 self.minAltInd = self.hval[0]
538 538
539 539 self.crosspairs = dataOut.pairsList
540 540 self.nPairs = len(self.crosspairs)
541 541 self.normFactor = dataOut.normFactor
542 542 self.nFFTPoints = dataOut.nFFTPoints
543 543 self.ippSeconds = dataOut.ippSeconds
544 544 self.currentTime = self.__initime
545 545 self.pairsArray = numpy.array(dataOut.pairsList)
546 546 self.factor_stdv = factor_stdv
547 print("CHANNELS: ",[x for x in self.channels])
547 #print("CHANNELS: ",[x for x in self.channels])
548 548
549 549 if n != None :
550 550 self.byProfiles = True
551 551 self.nIntProfiles = n
552 552 else:
553 553 self.__integrationtime = timeInterval
554 554
555 555 self.__dataReady = False
556 556 self.isConfig = True
557 557
558 558
559 559
560 560 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
561 561 #print (dataOut.utctime)
562 562 if not self.isConfig :
563 563 #print("Setting config")
564 564 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
565 565 #print("Config Done")
566 566 tini=dataOut.utctime
567 567
568 568 if self.byProfiles:
569 569 if self.__profIndex == self.nIntProfiles:
570 570 self.__dataReady = True
571 571 else:
572 572 if (tini - self.__initime) >= self.__integrationtime:
573 573 #print(tini - self.__initime,self.__profIndex)
574 574 self.__dataReady = True
575 575 self.__initime = tini
576 576
577 577 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
578 578
579 579 if self.__dataReady:
580 print("Data ready",self.__profIndex)
580 #print("Data ready",self.__profIndex)
581 581 self.__profIndex = 0
582 582 jspc = self.buffer
583 583 jcspc = self.buffer2
584 584 #jnoise = self.buffer3
585 585 self.buffer = dataOut.data_spc
586 586 self.buffer2 = dataOut.data_cspc
587 587 #self.buffer3 = dataOut.noise
588 588 self.currentTime = dataOut.utctime
589 589 if numpy.any(jspc) :
590 590 #print( jspc.shape, jcspc.shape)
591 591 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
592 592 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
593 593 self.__dataReady = False
594 594 #print( jspc.shape, jcspc.shape)
595 595 dataOut.flagNoData = False
596 596 else:
597 597 dataOut.flagNoData = True
598 598 self.__dataReady = False
599 599 return dataOut
600 600 else:
601 601 #print( len(self.buffer))
602 602 if numpy.any(self.buffer):
603 603 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
604 604 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
605 605 self.buffer3 += dataOut.data_dc
606 606 else:
607 607 self.buffer = dataOut.data_spc
608 608 self.buffer2 = dataOut.data_cspc
609 609 self.buffer3 = dataOut.data_dc
610 610 #print self.index, self.fint
611 611 #print self.buffer2.shape
612 612 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
613 613 self.__profIndex += 1
614 614 return dataOut ## NOTE: REV
615 615
616 616
617 617 #index = tini.tm_hour*12+tini.tm_min/5
618 618 '''REVISAR'''
619 619 # jspc = jspc/self.nFFTPoints/self.normFactor
620 620 # jcspc = jcspc/self.nFFTPoints/self.normFactor
621 621
622 622
623 #dataOut.data_spc,dataOut.data_cspc = self.CleanRayleigh(dataOut,jspc,jcspc,crosspairs,heights,channels,nProf,nHei,nChan,nPairs,nIncohInt,nBlocks=nBlocks)
624 #tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.min_hei,self.max_hei)
625 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
626 #jspectra = tmp_spectra*len(jspc[:,0,0,0])
627 #jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
628 623
624 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
629 625 dataOut.data_spc = tmp_spectra
630 626 dataOut.data_cspc = tmp_cspectra
627
628 #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
629
631 630 dataOut.data_dc = self.buffer3
632 631 dataOut.nIncohInt *= self.nIntProfiles
633 632 dataOut.utctime = self.currentTime #tiempo promediado
634 633 #print("Time: ",time.localtime(dataOut.utctime))
635 634 # dataOut.data_spc = sat_spectra
636 635 # dataOut.data_cspc = sat_cspectra
637 636 self.buffer = 0
638 637 self.buffer2 = 0
639 638 self.buffer3 = 0
640 639
641 640 return dataOut
642 641
643 642 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
644 print("OP cleanRayleigh")
643 #print("OP cleanRayleigh")
645 644 #import matplotlib.pyplot as plt
646 645 #for k in range(149):
647 646
648 647 rfunc = cspectra.copy() #self.bloques
649 val_spc = spectra*0.0 #self.bloque0*0.0
650 val_cspc = cspectra*0.0 #self.bloques*0.0
651 in_sat_spectra = spectra.copy() #self.bloque0
652 in_sat_cspectra = cspectra.copy() #self.bloques
648 #rfunc = cspectra
649 #val_spc = spectra*0.0 #self.bloque0*0.0
650 #val_cspc = cspectra*0.0 #self.bloques*0.0
651 #in_sat_spectra = spectra.copy() #self.bloque0
652 #in_sat_cspectra = cspectra.copy() #self.bloques
653 653
654 raxs = math.ceil(math.sqrt(self.nPairs))
655 caxs = math.ceil(self.nPairs/raxs)
654 #raxs = math.ceil(math.sqrt(self.nPairs))
655 #caxs = math.ceil(self.nPairs/raxs)
656 656
657 657 #print(self.hval)
658 658 #print numpy.absolute(rfunc[:,0,0,14])
659 gauss_fit, covariance = None, None
659 660 for ih in range(self.minAltInd,self.maxAltInd):
660 661 for ifreq in range(self.nFFTPoints):
661 662 # fig, axs = plt.subplots(raxs, caxs)
662 663 # fig2, axs2 = plt.subplots(raxs, caxs)
663 col_ax = 0
664 row_ax = 0
664 # col_ax = 0
665 # row_ax = 0
666 #print(len(self.nPairs))
665 667 for ii in range(self.nPairs): #PARES DE CANALES SELF y CROSS
666 668 #print("ii: ",ii)
667 if (col_ax%caxs==0 and col_ax!=0):
668 col_ax = 0
669 row_ax += 1
669 # if (col_ax%caxs==0 and col_ax!=0):
670 # col_ax = 0
671 # row_ax += 1
670 672 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
671 673 #print(func2clean.shape)
672 674 val = (numpy.isfinite(func2clean)==True).nonzero()
673 675
674 676 if len(val)>0: #limitador
675 677 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
676 678 if min_val <= -40 :
677 679 min_val = -40
678 680 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
679 681 if max_val >= 200 :
680 682 max_val = 200
681 683 #print min_val, max_val
682 684 step = 1
683 685 #print("Getting bins and the histogram")
684 686 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
685 687 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
686 688 #print(len(y_dist),len(binstep[:-1]))
687 689 #print(row_ax,col_ax, " ..")
688 690 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
689 691 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
690 692 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
691 693 parg = [numpy.amax(y_dist),mean,sigma]
692 gauss_fit, covariance = None, None
693 newY = None
694
695 #newY = None
696
694 697 try :
695 698 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
696 699 mode = gauss_fit[1]
697 700 stdv = gauss_fit[2]
698 701 #print(" FIT OK",gauss_fit)
699 702 '''
700 703 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
701 704 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
702 705 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
703 706 axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))'''
704 707 except:
705 708 mode = mean
706 709 stdv = sigma
707 710 #print("FIT FAIL")
708 711
709 712
710 713 #print(mode,stdv)
711 #Removing echoes greater than mode + 3*stdv
712 #factor_stdv = 2
714 #Removing echoes greater than mode + std_factor*stdv
713 715 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
714 716 #noval tiene los indices que se van a remover
715 717 #print("Pair ",ii," novals: ",len(noval[0]))
716 718 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
717 719 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
718 720 #print(novall)
719 721 #print(" ",self.pairsArray[ii])
720 722 cross_pairs = self.pairsArray[ii]
721 723 #Getting coherent echoes which are removed.
722 724 # if len(novall[0]) > 0:
723 725 #
724 726 # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
725 727 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
726 728 # val_cspc[novall[0],ii,ifreq,ih] = 1
727 729 #print("OUT NOVALL 1")
728 730 #Removing coherent from ISR data
729 731 chA = self.channels.index(cross_pairs[0])
730 732 chB = self.channels.index(cross_pairs[1])
731 733
732 734 new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0])
733 mean_cspc = numpy.mean(new_a)
735 cspectra[noval,ii,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
734 736 new_b = numpy.delete(spectra[:,chA,ifreq,ih], noval[0])
735 mean_spc0 = numpy.mean(new_b)
737 spectra[noval,chA,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
736 738 new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
737 mean_spc1 = numpy.mean(new_c)
738 spectra[noval,chA,ifreq,ih] = mean_spc0
739 spectra[noval,chB,ifreq,ih] = mean_spc1
740 cspectra[noval,ii,ifreq,ih] = mean_cspc
739 spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
740
741 741
742 742 '''
743 743 func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
744 744 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
745 745 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
746 746 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
747 747 axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
748 748 '''
749 749
750 col_ax += 1 #contador de ploteo columnas
750 #col_ax += 1 #contador de ploteo columnas
751 751 ##print(col_ax)
752 752 '''
753 753 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
754 754 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
755 755 fig.suptitle(title)
756 756 fig2.suptitle(title2)
757 757 plt.show()'''
758 758
759 759 ''' channels = channels
760 760 cross_pairs = cross_pairs
761 761 #print("OUT NOVALL 2")
762 762
763 763 vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
764 764 vcross1 = (cross_pairs[1] == channels[ii]).nonzero()
765 765 vcross = numpy.concatenate((vcross0,vcross1),axis=None)
766 766 #print('vcros =', vcross)
767 767
768 768 #Getting coherent echoes which are removed.
769 769 if len(novall) > 0:
770 770 #val_spc[novall,ii,ifreq,ih] = 1
771 771 val_spc[ii,ifreq,ih,novall] = 1
772 772 if len(vcross) > 0:
773 773 val_cspc[vcross,ifreq,ih,novall] = 1
774 774
775 775 #Removing coherent from ISR data.
776 776 self.bloque0[ii,ifreq,ih,noval] = numpy.nan
777 777 if len(vcross) > 0:
778 778 self.bloques[vcross,ifreq,ih,noval] = numpy.nan
779 779 '''
780 780
781 print("Getting average of the spectra and cross-spectra from incoherent echoes.")
781 #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
782 782 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
783 783 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
784 784 for ih in range(self.nHeights):
785 785 for ifreq in range(self.nFFTPoints):
786 786 for ich in range(self.nChan):
787 787 tmp = spectra[:,ich,ifreq,ih]
788 788 valid = (numpy.isfinite(tmp[:])==True).nonzero()
789 789 # if ich == 0 and ifreq == 0 and ih == 17 :
790 790 # print tmp
791 791 # print valid
792 792 # print len(valid[0])
793 793 #print('TMP',tmp)
794 794 if len(valid[0]) >0 :
795 795 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
796 796 #for icr in range(nPairs):
797 797 for icr in range(self.nPairs):
798 798 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
799 799 valid = (numpy.isfinite(tmp)==True).nonzero()
800 800 if len(valid[0]) > 0:
801 801 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
802 802 '''
803 803 # print('##########################################################')
804 804 print("Removing fake coherent echoes (at least 4 points around the point)")
805 805
806 806 val_spectra = numpy.sum(val_spc,0)
807 807 val_cspectra = numpy.sum(val_cspc,0)
808 808
809 809 val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
810 810 val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
811 811
812 812 for i in range(nChan):
813 813 for j in range(nProf):
814 814 for k in range(nHeights):
815 815 if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 :
816 816 val_spc[:,i,j,k] = 0.0
817 817 for i in range(nPairs):
818 818 for j in range(nProf):
819 819 for k in range(nHeights):
820 820 if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
821 821 val_cspc[:,i,j,k] = 0.0
822 822
823 823 # val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProf*nHeights*nChan))
824 824 # if numpy.isfinite(val_spectra)==str(True):
825 825 # noval = (val_spectra<1).nonzero()
826 826 # if len(noval) > 0:
827 827 # val_spc[:,noval] = 0.0
828 828 # val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHeights))
829 829
830 830 #val_cspc = numpy.reshape(val_spc, (149,nChan*nHeights*nProf))
831 831 #if numpy.isfinite(val_cspectra)==str(True):
832 832 # noval = (val_cspectra<1).nonzero()
833 833 # if len(noval) > 0:
834 834 # val_cspc[:,noval] = 0.0
835 835 # val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHeights))
836 836 tmp_sat_spectra = spectra.copy()
837 837 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
838 838 tmp_sat_cspectra = cspectra.copy()
839 839 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
840 840 '''
841 841 # fig = plt.figure(figsize=(6,5))
842 842 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
843 843 # ax = fig.add_axes([left, bottom, width, height])
844 844 # cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:])))
845 845 # ax.clabel(cp, inline=True,fontsize=10)
846 846 # plt.show()
847 847 '''
848 848 val = (val_spc > 0).nonzero()
849 849 if len(val[0]) > 0:
850 850 tmp_sat_spectra[val] = in_sat_spectra[val]
851 851 val = (val_cspc > 0).nonzero()
852 852 if len(val[0]) > 0:
853 853 tmp_sat_cspectra[val] = in_sat_cspectra[val]
854 854
855 855 print("Getting average of the spectra and cross-spectra from incoherent echoes 2")
856 856 sat_spectra = numpy.zeros((nChan,nProf,nHeights), dtype=float)
857 857 sat_cspectra = numpy.zeros((nPairs,nProf,nHeights), dtype=complex)
858 858 for ih in range(nHeights):
859 859 for ifreq in range(nProf):
860 860 for ich in range(nChan):
861 861 tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
862 862 valid = (numpy.isfinite(tmp)).nonzero()
863 863 if len(valid[0]) > 0:
864 864 sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
865 865
866 866 for icr in range(nPairs):
867 867 tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih])
868 868 valid = (numpy.isfinite(tmp)).nonzero()
869 869 if len(valid[0]) > 0:
870 870 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
871 871 '''
872 872 #self.__dataReady= True
873 873 #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra
874 874 #if not self.__dataReady:
875 875 #return None, None
876 876 #return out_spectra, out_cspectra ,sat_spectra,sat_cspectra
877 877 return out_spectra, out_cspectra
878 878
879 879 def REM_ISOLATED_POINTS(self,array,rth):
880 880 # import matplotlib.pyplot as plt
881 881 if rth == None :
882 882 rth = 4
883 883 print("REM ISO")
884 884 num_prof = len(array[0,:,0])
885 885 num_hei = len(array[0,0,:])
886 886 n2d = len(array[:,0,0])
887 887
888 888 for ii in range(n2d) :
889 889 #print ii,n2d
890 890 tmp = array[ii,:,:]
891 891 #print tmp.shape, array[ii,101,:],array[ii,102,:]
892 892
893 893 # fig = plt.figure(figsize=(6,5))
894 894 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
895 895 # ax = fig.add_axes([left, bottom, width, height])
896 896 # x = range(num_prof)
897 897 # y = range(num_hei)
898 898 # cp = ax.contour(y,x,tmp)
899 899 # ax.clabel(cp, inline=True,fontsize=10)
900 900 # plt.show()
901 901
902 902 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
903 903 tmp = numpy.reshape(tmp,num_prof*num_hei)
904 904 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
905 905 indxs2 = (tmp > 0).nonzero()
906 906
907 907 indxs1 = (indxs1[0])
908 908 indxs2 = indxs2[0]
909 909 #indxs1 = numpy.array(indxs1[0])
910 910 #indxs2 = numpy.array(indxs2[0])
911 911 indxs = None
912 912 #print indxs1 , indxs2
913 913 for iv in range(len(indxs2)):
914 914 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
915 915 #print len(indxs2), indv
916 916 if len(indv[0]) > 0 :
917 917 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
918 918 # print indxs
919 919 indxs = indxs[1:]
920 920 #print(indxs, len(indxs))
921 921 if len(indxs) < 4 :
922 922 array[ii,:,:] = 0.
923 923 return
924 924
925 925 xpos = numpy.mod(indxs ,num_hei)
926 926 ypos = (indxs / num_hei)
927 927 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
928 928 #print sx
929 929 xpos = xpos[sx]
930 930 ypos = ypos[sx]
931 931
932 932 # *********************************** Cleaning isolated points **********************************
933 933 ic = 0
934 934 while True :
935 935 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
936 936 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
937 937 #plt.plot(r)
938 938 #plt.show()
939 939 no_coh1 = (numpy.isfinite(r)==True).nonzero()
940 940 no_coh2 = (r <= rth).nonzero()
941 941 #print r, no_coh1, no_coh2
942 942 no_coh1 = numpy.array(no_coh1[0])
943 943 no_coh2 = numpy.array(no_coh2[0])
944 944 no_coh = None
945 945 #print valid1 , valid2
946 946 for iv in range(len(no_coh2)):
947 947 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
948 948 if len(indv[0]) > 0 :
949 949 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
950 950 no_coh = no_coh[1:]
951 951 #print len(no_coh), no_coh
952 952 if len(no_coh) < 4 :
953 953 #print xpos[ic], ypos[ic], ic
954 954 # plt.plot(r)
955 955 # plt.show()
956 956 xpos[ic] = numpy.nan
957 957 ypos[ic] = numpy.nan
958 958
959 959 ic = ic + 1
960 960 if (ic == len(indxs)) :
961 961 break
962 962 #print( xpos, ypos)
963 963
964 964 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
965 965 #print indxs[0]
966 966 if len(indxs[0]) < 4 :
967 967 array[ii,:,:] = 0.
968 968 return
969 969
970 970 xpos = xpos[indxs[0]]
971 971 ypos = ypos[indxs[0]]
972 972 for i in range(0,len(ypos)):
973 973 ypos[i]=int(ypos[i])
974 974 junk = tmp
975 975 tmp = junk*0.0
976 976
977 977 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
978 978 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
979 979
980 980 #print array.shape
981 981 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
982 982 #print tmp.shape
983 983
984 984 # fig = plt.figure(figsize=(6,5))
985 985 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
986 986 # ax = fig.add_axes([left, bottom, width, height])
987 987 # x = range(num_prof)
988 988 # y = range(num_hei)
989 989 # cp = ax.contour(y,x,array[ii,:,:])
990 990 # ax.clabel(cp, inline=True,fontsize=10)
991 991 # plt.show()
992 992 return array
993 993
994 994 class removeInterference(Operation):
995 995
996 996 def removeInterference2(self):
997 997
998 998 cspc = self.dataOut.data_cspc
999 999 spc = self.dataOut.data_spc
1000 1000 Heights = numpy.arange(cspc.shape[2])
1001 1001 realCspc = numpy.abs(cspc)
1002 1002
1003 1003 for i in range(cspc.shape[0]):
1004 1004 LinePower= numpy.sum(realCspc[i], axis=0)
1005 1005 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1006 1006 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1007 1007 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1008 1008 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1009 1009 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1010 1010
1011 1011
1012 1012 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1013 1013 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1014 1014 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1015 1015 cspc[i,InterferenceRange,:] = numpy.NaN
1016 1016
1017 1017 self.dataOut.data_cspc = cspc
1018 1018
1019 1019 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1020 1020
1021 1021 jspectra = self.dataOut.data_spc
1022 1022 jcspectra = self.dataOut.data_cspc
1023 1023 jnoise = self.dataOut.getNoise()
1024 1024 num_incoh = self.dataOut.nIncohInt
1025 1025
1026 1026 num_channel = jspectra.shape[0]
1027 1027 num_prof = jspectra.shape[1]
1028 1028 num_hei = jspectra.shape[2]
1029 1029
1030 1030 # hei_interf
1031 1031 if hei_interf is None:
1032 1032 count_hei = int(num_hei / 2)
1033 1033 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1034 1034 hei_interf = numpy.asarray(hei_interf)[0]
1035 1035 # nhei_interf
1036 1036 if (nhei_interf == None):
1037 1037 nhei_interf = 5
1038 1038 if (nhei_interf < 1):
1039 1039 nhei_interf = 1
1040 1040 if (nhei_interf > count_hei):
1041 1041 nhei_interf = count_hei
1042 1042 if (offhei_interf == None):
1043 1043 offhei_interf = 0
1044 1044
1045 1045 ind_hei = list(range(num_hei))
1046 1046 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1047 1047 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1048 1048 mask_prof = numpy.asarray(list(range(num_prof)))
1049 1049 num_mask_prof = mask_prof.size
1050 1050 comp_mask_prof = [0, num_prof / 2]
1051 1051
1052 1052 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1053 1053 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1054 1054 jnoise = numpy.nan
1055 1055 noise_exist = jnoise[0] < numpy.Inf
1056 1056
1057 1057 # Subrutina de Remocion de la Interferencia
1058 1058 for ich in range(num_channel):
1059 1059 # Se ordena los espectros segun su potencia (menor a mayor)
1060 1060 power = jspectra[ich, mask_prof, :]
1061 1061 power = power[:, hei_interf]
1062 1062 power = power.sum(axis=0)
1063 1063 psort = power.ravel().argsort()
1064 1064
1065 1065 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1066 1066 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1067 1067 offhei_interf, nhei_interf + offhei_interf))]]]
1068 1068
1069 1069 if noise_exist:
1070 1070 # tmp_noise = jnoise[ich] / num_prof
1071 1071 tmp_noise = jnoise[ich]
1072 1072 junkspc_interf = junkspc_interf - tmp_noise
1073 1073 #junkspc_interf[:,comp_mask_prof] = 0
1074 1074
1075 1075 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1076 1076 jspc_interf = jspc_interf.transpose()
1077 1077 # Calculando el espectro de interferencia promedio
1078 1078 noiseid = numpy.where(
1079 1079 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1080 1080 noiseid = noiseid[0]
1081 1081 cnoiseid = noiseid.size
1082 1082 interfid = numpy.where(
1083 1083 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1084 1084 interfid = interfid[0]
1085 1085 cinterfid = interfid.size
1086 1086
1087 1087 if (cnoiseid > 0):
1088 1088 jspc_interf[noiseid] = 0
1089 1089
1090 1090 # Expandiendo los perfiles a limpiar
1091 1091 if (cinterfid > 0):
1092 1092 new_interfid = (
1093 1093 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1094 1094 new_interfid = numpy.asarray(new_interfid)
1095 1095 new_interfid = {x for x in new_interfid}
1096 1096 new_interfid = numpy.array(list(new_interfid))
1097 1097 new_cinterfid = new_interfid.size
1098 1098 else:
1099 1099 new_cinterfid = 0
1100 1100
1101 1101 for ip in range(new_cinterfid):
1102 1102 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1103 1103 jspc_interf[new_interfid[ip]
1104 1104 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1105 1105
1106 1106 jspectra[ich, :, ind_hei] = jspectra[ich, :,
1107 1107 ind_hei] - jspc_interf # Corregir indices
1108 1108
1109 1109 # Removiendo la interferencia del punto de mayor interferencia
1110 1110 ListAux = jspc_interf[mask_prof].tolist()
1111 1111 maxid = ListAux.index(max(ListAux))
1112 1112
1113 1113 if cinterfid > 0:
1114 1114 for ip in range(cinterfid * (interf == 2) - 1):
1115 1115 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1116 1116 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1117 1117 cind = len(ind)
1118 1118
1119 1119 if (cind > 0):
1120 1120 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1121 1121 (1 + (numpy.random.uniform(cind) - 0.5) /
1122 1122 numpy.sqrt(num_incoh))
1123 1123
1124 1124 ind = numpy.array([-2, -1, 1, 2])
1125 1125 xx = numpy.zeros([4, 4])
1126 1126
1127 1127 for id1 in range(4):
1128 1128 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1129 1129
1130 1130 xx_inv = numpy.linalg.inv(xx)
1131 1131 xx = xx_inv[:, 0]
1132 1132 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1133 1133 yy = jspectra[ich, mask_prof[ind], :]
1134 1134 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
1135 1135 yy.transpose(), xx)
1136 1136
1137 1137 indAux = (jspectra[ich, :, :] < tmp_noise *
1138 1138 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1139 1139 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1140 1140 (1 - 1 / numpy.sqrt(num_incoh))
1141 1141
1142 1142 # Remocion de Interferencia en el Cross Spectra
1143 1143 if jcspectra is None:
1144 1144 return jspectra, jcspectra
1145 1145 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1146 1146 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1147 1147
1148 1148 for ip in range(num_pairs):
1149 1149
1150 1150 #-------------------------------------------
1151 1151
1152 1152 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1153 1153 cspower = cspower[:, hei_interf]
1154 1154 cspower = cspower.sum(axis=0)
1155 1155
1156 1156 cspsort = cspower.ravel().argsort()
1157 1157 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1158 1158 offhei_interf, nhei_interf + offhei_interf))]]]
1159 1159 junkcspc_interf = junkcspc_interf.transpose()
1160 1160 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1161 1161
1162 1162 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1163 1163
1164 1164 median_real = int(numpy.median(numpy.real(
1165 1165 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1166 1166 median_imag = int(numpy.median(numpy.imag(
1167 1167 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1168 1168 comp_mask_prof = [int(e) for e in comp_mask_prof]
1169 1169 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1170 1170 median_real, median_imag)
1171 1171
1172 1172 for iprof in range(num_prof):
1173 1173 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1174 1174 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1175 1175
1176 1176 # Removiendo la Interferencia
1177 1177 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1178 1178 :, ind_hei] - jcspc_interf
1179 1179
1180 1180 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1181 1181 maxid = ListAux.index(max(ListAux))
1182 1182
1183 1183 ind = numpy.array([-2, -1, 1, 2])
1184 1184 xx = numpy.zeros([4, 4])
1185 1185
1186 1186 for id1 in range(4):
1187 1187 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1188 1188
1189 1189 xx_inv = numpy.linalg.inv(xx)
1190 1190 xx = xx_inv[:, 0]
1191 1191
1192 1192 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1193 1193 yy = jcspectra[ip, mask_prof[ind], :]
1194 1194 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1195 1195
1196 1196 # Guardar Resultados
1197 1197 self.dataOut.data_spc = jspectra
1198 1198 self.dataOut.data_cspc = jcspectra
1199 1199
1200 1200 return 1
1201 1201
1202 1202 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
1203 1203
1204 1204 self.dataOut = dataOut
1205 1205
1206 1206 if mode == 1:
1207 1207 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1208 1208 elif mode == 2:
1209 1209 self.removeInterference2()
1210 1210
1211 1211 return self.dataOut
1212 1212
1213 1213
1214 1214 class IncohInt(Operation):
1215 1215
1216 1216 __profIndex = 0
1217 1217 __withOverapping = False
1218 1218
1219 1219 __byTime = False
1220 1220 __initime = None
1221 1221 __lastdatatime = None
1222 1222 __integrationtime = None
1223 1223
1224 1224 __buffer_spc = None
1225 1225 __buffer_cspc = None
1226 1226 __buffer_dc = None
1227 1227
1228 1228 __dataReady = False
1229 1229
1230 1230 __timeInterval = None
1231 1231
1232 1232 n = None
1233 1233
1234 1234 def __init__(self):
1235 1235
1236 1236 Operation.__init__(self)
1237 1237
1238 1238 def setup(self, n=None, timeInterval=None, overlapping=False):
1239 1239 """
1240 1240 Set the parameters of the integration class.
1241 1241
1242 1242 Inputs:
1243 1243
1244 1244 n : Number of coherent integrations
1245 1245 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1246 1246 overlapping :
1247 1247
1248 1248 """
1249 1249
1250 1250 self.__initime = None
1251 1251 self.__lastdatatime = 0
1252 1252
1253 1253 self.__buffer_spc = 0
1254 1254 self.__buffer_cspc = 0
1255 1255 self.__buffer_dc = 0
1256 1256
1257 1257 self.__profIndex = 0
1258 1258 self.__dataReady = False
1259 1259 self.__byTime = False
1260 1260
1261 1261 if n is None and timeInterval is None:
1262 1262 raise ValueError("n or timeInterval should be specified ...")
1263 1263
1264 1264 if n is not None:
1265 1265 self.n = int(n)
1266 1266 else:
1267 1267
1268 1268 self.__integrationtime = int(timeInterval)
1269 1269 self.n = None
1270 1270 self.__byTime = True
1271 1271
1272 1272 def putData(self, data_spc, data_cspc, data_dc):
1273 1273 """
1274 1274 Add a profile to the __buffer_spc and increase in one the __profileIndex
1275 1275
1276 1276 """
1277 1277
1278 1278 self.__buffer_spc += data_spc
1279 1279
1280 1280 if data_cspc is None:
1281 1281 self.__buffer_cspc = None
1282 1282 else:
1283 1283 self.__buffer_cspc += data_cspc
1284 1284
1285 1285 if data_dc is None:
1286 1286 self.__buffer_dc = None
1287 1287 else:
1288 1288 self.__buffer_dc += data_dc
1289 1289
1290 1290 self.__profIndex += 1
1291 1291
1292 1292 return
1293 1293
1294 1294 def pushData(self):
1295 1295 """
1296 1296 Return the sum of the last profiles and the profiles used in the sum.
1297 1297
1298 1298 Affected:
1299 1299
1300 1300 self.__profileIndex
1301 1301
1302 1302 """
1303 1303
1304 1304 data_spc = self.__buffer_spc
1305 1305 data_cspc = self.__buffer_cspc
1306 1306 data_dc = self.__buffer_dc
1307 1307 n = self.__profIndex
1308 1308
1309 1309 self.__buffer_spc = 0
1310 1310 self.__buffer_cspc = 0
1311 1311 self.__buffer_dc = 0
1312 1312 self.__profIndex = 0
1313 1313
1314 1314 return data_spc, data_cspc, data_dc, n
1315 1315
1316 1316 def byProfiles(self, *args):
1317 1317
1318 1318 self.__dataReady = False
1319 1319 avgdata_spc = None
1320 1320 avgdata_cspc = None
1321 1321 avgdata_dc = None
1322 1322
1323 1323 self.putData(*args)
1324 1324
1325 1325 if self.__profIndex == self.n:
1326 1326
1327 1327 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1328 1328 self.n = n
1329 1329 self.__dataReady = True
1330 1330
1331 1331 return avgdata_spc, avgdata_cspc, avgdata_dc
1332 1332
1333 1333 def byTime(self, datatime, *args):
1334 1334
1335 1335 self.__dataReady = False
1336 1336 avgdata_spc = None
1337 1337 avgdata_cspc = None
1338 1338 avgdata_dc = None
1339 1339
1340 1340 self.putData(*args)
1341 1341
1342 1342 if (datatime - self.__initime) >= self.__integrationtime:
1343 1343 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1344 1344 self.n = n
1345 1345 self.__dataReady = True
1346 1346
1347 1347 return avgdata_spc, avgdata_cspc, avgdata_dc
1348 1348
1349 1349 def integrate(self, datatime, *args):
1350 1350
1351 1351 if self.__profIndex == 0:
1352 1352 self.__initime = datatime
1353 1353
1354 1354 if self.__byTime:
1355 1355 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1356 1356 datatime, *args)
1357 1357 else:
1358 1358 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1359 1359
1360 1360 if not self.__dataReady:
1361 1361 return None, None, None, None
1362 1362
1363 1363 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1364 1364
1365 1365 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1366 1366 if n == 1:
1367 1367 return dataOut
1368 1368
1369 1369 dataOut.flagNoData = True
1370 1370
1371 1371 if not self.isConfig:
1372 1372 self.setup(n, timeInterval, overlapping)
1373 1373 self.isConfig = True
1374 1374
1375 1375 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1376 1376 dataOut.data_spc,
1377 1377 dataOut.data_cspc,
1378 1378 dataOut.data_dc)
1379 1379
1380 1380 if self.__dataReady:
1381 1381
1382 1382 dataOut.data_spc = avgdata_spc
1383 1383 dataOut.data_cspc = avgdata_cspc
1384 1384 dataOut.data_dc = avgdata_dc
1385 1385 dataOut.nIncohInt *= self.n
1386 1386 dataOut.utctime = avgdatatime
1387 1387 dataOut.flagNoData = False
1388 1388
1389 1389 return dataOut
1390 1390
1391 1391 class dopplerFlip(Operation):
1392 1392
1393 1393 def run(self, dataOut):
1394 1394 # arreglo 1: (num_chan, num_profiles, num_heights)
1395 1395 self.dataOut = dataOut
1396 1396 # JULIA-oblicua, indice 2
1397 1397 # arreglo 2: (num_profiles, num_heights)
1398 1398 jspectra = self.dataOut.data_spc[2]
1399 1399 jspectra_tmp = numpy.zeros(jspectra.shape)
1400 1400 num_profiles = jspectra.shape[0]
1401 1401 freq_dc = int(num_profiles / 2)
1402 1402 # Flip con for
1403 1403 for j in range(num_profiles):
1404 1404 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1405 1405 # Intercambio perfil de DC con perfil inmediato anterior
1406 1406 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1407 1407 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1408 1408 # canal modificado es re-escrito en el arreglo de canales
1409 1409 self.dataOut.data_spc[2] = jspectra_tmp
1410 1410
1411 1411 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now