##// END OF EJS Templates
schain
RSS
upgrade filterbyHeights, fix initial parameters. removeProfileSats2 improved debug
joabAM -
r1655:13a9a26f11d7
parent child
Show More
Side by Side Unified
Context file:
r1655:13a9a26f11d7 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,705 +1,707
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 if self.t_units == "h_m":
195 195 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
196 196 if self.t_units == "h":
197 197 return '{}'.format(self.getDateTime(x).strftime('%H'))
198 198
199 199 def __setup(self, **kwargs):
200 200 '''
201 201 Initialize variables
202 202 '''
203 203
204 204 self.figures = []
205 205 self.axes = []
206 206 self.cb_axes = []
207 207 self.pf_axes = []
208 208 self.localtime = kwargs.pop('localtime', True)
209 209 self.show = kwargs.get('show', True)
210 210 self.save = kwargs.get('save', False)
211 211 self.save_period = kwargs.get('save_period', 0)
212 212 self.colormap = kwargs.get('colormap', self.colormap)
213 213 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
214 214 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
215 215 self.colormaps = kwargs.get('colormaps', None)
216 216 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
217 217 self.showprofile = kwargs.get('showprofile', False)
218 218 self.title = kwargs.get('wintitle', self.CODE.upper())
219 219 self.cb_label = kwargs.get('cb_label', None)
220 220 self.cb_labels = kwargs.get('cb_labels', None)
221 221 self.labels = kwargs.get('labels', None)
222 222 self.xaxis = kwargs.get('xaxis', 'frequency')
223 223 self.zmin = kwargs.get('zmin', None)
224 224 self.zmax = kwargs.get('zmax', None)
225 225 self.zlimits = kwargs.get('zlimits', None)
226 226 self.xmin = kwargs.get('xmin', None)
227 227 self.xmax = kwargs.get('xmax', None)
228 228 self.xrange = kwargs.get('xrange', 12)
229 229 self.xscale = kwargs.get('xscale', None)
230 230 self.ymin = kwargs.get('ymin', None)
231 231 self.ymax = kwargs.get('ymax', None)
232 232 self.yscale = kwargs.get('yscale', None)
233 233 self.xlabel = kwargs.get('xlabel', None)
234 234 self.attr_time = kwargs.get('attr_time', 'utctime')
235 235 self.attr_data = kwargs.get('attr_data', 'data_param')
236 236 self.decimation = kwargs.get('decimation', None)
237 237 self.oneFigure = kwargs.get('oneFigure', True)
238 238 self.width = kwargs.get('width', None)
239 239 self.height = kwargs.get('height', None)
240 240 self.colorbar = kwargs.get('colorbar', True)
241 241 self.factors = kwargs.get('factors', range(18))
242 242 self.channels = kwargs.get('channels', None)
243 243 self.titles = kwargs.get('titles', [])
244 244 self.polar = False
245 245 self.type = kwargs.get('type', 'iq')
246 246 self.grid = kwargs.get('grid', False)
247 247 self.pause = kwargs.get('pause', False)
248 248 self.save_code = kwargs.get('save_code', self.CODE)
249 249 self.throttle = kwargs.get('throttle', 0)
250 250 self.exp_code = kwargs.get('exp_code', None)
251 251 self.server = kwargs.get('server', False)
252 252 self.sender_period = kwargs.get('sender_period', 60)
253 253 self.tag = kwargs.get('tag', '')
254 254 self.height_index = kwargs.get('height_index', [])
255 255 self.__throttle_plot = apply_throttle(self.throttle)
256 256 code = self.attr_data if self.attr_data else self.CODE
257 257 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
258 258 self.tmin = kwargs.get('tmin', None)
259 259 self.t_units = kwargs.get('t_units', "h_m")
260 260 self.selectedHeightsList = kwargs.get('selectedHeightsList', [])
261 261 if isinstance(self.selectedHeightsList, int):
262 262 self.selectedHeightsList = [self.selectedHeightsList]
263 263
264 264 if self.server:
265 265 if not self.server.startswith('tcp://'):
266 266 self.server = 'tcp://{}'.format(self.server)
267 267 log.success(
268 268 'Sending to server: {}'.format(self.server),
269 269 self.name
270 270 )
271 271
272 272 if isinstance(self.attr_data, str):
273 273 self.attr_data = [self.attr_data]
274 274
275 275 def __setup_plot(self):
276 276 '''
277 277 Common setup for all figures, here figures and axes are created
278 278 '''
279 279
280 280 self.setup()
281 281
282 282 self.time_label = 'LT' if self.localtime else 'UTC'
283 283
284 284 if self.width is None:
285 285 self.width = 8
286 286
287 287 self.figures = []
288 288 self.axes = []
289 289 self.cb_axes = []
290 290 self.pf_axes = []
291 291 self.cmaps = []
292 292
293 293 size = '15%' if self.ncols == 1 else '30%'
294 294 pad = '4%' if self.ncols == 1 else '8%'
295 295
296 296 if self.oneFigure:
297 297 if self.height is None:
298 298 self.height = 1.4 * self.nrows + 1
299 299 fig = plt.figure(figsize=(self.width, self.height),
300 300 edgecolor='k',
301 301 facecolor='w')
302 302 self.figures.append(fig)
303 303 for n in range(self.nplots):
304 304 ax = fig.add_subplot(self.nrows, self.ncols,
305 305 n + 1, polar=self.polar)
306 306 ax.tick_params(labelsize=8)
307 307 ax.firsttime = True
308 308 ax.index = 0
309 309 ax.press = None
310 310 self.axes.append(ax)
311 311 if self.showprofile:
312 312 cax = self.__add_axes(ax, size=size, pad=pad)
313 313 cax.tick_params(labelsize=8)
314 314 self.pf_axes.append(cax)
315 315 else:
316 316 if self.height is None:
317 317 self.height = 3
318 318 for n in range(self.nplots):
319 319 fig = plt.figure(figsize=(self.width, self.height),
320 320 edgecolor='k',
321 321 facecolor='w')
322 322 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
323 323 ax.tick_params(labelsize=8)
324 324 ax.firsttime = True
325 325 ax.index = 0
326 326 ax.press = None
327 327 self.figures.append(fig)
328 328 self.axes.append(ax)
329 329 if self.showprofile:
330 330 cax = self.__add_axes(ax, size=size, pad=pad)
331 331 cax.tick_params(labelsize=8)
332 332 self.pf_axes.append(cax)
333 333
334 334 for n in range(self.nrows):
335 335 if self.colormaps is not None:
336 336 cmap = plt.get_cmap(self.colormaps[n])
337 337 else:
338 338 cmap = plt.get_cmap(self.colormap)
339 339 cmap.set_bad(self.bgcolor, 1.)
340 340 self.cmaps.append(cmap)
341 341
342 342 def __add_axes(self, ax, size='30%', pad='8%'):
343 343 '''
344 344 Add new axes to the given figure
345 345 '''
346 346 divider = make_axes_locatable(ax)
347 347 nax = divider.new_horizontal(size=size, pad=pad)
348 348 ax.figure.add_axes(nax)
349 349 return nax
350 350
351 351 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
352 352 '''
353 353 Create a masked array for missing data
354 354 '''
355 355 if x_buffer.shape[0] < 2:
356 356 return x_buffer, y_buffer, z_buffer
357 357
358 358 deltas = x_buffer[1:] - x_buffer[0:-1]
359 359 x_median = numpy.median(deltas)
360 360
361 361 index = numpy.where(deltas > 5 * x_median)
362 362
363 363 if len(index[0]) != 0:
364 364 z_buffer[::, index[0], ::] = self.__missing
365 365 z_buffer = numpy.ma.masked_inside(z_buffer,
366 366 0.99 * self.__missing,
367 367 1.01 * self.__missing)
368 368
369 369 return x_buffer, y_buffer, z_buffer
370 370
371 371 def decimate(self):
372 372
373 373 # dx = int(len(self.x)/self.__MAXNUMX) + 1
374 374 dy = int(len(self.y) / self.decimation) + 1
375 375
376 376 # x = self.x[::dx]
377 377 x = self.x
378 378 y = self.y[::dy]
379 379 z = self.z[::, ::, ::dy]
380 380
381 381 return x, y, z
382 382
383 383 def format(self):
384 384 '''
385 385 Set min and max values, labels, ticks and titles
386 386 '''
387 387
388 388 for n, ax in enumerate(self.axes):
389 389 if ax.firsttime:
390 390 if self.xaxis != 'time':
391 391 xmin = self.xmin
392 392 xmax = self.xmax
393 393 else:
394 394 xmin = self.tmin
395 395 xmax = self.tmin + self.xrange*60*60
396 396 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
397 397 if self.t_units == "h_m":
398 398 ax.xaxis.set_major_locator(LinearLocator(9))
399 399 if self.t_units == "h":
400 400 ax.xaxis.set_major_locator(LinearLocator(int((xmax-xmin)/3600)+1))
401 401 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
402 402 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
403 403 ax.set_facecolor(self.bgcolor)
404 404 if self.xscale:
405 405 ax.xaxis.set_major_formatter(FuncFormatter(
406 406 lambda x, pos: '{0:g}'.format(x*self.xscale)))
407 407 if self.yscale:
408 408 ax.yaxis.set_major_formatter(FuncFormatter(
409 409 lambda x, pos: '{0:g}'.format(x*self.yscale)))
410 410 if self.xlabel is not None:
411 411 ax.set_xlabel(self.xlabel)
412 412 if self.ylabel is not None:
413 413 ax.set_ylabel(self.ylabel)
414 414 if self.showprofile:
415 415 self.pf_axes[n].set_ylim(ymin, ymax)
416 416 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
417 417 self.pf_axes[n].set_xlabel('dB')
418 418 self.pf_axes[n].grid(b=True, axis='x')
419 419 [tick.set_visible(False)
420 420 for tick in self.pf_axes[n].get_yticklabels()]
421 421 if self.colorbar:
422 422 ax.cbar = plt.colorbar(
423 423 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
424 424 ax.cbar.ax.tick_params(labelsize=8)
425 425 ax.cbar.ax.press = None
426 426 if self.cb_label:
427 427 ax.cbar.set_label(self.cb_label, size=8)
428 428 elif self.cb_labels:
429 429 ax.cbar.set_label(self.cb_labels[n], size=8)
430 430 else:
431 431 ax.cbar = None
432 432 ax.set_xlim(xmin, xmax)
433 433 ax.set_ylim(ymin, ymax)
434 434 ax.firsttime = False
435 435 if self.grid:
436 436 ax.grid(True)
437 437
438 438 if not self.polar:
439 439 ax.set_title('{} {} {}'.format(
440 440 self.titles[n],
441 441 self.getDateTime(self.data.max_time).strftime(
442 442 '%Y-%m-%d %H:%M:%S'),
443 443 self.time_label),
444 444 size=8)
445 445 else:
446 446
447 447 ax.set_title('{}'.format(self.titles[n]), size=8)
448 448 ax.set_ylim(0, 90)
449 449 ax.set_yticks(numpy.arange(0, 90, 20))
450 450 ax.yaxis.labelpad = 40
451 451
452 452 if self.firsttime:
453 453 for n, fig in enumerate(self.figures):
454 454 fig.subplots_adjust(**self.plots_adjust)
455 455 self.firsttime = False
456 456
457 457 def clear_figures(self):
458 458 '''
459 459 Reset axes for redraw plots
460 460 '''
461 461
462 462 for ax in self.axes+self.pf_axes+self.cb_axes:
463 463 ax.clear()
464 464 ax.firsttime = True
465 465 if hasattr(ax, 'cbar') and ax.cbar:
466 466 ax.cbar.remove()
467 467
468 468 def __plot(self):
469 469 '''
470 470 Main function to plot, format and save figures
471 471 '''
472 472
473 473 self.plot()
474 474 self.format()
475 475
476 476 for n, fig in enumerate(self.figures):
477 477 if self.nrows == 0 or self.nplots == 0:
478 478 log.warning('No data', self.name)
479 479 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
480 480 fig.canvas.manager.set_window_title(self.CODE)
481 481 continue
482 482
483 483 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
484 484 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
485 485
486 486 fig.canvas.draw()
487 487 if self.show:
488 488 fig.show()
489 489 figpause(0.01)
490 490
491 491 if self.save:
492 492 self.save_figure(n)
493 493
494 494 if self.server:
495 495 self.send_to_server()
496 496
497 497 def __update(self, dataOut, timestamp):
498 498 '''
499 499 '''
500 500
501 501 metadata = {
502 502 'yrange': dataOut.heightList,
503 503 'interval': dataOut.timeInterval,
504 504 'channels': dataOut.channelList
505 505 }
506 506 data, meta = self.update(dataOut)
507 507 metadata.update(meta)
508 508 self.data.update(data, timestamp, metadata)
509 509
510 510 def save_figure(self, n):
511 511 '''
512 512 '''
513 513
514 514 if (self.data.max_time - self.save_time) <= self.save_period:
515 515 return
516 516
517 517 self.save_time = self.data.max_time
518 518
519 519 fig = self.figures[n]
520 520
521 521 if self.throttle == 0:
522 522 figname = os.path.join(
523 523 self.save,
524 524 self.save_code,
525 525 '{}_{}.png'.format(
526 526 self.save_code,
527 527 self.getDateTime(self.data.max_time).strftime(
528 528 '%Y%m%d_%H%M%S'
529 529 ),
530 530 )
531 531 )
532 532 log.log('Saving figure: {}'.format(figname), self.name)
533 533 if not os.path.isdir(os.path.dirname(figname)):
534 534 os.makedirs(os.path.dirname(figname))
535 535 fig.savefig(figname)
536 536
537 537 figname = os.path.join(
538 538 self.save,
539 539 '{}_{}.png'.format(
540 540 self.save_code,
541 541 self.getDateTime(self.data.min_time).strftime(
542 542 '%Y%m%d'
543 543 ),
544 544 )
545 545 )
546 546
547 547 log.log('Saving figure: {}'.format(figname), self.name)
548 548 if not os.path.isdir(os.path.dirname(figname)):
549 549 os.makedirs(os.path.dirname(figname))
550 550 fig.savefig(figname)
551 551
552 552 def send_to_server(self):
553 553 '''
554 554 '''
555 555
556 556 if self.exp_code == None:
557 557 log.warning('Missing `exp_code` skipping sending to server...')
558 558
559 559 last_time = self.data.max_time
560 560 interval = last_time - self.sender_time
561 561 if interval < self.sender_period:
562 562 return
563 563
564 564 self.sender_time = last_time
565 565
566 566 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
567 567 for attr in attrs:
568 568 value = getattr(self, attr)
569 569 if value:
570 570 if isinstance(value, (numpy.float32, numpy.float64)):
571 571 value = round(float(value), 2)
572 572 self.data.meta[attr] = value
573 573 if self.colormap == 'jet':
574 574 self.data.meta['colormap'] = 'Jet'
575 575 elif 'RdBu' in self.colormap:
576 576 self.data.meta['colormap'] = 'RdBu'
577 577 else:
578 578 self.data.meta['colormap'] = 'Viridis'
579 579 self.data.meta['interval'] = int(interval)
580 580
581 581 self.sender_queue.append(last_time)
582 582
583 583 while 1:
584 584 try:
585 585 tm = self.sender_queue.popleft()
586 586 except IndexError:
587 587 break
588 588 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
589 589 self.socket.send_string(msg)
590 590 socks = dict(self.poll.poll(2000))
591 591 if socks.get(self.socket) == zmq.POLLIN:
592 592 reply = self.socket.recv_string()
593 593 if reply == 'ok':
594 594 log.log("Response from server ok", self.name)
595 595 time.sleep(0.1)
596 596 continue
597 597 else:
598 598 log.warning(
599 599 "Malformed reply from server: {}".format(reply), self.name)
600 600 else:
601 601 log.warning(
602 602 "No response from server, retrying...", self.name)
603 603 self.sender_queue.appendleft(tm)
604 604 self.socket.setsockopt(zmq.LINGER, 0)
605 605 self.socket.close()
606 606 self.poll.unregister(self.socket)
607 607 self.socket = self.context.socket(zmq.REQ)
608 608 self.socket.connect(self.server)
609 609 self.poll.register(self.socket, zmq.POLLIN)
610 610 break
611 611
612 612 def setup(self):
613 613 '''
614 614 This method should be implemented in the child class, the following
615 615 attributes should be set:
616 616
617 617 self.nrows: number of rows
618 618 self.ncols: number of cols
619 619 self.nplots: number of plots (channels or pairs)
620 620 self.ylabel: label for Y axes
621 621 self.titles: list of axes title
622 622
623 623 '''
624 624 raise NotImplementedError
625 625
626 626 def plot(self):
627 627 '''
628 628 Must be defined in the child class, the actual plotting method
629 629 '''
630 630 raise NotImplementedError
631 631
632 632 def update(self, dataOut):
633 633 '''
634 634 Must be defined in the child class, update self.data with new data
635 635 '''
636 636
637 637 data = {
638 638 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
639 639 }
640 640 meta = {}
641 641
642 642 return data, meta
643 643
644 644 def run(self, dataOut, **kwargs):
645 645 '''
646 646 Main plotting routine
647 647 '''
648 648 if self.isConfig is False:
649 649 self.__setup(**kwargs)
650 650
651 651 if self.localtime:
652 652 self.getDateTime = datetime.datetime.fromtimestamp
653 653 else:
654 654 self.getDateTime = datetime.datetime.utcfromtimestamp
655 655
656 656 self.data.setup()
657 657 self.isConfig = True
658 658 if self.server:
659 659 self.context = zmq.Context()
660 660 self.socket = self.context.socket(zmq.REQ)
661 661 self.socket.connect(self.server)
662 662 self.poll = zmq.Poller()
663 663 self.poll.register(self.socket, zmq.POLLIN)
664 664
665 665 tm = getattr(dataOut, self.attr_time)
666 666
667 667 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
668 self.clear_figures()
668 669 self.save_time = tm
670 self.__plot()
669 671 self.tmin += self.xrange*60*60
670 672 self.data.setup()
671 self.clear_figures()
672 self.__plot()
673 #self.clear_figures()
674
673 675
674 676 self.__update(dataOut, tm)
675 677
676 678 if self.isPlotConfig is False:
677 679 self.__setup_plot()
678 680 self.isPlotConfig = True
679 681 if self.xaxis == 'time':
680 682 dt = self.getDateTime(tm)
681 683 if self.xmin is None:
682 684 self.tmin = tm
683 685 self.xmin = dt.hour
684 686 minutes = (self.xmin-int(self.xmin)) * 60
685 687 seconds = (minutes - int(minutes)) * 60
686 688 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
687 689 datetime.datetime(1970, 1, 1)).total_seconds()
688 690 if self.localtime:
689 691 self.tmin += time.timezone
690 692
691 693 if self.xmin is not None and self.xmax is not None:
692 694 self.xrange = self.xmax - self.xmin
693 695
694 696 if self.throttle == 0:
695 697 self.__plot()
696 698 else:
697 699 self.__throttle_plot(self.__plot)#, coerce=coerce)
698 700
699 701 def close(self):
700 702
701 703 if self.data and not self.data.flagNoData:
702 704 self.save_time = 0
703 705 self.__plot()
704 706 if self.data and not self.data.flagNoData and self.pause:
705 707 figpause(10)
Show file before | Show file after | Hide comments
@@ -1,3787 +1,3813
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from schainpy.model.io.utilsIO import getHei_index
8 8 from time import time
9 9 import datetime
10 10 import numpy
11 11 #import copy
12 12 from schainpy.model.data import _noise
13 13
14 14 from matplotlib import pyplot as plt
15 15
16 16 class VoltageProc(ProcessingUnit):
17 17
18 18 def __init__(self):
19 19
20 20 ProcessingUnit.__init__(self)
21 21
22 22 self.dataOut = Voltage()
23 23 self.flip = 1
24 24 self.setupReq = False
25 25
26 26 def run(self):
27 27 #print("running volt proc")
28 28
29 29 if self.dataIn.type == 'AMISR':
30 30 self.__updateObjFromAmisrInput()
31 31
32 32 if self.dataOut.buffer_empty:
33 33 if self.dataIn.type == 'Voltage':
34 34 self.dataOut.copy(self.dataIn)
35 35 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
36 36 self.dataOut.ippSeconds = self.dataIn.ippSeconds
37 37 self.dataOut.ipp = self.dataIn.ipp
38 38
39 39 #update Processing Header:
40 40 self.dataOut.processingHeaderObj.heightList = self.dataOut.heightList
41 41 self.dataOut.processingHeaderObj.ipp = self.dataOut.ipp
42 42 self.dataOut.processingHeaderObj.nCohInt = self.dataOut.nCohInt
43 43 self.dataOut.processingHeaderObj.dtype = self.dataOut.type
44 44 self.dataOut.processingHeaderObj.channelList = self.dataOut.channelList
45 45 self.dataOut.processingHeaderObj.azimuthList = self.dataOut.azimuthList
46 46 self.dataOut.processingHeaderObj.elevationList = self.dataOut.elevationList
47 47 self.dataOut.processingHeaderObj.codeList = self.dataOut.nChannels
48 48 self.dataOut.processingHeaderObj.heightList = self.dataOut.heightList
49 49 self.dataOut.processingHeaderObj.heightResolution = self.dataOut.heightList[1] - self.dataOut.heightList[0]
50 50
51 51
52 52
53 53 def __updateObjFromAmisrInput(self):
54 54
55 55 self.dataOut.timeZone = self.dataIn.timeZone
56 56 self.dataOut.dstFlag = self.dataIn.dstFlag
57 57 self.dataOut.errorCount = self.dataIn.errorCount
58 58 self.dataOut.useLocalTime = self.dataIn.useLocalTime
59 59
60 60 self.dataOut.flagNoData = self.dataIn.flagNoData
61 61 self.dataOut.data = self.dataIn.data
62 62 self.dataOut.utctime = self.dataIn.utctime
63 63 self.dataOut.channelList = self.dataIn.channelList
64 64 #self.dataOut.timeInterval = self.dataIn.timeInterval
65 65 self.dataOut.heightList = self.dataIn.heightList
66 66 self.dataOut.nProfiles = self.dataIn.nProfiles
67 67
68 68 self.dataOut.nCohInt = self.dataIn.nCohInt
69 69 self.dataOut.ippSeconds = self.dataIn.ippSeconds
70 70 self.dataOut.frequency = self.dataIn.frequency
71 71
72 72 self.dataOut.azimuth = self.dataIn.azimuth
73 73 self.dataOut.zenith = self.dataIn.zenith
74 74
75 75 self.dataOut.beam.codeList = self.dataIn.beam.codeList
76 76 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
77 77 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
78 78
79 79
80 80 class selectChannels(Operation):
81 81
82 82 def run(self, dataOut, channelList=None):
83 83 self.channelList = channelList
84 84 if self.channelList == None:
85 85 print("Missing channelList")
86 86 return dataOut
87 87 channelIndexList = []
88 88 if not dataOut.buffer_empty: # cuando se usa proc volts como buffer de datos
89 89 return dataOut
90 90 #print("channel List: ", dataOut.channelList)
91 91 if type(dataOut.channelList) is not list: #leer array desde HDF5
92 92 try:
93 93 dataOut.channelList = dataOut.channelList.tolist()
94 94 except Exception as e:
95 95 print("Select Channels: ",e)
96 96 for channel in self.channelList:
97 97 if channel not in dataOut.channelList:
98 98 raise ValueError("Channel %d is not in %s" %(channel, str(dataOut.channelList)))
99 99
100 100 index = dataOut.channelList.index(channel)
101 101 channelIndexList.append(index)
102 102 dataOut = self.selectChannelsByIndex(dataOut,channelIndexList)
103 103
104 104 #update Processing Header:
105 105 dataOut.processingHeaderObj.channelList = dataOut.channelList
106 106 dataOut.processingHeaderObj.elevationList = dataOut.elevationList
107 107 dataOut.processingHeaderObj.azimuthList = dataOut.azimuthList
108 108 dataOut.processingHeaderObj.codeList = dataOut.codeList
109 109 dataOut.processingHeaderObj.nChannels = len(dataOut.channelList)
110 110
111 111 return dataOut
112 112
113 113 def selectChannelsByIndex(self, dataOut, channelIndexList):
114 114 """
115 115 Selecciona un bloque de datos en base a canales segun el channelIndexList
116 116
117 117 Input:
118 118 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
119 119
120 120 Affected:
121 121 dataOut.data
122 122 dataOut.channelIndexList
123 123 dataOut.nChannels
124 124 dataOut.m_ProcessingHeader.totalSpectra
125 125 dataOut.systemHeaderObj.numChannels
126 126 dataOut.m_ProcessingHeader.blockSize
127 127
128 128 Return:
129 129 None
130 130 """
131 131 #print("selectChannelsByIndex")
132 132 # for channelIndex in channelIndexList:
133 133 # if channelIndex not in dataOut.channelIndexList:
134 134 # raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
135 135
136 136 if dataOut.type == 'Voltage':
137 137 if dataOut.flagDataAsBlock:
138 138 """
139 139 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
140 140 """
141 141 data = dataOut.data[channelIndexList,:,:]
142 142 else:
143 143 data = dataOut.data[channelIndexList,:]
144 144
145 145 dataOut.data = data
146 146 # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
147 147 dataOut.channelList = [n for n in range(len(channelIndexList))]
148 148
149 149 elif dataOut.type == 'Spectra':
150 150 if hasattr(dataOut, 'data_spc'):
151 151 if dataOut.data_spc is None:
152 152 raise ValueError("data_spc is None")
153 153 return dataOut
154 154 else:
155 155 data_spc = dataOut.data_spc[channelIndexList, :]
156 156 dataOut.data_spc = data_spc
157 157
158 158 # if hasattr(dataOut, 'data_dc') :# and
159 159 # if dataOut.data_dc is None:
160 160 # raise ValueError("data_dc is None")
161 161 # return dataOut
162 162 # else:
163 163 # data_dc = dataOut.data_dc[channelIndexList, :]
164 164 # dataOut.data_dc = data_dc
165 165 # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
166 166 dataOut.channelList = channelIndexList
167 167 dataOut = self.__selectPairsByChannel(dataOut,channelIndexList)
168 168 if len(dataOut.elevationList>0):
169 169 dataOut.elevationList = dataOut.elevationList[channelIndexList]
170 170 dataOut.azimuthList = dataOut.azimuthList[channelIndexList]
171 171 dataOut.codeList = dataOut.codeList[channelIndexList]
172 172 return dataOut
173 173
174 174 def __selectPairsByChannel(self, dataOut, channelList=None):
175 175 #print("__selectPairsByChannel")
176 176 if channelList == None:
177 177 return
178 178
179 179 pairsIndexListSelected = []
180 180 for pairIndex in dataOut.pairsIndexList:
181 181 # First pair
182 182 if dataOut.pairsList[pairIndex][0] not in channelList:
183 183 continue
184 184 # Second pair
185 185 if dataOut.pairsList[pairIndex][1] not in channelList:
186 186 continue
187 187
188 188 pairsIndexListSelected.append(pairIndex)
189 189 if not pairsIndexListSelected:
190 190 dataOut.data_cspc = None
191 191 dataOut.pairsList = []
192 192 return
193 193
194 194 dataOut.data_cspc = dataOut.data_cspc[pairsIndexListSelected]
195 195 dataOut.pairsList = [dataOut.pairsList[i]
196 196 for i in pairsIndexListSelected]
197 197
198 198 return dataOut
199 199
200 200 class selectHeights(Operation):
201 201
202 202 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
203 203 """
204 204 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
205 205 minHei <= height <= maxHei
206 206
207 207 Input:
208 208 minHei : valor minimo de altura a considerar
209 209 maxHei : valor maximo de altura a considerar
210 210
211 211 Affected:
212 212 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
213 213
214 214 Return:
215 215 1 si el metodo se ejecuto con exito caso contrario devuelve 0
216 216 """
217 217
218 218 self.dataOut = dataOut
219 219
220 220 if minHei and maxHei:
221 221
222 222 if (minHei < dataOut.heightList[0]):
223 223 minHei = dataOut.heightList[0]
224 224
225 225 if (maxHei > dataOut.heightList[-1]):
226 226 maxHei = dataOut.heightList[-1]
227 227
228 228 minIndex = 0
229 229 maxIndex = 0
230 230 heights = dataOut.heightList
231 231
232 232 inda = numpy.where(heights >= minHei)
233 233 indb = numpy.where(heights <= maxHei)
234 234
235 235 try:
236 236 minIndex = inda[0][0]
237 237 except:
238 238 minIndex = 0
239 239
240 240 try:
241 241 maxIndex = indb[0][-1]
242 242 except:
243 243 maxIndex = len(heights)
244 244
245 245 self.selectHeightsByIndex(minIndex, maxIndex)
246 246
247 247 #update Processing Header:
248 248 dataOut.processingHeaderObj.heightList = dataOut.heightList
249 249
250 250
251 251
252 252 return dataOut
253 253
254 254 def selectHeightsByIndex(self, minIndex, maxIndex):
255 255 """
256 256 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
257 257 minIndex <= index <= maxIndex
258 258
259 259 Input:
260 260 minIndex : valor de indice minimo de altura a considerar
261 261 maxIndex : valor de indice maximo de altura a considerar
262 262
263 263 Affected:
264 264 self.dataOut.data
265 265 self.dataOut.heightList
266 266
267 267 Return:
268 268 1 si el metodo se ejecuto con exito caso contrario devuelve 0
269 269 """
270 270
271 271 if self.dataOut.type == 'Voltage':
272 272 if (minIndex < 0) or (minIndex > maxIndex):
273 273 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
274 274
275 275 if (maxIndex >= self.dataOut.nHeights):
276 276 maxIndex = self.dataOut.nHeights
277 277
278 278 #voltage
279 279 if self.dataOut.flagDataAsBlock:
280 280 """
281 281 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
282 282 """
283 283 data = self.dataOut.data[:,:, minIndex:maxIndex]
284 284 else:
285 285 data = self.dataOut.data[:, minIndex:maxIndex]
286 286
287 287 # firstHeight = self.dataOut.heightList[minIndex]
288 288
289 289 self.dataOut.data = data
290 290 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
291 291
292 292 if self.dataOut.nHeights <= 1:
293 293 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
294 294 elif self.dataOut.type == 'Spectra':
295 295 if (minIndex < 0) or (minIndex > maxIndex):
296 296 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
297 297 minIndex, maxIndex))
298 298
299 299 if (maxIndex >= self.dataOut.nHeights):
300 300 maxIndex = self.dataOut.nHeights - 1
301 301
302 302 # Spectra
303 303 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
304 304
305 305 data_cspc = None
306 306 if self.dataOut.data_cspc is not None:
307 307 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
308 308
309 309 data_dc = None
310 310 if self.dataOut.data_dc is not None:
311 311 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
312 312
313 313 self.dataOut.data_spc = data_spc
314 314 self.dataOut.data_cspc = data_cspc
315 315 self.dataOut.data_dc = data_dc
316 316
317 317 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
318 318
319 319 return 1
320 320
321 321
322 322 class filterByHeights(Operation):
323
323 ifConfig=False
324 deltaHeight = None
325 newdelta=None
326 newheights=None
327 r=None
328 h0=None
329 nHeights=None
324 330 def run(self, dataOut, window):
325
326 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
327
331
332
333 # print("1",dataOut.data.shape)
334 # print(dataOut.nHeights)
328 335 if window == None:
329 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
330
331 newdelta = deltaHeight * window
332 r = dataOut.nHeights % window
333 newheights = (dataOut.nHeights-r)/window
334
335 if newheights <= 1:
336 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
336 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / self.deltaHeight
337
338 if not self.ifConfig: #and dataOut.useInputBuffer:
339 self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
340 self.ifConfig = True
341 self.newdelta = self.deltaHeight * window
342 self.r = dataOut.nHeights % window
343 self.newheights = (dataOut.nHeights-self.r)/window
344 self.h0 = dataOut.heightList[0]
345 self.nHeights = dataOut.nHeights
346 if self.newheights <= 1:
347 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
337 348
338 349 if dataOut.flagDataAsBlock:
339 350 """
340 351 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
341 352 """
342 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
343 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
353 buffer = dataOut.data[:, :, 0:int(self.nHeights-self.r)]
354 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(self.nHeights/window), window)
344 355 buffer = numpy.sum(buffer,3)
345 356
346 357 else:
347 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
348 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
358 buffer = dataOut.data[:,0:int(self.nHeights-self.r)]
359 buffer = buffer.reshape(dataOut.nChannels,int(self.nHeights/window),int(window))
349 360 buffer = numpy.sum(buffer,2)
350 361
351 362 dataOut.data = buffer
352 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
363 dataOut.heightList = self.h0 + numpy.arange( self.newheights )*self.newdelta
353 364 dataOut.windowOfFilter = window
354 365
355 366 #update Processing Header:
356 367 dataOut.processingHeaderObj.heightList = dataOut.heightList
357 368 dataOut.processingHeaderObj.nWindows = window
358
369
359 370 return dataOut
360 371
361 372
373
362 374 class setH0(Operation):
363 375
364 376 def run(self, dataOut, h0, deltaHeight = None):
365 377
366 378 if not deltaHeight:
367 379 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
368 380
369 381 nHeights = dataOut.nHeights
370 382
371 383 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
372 384
373 385 dataOut.heightList = newHeiRange
374 386
375 387 #update Processing Header:
376 388 dataOut.processingHeaderObj.heightList = dataOut.heightList
377 389
378 390 return dataOut
379 391
380 392
381 393 class deFlip(Operation):
382 394
383 395 def run(self, dataOut, channelList = []):
384 396
385 397 data = dataOut.data.copy()
386 398
387 399 if dataOut.flagDataAsBlock:
388 400 flip = self.flip
389 401 profileList = list(range(dataOut.nProfiles))
390 402
391 403 if not channelList:
392 404 for thisProfile in profileList:
393 405 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
394 406 flip *= -1.0
395 407 else:
396 408 for thisChannel in channelList:
397 409 if thisChannel not in dataOut.channelList:
398 410 continue
399 411
400 412 for thisProfile in profileList:
401 413 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
402 414 flip *= -1.0
403 415
404 416 self.flip = flip
405 417
406 418 else:
407 419 if not channelList:
408 420 data[:,:] = data[:,:]*self.flip
409 421 else:
410 422 for thisChannel in channelList:
411 423 if thisChannel not in dataOut.channelList:
412 424 continue
413 425
414 426 data[thisChannel,:] = data[thisChannel,:]*self.flip
415 427
416 428 self.flip *= -1.
417 429
418 430 dataOut.data = data
419 431
420 432 return dataOut
421 433
422 434
423 435 class setAttribute(Operation):
424 436 '''
425 437 Set an arbitrary attribute(s) to dataOut
426 438 '''
427 439
428 440 def __init__(self):
429 441
430 442 Operation.__init__(self)
431 443 self._ready = False
432 444
433 445 def run(self, dataOut, **kwargs):
434 446
435 447 for key, value in kwargs.items():
436 448 setattr(dataOut, key, value)
437 449
438 450 return dataOut
439 451
440 452
441 453 @MPDecorator
442 454 class printAttribute(Operation):
443 455 '''
444 456 Print an arbitrary attribute of dataOut
445 457 '''
446 458
447 459 def __init__(self):
448 460
449 461 Operation.__init__(self)
450 462
451 463 def run(self, dataOut, attributes):
452 464
453 465 if isinstance(attributes, str):
454 466 attributes = [attributes]
455 467 for attr in attributes:
456 468 if hasattr(dataOut, attr):
457 469 log.log(getattr(dataOut, attr), attr)
458 470
459 471 class cleanHeightsInterf(Operation):
460 472 __slots__ =('heights_indx', 'repeats', 'step', 'factor', 'idate', 'idxs','config','wMask')
461 473 def __init__(self):
462 474 self.repeats = 0
463 475 self.factor=1
464 476 self.wMask = None
465 477 self.config = False
466 478 self.idxs = None
467 479 self.heights_indx = None
468 480
469 481 def run(self, dataOut, heightsList, repeats=0, step=0, factor=1, idate=None, startH=None, endH=None):
470 482
471 483 #print(dataOut.data.shape)
472 484
473 485 startTime = datetime.datetime.combine(idate,startH)
474 486 endTime = datetime.datetime.combine(idate,endH)
475 487 currentTime = datetime.datetime.fromtimestamp(dataOut.utctime)
476 488
477 489 if currentTime < startTime or currentTime > endTime:
478 490 return dataOut
479 491 if not self.config:
480 492
481 493 #print(wMask)
482 494 heights = [float(hei) for hei in heightsList]
483 495 for r in range(repeats):
484 496 heights += [ (h+(step*(r+1))) for h in heights]
485 497 #print(heights)
486 498 heiList = dataOut.heightList
487 499 self.heights_indx = [getHei_index(h,h,heiList)[0] for h in heights]
488 500
489 501 self.wMask = numpy.asarray(factor)
490 502 self.wMask = numpy.tile(self.wMask,(repeats+2))
491 503 self.config = True
492 504
493 505 """
494 506 getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None)
495 507 """
496 508 #print(self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref)))
497 509
498 510
499 511 for ch in range(dataOut.data.shape[0]):
500 512 i = 0
501 513
502 514
503 515 for hei in self.heights_indx:
504 516 h = hei - 1
505 517
506 518
507 519 if dataOut.data.ndim < 3:
508 520 module = numpy.absolute(dataOut.data[ch,h])
509 521 prev_h1 = numpy.absolute(dataOut.data[ch,h-1])
510 522 dataOut.data[ch,h] = (dataOut.data[ch,h])/module * prev_h1
511 523
512 524 #dataOut.data[ch,hei-1] = (dataOut.data[ch,hei-1])*self.wMask[i]
513 525 else:
514 526 module = numpy.absolute(dataOut.data[ch,:,h])
515 527 prev_h1 = numpy.absolute(dataOut.data[ch,:,h-1])
516 528 dataOut.data[ch,:,h] = (dataOut.data[ch,:,h])/module * prev_h1
517 529 #dataOut.data[ch,:,hei-1] = (dataOut.data[ch,:,hei-1])*self.wMask[i]
518 530 #print("done")
519 531 i += 1
520 532
521 533
522 534 return dataOut
523 535
524 536
525 537
526 538 class interpolateHeights(Operation):
527 539
528 540 def run(self, dataOut, topLim, botLim):
529 541 #69 al 72 para julia
530 542 #82-84 para meteoros
531 543 if len(numpy.shape(dataOut.data))==2:
532 544 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
533 545 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
534 546 #dataOut.data[:,botLim:limSup+1] = sampInterp
535 547 dataOut.data[:,botLim:topLim+1] = sampInterp
536 548 else:
537 549 nHeights = dataOut.data.shape[2]
538 550 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
539 551 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
540 552 f = interpolate.interp1d(x, y, axis = 2)
541 553 xnew = numpy.arange(botLim,topLim+1)
542 554 ynew = f(xnew)
543 555 dataOut.data[:,:,botLim:topLim+1] = ynew
544 556
545 557 return dataOut
546 558
547 559
548 560 class CohInt(Operation):
549 561
550 562 isConfig = False
551 563 __profIndex = 0
552 564 __byTime = False
553 565 __initime = None
554 566 __lastdatatime = None
555 567 __integrationtime = None
556 568 __buffer = None
557 569 __bufferStride = []
558 570 __dataReady = False
559 571 __profIndexStride = 0
560 572 __dataToPutStride = False
561 573 n = None
562 574
563 575 def __init__(self, **kwargs):
564 576
565 577 Operation.__init__(self, **kwargs)
566 578
567 579 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
568 580 """
569 581 Set the parameters of the integration class.
570 582
571 583 Inputs:
572 584
573 585 n : Number of coherent integrations
574 586 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
575 587 overlapping :
576 588 """
577 589
578 590 self.__initime = None
579 591 self.__lastdatatime = 0
580 592 self.__buffer = None
581 593 self.__dataReady = False
582 594 self.byblock = byblock
583 595 self.stride = stride
584 596
585 597 if n == None and timeInterval == None:
586 598 raise ValueError("n or timeInterval should be specified ...")
587 599
588 600 if n != None:
589 601 self.n = n
590 602 self.__byTime = False
591 603 else:
592 604 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
593 605 self.n = 9999
594 606 self.__byTime = True
595 607
596 608 if overlapping:
597 609 self.__withOverlapping = True
598 610 self.__buffer = None
599 611 else:
600 612 self.__withOverlapping = False
601 613 self.__buffer = 0
602 614
603 615 self.__profIndex = 0
604 616
605 617 def putData(self, data):
606 618
607 619 """
608 620 Add a profile to the __buffer and increase in one the __profileIndex
609 621
610 622 """
611 623
612 624 if not self.__withOverlapping:
613 625 self.__buffer += data.copy()
614 626 self.__profIndex += 1
615 627 return
616 628
617 629 #Overlapping data
618 630 nChannels, nHeis = data.shape
619 631 data = numpy.reshape(data, (1, nChannels, nHeis))
620 632
621 633 #If the buffer is empty then it takes the data value
622 634 if self.__buffer is None:
623 635 self.__buffer = data
624 636 self.__profIndex += 1
625 637 return
626 638
627 639 #If the buffer length is lower than n then stakcing the data value
628 640 if self.__profIndex < self.n:
629 641 self.__buffer = numpy.vstack((self.__buffer, data))
630 642 self.__profIndex += 1
631 643 return
632 644
633 645 #If the buffer length is equal to n then replacing the last buffer value with the data value
634 646 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
635 647 self.__buffer[self.n-1] = data
636 648 self.__profIndex = self.n
637 649 return
638 650
639 651
640 652 def pushData(self):
641 653 """
642 654 Return the sum of the last profiles and the profiles used in the sum.
643 655
644 656 Affected:
645 657
646 658 self.__profileIndex
647 659
648 660 """
649 661
650 662 if not self.__withOverlapping:
651 663 data = self.__buffer
652 664 n = self.__profIndex
653 665
654 666 self.__buffer = 0
655 667 self.__profIndex = 0
656 668
657 669 return data, n
658 670
659 671 #Integration with Overlapping
660 672 data = numpy.sum(self.__buffer, axis=0)
661 673 # print data
662 674 # raise
663 675 n = self.__profIndex
664 676
665 677 return data, n
666 678
667 679 def byProfiles(self, data):
668 680
669 681 self.__dataReady = False
670 682 avgdata = None
671 683 # n = None
672 684 # print data
673 685 # raise
674 686 self.putData(data)
675 687
676 688 if self.__profIndex == self.n:
677 689 avgdata, n = self.pushData()
678 690 self.__dataReady = True
679 691
680 692 return avgdata
681 693
682 694 def byTime(self, data, datatime):
683 695
684 696 self.__dataReady = False
685 697 avgdata = None
686 698 n = None
687 699
688 700 self.putData(data)
689 701
690 702 if (datatime - self.__initime) >= self.__integrationtime:
691 703 avgdata, n = self.pushData()
692 704 self.n = n
693 705 self.__dataReady = True
694 706
695 707 return avgdata
696 708
697 709 def integrateByStride(self, data, datatime):
698 710 # print data
699 711 if self.__profIndex == 0:
700 712 self.__buffer = [[data.copy(), datatime]]
701 713 else:
702 714 self.__buffer.append([data.copy(),datatime])
703 715 self.__profIndex += 1
704 716 self.__dataReady = False
705 717
706 718 if self.__profIndex == self.n * self.stride :
707 719 self.__dataToPutStride = True
708 720 self.__profIndexStride = 0
709 721 self.__profIndex = 0
710 722 self.__bufferStride = []
711 723 for i in range(self.stride):
712 724 current = self.__buffer[i::self.stride]
713 725 data = numpy.sum([t[0] for t in current], axis=0)
714 726 avgdatatime = numpy.average([t[1] for t in current])
715 727 # print data
716 728 self.__bufferStride.append((data, avgdatatime))
717 729
718 730 if self.__dataToPutStride:
719 731 self.__dataReady = True
720 732 self.__profIndexStride += 1
721 733 if self.__profIndexStride == self.stride:
722 734 self.__dataToPutStride = False
723 735 # print self.__bufferStride[self.__profIndexStride - 1]
724 736 # raise
725 737 return self.__bufferStride[self.__profIndexStride - 1]
726 738
727 739
728 740 return None, None
729 741
730 742 def integrate(self, data, datatime=None):
731 743
732 744 if self.__initime == None:
733 745 self.__initime = datatime
734 746
735 747 if self.__byTime:
736 748 avgdata = self.byTime(data, datatime)
737 749 else:
738 750 avgdata = self.byProfiles(data)
739 751
740 752
741 753 self.__lastdatatime = datatime
742 754
743 755 if avgdata is None:
744 756 return None, None
745 757
746 758 avgdatatime = self.__initime
747 759
748 760 deltatime = datatime - self.__lastdatatime
749 761
750 762 if not self.__withOverlapping:
751 763 self.__initime = datatime
752 764 else:
753 765 self.__initime += deltatime
754 766
755 767 return avgdata, avgdatatime
756 768
757 769 def integrateByBlock(self, dataOut):
758 770
759 771 times = int(dataOut.data.shape[1]/self.n)
760 772 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
761 773
762 774 id_min = 0
763 775 id_max = self.n
764 776
765 777 for i in range(times):
766 778 junk = dataOut.data[:,id_min:id_max,:]
767 779 avgdata[:,i,:] = junk.sum(axis=1)
768 780 id_min += self.n
769 781 id_max += self.n
770 782
771 783 timeInterval = dataOut.ippSeconds*self.n
772 784 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
773 785 self.__dataReady = True
774 786 return avgdata, avgdatatime
775 787
776 788 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
777 789
778 790 if not self.isConfig:
779 791 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
780 792 self.isConfig = True
781 793
782 794 if dataOut.flagDataAsBlock:
783 795 """
784 796 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
785 797 """
786 798 avgdata, avgdatatime = self.integrateByBlock(dataOut)
787 799 dataOut.nProfiles /= self.n
788 800 else:
789 801 if stride is None:
790 802 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
791 803 else:
792 804 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
793 805
794 806
795 807 # dataOut.timeInterval *= n
796 808 dataOut.flagNoData = True
797 809
798 810 if self.__dataReady:
799 811 dataOut.data = avgdata
800 812 if not dataOut.flagCohInt:
801 813 dataOut.nCohInt *= self.n
802 814 dataOut.flagCohInt = True
803 815 dataOut.utctime = avgdatatime
804 816 # print avgdata, avgdatatime
805 817 # raise
806 818 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
807 819 dataOut.flagNoData = False
808 820
809 821 #update Processing Header:
810 822 dataOut.processingHeaderObj.nCohInt = dataOut.nCohInt
811 823
812 824
813 825 return dataOut
814 826
815 827 class Decoder(Operation):
816 828
817 829 isConfig = False
818 830 __profIndex = 0
819 831
820 832 code = None
821 833
822 834 nCode = None
823 835 nBaud = None
824 836
825 837 def __init__(self, **kwargs):
826 838
827 839 Operation.__init__(self, **kwargs)
828 840
829 841 self.times = None
830 842 self.osamp = None
831 843 # self.__setValues = False
832 844 self.isConfig = False
833 845 self.setupReq = False
834 846 def setup(self, code, osamp, dataOut):
835 847
836 848 self.__profIndex = 0
837 849
838 850 self.code = code
839 851
840 852 self.nCode = len(code)
841 853 self.nBaud = len(code[0])
842 854 if (osamp != None) and (osamp >1):
843 855 self.osamp = osamp
844 856 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
845 857 self.nBaud = self.nBaud*self.osamp
846 858
847 859 self.__nChannels = dataOut.nChannels
848 860 self.__nProfiles = dataOut.nProfiles
849 861 self.__nHeis = dataOut.nHeights
850 862
851 863 if self.__nHeis < self.nBaud:
852 864 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
853 865
854 866 #Frequency
855 867 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
856 868
857 869 __codeBuffer[:,0:self.nBaud] = self.code
858 870
859 871 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
860 872
861 873 if dataOut.flagDataAsBlock:
862 874
863 875 self.ndatadec = self.__nHeis #- self.nBaud + 1
864 876
865 877 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
866 878
867 879 else:
868 880
869 881 #Time
870 882 self.ndatadec = self.__nHeis #- self.nBaud + 1
871 883
872 884 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
873 885
874 886 def __convolutionInFreq(self, data):
875 887
876 888 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
877 889
878 890 fft_data = numpy.fft.fft(data, axis=1)
879 891
880 892 conv = fft_data*fft_code
881 893
882 894 data = numpy.fft.ifft(conv,axis=1)
883 895
884 896 return data
885 897
886 898 def __convolutionInFreqOpt(self, data):
887 899
888 900 raise NotImplementedError
889 901
890 902 def __convolutionInTime(self, data):
891 903
892 904 code = self.code[self.__profIndex]
893 905 for i in range(self.__nChannels):
894 906 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
895 907
896 908 return self.datadecTime
897 909
898 910 def __convolutionByBlockInTime(self, data):
899 911
900 912 repetitions = int(self.__nProfiles / self.nCode)
901 913 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
902 914 junk = junk.flatten()
903 915 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
904 916 profilesList = range(self.__nProfiles)
905 917
906 918 for i in range(self.__nChannels):
907 919 for j in profilesList:
908 920 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
909 921 return self.datadecTime
910 922
911 923 def __convolutionByBlockInFreq(self, data):
912 924
913 925 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
914 926
915 927
916 928 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
917 929
918 930 fft_data = numpy.fft.fft(data, axis=2)
919 931
920 932 conv = fft_data*fft_code
921 933
922 934 data = numpy.fft.ifft(conv,axis=2)
923 935
924 936 return data
925 937
926 938
927 939 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
928 940
929 941 if dataOut.flagDecodeData:
930 942 print("This data is already decoded, recoding again ...")
931 943
932 944 if not self.isConfig:
933 945
934 946 if code is None:
935 947 if dataOut.code is None:
936 948 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
937 949
938 950 code = dataOut.code
939 951 else:
940 952 code = numpy.array(code).reshape(nCode,nBaud)
941 953 self.setup(code, osamp, dataOut)
942 954
943 955 self.isConfig = True
944 956
945 957 if mode == 3:
946 958 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
947 959
948 960 if times != None:
949 961 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
950 962
951 963 if self.code is None:
952 964 print("Fail decoding: Code is not defined.")
953 965 return
954 966
955 967 self.__nProfiles = dataOut.nProfiles
956 968 datadec = None
957 969
958 970 if mode == 3:
959 971 mode = 0
960 972
961 973 if dataOut.flagDataAsBlock:
962 974 """
963 975 Decoding when data have been read as block,
964 976 """
965 977
966 978 if mode == 0:
967 979 datadec = self.__convolutionByBlockInTime(dataOut.data)
968 980 if mode == 1:
969 981 datadec = self.__convolutionByBlockInFreq(dataOut.data)
970 982 else:
971 983 """
972 984 Decoding when data have been read profile by profile
973 985 """
974 986 if mode == 0:
975 987 datadec = self.__convolutionInTime(dataOut.data)
976 988
977 989 if mode == 1:
978 990 datadec = self.__convolutionInFreq(dataOut.data)
979 991
980 992 if mode == 2:
981 993 datadec = self.__convolutionInFreqOpt(dataOut.data)
982 994
983 995 if datadec is None:
984 996 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
985 997
986 998 dataOut.code = self.code
987 999 dataOut.nCode = self.nCode
988 1000 dataOut.nBaud = self.nBaud
989 1001
990 1002 dataOut.data = datadec
991 1003 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
992 1004 dataOut.flagDecodeData = True #asumo q la data esta decodificada
993 1005
994 1006
995 1007 #update Processing Header:
996 1008 dataOut.radarControllerHeaderObj.code = self.code
997 1009 dataOut.radarControllerHeaderObj.nCode = self.nCode
998 1010 dataOut.radarControllerHeaderObj.nBaud = self.nBaud
999 1011 dataOut.radarControllerHeaderObj.nOsamp = osamp
1000 1012 #update Processing Header:
1001 1013 dataOut.processingHeaderObj.heightList = dataOut.heightList
1002 1014 dataOut.processingHeaderObj.heightResolution = dataOut.heightList[1]-dataOut.heightList[0]
1003 1015
1004 1016 if self.__profIndex == self.nCode-1:
1005 1017 self.__profIndex = 0
1006 1018 return dataOut
1007 1019
1008 1020 self.__profIndex += 1
1009 1021
1010 1022 return dataOut
1011 1023
1012 1024 class ProfileConcat(Operation):
1013 1025
1014 1026 isConfig = False
1015 1027 buffer = None
1016 1028
1017 1029 def __init__(self, **kwargs):
1018 1030
1019 1031 Operation.__init__(self, **kwargs)
1020 1032 self.profileIndex = 0
1021 1033
1022 1034 def reset(self):
1023 1035 self.buffer = numpy.zeros_like(self.buffer)
1024 1036 self.start_index = 0
1025 1037 self.times = 1
1026 1038
1027 1039 def setup(self, data, m, n=1):
1028 1040 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
1029 1041 self.nHeights = data.shape[1]#.nHeights
1030 1042 self.start_index = 0
1031 1043 self.times = 1
1032 1044
1033 1045 def concat(self, data):
1034 1046
1035 1047 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
1036 1048 self.start_index = self.start_index + self.nHeights
1037 1049
1038 1050 def run(self, dataOut, m):
1039 1051 dataOut.flagNoData = True
1040 1052
1041 1053 if not self.isConfig:
1042 1054 self.setup(dataOut.data, m, 1)
1043 1055 self.isConfig = True
1044 1056
1045 1057 if dataOut.flagDataAsBlock:
1046 1058 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
1047 1059
1048 1060 else:
1049 1061 self.concat(dataOut.data)
1050 1062 self.times += 1
1051 1063 if self.times > m:
1052 1064 dataOut.data = self.buffer
1053 1065 self.reset()
1054 1066 dataOut.flagNoData = False
1055 1067 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
1056 1068 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1057 1069 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
1058 1070 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
1059 1071 dataOut.ippSeconds *= m
1060 1072
1061 1073 #update Processing Header:
1062 1074 dataOut.processingHeaderObj.heightList = dataOut.heightList
1063 1075 dataOut.processingHeaderObj.ipp = dataOut.ippSeconds
1064 1076
1065 1077 return dataOut
1066 1078
1067 1079 class ProfileSelector(Operation):
1068 1080
1069 1081 profileIndex = None
1070 1082 # Tamanho total de los perfiles
1071 1083 nProfiles = None
1072 1084
1073 1085 def __init__(self, **kwargs):
1074 1086
1075 1087 Operation.__init__(self, **kwargs)
1076 1088 self.profileIndex = 0
1077 1089
1078 1090 def incProfileIndex(self):
1079 1091
1080 1092 self.profileIndex += 1
1081 1093
1082 1094 if self.profileIndex >= self.nProfiles:
1083 1095 self.profileIndex = 0
1084 1096
1085 1097 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
1086 1098
1087 1099 if profileIndex < minIndex:
1088 1100 return False
1089 1101
1090 1102 if profileIndex > maxIndex:
1091 1103 return False
1092 1104
1093 1105 return True
1094 1106
1095 1107 def isThisProfileInList(self, profileIndex, profileList):
1096 1108
1097 1109 if profileIndex not in profileList:
1098 1110 return False
1099 1111
1100 1112 return True
1101 1113
1102 1114 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
1103 1115
1104 1116 """
1105 1117 ProfileSelector:
1106 1118
1107 1119 Inputs:
1108 1120 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
1109 1121
1110 1122 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
1111 1123
1112 1124 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
1113 1125
1114 1126 """
1115 1127
1116 1128 if rangeList is not None:
1117 1129 if type(rangeList[0]) not in (tuple, list):
1118 1130 rangeList = [rangeList]
1119 1131
1120 1132 dataOut.flagNoData = True
1121 1133
1122 1134 if dataOut.flagDataAsBlock:
1123 1135 """
1124 1136 data dimension = [nChannels, nProfiles, nHeis]
1125 1137 """
1126 1138 if profileList != None:
1127 1139 dataOut.data = dataOut.data[:,profileList,:]
1128 1140
1129 1141 if profileRangeList != None:
1130 1142 minIndex = profileRangeList[0]
1131 1143 maxIndex = profileRangeList[1]
1132 1144 profileList = list(range(minIndex, maxIndex+1))
1133 1145
1134 1146 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
1135 1147
1136 1148 if rangeList != None:
1137 1149
1138 1150 profileList = []
1139 1151
1140 1152 for thisRange in rangeList:
1141 1153 minIndex = thisRange[0]
1142 1154 maxIndex = thisRange[1]
1143 1155
1144 1156 profileList.extend(list(range(minIndex, maxIndex+1)))
1145 1157
1146 1158 dataOut.data = dataOut.data[:,profileList,:]
1147 1159
1148 1160 dataOut.nProfiles = len(profileList)
1149 1161 dataOut.profileIndex = dataOut.nProfiles - 1
1150 1162 dataOut.flagNoData = False
1151 1163
1152 1164 return dataOut
1153 1165
1154 1166 """
1155 1167 data dimension = [nChannels, nHeis]
1156 1168 """
1157 1169
1158 1170 if profileList != None:
1159 1171
1160 1172 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1161 1173
1162 1174 self.nProfiles = len(profileList)
1163 1175 dataOut.nProfiles = self.nProfiles
1164 1176 dataOut.profileIndex = self.profileIndex
1165 1177 dataOut.flagNoData = False
1166 1178
1167 1179 self.incProfileIndex()
1168 1180 return dataOut
1169 1181
1170 1182 if profileRangeList != None:
1171 1183
1172 1184 minIndex = profileRangeList[0]
1173 1185 maxIndex = profileRangeList[1]
1174 1186
1175 1187 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1176 1188
1177 1189 self.nProfiles = maxIndex - minIndex + 1
1178 1190 dataOut.nProfiles = self.nProfiles
1179 1191 dataOut.profileIndex = self.profileIndex
1180 1192 dataOut.flagNoData = False
1181 1193
1182 1194 self.incProfileIndex()
1183 1195 return dataOut
1184 1196
1185 1197 if rangeList != None:
1186 1198
1187 1199 nProfiles = 0
1188 1200
1189 1201 for thisRange in rangeList:
1190 1202 minIndex = thisRange[0]
1191 1203 maxIndex = thisRange[1]
1192 1204
1193 1205 nProfiles += maxIndex - minIndex + 1
1194 1206
1195 1207 for thisRange in rangeList:
1196 1208
1197 1209 minIndex = thisRange[0]
1198 1210 maxIndex = thisRange[1]
1199 1211
1200 1212 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1201 1213
1202 1214 self.nProfiles = nProfiles
1203 1215 dataOut.nProfiles = self.nProfiles
1204 1216 dataOut.profileIndex = self.profileIndex
1205 1217 dataOut.flagNoData = False
1206 1218
1207 1219 self.incProfileIndex()
1208 1220
1209 1221 break
1210 1222
1211 1223 return dataOut
1212 1224
1213 1225
1214 1226 if beam != None: #beam is only for AMISR data
1215 1227 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1216 1228 dataOut.flagNoData = False
1217 1229 dataOut.profileIndex = self.profileIndex
1218 1230
1219 1231 self.incProfileIndex()
1220 1232
1221 1233 return dataOut
1222 1234
1223 1235 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1224 1236
1225 1237
1226 1238 class Reshaper(Operation):
1227 1239
1228 1240 def __init__(self, **kwargs):
1229 1241
1230 1242 Operation.__init__(self, **kwargs)
1231 1243
1232 1244 self.__buffer = None
1233 1245 self.__nitems = 0
1234 1246
1235 1247 def __appendProfile(self, dataOut, nTxs):
1236 1248
1237 1249 if self.__buffer is None:
1238 1250 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1239 1251 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1240 1252
1241 1253 ini = dataOut.nHeights * self.__nitems
1242 1254 end = ini + dataOut.nHeights
1243 1255
1244 1256 self.__buffer[:, ini:end] = dataOut.data
1245 1257
1246 1258 self.__nitems += 1
1247 1259
1248 1260 return int(self.__nitems*nTxs)
1249 1261
1250 1262 def __getBuffer(self):
1251 1263
1252 1264 if self.__nitems == int(1./self.__nTxs):
1253 1265
1254 1266 self.__nitems = 0
1255 1267
1256 1268 return self.__buffer.copy()
1257 1269
1258 1270 return None
1259 1271
1260 1272 def __checkInputs(self, dataOut, shape, nTxs):
1261 1273
1262 1274 if shape is None and nTxs is None:
1263 1275 raise ValueError("Reshaper: shape of factor should be defined")
1264 1276
1265 1277 if nTxs:
1266 1278 if nTxs < 0:
1267 1279 raise ValueError("nTxs should be greater than 0")
1268 1280
1269 1281 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1270 1282 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1271 1283
1272 1284 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1273 1285
1274 1286 return shape, nTxs
1275 1287
1276 1288 if len(shape) != 2 and len(shape) != 3:
1277 1289 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1278 1290
1279 1291 if len(shape) == 2:
1280 1292 shape_tuple = [dataOut.nChannels]
1281 1293 shape_tuple.extend(shape)
1282 1294 else:
1283 1295 shape_tuple = list(shape)
1284 1296
1285 1297 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1286 1298
1287 1299 return shape_tuple, nTxs
1288 1300
1289 1301 def run(self, dataOut, shape=None, nTxs=None):
1290 1302
1291 1303 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1292 1304
1293 1305 dataOut.flagNoData = True
1294 1306 profileIndex = None
1295 1307
1296 1308 if dataOut.flagDataAsBlock:
1297 1309
1298 1310 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1299 1311 dataOut.flagNoData = False
1300 1312
1301 1313 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1302 1314
1303 1315 else:
1304 1316
1305 1317 if self.__nTxs < 1:
1306 1318
1307 1319 self.__appendProfile(dataOut, self.__nTxs)
1308 1320 new_data = self.__getBuffer()
1309 1321
1310 1322 if new_data is not None:
1311 1323 dataOut.data = new_data
1312 1324 dataOut.flagNoData = False
1313 1325
1314 1326 profileIndex = dataOut.profileIndex*nTxs
1315 1327
1316 1328 else:
1317 1329 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1318 1330
1319 1331 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1320 1332
1321 1333 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1322 1334
1323 1335 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1324 1336
1325 1337 dataOut.profileIndex = profileIndex
1326 1338
1327 1339 dataOut.ippSeconds /= self.__nTxs
1328 1340
1329 1341 return dataOut
1330 1342
1331 1343 class SplitProfiles(Operation):
1332 1344
1333 1345 def __init__(self, **kwargs):
1334 1346
1335 1347 Operation.__init__(self, **kwargs)
1336 1348
1337 1349 def run(self, dataOut, n):
1338 1350
1339 1351 dataOut.flagNoData = True
1340 1352 profileIndex = None
1341 1353
1342 1354 if dataOut.flagDataAsBlock:
1343 1355
1344 1356 #nchannels, nprofiles, nsamples
1345 1357 shape = dataOut.data.shape
1346 1358
1347 1359 if shape[2] % n != 0:
1348 1360 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1349 1361
1350 1362 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1351 1363
1352 1364 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1353 1365 dataOut.flagNoData = False
1354 1366
1355 1367 profileIndex = int(dataOut.nProfiles/n) - 1
1356 1368
1357 1369 else:
1358 1370
1359 1371 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1360 1372
1361 1373 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1362 1374
1363 1375 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1364 1376
1365 1377 dataOut.nProfiles = int(dataOut.nProfiles*n)
1366 1378
1367 1379 dataOut.profileIndex = profileIndex
1368 1380
1369 1381 dataOut.ippSeconds /= n
1370 1382
1371 1383 return dataOut
1372 1384
1373 1385 class CombineProfiles(Operation):
1374 1386 def __init__(self, **kwargs):
1375 1387
1376 1388 Operation.__init__(self, **kwargs)
1377 1389
1378 1390 self.__remData = None
1379 1391 self.__profileIndex = 0
1380 1392
1381 1393 def run(self, dataOut, n):
1382 1394
1383 1395 dataOut.flagNoData = True
1384 1396 profileIndex = None
1385 1397
1386 1398 if dataOut.flagDataAsBlock:
1387 1399
1388 1400 #nchannels, nprofiles, nsamples
1389 1401 shape = dataOut.data.shape
1390 1402 new_shape = shape[0], shape[1]/n, shape[2]*n
1391 1403
1392 1404 if shape[1] % n != 0:
1393 1405 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1394 1406
1395 1407 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1396 1408 dataOut.flagNoData = False
1397 1409
1398 1410 profileIndex = int(dataOut.nProfiles*n) - 1
1399 1411
1400 1412 else:
1401 1413
1402 1414 #nchannels, nsamples
1403 1415 if self.__remData is None:
1404 1416 newData = dataOut.data
1405 1417 else:
1406 1418 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1407 1419
1408 1420 self.__profileIndex += 1
1409 1421
1410 1422 if self.__profileIndex < n:
1411 1423 self.__remData = newData
1412 1424 #continue
1413 1425 return
1414 1426
1415 1427 self.__profileIndex = 0
1416 1428 self.__remData = None
1417 1429
1418 1430 dataOut.data = newData
1419 1431 dataOut.flagNoData = False
1420 1432
1421 1433 profileIndex = dataOut.profileIndex/n
1422 1434
1423 1435
1424 1436 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1425 1437
1426 1438 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1427 1439
1428 1440 dataOut.nProfiles = int(dataOut.nProfiles/n)
1429 1441
1430 1442 dataOut.profileIndex = profileIndex
1431 1443
1432 1444 dataOut.ippSeconds *= n
1433 1445
1434 1446 return dataOut
1435 1447
1436 1448 class PulsePairVoltage(Operation):
1437 1449 '''
1438 1450 Function PulsePair(Signal Power, Velocity)
1439 1451 The real component of Lag[0] provides Intensity Information
1440 1452 The imag component of Lag[1] Phase provides Velocity Information
1441 1453
1442 1454 Configuration Parameters:
1443 1455 nPRF = Number of Several PRF
1444 1456 theta = Degree Azimuth angel Boundaries
1445 1457
1446 1458 Input:
1447 1459 self.dataOut
1448 1460 lag[N]
1449 1461 Affected:
1450 1462 self.dataOut.spc
1451 1463 '''
1452 1464 isConfig = False
1453 1465 __profIndex = 0
1454 1466 __initime = None
1455 1467 __lastdatatime = None
1456 1468 __buffer = None
1457 1469 noise = None
1458 1470 __dataReady = False
1459 1471 n = None
1460 1472 __nch = 0
1461 1473 __nHeis = 0
1462 1474 removeDC = False
1463 1475 ipp = None
1464 1476 lambda_ = 0
1465 1477
1466 1478 def __init__(self,**kwargs):
1467 1479 Operation.__init__(self,**kwargs)
1468 1480
1469 1481 def setup(self, dataOut, n = None, removeDC=False):
1470 1482 '''
1471 1483 n= Numero de PRF's de entrada
1472 1484 '''
1473 1485 self.__initime = None
1474 1486 self.__lastdatatime = 0
1475 1487 self.__dataReady = False
1476 1488 self.__buffer = 0
1477 1489 self.__profIndex = 0
1478 1490 self.noise = None
1479 1491 self.__nch = dataOut.nChannels
1480 1492 self.__nHeis = dataOut.nHeights
1481 1493 self.removeDC = removeDC
1482 1494 self.lambda_ = 3.0e8/(9345.0e6)
1483 1495 self.ippSec = dataOut.ippSeconds
1484 1496 self.nCohInt = dataOut.nCohInt
1485 1497
1486 1498 if n == None:
1487 1499 raise ValueError("n should be specified.")
1488 1500
1489 1501 if n != None:
1490 1502 if n<2:
1491 1503 raise ValueError("n should be greater than 2")
1492 1504
1493 1505 self.n = n
1494 1506 self.__nProf = n
1495 1507
1496 1508 self.__buffer = numpy.zeros((dataOut.nChannels,
1497 1509 n,
1498 1510 dataOut.nHeights),
1499 1511 dtype='complex')
1500 1512
1501 1513 def putData(self,data):
1502 1514 '''
1503 1515 Add a profile to he __buffer and increase in one the __profiel Index
1504 1516 '''
1505 1517 self.__buffer[:,self.__profIndex,:]= data
1506 1518 self.__profIndex += 1
1507 1519 return
1508 1520
1509 1521 def pushData(self,dataOut):
1510 1522 '''
1511 1523 Return the PULSEPAIR and the profiles used in the operation
1512 1524 Affected : self.__profileIndex
1513 1525 '''
1514 1526 #----------------- Remove DC-----------------------------------
1515 1527 if self.removeDC==True:
1516 1528 mean = numpy.mean(self.__buffer,1)
1517 1529 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1518 1530 dc= numpy.tile(tmp,[1,self.__nProf,1])
1519 1531 self.__buffer = self.__buffer - dc
1520 1532 #------------------Calculo de Potencia ------------------------
1521 1533 pair0 = self.__buffer*numpy.conj(self.__buffer)
1522 1534 pair0 = pair0.real
1523 1535 lag_0 = numpy.sum(pair0,1)
1524 1536 #------------------Calculo de Ruido x canal--------------------
1525 1537 self.noise = numpy.zeros(self.__nch)
1526 1538 for i in range(self.__nch):
1527 1539 daux = numpy.sort(pair0[i,:,:],axis= None)
1528 1540 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1529 1541
1530 1542 self.noise = self.noise.reshape(self.__nch,1)
1531 1543 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1532 1544 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1533 1545 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1534 1546 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1535 1547 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1536 1548 #-------------------- Power --------------------------------------------------
1537 1549 data_power = lag_0/(self.n*self.nCohInt)
1538 1550 #------------------ Senal ---------------------------------------------------
1539 1551 data_intensity = pair0 - noise_buffer
1540 1552 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1541 1553 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1542 1554 for i in range(self.__nch):
1543 1555 for j in range(self.__nHeis):
1544 1556 if data_intensity[i][j] < 0:
1545 1557 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1546 1558
1547 1559 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1548 1560 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1549 1561 lag_1 = numpy.sum(pair1,1)
1550 1562 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1551 1563 data_velocity = (self.lambda_/2.0)*data_freq
1552 1564
1553 1565 #---------------- Potencia promedio estimada de la Senal-----------
1554 1566 lag_0 = lag_0/self.n
1555 1567 S = lag_0-self.noise
1556 1568
1557 1569 #---------------- Frecuencia Doppler promedio ---------------------
1558 1570 lag_1 = lag_1/(self.n-1)
1559 1571 R1 = numpy.abs(lag_1)
1560 1572
1561 1573 #---------------- Calculo del SNR----------------------------------
1562 1574 data_snrPP = S/self.noise
1563 1575 for i in range(self.__nch):
1564 1576 for j in range(self.__nHeis):
1565 1577 if data_snrPP[i][j] < 1.e-20:
1566 1578 data_snrPP[i][j] = 1.e-20
1567 1579
1568 1580 #----------------- Calculo del ancho espectral ----------------------
1569 1581 L = S/R1
1570 1582 L = numpy.where(L<0,1,L)
1571 1583 L = numpy.log(L)
1572 1584 tmp = numpy.sqrt(numpy.absolute(L))
1573 1585 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1574 1586 n = self.__profIndex
1575 1587
1576 1588 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1577 1589 self.__profIndex = 0
1578 1590 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1579 1591
1580 1592
1581 1593 def pulsePairbyProfiles(self,dataOut):
1582 1594
1583 1595 self.__dataReady = False
1584 1596 data_power = None
1585 1597 data_intensity = None
1586 1598 data_velocity = None
1587 1599 data_specwidth = None
1588 1600 data_snrPP = None
1589 1601 self.putData(data=dataOut.data)
1590 1602 if self.__profIndex == self.n:
1591 1603 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1592 1604 self.__dataReady = True
1593 1605
1594 1606 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1595 1607
1596 1608
1597 1609 def pulsePairOp(self, dataOut, datatime= None):
1598 1610
1599 1611 if self.__initime == None:
1600 1612 self.__initime = datatime
1601 1613 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1602 1614 self.__lastdatatime = datatime
1603 1615
1604 1616 if data_power is None:
1605 1617 return None, None, None,None,None,None
1606 1618
1607 1619 avgdatatime = self.__initime
1608 1620 deltatime = datatime - self.__lastdatatime
1609 1621 self.__initime = datatime
1610 1622
1611 1623 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1612 1624
1613 1625 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1614 1626
1615 1627 if not self.isConfig:
1616 1628 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1617 1629 self.isConfig = True
1618 1630 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1619 1631 dataOut.flagNoData = True
1620 1632
1621 1633 if self.__dataReady:
1622 1634 dataOut.nCohInt *= self.n
1623 1635 dataOut.dataPP_POW = data_intensity # S
1624 1636 dataOut.dataPP_POWER = data_power # P
1625 1637 dataOut.dataPP_DOP = data_velocity
1626 1638 dataOut.dataPP_SNR = data_snrPP
1627 1639 dataOut.dataPP_WIDTH = data_specwidth
1628 1640 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1629 1641 dataOut.utctime = avgdatatime
1630 1642 dataOut.flagNoData = False
1631 1643 return dataOut
1632 1644
1633 1645
1634 1646
1635 1647 # import collections
1636 1648 # from scipy.stats import mode
1637 1649 #
1638 1650 # class Synchronize(Operation):
1639 1651 #
1640 1652 # isConfig = False
1641 1653 # __profIndex = 0
1642 1654 #
1643 1655 # def __init__(self, **kwargs):
1644 1656 #
1645 1657 # Operation.__init__(self, **kwargs)
1646 1658 # # self.isConfig = False
1647 1659 # self.__powBuffer = None
1648 1660 # self.__startIndex = 0
1649 1661 # self.__pulseFound = False
1650 1662 #
1651 1663 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1652 1664 #
1653 1665 # #Read data
1654 1666 #
1655 1667 # powerdB = dataOut.getPower(channel = channel)
1656 1668 # noisedB = dataOut.getNoise(channel = channel)[0]
1657 1669 #
1658 1670 # self.__powBuffer.extend(powerdB.flatten())
1659 1671 #
1660 1672 # dataArray = numpy.array(self.__powBuffer)
1661 1673 #
1662 1674 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1663 1675 #
1664 1676 # maxValue = numpy.nanmax(filteredPower)
1665 1677 #
1666 1678 # if maxValue < noisedB + 10:
1667 1679 # #No se encuentra ningun pulso de transmision
1668 1680 # return None
1669 1681 #
1670 1682 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1671 1683 #
1672 1684 # if len(maxValuesIndex) < 2:
1673 1685 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1674 1686 # return None
1675 1687 #
1676 1688 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1677 1689 #
1678 1690 # #Seleccionar solo valores con un espaciamiento de nSamples
1679 1691 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1680 1692 #
1681 1693 # if len(pulseIndex) < 2:
1682 1694 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1683 1695 # return None
1684 1696 #
1685 1697 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1686 1698 #
1687 1699 # #remover senales que se distancien menos de 10 unidades o muestras
1688 1700 # #(No deberian existir IPP menor a 10 unidades)
1689 1701 #
1690 1702 # realIndex = numpy.where(spacing > 10 )[0]
1691 1703 #
1692 1704 # if len(realIndex) < 2:
1693 1705 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1694 1706 # return None
1695 1707 #
1696 1708 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1697 1709 # realPulseIndex = pulseIndex[realIndex]
1698 1710 #
1699 1711 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1700 1712 #
1701 1713 # print "IPP = %d samples" %period
1702 1714 #
1703 1715 # self.__newNSamples = dataOut.nHeights #int(period)
1704 1716 # self.__startIndex = int(realPulseIndex[0])
1705 1717 #
1706 1718 # return 1
1707 1719 #
1708 1720 #
1709 1721 # def setup(self, nSamples, nChannels, buffer_size = 4):
1710 1722 #
1711 1723 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1712 1724 # maxlen = buffer_size*nSamples)
1713 1725 #
1714 1726 # bufferList = []
1715 1727 #
1716 1728 # for i in range(nChannels):
1717 1729 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1718 1730 # maxlen = buffer_size*nSamples)
1719 1731 #
1720 1732 # bufferList.append(bufferByChannel)
1721 1733 #
1722 1734 # self.__nSamples = nSamples
1723 1735 # self.__nChannels = nChannels
1724 1736 # self.__bufferList = bufferList
1725 1737 #
1726 1738 # def run(self, dataOut, channel = 0):
1727 1739 #
1728 1740 # if not self.isConfig:
1729 1741 # nSamples = dataOut.nHeights
1730 1742 # nChannels = dataOut.nChannels
1731 1743 # self.setup(nSamples, nChannels)
1732 1744 # self.isConfig = True
1733 1745 #
1734 1746 # #Append new data to internal buffer
1735 1747 # for thisChannel in range(self.__nChannels):
1736 1748 # bufferByChannel = self.__bufferList[thisChannel]
1737 1749 # bufferByChannel.extend(dataOut.data[thisChannel])
1738 1750 #
1739 1751 # if self.__pulseFound:
1740 1752 # self.__startIndex -= self.__nSamples
1741 1753 #
1742 1754 # #Finding Tx Pulse
1743 1755 # if not self.__pulseFound:
1744 1756 # indexFound = self.__findTxPulse(dataOut, channel)
1745 1757 #
1746 1758 # if indexFound == None:
1747 1759 # dataOut.flagNoData = True
1748 1760 # return
1749 1761 #
1750 1762 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1751 1763 # self.__pulseFound = True
1752 1764 # self.__startIndex = indexFound
1753 1765 #
1754 1766 # #If pulse was found ...
1755 1767 # for thisChannel in range(self.__nChannels):
1756 1768 # bufferByChannel = self.__bufferList[thisChannel]
1757 1769 # #print self.__startIndex
1758 1770 # x = numpy.array(bufferByChannel)
1759 1771 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1760 1772 #
1761 1773 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1762 1774 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1763 1775 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1764 1776 #
1765 1777 # dataOut.data = self.__arrayBuffer
1766 1778 #
1767 1779 # self.__startIndex += self.__newNSamples
1768 1780 #
1769 1781 # return
1770 1782 class SSheightProfiles(Operation):
1771 1783
1772 1784 step = None
1773 1785 nsamples = None
1774 1786 bufferShape = None
1775 1787 profileShape = None
1776 1788 sshProfiles = None
1777 1789 profileIndex = None
1778 1790
1779 1791 def __init__(self, **kwargs):
1780 1792
1781 1793 Operation.__init__(self, **kwargs)
1782 1794 self.isConfig = False
1783 1795
1784 1796 def setup(self,dataOut ,step = None , nsamples = None):
1785 1797
1786 1798 if step == None and nsamples == None:
1787 1799 raise ValueError("step or nheights should be specified ...")
1788 1800
1789 1801 self.step = step
1790 1802 self.nsamples = nsamples
1791 1803 self.__nChannels = dataOut.nChannels
1792 1804 self.__nProfiles = dataOut.nProfiles
1793 1805 self.__nHeis = dataOut.nHeights
1794 1806 shape = dataOut.data.shape #nchannels, nprofiles, nsamples
1795 1807
1796 1808 residue = (shape[1] - self.nsamples) % self.step
1797 1809 if residue != 0:
1798 1810 print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
1799 1811
1800 1812 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1801 1813 numberProfile = self.nsamples
1802 1814 numberSamples = (shape[1] - self.nsamples)/self.step
1803 1815
1804 1816 self.bufferShape = int(shape[0]), int(numberSamples), int(numberProfile) # nchannels, nsamples , nprofiles
1805 1817 self.profileShape = int(shape[0]), int(numberProfile), int(numberSamples) # nchannels, nprofiles, nsamples
1806 1818
1807 1819 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1808 1820 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
1809 1821
1810 1822 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1811 1823 dataOut.flagNoData = True
1812 1824
1813 1825 profileIndex = None
1814 1826 #print("nProfiles, nHeights ",dataOut.nProfiles, dataOut.nHeights)
1815 1827 #print(dataOut.getFreqRange(1)/1000.)
1816 1828 #exit(1)
1817 1829 if dataOut.flagDataAsBlock:
1818 1830 dataOut.data = numpy.average(dataOut.data,axis=1)
1819 1831 #print("jee")
1820 1832 dataOut.flagDataAsBlock = False
1821 1833 if not self.isConfig:
1822 1834 self.setup(dataOut, step=step , nsamples=nsamples)
1823 1835 #print("Setup done")
1824 1836 self.isConfig = True
1825 1837
1826 1838
1827 1839 if code is not None:
1828 1840 code = numpy.array(code)
1829 1841 code_block = code
1830 1842
1831 1843 if repeat is not None:
1832 1844 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1833 1845 #print(code_block.shape)
1834 1846 for i in range(self.buffer.shape[1]):
1835 1847
1836 1848 if code is not None:
1837 1849 self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]*code_block
1838 1850
1839 1851 else:
1840 1852
1841 1853 self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1842 1854
1843 1855 #self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights])
1844 1856
1845 1857 for j in range(self.buffer.shape[0]):
1846 1858 self.sshProfiles[j] = numpy.transpose(self.buffer[j])
1847 1859
1848 1860 profileIndex = self.nsamples
1849 1861 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1850 1862 ippSeconds = (deltaHeight*1.0e-6)/(0.15)
1851 1863 #print("ippSeconds, dH: ",ippSeconds,deltaHeight)
1852 1864 try:
1853 1865 if dataOut.concat_m is not None:
1854 1866 ippSeconds= ippSeconds/float(dataOut.concat_m)
1855 1867 #print "Profile concat %d"%dataOut.concat_m
1856 1868 except:
1857 1869 pass
1858 1870
1859 1871 dataOut.data = self.sshProfiles
1860 1872 dataOut.flagNoData = False
1861 1873 dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0]
1862 1874 dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
1863 1875
1864 1876 dataOut.profileIndex = profileIndex
1865 1877 dataOut.flagDataAsBlock = True
1866 1878 dataOut.ippSeconds = ippSeconds
1867 1879 dataOut.step = self.step
1868 1880 #print(numpy.shape(dataOut.data))
1869 1881 #exit(1)
1870 1882 #print("new data shape and time:", dataOut.data.shape, dataOut.utctime)
1871 1883
1872 1884 return dataOut
1873 1885 ################################################################################3############################3
1874 1886 ################################################################################3############################3
1875 1887 ################################################################################3############################3
1876 1888 ################################################################################3############################3
1877 1889
1878 1890 class SSheightProfiles2(Operation):
1879 1891 '''
1880 1892 Procesa por perfiles y por bloques
1881 1893 Versión corregida y actualizada para trabajar con RemoveProfileSats2
1882 1894 Usar esto
1883 1895 '''
1884 1896
1885 1897
1886 1898 bufferShape = None
1887 1899 profileShape = None
1888 1900 sshProfiles = None
1889 1901 profileIndex = None
1890 1902 #nsamples = None
1891 1903 #step = None
1892 1904 #deltaHeight = None
1893 1905 #init_range = None
1894 1906 __slots__ = ('step', 'nsamples', 'deltaHeight', 'init_range', 'isConfig', '__nChannels',
1895 1907 '__nProfiles', '__nHeis', 'deltaHeight', 'new_nHeights')
1896 1908
1897 1909 def __init__(self, **kwargs):
1898 1910
1899 1911 Operation.__init__(self, **kwargs)
1900 1912 self.isConfig = False
1901 1913
1902 1914 def setup(self,dataOut ,step = None , nsamples = None):
1903 1915
1904 1916 if step == None and nsamples == None:
1905 1917 raise ValueError("step or nheights should be specified ...")
1906 1918
1907 1919 self.step = step
1908 1920 self.nsamples = nsamples
1909 1921 self.__nChannels = int(dataOut.nChannels)
1910 1922 self.__nProfiles = int(dataOut.nProfiles)
1911 1923 self.__nHeis = int(dataOut.nHeights)
1912 1924
1913 1925 residue = (self.__nHeis - self.nsamples) % self.step
1914 1926 if residue != 0:
1915 1927 print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,self.__nProfiles - self.nsamples,residue))
1916 1928
1917 1929 self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1918 1930 self.init_range = dataOut.heightList[0]
1919 1931 #numberProfile = self.nsamples
1920 1932 numberSamples = (self.__nHeis - self.nsamples)/self.step
1921 1933
1922 1934 self.new_nHeights = numberSamples
1923 1935
1924 1936 self.bufferShape = int(self.__nChannels), int(numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles
1925 1937 self.profileShape = int(self.__nChannels), int(self.nsamples), int(numberSamples) # nchannels, nprofiles, nsamples
1926 1938
1927 1939 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1928 1940 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
1929 1941
1930 1942 def getNewProfiles(self, data, code=None, repeat=None):
1931 1943
1932 1944 if code is not None:
1933 1945 code = numpy.array(code)
1934 1946 code_block = code
1935 1947
1936 1948 if repeat is not None:
1937 1949 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1938 1950 if data.ndim < 3:
1939 1951 data = data.reshape(self.__nChannels,1,self.__nHeis )
1940 1952 #print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape, code_block.shape)
1941 1953 for ch in range(self.__nChannels):
1942 1954 for i in range(int(self.new_nHeights)): #nuevas alturas
1943 1955 if code is not None:
1944 1956 self.buffer[ch,i,:] = data[ch,:,i*self.step:i*self.step + self.nsamples]*code_block
1945 1957 else:
1946 1958 self.buffer[ch,i,:] = data[ch,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1947 1959
1948 1960 for j in range(self.__nChannels): #en los cananles
1949 1961 self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:])
1950 1962 #print("new profs Done")
1951 1963
1952 1964
1953 1965
1954 1966 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1955 1967 # print("running")
1956 1968 if dataOut.flagNoData == True:
1957 1969 return dataOut
1958 1970 dataOut.flagNoData = True
1959 1971 #print("init data shape:", dataOut.data.shape)
1960 1972 #print("ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
1961 1973 # int(dataOut.nProfiles),int(dataOut.nHeights)))
1962 1974
1963 1975 profileIndex = None
1964 1976 # if not dataOut.flagDataAsBlock:
1965 1977 # dataOut.nProfiles = 1
1966 1978
1967 1979 if not self.isConfig:
1968 1980 self.setup(dataOut, step=step , nsamples=nsamples)
1969 1981 #print("Setup done")
1970 1982 self.isConfig = True
1971 1983
1972 1984 dataBlock = None
1973 1985
1974 1986 nprof = 1
1975 1987 if dataOut.flagDataAsBlock:
1976 1988 nprof = int(dataOut.nProfiles)
1977 1989
1978 1990 #print("dataOut nProfiles:", dataOut.nProfiles)
1979 1991 for profile in range(nprof):
1980 1992 if dataOut.flagDataAsBlock:
1981 1993 #print("read blocks")
1982 1994 self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat)
1983 1995 else:
1984 1996 #print("read profiles")
1985 1997 self.getNewProfiles(dataOut.data, code=code, repeat=repeat) #only one channe
1986 1998 if profile == 0:
1987 1999 dataBlock = self.sshProfiles.copy()
1988 2000 else: #by blocks
1989 2001 dataBlock = numpy.concatenate((dataBlock,self.sshProfiles), axis=1) #profile axis
1990 2002 #print("by blocks: ",dataBlock.shape, self.sshProfiles.shape)
1991 2003
1992 2004 profileIndex = self.nsamples
1993 2005 #deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1994 2006 ippSeconds = (self.deltaHeight*1.0e-6)/(0.15)
1995 2007
1996 2008
1997 2009 dataOut.data = dataBlock
1998 2010 #print("show me: ",self.step,self.deltaHeight, dataOut.heightList, self.new_nHeights)
1999 2011 dataOut.heightList = numpy.arange(int(self.new_nHeights)) *self.step*self.deltaHeight + self.init_range
2000 2012 dataOut.sampled_heightsFFT = self.nsamples
2001 2013 dataOut.ippSeconds = ippSeconds
2002 2014 dataOut.step = self.step
2003 2015 dataOut.deltaHeight = self.step*self.deltaHeight
2004 2016 dataOut.flagNoData = False
2005 2017 if dataOut.flagDataAsBlock:
2006 2018 dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
2007 2019
2008 2020 else:
2009 2021 dataOut.nProfiles = int(self.nsamples)
2010 2022 dataOut.profileIndex = dataOut.nProfiles
2011 2023 dataOut.flagDataAsBlock = True
2012 2024
2013 2025 dataBlock = None
2014 2026
2015 2027 #print("new data shape:", dataOut.data.shape, dataOut.utctime)
2016 2028
2017 2029 #update Processing Header:
2018 2030 dataOut.processingHeaderObj.heightList = dataOut.heightList
2019 2031 dataOut.processingHeaderObj.ipp = ippSeconds
2020 2032 dataOut.processingHeaderObj.heightResolution = dataOut.deltaHeight
2021 2033 #dataOut.processingHeaderObj.profilesPerBlock = nProfiles
2022 2034
2023 2035 # # dataOut.data = CH, PROFILES, HEIGHTS
2024 2036 #print(dataOut.data .shape)
2025 2037 if dataOut.flagProfilesByRange:
2026 2038 # #assuming the same remotion for all channels
2027 2039 aux = [ self.nsamples - numpy.count_nonzero(dataOut.data[0, :, h]==0) for h in range(len(dataOut.heightList))]
2028 2040 dataOut.nProfilesByRange = (numpy.asarray(aux)).reshape((1,len(dataOut.heightList) ))
2029 2041 #print(dataOut.nProfilesByRange.shape)
2030 2042 else:
2031 2043 dataOut.nProfilesByRange = numpy.ones((1, len(dataOut.heightList)))*dataOut.nProfiles
2032 2044 return dataOut
2033 2045
2034 2046
2035 2047
2036 2048
2037 2049
2038 2050 class removeProfileByFaradayHS(Operation):
2039 2051 '''
2040 2052
2041 2053 '''
2042 2054
2043 2055 __buffer_data = []
2044 2056 __buffer_times = []
2045 2057
2046 2058 buffer = None
2047 2059
2048 2060 outliers_IDs_list = []
2049 2061
2050 2062
2051 2063 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
2052 2064 '__dh','first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
2053 2065 '__count_exec','__initime','__dataReady','__ipp')
2054 2066 def __init__(self, **kwargs):
2055 2067
2056 2068 Operation.__init__(self, **kwargs)
2057 2069 self.isConfig = False
2058 2070
2059 2071 def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=3, minHei=None, maxHei=None):
2060 2072
2061 2073 if n == None and timeInterval == None:
2062 2074 raise ValueError("nprofiles or timeInterval should be specified ...")
2063 2075
2064 2076 if n != None:
2065 2077 self.n = n
2066 2078
2067 2079 self.navg = navg
2068 2080 self.profileMargin = profileMargin
2069 2081 self.thHistOutlier = thHistOutlier
2070 2082 self.__profIndex = 0
2071 2083 self.buffer = None
2072 2084 self._ipp = dataOut.ippSeconds
2073 2085 self.n_prof_released = 0
2074 2086 self.heightList = dataOut.heightList
2075 2087 self.init_prof = 0
2076 2088 self.end_prof = 0
2077 2089 self.__count_exec = 0
2078 2090 self.__profIndex = 0
2079 2091 self.first_utcBlock = None
2080 2092 self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
2081 2093 minHei = minHei
2082 2094 maxHei = maxHei
2083 2095 if minHei==None :
2084 2096 minHei = dataOut.heightList[0]
2085 2097 if maxHei==None :
2086 2098 maxHei = dataOut.heightList[-1]
2087 2099 self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
2088 2100
2089 2101 self.nChannels = dataOut.nChannels
2090 2102 self.nHeights = dataOut.nHeights
2091 2103 self.test_counter = 0
2092 2104
2093 2105 def filterSatsProfiles(self):
2094 2106 data = self.__buffer_data
2095 2107 #print(data.shape)
2096 2108 nChannels, profiles, heights = data.shape
2097 2109 indexes=[]
2098 2110 outliers_IDs=[]
2099 2111 for c in range(nChannels):
2100 2112 for h in range(self.minHei_idx, self.maxHei_idx):
2101 2113 power = 10* numpy.log10((data[c,:,h] * numpy.conjugate(data[c,:,h])).real)
2102 2114 #power = power.real
2103 2115 #power = (numpy.abs(data[c,:,h].real))
2104 2116 sortdata = numpy.sort(power, axis=None)
2105 2117 sortID=power.argsort()
2106 2118 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
2107 2119
2108 2120 indexes.append(index)
2109 2121 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
2110 2122
2111 2123 # print(sortdata.min(), sortdata.max(), sortdata.mean())
2112 2124 # fig,ax = plt.subplots()
2113 2125 # #ax.set_title(str(k)+" "+str(j))
2114 2126 # x=range(len(sortdata))
2115 2127 # ax.scatter(x,sortdata)
2116 2128 # ax.axvline(index)
2117 2129 # plt.grid()
2118 2130 # plt.show()
2119 2131
2120 2132
2121 2133
2122 2134
2123 2135 outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
2124 2136 outliers_IDs = numpy.unique(outliers_IDs)
2125 2137 outs_lines = numpy.sort(outliers_IDs)
2126 2138 # #print("outliers Ids: ", outs_lines, outs_lines.shape)
2127 2139 #hist, bin_edges = numpy.histogram(outs_lines, bins=10, density=True)
2128 2140
2129 2141
2130 2142 #Agrupando el histograma de outliers,
2131 2143 my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=False)
2132 2144 #my_bins = numpy.linspace(0,1600, 96, endpoint=False)
2133 2145
2134 2146 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
2135 2147 hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
2136 2148 #print(hist_outliers_indexes[0])
2137 2149 bins_outliers_indexes = [int(i) for i in bins[hist_outliers_indexes]] #
2138 2150 #print(bins_outliers_indexes)
2139 2151 outlier_loc_index = []
2140 2152
2141 2153
2142 2154 # for n in range(len(bins_outliers_indexes)-1):
2143 2155 # for k in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin):
2144 2156 # outlier_loc_index.append(k)
2145 2157
2146 2158 outlier_loc_index = [e for n in range(len(bins_outliers_indexes)-1) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin) ]
2147 2159
2148 2160 outlier_loc_index = numpy.asarray(outlier_loc_index)
2149 2161 #print(len(numpy.unique(outlier_loc_index)), numpy.unique(outlier_loc_index))
2150 2162
2151 2163
2152 2164
2153 2165 x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
2154 2166 fig, ax = plt.subplots(1,2,figsize=(8, 6))
2155 2167
2156 2168 dat = data[0,:,:].real
2157 2169 dat = 10* numpy.log10((data[0,:,:] * numpy.conjugate(data[0,:,:])).real)
2158 2170 m = numpy.nanmean(dat)
2159 2171 o = numpy.nanstd(dat)
2160 2172 #print(m, o, x.shape, y.shape)
2161 2173 c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
2162 2174 ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
2163 2175 fig.colorbar(c)
2164 2176 ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
2165 2177 ax[1].hist(outs_lines,bins=my_bins)
2166 2178 plt.show()
2167 2179
2168 2180
2169 2181 self.outliers_IDs_list = numpy.unique(outlier_loc_index)
2170 2182 print("outs list: ", self.outliers_IDs_list)
2171 2183 return data
2172 2184
2173 2185 def filterSatsProfiles2(self):
2174 2186 data = self.__buffer_data
2175 2187 #print(data.shape)
2176 2188 nChannels, profiles, heights = data.shape
2177 2189 indexes=numpy.zeros([], dtype=int)
2178 2190 outliers_IDs=[]
2179 2191 for c in range(nChannels):
2180 2192 noise_ref =10* numpy.log10((data[c,:,550:600] * numpy.conjugate(data[c,:,550:600])).real)
2181 2193 print("Noise ",noise_ref.mean())
2182 2194 for h in range(self.minHei_idx, self.maxHei_idx):
2183 2195 power = 10* numpy.log10((data[c,:,h] * numpy.conjugate(data[c,:,h])).real)
2184 2196 #power = power.real
2185 2197 #power = (numpy.abs(data[c,:,h].real))
2186 2198 #sortdata = numpy.sort(power, axis=None)
2187 2199 #sortID=power.argsort()
2188 2200 #print(sortID)
2189 2201 th = 60 + 10
2190 2202 index = numpy.where(power > th )
2191 2203 if index[0].size > 10 and index[0].size < int(0.8*profiles):
2192 2204 indexes = numpy.append(indexes, index[0])
2193 2205 #print(index[0])
2194 2206 #print(index[0])
2195 2207
2196 2208 # fig,ax = plt.subplots()
2197 2209 # #ax.set_title(str(k)+" "+str(j))
2198 2210 # x=range(len(power))
2199 2211 # ax.scatter(x,power)
2200 2212 # #ax.axvline(index)
2201 2213 # plt.grid()
2202 2214 # plt.show()
2203 2215 #print(indexes)
2204 2216
2205 2217 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
2206 2218 #outliers_IDs = numpy.unique(outliers_IDs)
2207 2219
2208 2220 outs_lines = numpy.unique(indexes)
2209 2221 print("outliers Ids: ", outs_lines, outs_lines.shape)
2210 2222 #hist, bin_edges = numpy.histogram(outs_lines, bins=10, density=True)
2211 2223
2212 2224
2213 2225 #Agrupando el histograma de outliers,
2214 2226 my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=False)
2215 2227 #my_bins = numpy.linspace(0,1600, 96, endpoint=False)
2216 2228
2217 2229 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
2218 2230 hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
2219 2231 #print(hist_outliers_indexes[0])
2220 2232 bins_outliers_indexes = [int(i) for i in bins[hist_outliers_indexes]] #
2221 2233 #print(bins_outliers_indexes)
2222 2234 outlier_loc_index = []
2223 2235
2224 2236
2225 2237
2226 2238 outlier_loc_index = [e for n in range(len(bins_outliers_indexes)-1) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin) ]
2227 2239
2228 2240 outlier_loc_index = numpy.asarray(outlier_loc_index)
2229 2241 outlier_loc_index = outlier_loc_index[~numpy.all(outlier_loc_index < 0)]
2230 2242
2231 2243 print("outliers final: ", outlier_loc_index)
2232 2244
2233 2245 x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
2234 2246 fig, ax = plt.subplots(1,2,figsize=(8, 6))
2235 2247
2236 2248 dat = data[0,:,:].real
2237 2249 dat = 10* numpy.log10((data[0,:,:] * numpy.conjugate(data[0,:,:])).real)
2238 2250 m = numpy.nanmean(dat)
2239 2251 o = numpy.nanstd(dat)
2240 2252 #print(m, o, x.shape, y.shape)
2241 2253 c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
2242 2254 ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
2243 2255 fig.colorbar(c)
2244 2256 ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
2245 2257 ax[1].hist(outs_lines,bins=my_bins)
2246 2258 plt.show()
2247 2259
2248 2260
2249 2261 self.outliers_IDs_list = numpy.unique(outlier_loc_index)
2250 2262 print("outs list: ", self.outliers_IDs_list)
2251 2263 return data
2252 2264
2253 2265 def cleanSpikesFFT2D(self):
2254 2266 incoh_int = 10
2255 2267 norm_img = 75
2256 2268 import matplotlib.pyplot as plt
2257 2269 import datetime
2258 2270 import cv2
2259 2271 data = self.__buffer_data
2260 2272 print("cleaning shape inpt: ",data.shape)
2261 2273 self.__buffer_data = []
2262 2274
2263 2275
2264 2276 channels , profiles, heights = data.shape
2265 2277 len_split_prof = profiles / incoh_int
2266 2278
2267 2279
2268 2280 for ch in range(channels):
2269 2281 data_10 = numpy.split(data[ch, :, self.minHei_idx:], incoh_int, axis=0) # división de los perfiles
2270 2282 print("splited data: ",len(data_10)," -> ", data_10[0].shape)
2271 2283 int_img = None
2272 2284 i_count = 0
2273 2285 n_x, n_y = data_10[0].shape
2274 2286 for s_data in data_10: #porciones de espectro
2275 2287 spectrum = numpy.fft.fft2(s_data, axes=(0,1))
2276 2288 z = numpy.abs(spectrum)
2277 2289 mg = z[2:n_y,:] #omitir dc y adjunto
2278 2290 dat = numpy.log10(mg.T)
2279 2291 i_count += 1
2280 2292 if i_count == 1:
2281 2293 int_img = dat
2282 2294 else:
2283 2295 int_img += dat
2284 2296 #print(i_count)
2285 2297
2286 2298 min, max = int_img.min(), int_img.max()
2287 2299 int_img = ((int_img-min)*255/(max-min)).astype(numpy.uint8)
2288 2300
2289 2301 cv2.imshow('integrated image', int_img) #numpy.fft.fftshift(img))
2290 2302 cv2.waitKey(0)
2291 2303 #####################################################################
2292 2304 kernel_h = numpy.zeros((3,3)) #
2293 2305 kernel_h[0, :] = -2
2294 2306 kernel_h[1, :] = 3
2295 2307 kernel_h[2, :] = -2
2296 2308
2297 2309
2298 2310 kernel_5h = numpy.zeros((5,5)) #
2299 2311 kernel_5h[0, :] = -2
2300 2312 kernel_5h[1, :] = -1
2301 2313 kernel_5h[2, :] = 5
2302 2314 kernel_5h[3, :] = -1
2303 2315 kernel_5h[4, :] = -2
2304 2316
2305 2317 #####################################################################
2306 2318 sharp_img = cv2.filter2D(src=int_img, ddepth=-1, kernel=kernel_5h)
2307 2319 # cv2.imshow('sharp image h ', sharp_img)
2308 2320 # cv2.waitKey(0)
2309 2321 #####################################################################
2310 2322 horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5,1)) #11
2311 2323 #####################################################################
2312 2324 detected_lines_h = cv2.morphologyEx(sharp_img, cv2.MORPH_OPEN, horizontal_kernel, iterations=1)
2313 2325 # cv2.imshow('lines horizontal', detected_lines_h) #numpy.fft.fftshift(detected_lines_h))
2314 2326 # cv2.waitKey(0)
2315 2327 #####################################################################
2316 2328 ret, detected_lines_h = cv2.threshold(detected_lines_h, 200, 255, cv2.THRESH_BINARY)#
2317 2329 cv2.imshow('binary img', detected_lines_h) #numpy.fft.fftshift(detected_lines_h))
2318 2330 cv2.waitKey(0)
2319 2331 #####################################################################
2320 2332 cnts_h, h0 = cv2.findContours(detected_lines_h, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2321 2333 #####################################################################
2322 2334 h_line_index = []
2323 2335 v_line_index = []
2324 2336
2325 2337 #cnts_h += cnts_h_s #combine large and small lines
2326 2338
2327 2339 # line indexes x1, x2, y
2328 2340 for c in cnts_h:
2329 2341 #print(c)
2330 2342 if len(c) < 3: #contorno linea
2331 2343 x1 = c[0][0][0]
2332 2344 x2 = c[1][0][0]
2333 2345 if x1 > 5 and x2 < (n_x-5) :
2334 2346 start = incoh_int + (x1 * incoh_int)
2335 2347 end = incoh_int + (x2 * incoh_int)
2336 2348 h_line_index.append( [start, end, c[0][0][1]] )
2337 2349
2338 2350 #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1])
2339 2351 else: #contorno poligono
2340 2352 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2341 2353 y = numpy.unique(pairs[:,1])
2342 2354 x = numpy.unique(pairs[:,0])
2343 2355 #print(x)
2344 2356 for yk in y:
2345 2357 x0 = x[0]
2346 2358 if x0 < 8:
2347 2359 x0 = 10
2348 2360 #print(x[0], x[-1], yk)
2349 2361 h_line_index.append( [x0, x[-1], yk])
2350 2362 #print("x1, x2, y ->p ", x[0], x[-1], yk)
2351 2363 ###################################################################
2352 2364 #print("Cleaning")
2353 2365 # # clean Spectrum
2354 2366 spectrum = numpy.fft.fft2(data[ch,:,self.minHei_idx:], axes=(0,1))
2355 2367 z = numpy.abs(spectrum)
2356 2368 phase = numpy.angle(spectrum)
2357 2369 print("Total Horizontal", len(h_line_index))
2358 2370 if len(h_line_index) < 75 :
2359 2371 for x1, x2, y in h_line_index:
2360 2372 print(x1, x2, y)
2361 2373 z[x1:x2,y] = 0
2362 2374
2363 2375
2364 2376 spcCleaned = z * numpy.exp(1j*phase)
2365 2377
2366 2378 dat2 = numpy.log10(z[1:-1,:].T)
2367 2379 min, max =dat2.min(), dat2.max()
2368 2380 print(min, max)
2369 2381 img2 = ((dat2-min)*255/(max-min)).astype(numpy.uint8)
2370 2382 cv2.imshow('cleaned', img2) #numpy.fft.fftshift(img_cleaned))
2371 2383 cv2.waitKey(0)
2372 2384 cv2.destroyAllWindows()
2373 2385
2374 2386 data[ch,:,self.minHei_idx:] = numpy.fft.ifft2(spcCleaned, axes=(0,1))
2375 2387
2376 2388
2377 2389 #print("cleanOutliersByBlock Done", data.shape)
2378 2390 self.__buffer_data = data
2379 2391 return data
2380 2392
2381 2393
2382 2394
2383 2395
2384 2396 def cleanOutliersByBlock(self):
2385 2397 import matplotlib.pyplot as plt
2386 2398 import datetime
2387 2399 import cv2
2388 2400 #print(self.__buffer_data[0].shape)
2389 2401 data = self.__buffer_data#.copy()
2390 2402 print("cleaning shape inpt: ",data.shape)
2391 2403 self.__buffer_data = []
2392 2404
2393 2405
2394 2406 spectrum = numpy.fft.fft2(data[:,:,self.minHei_idx:], axes=(1,2))
2395 2407 print("spc : ",spectrum.shape)
2396 2408 (nch,nsamples, nh) = spectrum.shape
2397 2409 data2 = None
2398 2410 #print(data.shape)
2399 2411 cleanedBlock = None
2400 2412 spectrum2 = spectrum.copy()
2401 2413 for ch in range(nch):
2402 2414 dh = self.__dh
2403 2415 dt1 = (dh*1.0e-6)/(0.15)
2404 2416 dt2 = self.__buffer_times[1]-self.__buffer_times[0]
2405 2417
2406 2418 freqv = numpy.fft.fftfreq(nh, d=dt1)
2407 2419 freqh = numpy.fft.fftfreq(self.n, d=dt2)
2408 2420
2409 2421 z = numpy.abs(spectrum[ch,:,:])
2410 2422 phase = numpy.angle(spectrum[ch,:,:])
2411 2423 z1 = z[0,:]
2412 2424 #print("shape z: ", z.shape, nsamples)
2413 2425
2414 2426 dat = numpy.log10(z[1:nsamples,:].T)
2415 2427
2416 2428 pdat = numpy.log10(phase.T)
2417 2429 #print("dat mean",dat.mean())
2418 2430
2419 2431 mean, min, max = dat.mean(), dat.min(), dat.max()
2420 2432 img = ((dat-min)*200/(max-min)).astype(numpy.uint8)
2421 2433
2422 2434 # print(img.shape)
2423 2435 cv2.imshow('image', img) #numpy.fft.fftshift(img))
2424 2436 cv2.waitKey(0)
2425 2437
2426 2438
2427 2439 ''' #FUNCIONA LINEAS PEQUEÑAS
2428 2440 kernel_5h = numpy.zeros((5,3)) #
2429 2441 kernel_5h[0, :] = 2
2430 2442 kernel_5h[1, :] = 1
2431 2443 kernel_5h[2, :] = 0
2432 2444 kernel_5h[3, :] = -1
2433 2445 kernel_5h[4, :] = -2
2434 2446
2435 2447 sharp_imgh = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_5h)
2436 2448 cv2.imshow('sharp image h',sharp_imgh)
2437 2449 cv2.waitKey(0)
2438 2450 horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (20,1))
2439 2451
2440 2452 detected_lines_h = cv2.morphologyEx(sharp_imgh, cv2.MORPH_OPEN, horizontal_kernel, iterations=1)
2441 2453 #detected_lines_h = cv2.medianBlur(detected_lines_h, 3)
2442 2454 #detected_lines_h = cv2.filter2D(src=img, ddepth=-1, kernel=kernel)
2443 2455 cv2.imshow('lines h gray', detected_lines_h)
2444 2456 cv2.waitKey(0)
2445 2457 reth, detected_lines_h = cv2.threshold(detected_lines_h, 90, 255, cv2.THRESH_BINARY)
2446 2458 cv2.imshow('lines h ', detected_lines_h)
2447 2459 cv2.waitKey(0)
2448 2460 '''
2449 2461
2450 2462
2451 2463 '''
2452 2464 kernel_3h = numpy.zeros((3,10)) #10
2453 2465 kernel_3h[0, :] = -1
2454 2466 kernel_3h[1, :] = 2
2455 2467 kernel_3h[2, :] = -1
2456 2468
2457 2469
2458 2470 kernel_h = numpy.zeros((3,20)) #20
2459 2471 kernel_h[0, :] = -1
2460 2472 kernel_h[1, :] = 2
2461 2473 kernel_h[2, :] = -1
2462 2474
2463 2475 kernel_v = numpy.zeros((30,3)) #30
2464 2476 kernel_v[:, 0] = -1
2465 2477 kernel_v[:, 1] = 2
2466 2478 kernel_v[:, 2] = -1
2467 2479
2468 2480 kernel_4h = numpy.zeros((4,20)) #
2469 2481 kernel_4h[0, :] = 1
2470 2482 kernel_4h[1, :] = 0
2471 2483 kernel_4h[2, :] = 0
2472 2484 kernel_4h[3, :] = -1
2473 2485
2474 2486 kernel_5h = numpy.zeros((5,30)) #
2475 2487 kernel_5h[0, :] = 2
2476 2488 kernel_5h[1, :] = 1
2477 2489 kernel_5h[2, :] = 0
2478 2490 kernel_5h[3, :] = -1
2479 2491 kernel_5h[4, :] = -2
2480 2492
2481 2493
2482 2494 sharp_img0 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_3h)
2483 2495 # cv2.imshow('sharp image small h',sharp_img0) # numpy.fft.fftshift(sharp_img1))
2484 2496 # cv2.waitKey(0)
2485 2497
2486 2498 sharp_img1 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_h)
2487 2499 # cv2.imshow('sharp image h',sharp_img1) # numpy.fft.fftshift(sharp_img1))
2488 2500 # cv2.waitKey(0)
2489 2501
2490 2502 sharp_img2 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_v)
2491 2503 # cv2.imshow('sharp image v', sharp_img2) #numpy.fft.fftshift(sharp_img2))
2492 2504 # cv2.waitKey(0)
2493 2505
2494 2506 sharp_imgw = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_4h)
2495 2507 # cv2.imshow('sharp image h wide', sharp_imgw) #numpy.fft.fftshift(sharp_img2))
2496 2508 # cv2.waitKey(0)
2497 2509
2498 2510 sharp_imgwl = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_5h, borderType = cv2.BORDER_ISOLATED)
2499 2511 cv2.imshow('sharp image h long wide', sharp_imgwl) #numpy.fft.fftshift(sharp_img2))
2500 2512 cv2.waitKey(0)
2501 2513
2502 2514 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/sharp_h.jpg', sharp_img1)
2503 2515 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/sharp_v.jpg', sharp_img2)
2504 2516 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/input_img.jpg', img)
2505 2517
2506 2518 ########################small horizontal
2507 2519 horizontal_kernel_s = cv2.getStructuringElement(cv2.MORPH_RECT, (11,1)) #11
2508 2520 ######################## horizontal
2509 2521 horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (30,1)) #30
2510 2522 ######################## vertical
2511 2523 vertical_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (1,50)) #50
2512 2524 ######################## horizontal wide
2513 2525 horizontal_kernel_w = cv2.getStructuringElement(cv2.MORPH_RECT, (30,1)) # 30
2514 2526
2515 2527 horizontal_kernel_expand = cv2.getStructuringElement(cv2.MORPH_RECT, (3,3)) #
2516 2528
2517 2529 horizontal_kernel_wl = cv2.getStructuringElement(cv2.MORPH_RECT, (50,1)) #
2518 2530
2519 2531 detected_lines_h_s = cv2.morphologyEx(sharp_img0, cv2.MORPH_OPEN, horizontal_kernel_s, iterations=7) #7
2520 2532 detected_lines_h = cv2.morphologyEx(sharp_img1, cv2.MORPH_OPEN, horizontal_kernel, iterations=7) #7
2521 2533 detected_lines_v = cv2.morphologyEx(sharp_img2, cv2.MORPH_OPEN, vertical_kernel, iterations=7) #7
2522 2534 detected_lines_h_w = cv2.morphologyEx(sharp_imgw, cv2.MORPH_OPEN, horizontal_kernel_w, iterations=5) #5
2523 2535
2524 2536 detected_lines_h_wl = cv2.morphologyEx(sharp_imgwl, cv2.MORPH_OPEN, horizontal_kernel_wl, iterations=5) #
2525 2537 detected_lines_h_wl = cv2.filter2D(src=detected_lines_h_wl, ddepth=-1, kernel=horizontal_kernel_expand)
2526 2538
2527 2539 # cv2.imshow('lines h small gray', detected_lines_h_s) #numpy.fft.fftshift(detected_lines_h))
2528 2540 # cv2.waitKey(0)
2529 2541 # cv2.imshow('lines h gray', detected_lines_h) #numpy.fft.fftshift(detected_lines_h))
2530 2542 # cv2.waitKey(0)
2531 2543 # cv2.imshow('lines v gray', detected_lines_v) #numpy.fft.fftshift(detected_lines_h))
2532 2544 # cv2.waitKey(0)
2533 2545 # cv2.imshow('lines h wide gray', detected_lines_h_w) #numpy.fft.fftshift(detected_lines_h))
2534 2546 # cv2.waitKey(0)
2535 2547 cv2.imshow('lines h long wide gray', detected_lines_h_wl) #numpy.fft.fftshift(detected_lines_h))
2536 2548 cv2.waitKey(0)
2537 2549
2538 2550 reth_s, detected_lines_h_s = cv2.threshold(detected_lines_h_s, 85, 255, cv2.THRESH_BINARY)# 85
2539 2551 reth, detected_lines_h = cv2.threshold(detected_lines_h, 30, 255, cv2.THRESH_BINARY) #30
2540 2552 retv, detected_lines_v = cv2.threshold(detected_lines_v, 30, 255, cv2.THRESH_BINARY) #30
2541 2553 reth_w, detected_lines_h_w = cv2.threshold(detected_lines_h_w, 35, 255, cv2.THRESH_BINARY)#
2542 2554 reth_wl, detected_lines_h_wl = cv2.threshold(detected_lines_h_wl, 200, 255, cv2.THRESH_BINARY)#
2543 2555
2544 2556 cnts_h_s, h0 = cv2.findContours(detected_lines_h_s, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2545 2557 cnts_h, h1 = cv2.findContours(detected_lines_h, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2546 2558 cnts_v, h2 = cv2.findContours(detected_lines_v, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2547 2559 cnts_h_w, h3 = cv2.findContours(detected_lines_h_w, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2548 2560 cnts_h_wl, h4 = cv2.findContours(detected_lines_h_wl, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2549 2561 #print("horizontal ", cnts_h)
2550 2562 #print("vertical ", cnts_v)
2551 2563 # cnts, h = cv2.findContours(detected_lines_h, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
2552 2564 # print(cnts)
2553 2565 # cv2.imshow('lines h wide', detected_lines_h_w) #numpy.fft.fftshift(detected_lines_h))
2554 2566 # cv2.waitKey(0)
2555 2567 cv2.imshow('lines h wide long ', detected_lines_h_wl) #numpy.fft.fftshift(detected_lines_v))
2556 2568 cv2.waitKey(0)
2557 2569 # cv2.imshow('lines h small', detected_lines_h_s) #numpy.fft.fftshift(detected_lines_h))
2558 2570 # cv2.waitKey(0)
2559 2571 # cv2.imshow('lines h ', detected_lines_h) #numpy.fft.fftshift(detected_lines_h))
2560 2572 # cv2.waitKey(0)
2561 2573 # cv2.imshow('lines v ', detected_lines_v) #numpy.fft.fftshift(detected_lines_v))
2562 2574 # cv2.waitKey(0)
2563 2575
2564 2576 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/lines_h.jpg', detected_lines_h)
2565 2577 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/lines_v.jpg', detected_lines_v)
2566 2578
2567 2579 #cnts = cnts[0] if len(cnts) == 2 else cnts[1]
2568 2580 #y_line_index = numpy.asarray([ [c[0][0][0],c[1][0][0], c[0][0][1]] for c in cnts_v ])
2569 2581 h_line_index = []
2570 2582 v_line_index = []
2571 2583
2572 2584 cnts_h += cnts_h_s #combine large and small lines
2573 2585
2574 2586 # line indexes x1, x2, y
2575 2587 for c in cnts_h:
2576 2588 #print(c)
2577 2589 if len(c) < 3: #contorno linea
2578 2590 x1 = c[0][0][0]
2579 2591 if x1 < 8:
2580 2592 x1 = 10
2581 2593 h_line_index.append( [x1,c[1][0][0], c[0][0][1]] )
2582 2594 #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1])
2583 2595 else: #contorno poligono
2584 2596 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2585 2597 y = numpy.unique(pairs[:,1])
2586 2598 x = numpy.unique(pairs[:,0])
2587 2599 #print(x)
2588 2600 for yk in y:
2589 2601 x0 = x[0]
2590 2602 if x0 < 8:
2591 2603 x0 = 10
2592 2604 #print(x[0], x[-1], yk)
2593 2605 h_line_index.append( [x0, x[-1], yk])
2594 2606 #print("x1, x2, y ->p ", x[0], x[-1], yk)
2595 2607 for c in cnts_h_w:
2596 2608 #print(c)
2597 2609 if len(c) < 3: #contorno linea
2598 2610 x1 = c[0][0][0]
2599 2611 if x1 < 8:
2600 2612 x1 = 10
2601 2613 y = c[0][0][1] - 2 # se incrementa 2 líneas x el filtro
2602 2614 h_line_index.append( [x1,c[1][0][0],y] )
2603 2615 #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1])
2604 2616 else: #contorno poligono
2605 2617 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2606 2618 y = numpy.unique(pairs[:,1])
2607 2619 x = numpy.unique(pairs[:,0])
2608 2620 #print(x)
2609 2621 for yk in y:
2610 2622
2611 2623 x0 = x[0]
2612 2624 if x0 < 8:
2613 2625 x0 = 10
2614 2626 h_line_index.append( [x0, x[-1], yk-2])
2615 2627
2616 2628 for c in cnts_h_wl: # # revisar
2617 2629 #print(c)
2618 2630 if len(c) < 3: #contorno linea
2619 2631 x1 = c[0][0][0]
2620 2632 if x1 < 8:
2621 2633 x1 = 10
2622 2634 y = c[0][0][1] - 2 # se incrementa 2 líneas x el filtro
2623 2635 h_line_index.append( [x1,c[1][0][0],y] )
2624 2636 #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1])
2625 2637 else: #contorno poligono
2626 2638 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2627 2639 y = numpy.unique(pairs[:,1])
2628 2640 x = numpy.unique(pairs[:,0])
2629 2641 for yk in range(y[-1]-y[0]):
2630 2642 y_k = yk +y[0]
2631 2643
2632 2644 x0 = x[0]
2633 2645 if x0 < 8:
2634 2646 x0 = 10
2635 2647 h_line_index.append( [x0, x[-1], y_k-2])
2636 2648
2637 2649 print([[c[0][0][1],c[1][0][1], c[0][0][0] ] for c in cnts_v])
2638 2650 # line indexes y1, y2, x
2639 2651 for c in cnts_v:
2640 2652 if len(c) < 3: #contorno linea
2641 2653 v_line_index.append( [c[0][0][1],c[1][0][1], c[0][0][0] ] )
2642 2654 else: #contorno poligono
2643 2655 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2644 2656 #print(pairs)
2645 2657 y = numpy.unique(pairs[:,1])
2646 2658 x = numpy.unique(pairs[:,0])
2647 2659
2648 2660 for xk in x:
2649 2661 #print(x[0], x[-1], yk)
2650 2662 v_line_index.append( [y[0],y[-1], xk])
2651 2663
2652 2664 ###################################################################
2653 2665 # # clean Horizontal
2654 2666 print("Total Horizontal", len(h_line_index))
2655 2667 if len(h_line_index) < 75 :
2656 2668 for x1, x2, y in h_line_index:
2657 2669 #print("cleaning ",x1, x2, y)
2658 2670 len_line = x2 - x1
2659 2671 if y > 10 and y < (nh -10):
2660 2672 # if y != (nh-1):
2661 2673 # list = [ ((z[n, y-1] + z[n,y+1])/2) for n in range(len_line)]
2662 2674 # else:
2663 2675 # list = [ ((z[n, y-1] + z[n,0])/2) for n in range(len_line)]
2664 2676 #
2665 2677 # z[x1:x2,y] = numpy.asarray(list)
2666 2678 z[x1-5:x2+5,y:y+1] = 0
2667 2679
2668 2680 # clean vertical
2669 2681 for y1, y2, x in v_line_index:
2670 2682 len_line = y2 - y1
2671 2683 #print(x)
2672 2684 if x > 0 and x < (nsamples-2):
2673 2685 # if x != (nsamples-1):
2674 2686 # list = [ ((z[x-2, n] + z[x+2,n])/2) for n in range(len_line)]
2675 2687 # else:
2676 2688 # list = [ ((z[x-2, n] + z[1,n])/2) for n in range(len_line)]
2677 2689 #
2678 2690 # #z[x-1:x+1,y1:y2] = numpy.asarray(list)
2679 2691 #
2680 2692 z[x+1,y1:y2] = 0
2681 2693
2682 2694 '''
2683 2695 #z[: ,[215, 217, 221, 223, 225, 340, 342, 346, 348, 350, 465, 467, 471, 473, 475]]=0
2684 2696 z[1: ,[112, 114, 118, 120, 122, 237, 239, 245, 247, 249, 362, 364, 368, 370, 372]]=0
2685 2697 # z[: ,217]=0
2686 2698 # z[: ,221]=0
2687 2699 # z[: ,223]=0
2688 2700 # z[: ,225]=0
2689 2701
2690 2702 dat2 = numpy.log10(z.T)
2691 2703 #print(dat2)
2692 2704 max = dat2.max()
2693 2705 #print(" min, max ", max, min)
2694 2706 img2 = ((dat2-min)*255/(max-min)).astype(numpy.uint8)
2695 2707 #img_cleaned = img2.copy()
2696 2708 #cv2.drawContours(img2, cnts_h, -1, (255,255,255), 1)
2697 2709 #cv2.drawContours(img2, cnts_v, -1, (255,255,255), 1)
2698 2710
2699 2711
2700 2712 spcCleaned = z * numpy.exp(1j*phase)
2701 2713 #print(spcCleaned)
2702 2714
2703 2715
2704 2716 # cv2.imshow('image contours', img2) #numpy.fft.fftshift(img))
2705 2717 # cv2.waitKey(0)
2706 2718
2707 2719 cv2.imshow('cleaned', img2) #numpy.fft.fftshift(img_cleaned))
2708 2720 cv2.waitKey(0)
2709 2721 # # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/cleaned_{}.jpg'.format(self.test_counter), img2)
2710 2722 cv2.destroyAllWindows()
2711 2723 # self.test_counter += 1
2712 2724
2713 2725
2714 2726 #print("DC difference " ,z1 - z[0,:])
2715 2727
2716 2728 # m = numpy.mean(dat)
2717 2729 # o = numpy.std(dat)
2718 2730 # print("mean ", m, " std ", o)
2719 2731 # fig, ax = plt.subplots(1,2,figsize=(12, 6))
2720 2732 # #X, Y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
2721 2733 # X, Y = numpy.meshgrid(numpy.fft.fftshift(freqh),numpy.fft.fftshift(freqv))
2722 2734 #
2723 2735 # colormap = 'jet'
2724 2736 # #c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o))
2725 2737 # #c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = 6.5, vmax = 6.8)
2726 2738 # c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = (m-2*o), vmax = (m+1.5*o))
2727 2739 # fig.colorbar(c, ax=ax[0])
2728 2740 #
2729 2741 #
2730 2742 # #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o))
2731 2743 # #date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
2732 2744 # #strftime('%Y-%m-%d %H:%M:%S')
2733 2745 # #ax[0].set_title('Spectrum magnitude '+date_time)
2734 2746 # #fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
2735 2747 # #print("aqui estoy2",dat2[:,:,0].shape)
2736 2748 # #c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(pdat), cmap =colormap, vmin = 4.2, vmax = 5.0)
2737 2749 # c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat2), cmap =colormap, vmin = (m-2*o), vmax = (m+1.5*o))
2738 2750 # #c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(pdat), cmap =colormap ) #, vmin = 0.0, vmax = 0.5)
2739 2751 # #c = ax[1].pcolormesh(x, y, dat2[:,:,0], cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o)-1)
2740 2752 # #print("aqui estoy3")
2741 2753 # fig.colorbar(c, ax=ax[1])
2742 2754 # plt.show()
2743 2755
2744 2756 spectrum[ch,:,:] = spcCleaned
2745 2757
2746 2758 #print(data2.shape)
2747 2759
2748 2760
2749 2761
2750 2762 data[:,:,self.minHei_idx:] = numpy.fft.ifft2(spectrum, axes=(1,2))
2751 2763
2752 2764 #print("cleanOutliersByBlock Done", data.shape)
2753 2765 self.__buffer_data = data
2754 2766 return data
2755 2767
2756 2768
2757 2769
2758 2770 def fillBuffer(self, data, datatime):
2759 2771
2760 2772 if self.__profIndex == 0:
2761 2773 self.__buffer_data = data.copy()
2762 2774
2763 2775 else:
2764 2776 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
2765 2777 self.__profIndex += 1
2766 2778 self.__buffer_times.append(datatime)
2767 2779
2768 2780 def getData(self, data, datatime=None):
2769 2781
2770 2782 if self.__profIndex == 0:
2771 2783 self.__initime = datatime
2772 2784
2773 2785
2774 2786 self.__dataReady = False
2775 2787
2776 2788 self.fillBuffer(data, datatime)
2777 2789 dataBlock = None
2778 2790
2779 2791 if self.__profIndex == self.n:
2780 2792 #print("apnd : ",data)
2781 2793 dataBlock = self.cleanOutliersByBlock()
2782 2794 #dataBlock = self.cleanSpikesFFT2D()
2783 2795 #dataBlock = self.filterSatsProfiles2()
2784 2796 self.__dataReady = True
2785 2797
2786 2798 return dataBlock
2787 2799
2788 2800 if dataBlock is None:
2789 2801 return None, None
2790 2802
2791 2803
2792 2804
2793 2805 return dataBlock
2794 2806
2795 2807 def releaseBlock(self):
2796 2808
2797 2809 if self.n % self.lenProfileOut != 0:
2798 2810 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
2799 2811 return None
2800 2812
2801 2813 data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2802 2814
2803 2815 self.init_prof = self.end_prof
2804 2816 self.end_prof += self.lenProfileOut
2805 2817 #print("data release shape: ",dataOut.data.shape, self.end_prof)
2806 2818 self.n_prof_released += 1
2807 2819
2808 2820
2809 2821 #print("f_no_data ", dataOut.flagNoData)
2810 2822 return data
2811 2823
2812 2824 def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,th_hist_outlier=3,minHei=None, maxHei=None):
2813 2825 #print("run op buffer 2D",dataOut.ippSeconds)
2814 2826 # self.nChannels = dataOut.nChannels
2815 2827 # self.nHeights = dataOut.nHeights
2816 2828
2817 2829 if not self.isConfig:
2818 2830 #print("init p idx: ", dataOut.profileIndex )
2819 2831 self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,
2820 2832 thHistOutlier=th_hist_outlier,minHei=minHei, maxHei=maxHei)
2821 2833 self.isConfig = True
2822 2834
2823 2835 dataBlock = None
2824 2836
2825 2837 if not dataOut.buffer_empty: #hay datos acumulados
2826 2838
2827 2839 if self.init_prof == 0:
2828 2840 self.n_prof_released = 0
2829 2841 self.lenProfileOut = nProfilesOut
2830 2842 dataOut.flagNoData = False
2831 2843 #print("tp 2 ",dataOut.data.shape)
2832 2844
2833 2845 self.init_prof = 0
2834 2846 self.end_prof = self.lenProfileOut
2835 2847
2836 2848 dataOut.nProfiles = self.lenProfileOut
2837 2849 if nProfilesOut == 1:
2838 2850 dataOut.flagDataAsBlock = False
2839 2851 else:
2840 2852 dataOut.flagDataAsBlock = True
2841 2853 #print("prof: ",self.init_prof)
2842 2854 dataOut.flagNoData = False
2843 2855 if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
2844 2856 print("omitting: ", self.n_prof_released)
2845 2857 dataOut.flagNoData = True
2846 2858 dataOut.ippSeconds = self._ipp
2847 2859 dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
2848 2860 # print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
2849 2861 #dataOut.data = self.releaseBlock()
2850 2862 #########################################################3
2851 2863 if self.n % self.lenProfileOut != 0:
2852 2864 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
2853 2865 return None
2854 2866
2855 2867 dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2856 2868
2857 2869 self.init_prof = self.end_prof
2858 2870 self.end_prof += self.lenProfileOut
2859 2871 #print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)
2860 2872 self.n_prof_released += 1
2861 2873
2862 2874 if self.end_prof >= (self.n +self.lenProfileOut):
2863 2875
2864 2876 self.init_prof = 0
2865 2877 self.__profIndex = 0
2866 2878 self.buffer = None
2867 2879 dataOut.buffer_empty = True
2868 2880 self.outliers_IDs_list = []
2869 2881 self.n_prof_released = 0
2870 2882 dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
2871 2883 #print("cleaning...", dataOut.buffer_empty)
2872 2884 dataOut.profileIndex = 0 #self.lenProfileOut
2873 2885 ####################################################################
2874 2886 return dataOut
2875 2887
2876 2888
2877 2889 #print("tp 223 ",dataOut.data.shape)
2878 2890 dataOut.flagNoData = True
2879 2891
2880 2892
2881 2893
2882 2894 try:
2883 2895 #dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
2884 2896 dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
2885 2897 self.__count_exec +=1
2886 2898 except Exception as e:
2887 2899 print("Error getting profiles data",self.__count_exec )
2888 2900 print(e)
2889 2901 sys.exit()
2890 2902
2891 2903 if self.__dataReady:
2892 2904 #print("omitting: ", len(self.outliers_IDs_list))
2893 2905 self.__count_exec = 0
2894 2906 #dataOut.data =
2895 2907 #self.buffer = numpy.flip(dataBlock, axis=1)
2896 2908 self.buffer = dataBlock
2897 2909 self.first_utcBlock = self.__initime
2898 2910 dataOut.utctime = self.__initime
2899 2911 dataOut.nProfiles = self.__profIndex
2900 2912 #dataOut.flagNoData = False
2901 2913 self.init_prof = 0
2902 2914 self.__profIndex = 0
2903 2915 self.__initime = None
2904 2916 dataBlock = None
2905 2917 self.__buffer_times = []
2906 2918 dataOut.error = False
2907 2919 dataOut.useInputBuffer = True
2908 2920 dataOut.buffer_empty = False
2909 2921 #print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
2910 2922
2911 2923
2912 2924
2913 2925 #print(self.__count_exec)
2914 2926
2915 2927 return dataOut
2916 2928
2917 2929
2918 2930 class RemoveProfileSats(Operation):
2919 2931 '''
2920 2932 Escrito: Joab Apaza
2921 2933
2922 2934 Omite los perfiles contaminados con señal de satélites, usando una altura de referencia
2923 2935 In: minHei = min_sat_range
2924 2936 max_sat_range
2925 2937 min_hei_ref
2926 2938 max_hei_ref
2927 2939 th = diference between profiles mean, ref and sats
2928 2940 Out:
2929 2941 profile clean
2930 2942 '''
2931 2943
2932 2944
2933 2945 __buffer_data = []
2934 2946 __buffer_times = []
2935 2947
2936 2948 buffer = None
2937 2949
2938 2950 outliers_IDs_list = []
2939 2951
2940 2952
2941 2953 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
2942 2954 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
2943 2955 '__count_exec','__initime','__dataReady','__ipp', 'minRef', 'maxRef', 'thdB')
2944 2956 def __init__(self, **kwargs):
2945 2957
2946 2958 Operation.__init__(self, **kwargs)
2947 2959 self.isConfig = False
2948 2960
2949 2961 def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=15,
2950 2962 minHei=None, maxHei=None, minRef=None, maxRef=None, thdB=10):
2951 2963
2952 2964 if n == None and timeInterval == None:
2953 2965 raise ValueError("nprofiles or timeInterval should be specified ...")
2954 2966
2955 2967 if n != None:
2956 2968 self.n = n
2957 2969
2958 2970 self.navg = navg
2959 2971 self.profileMargin = profileMargin
2960 2972 self.thHistOutlier = thHistOutlier
2961 2973 self.__profIndex = 0
2962 2974 self.buffer = None
2963 2975 self._ipp = dataOut.ippSeconds
2964 2976 self.n_prof_released = 0
2965 2977 self.heightList = dataOut.heightList
2966 2978 self.init_prof = 0
2967 2979 self.end_prof = 0
2968 2980 self.__count_exec = 0
2969 2981 self.__profIndex = 0
2970 2982 self.first_utcBlock = None
2971 2983 #self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
2972 2984 minHei = minHei
2973 2985 maxHei = maxHei
2974 2986 if minHei==None :
2975 2987 minHei = dataOut.heightList[0]
2976 2988 if maxHei==None :
2977 2989 maxHei = dataOut.heightList[-1]
2978 2990 self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
2979 2991 self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)
2980 2992 self.nChannels = dataOut.nChannels
2981 2993 self.nHeights = dataOut.nHeights
2982 2994 self.test_counter = 0
2983 2995 self.thdB = thdB
2984 2996
2985 2997 def filterSatsProfiles(self):
2986 2998 data = self.__buffer_data
2987 2999 #print(data.shape)
2988 3000 nChannels, profiles, heights = data.shape
2989 3001 indexes=numpy.zeros([], dtype=int)
2990 3002 outliers_IDs=[]
2991 3003 for c in range(nChannels):
2992 3004 #print(self.min_ref,self.max_ref)
2993 3005 noise_ref = 10* numpy.log10((data[c,:,self.min_ref:self.max_ref] * numpy.conjugate(data[c,:,self.min_ref:self.max_ref])).real)
2994 3006 #print("Noise ",numpy.percentile(noise_ref,95))
2995 3007 p95 = numpy.percentile(noise_ref,95)
2996 3008 noise_ref = noise_ref.mean()
2997 3009 #print("Noise ",noise_ref
2998 3010
2999 3011
3000 3012 for h in range(self.minHei_idx, self.maxHei_idx):
3001 3013 power = 10* numpy.log10((data[c,:,h] * numpy.conjugate(data[c,:,h])).real)
3002 3014 #th = noise_ref + self.thdB
3003 3015 th = noise_ref + 1.5*(p95-noise_ref)
3004 3016 index = numpy.where(power > th )
3005 3017 if index[0].size > 10 and index[0].size < int(self.navg*profiles):
3006 3018 indexes = numpy.append(indexes, index[0])
3007 3019 #print(index[0])
3008 3020 #print(index[0])
3009 3021
3010 3022 # fig,ax = plt.subplots()
3011 3023 # #ax.set_title(str(k)+" "+str(j))
3012 3024 # x=range(len(power))
3013 3025 # ax.scatter(x,power)
3014 3026 # #ax.axvline(index)
3015 3027 # plt.grid()
3016 3028 # plt.show()
3017 3029 #print(indexes)
3018 3030
3019 3031 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
3020 3032 #outliers_IDs = numpy.unique(outliers_IDs)
3021 3033
3022 3034 outs_lines = numpy.unique(indexes)
3023 3035
3024 3036
3025 3037 #Agrupando el histograma de outliers,
3026 3038 my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=True)
3027 3039
3028 3040
3029 3041 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
3030 3042 hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
3031 3043 hist_outliers_indexes = hist_outliers_indexes[0]
3032 3044 # if len(hist_outliers_indexes>0):
3033 3045 # hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1)
3034 3046 #print(hist_outliers_indexes)
3035 3047 #print(bins, hist_outliers_indexes)
3036 3048 bins_outliers_indexes = [int(i) for i in (bins[hist_outliers_indexes])] #
3037 3049 outlier_loc_index = []
3038 3050 # for n in range(len(bins_outliers_indexes)):
3039 3051 # for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n]+ self.profileMargin):
3040 3052 # outlier_loc_index.append(e)
3041 3053 outlier_loc_index = [e for n in range(len(bins_outliers_indexes)) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n]+ profiles//100 + self.profileMargin) ]
3042 3054 outlier_loc_index = numpy.asarray(outlier_loc_index)
3043 3055
3044 3056
3045 3057
3046 3058
3047 3059 #print("outliers Ids: ", outlier_loc_index, outlier_loc_index.shape)
3048 3060 outlier_loc_index = outlier_loc_index[ (outlier_loc_index >= 0) & (outlier_loc_index<profiles)]
3049 3061 #print("outliers final: ", outlier_loc_index)
3050 3062
3051 3063 from matplotlib import pyplot as plt
3052 3064 x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
3053 3065 fig, ax = plt.subplots(1,2,figsize=(8, 6))
3054 3066 dat = data[0,:,:].real
3055 3067 dat = 10* numpy.log10((data[0,:,:] * numpy.conjugate(data[0,:,:])).real)
3056 3068 m = numpy.nanmean(dat)
3057 3069 o = numpy.nanstd(dat)
3058 3070 #print(m, o, x.shape, y.shape)
3059 3071 #c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
3060 3072 c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = 50, vmax = 75)
3061 3073 ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
3062 3074 fig.colorbar(c)
3063 3075 ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
3064 3076 ax[1].hist(outs_lines,bins=my_bins)
3065 3077 plt.show()
3066 3078
3067 3079
3068 3080 self.outliers_IDs_list = outlier_loc_index
3069 3081 #print("outs list: ", self.outliers_IDs_list)
3070 3082 return data
3071 3083
3072 3084
3073 3085
3074 3086 def fillBuffer(self, data, datatime):
3075 3087
3076 3088 if self.__profIndex == 0:
3077 3089 self.__buffer_data = data.copy()
3078 3090
3079 3091 else:
3080 3092 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
3081 3093 self.__profIndex += 1
3082 3094 self.__buffer_times.append(datatime)
3083 3095
3084 3096 def getData(self, data, datatime=None):
3085 3097
3086 3098 if self.__profIndex == 0:
3087 3099 self.__initime = datatime
3088 3100
3089 3101
3090 3102 self.__dataReady = False
3091 3103
3092 3104 self.fillBuffer(data, datatime)
3093 3105 dataBlock = None
3094 3106
3095 3107 if self.__profIndex == self.n:
3096 3108 #print("apnd : ",data)
3097 3109 dataBlock = self.filterSatsProfiles()
3098 3110 self.__dataReady = True
3099 3111
3100 3112 return dataBlock
3101 3113
3102 3114 if dataBlock is None:
3103 3115 return None, None
3104 3116
3105 3117
3106 3118
3107 3119 return dataBlock
3108 3120
3109 3121 def releaseBlock(self):
3110 3122
3111 3123 if self.n % self.lenProfileOut != 0:
3112 3124 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
3113 3125 return None
3114 3126
3115 3127 data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
3116 3128
3117 3129 self.init_prof = self.end_prof
3118 3130 self.end_prof += self.lenProfileOut
3119 3131 #print("data release shape: ",dataOut.data.shape, self.end_prof)
3120 3132 self.n_prof_released += 1
3121 3133
3122 3134 return data
3123 3135
3124 3136 def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,
3125 3137 th_hist_outlier=15,minHei=None, maxHei=None, minRef=None, maxRef=None, thdB=10):
3126 3138
3127 3139 if not self.isConfig:
3128 3140 #print("init p idx: ", dataOut.profileIndex )
3129 3141 self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,thHistOutlier=th_hist_outlier,
3130 3142 minHei=minHei, maxHei=maxHei, minRef=minRef, maxRef=maxRef, thdB=thdB)
3131 3143 self.isConfig = True
3132 3144
3133 3145 dataBlock = None
3134 3146
3135 3147 if not dataOut.buffer_empty: #hay datos acumulados
3136 3148
3137 3149 if self.init_prof == 0:
3138 3150 self.n_prof_released = 0
3139 3151 self.lenProfileOut = nProfilesOut
3140 3152 dataOut.flagNoData = False
3141 3153 #print("tp 2 ",dataOut.data.shape)
3142 3154
3143 3155 self.init_prof = 0
3144 3156 self.end_prof = self.lenProfileOut
3145 3157
3146 3158 dataOut.nProfiles = self.lenProfileOut
3147 3159 if nProfilesOut == 1:
3148 3160 dataOut.flagDataAsBlock = False
3149 3161 else:
3150 3162 dataOut.flagDataAsBlock = True
3151 3163 #print("prof: ",self.init_prof)
3152 3164 dataOut.flagNoData = False
3153 3165 if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
3154 3166 #print("omitting: ", self.n_prof_released)
3155 3167 dataOut.flagNoData = True
3156 3168 dataOut.ippSeconds = self._ipp
3157 3169 dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
3158 3170 # print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
3159 3171 #dataOut.data = self.releaseBlock()
3160 3172 #########################################################3
3161 3173 if self.n % self.lenProfileOut != 0:
3162 3174 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
3163 3175 return None
3164 3176
3165 3177 dataOut.data = None
3166 3178
3167 3179 if nProfilesOut == 1:
3168 3180 dataOut.data = self.buffer[:,self.end_prof-1,:] #ch, prof, alt
3169 3181 else:
3170 3182 dataOut.data = self.buffer[:,self.init_prof:self.end_prof,:] #ch, prof, alt
3171 3183
3172 3184 self.init_prof = self.end_prof
3173 3185 self.end_prof += self.lenProfileOut
3174 3186 #print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)
3175 3187 self.n_prof_released += 1
3176 3188
3177 3189 if self.end_prof >= (self.n +self.lenProfileOut):
3178 3190
3179 3191 self.init_prof = 0
3180 3192 self.__profIndex = 0
3181 3193 self.buffer = None
3182 3194 dataOut.buffer_empty = True
3183 3195 self.outliers_IDs_list = []
3184 3196 self.n_prof_released = 0
3185 3197 dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
3186 3198 #print("cleaning...", dataOut.buffer_empty)
3187 3199 dataOut.profileIndex = 0 #self.lenProfileOut
3188 3200 ####################################################################
3189 3201 return dataOut
3190 3202
3191 3203
3192 3204 #print("tp 223 ",dataOut.data.shape)
3193 3205 dataOut.flagNoData = True
3194 3206
3195 3207
3196 3208
3197 3209 try:
3198 3210 #dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
3199 3211 dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
3200 3212 self.__count_exec +=1
3201 3213 except Exception as e:
3202 3214 print("Error getting profiles data",self.__count_exec )
3203 3215 print(e)
3204 3216 sys.exit()
3205 3217
3206 3218 if self.__dataReady:
3207 3219 #print("omitting: ", len(self.outliers_IDs_list))
3208 3220 self.__count_exec = 0
3209 3221 #dataOut.data =
3210 3222 #self.buffer = numpy.flip(dataBlock, axis=1)
3211 3223 self.buffer = dataBlock
3212 3224 self.first_utcBlock = self.__initime
3213 3225 dataOut.utctime = self.__initime
3214 3226 dataOut.nProfiles = self.__profIndex
3215 3227 #dataOut.flagNoData = False
3216 3228 self.init_prof = 0
3217 3229 self.__profIndex = 0
3218 3230 self.__initime = None
3219 3231 dataBlock = None
3220 3232 self.__buffer_times = []
3221 3233 dataOut.error = False
3222 3234 dataOut.useInputBuffer = True
3223 3235 dataOut.buffer_empty = False
3224 3236 #print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
3225 3237
3226 3238
3227 3239
3228 3240 #print(self.__count_exec)
3229 3241
3230 3242 return dataOut
3231 3243
3232 3244
3233 3245 class RemoveProfileSats2(Operation):
3234 3246 '''
3235 3247 Escrito: Joab Apaza
3236 3248
3237 3249 Omite los perfiles contaminados con señal de satélites, usando una altura de referencia
3238 3250 promedia todas las alturas para los cálculos
3239 3251 In:
3240 3252 n = Cantidad de perfiles que se acumularan, usualmente 10 segundos
3241 3253 navg = Porcentaje de perfiles que puede considerarse como satélite, máximo 90%
3242 3254 minHei =
3243 3255 minRef =
3244 3256 maxRef =
3245 3257 nBins =
3246 3258 profile_margin =
3247 3259 th_hist_outlier =
3248 3260 nProfilesOut =
3249 3261
3250 3262 Pensado para remover interferencias de las YAGI, se puede adaptar a otras interferencias
3251 3263
3252 3264 remYagi = Activa la funcion de remoción de interferencias de la YAGI
3253 3265 nProfYagi = Cantidad de perfiles que son afectados, acorde NTX de la YAGI
3254 3266 offYagi =
3255 3267 minHJULIA = Altura mínima donde aparece la señal referencia de JULIA (-50)
3256 3268 maxHJULIA = Altura máxima donde aparece la señal referencia de JULIA (-15)
3257 3269
3258 3270 debug = Activa los gráficos, recomendable ejecutar para ajustar los parámetros
3259 3271 para un experimento en específico.
3260 3272
3261 3273 ** se modifica para remover interferencias puntuales, es decir, desde otros radares.
3262 3274 Inicialmente se ha configurado para omitir también los perfiles de la YAGI en los datos
3263 3275 de AMISR-ISR.
3264 3276
3265 3277 Out:
3266 3278 profile clean
3267 3279 '''
3268 3280
3269 3281
3270 3282 __buffer_data = []
3271 3283 __buffer_times = []
3272 3284
3273 3285 buffer = None
3274 3286
3275 3287 outliers_IDs_list = []
3276 3288
3277 3289
3278 3290 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
3279 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
3291 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels','cohFactor',
3280 3292 '__count_exec','__initime','__dataReady','__ipp', 'minRef', 'maxRef', 'debug','prev_pnoise')
3281 3293 def __init__(self, **kwargs):
3282 3294
3283 3295 Operation.__init__(self, **kwargs)
3284 3296 self.isConfig = False
3285 3297 self.currentTime = None
3286 3298
3287 def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=15,minHei=None, maxHei=None, nBins=10,
3299 def setup(self,dataOut, n=None , navg=0.9, profileMargin=50,thHistOutlier=15,minHei=None, maxHei=None, nBins=10,
3288 3300 minRef=None, maxRef=None, debug=False, remYagi=False, nProfYagi = 0, offYagi=0, minHJULIA=None, maxHJULIA=None,
3289 3301 idate=None,startH=None,endH=None ):
3290 3302
3291 3303 if n == None and timeInterval == None:
3292 3304 raise ValueError("nprofiles or timeInterval should be specified ...")
3293 3305
3294 3306 if n != None:
3295 3307 self.n = n
3296 3308
3297 3309 self.navg = navg
3298 3310 self.profileMargin = profileMargin
3299 3311 self.thHistOutlier = thHistOutlier
3300 3312 self.__profIndex = 0
3301 3313 self.buffer = None
3302 3314 self._ipp = dataOut.ippSeconds
3303 3315 self.n_prof_released = 0
3304 3316 self.heightList = dataOut.heightList
3305 3317 self.init_prof = 0
3306 3318 self.end_prof = 0
3307 3319 self.__count_exec = 0
3308 3320 self.__profIndex = 0
3309 3321 self.first_utcBlock = None
3310 3322 self.prev_pnoise = None
3311 3323 self.nBins = nBins
3312 3324 #self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
3313 3325 minHei = minHei
3314 3326 maxHei = maxHei
3315 3327 if minHei==None :
3316 3328 minHei = dataOut.heightList[0]
3317 3329 if maxHei==None :
3318 3330 maxHei = dataOut.heightList[-1]
3319 3331 self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
3320 3332 self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)
3321 3333 self.nChannels = dataOut.nChannels
3322 3334 self.nHeights = dataOut.nHeights
3323 3335 self.test_counter = 0
3324 3336 self.debug = debug
3325 3337 self.remYagi = remYagi
3326
3338 self.cohFactor = dataOut.nCohInt
3327 3339 if self.remYagi :
3328 3340 if minHJULIA==None or maxHJULIA==None:
3329 3341 raise ValueError("Parameters minHYagi and minHYagi are necessary!")
3330 3342 return
3331 3343 if idate==None or startH==None or endH==None:
3332 3344 raise ValueError("Date and hour parameters are necessary!")
3333 3345 return
3334 3346 self.minHJULIA_idx,self.maxHJULIA_idx = getHei_index(minHJULIA, maxHJULIA, dataOut.heightList)
3335 3347 self.offYagi = offYagi
3336 3348 self.nTxYagi = nProfYagi
3337 3349
3338 3350 self.startTime = datetime.datetime.combine(idate,startH)
3339 3351 self.endTime = datetime.datetime.combine(idate,endH)
3340 3352
3341 log.warning("Be careful with the selection of parameters for sat removal! Is recommended to \
3353 log.warning("Be careful with the selection of parameters for sats removal! It is avisable to \
3342 3354 activate the debug parameter in this operation for calibration", self.name)
3343 3355
3344 3356
3345 3357 def filterSatsProfiles(self):
3346 3358
3347 3359 data = self.__buffer_data.copy()
3348 3360 #print(data.shape)
3349 3361 nChannels, profiles, heights = data.shape
3350 3362 indexes=numpy.zeros([], dtype=int)
3351 3363 indexes = numpy.delete(indexes,0)
3352 3364
3353 3365 indexesYagi=numpy.zeros([], dtype=int)
3354 3366 indexesYagi = numpy.delete(indexesYagi,0)
3355 3367
3356 3368 indexesYagi_up=numpy.zeros([], dtype=int)
3357 3369 indexesYagi_up = numpy.delete(indexesYagi_up,0)
3358 3370 indexesYagi_down=numpy.zeros([], dtype=int)
3359 3371 indexesYagi_down = numpy.delete(indexesYagi_down,0)
3360 3372
3361 3373
3362 3374 indexesJULIA=numpy.zeros([], dtype=int)
3363 3375 indexesJULIA = numpy.delete(indexesJULIA,0)
3364 3376
3365 3377 outliers_IDs=[]
3366 3378
3367 3379 div = profiles//self.nBins
3368 3380
3369 3381 for c in range(nChannels):
3370 3382 #print(self.min_ref,self.max_ref)
3371 3383
3372 3384 import scipy.signal
3373 b, a = scipy.signal.butter(3, 0.1)
3374 noise_ref = (data[c,:,self.min_ref:self.max_ref] * numpy.conjugate(data[c,:,self.min_ref:self.max_ref])).real
3385 b, a = scipy.signal.butter(3, 0.5)
3386 #noise_ref = (data[c,:,self.min_ref:self.max_ref] * numpy.conjugate(data[c,:,self.min_ref:self.max_ref]))
3387 noise_ref = numpy.abs(data[c,:,self.min_ref:self.max_ref])
3388 lnoise = len(noise_ref[0,:])
3389 #print(noise_ref.shape)
3375 3390 noise_ref = noise_ref.mean(axis=1)
3391 #fnoise = noise_ref
3376 3392 fnoise = scipy.signal.filtfilt(b, a, noise_ref)
3377 3393 #noise_refdB = 10* numpy.log10(noise_ref)
3378 3394 #print("Noise ",numpy.percentile(noise_ref,95))
3379 p85 = numpy.percentile(fnoise,83)
3395 p95 = numpy.percentile(fnoise,90)
3380 3396 mean_noise = fnoise.mean()
3397
3381 3398 if self.prev_pnoise != None:
3382 3399 if mean_noise < (1.5 * self.prev_pnoise) :
3383 3400 self.prev_pnoise = mean_noise
3384 3401 else:
3385 3402 mean_noise = self.prev_pnoise
3386 3403 else:
3387 3404 self.prev_pnoise = mean_noise
3388 3405
3389 3406 std = fnoise.std()+ fnoise.mean()
3390 3407
3391 3408
3392 3409
3393 power = (data[c,:,self.minHei_idx:self.maxHei_idx] * numpy.conjugate(data[c,:,self.minHei_idx:self.maxHei_idx])).real
3394 heis = len(power[0,:])
3395 power = power.mean(axis=1)
3410 #power = (data[c,:,self.minHei_idx:self.maxHei_idx] * numpy.conjugate(data[c,:,self.minHei_idx:self.maxHei_idx]))
3411 power = numpy.abs(data[c,:,self.minHei_idx:self.maxHei_idx])
3412 npower = len(power[0,:])
3413 #print(power.shape)
3414 power = power.mean(axis=1)
3396 3415
3397 3416 fpower = scipy.signal.filtfilt(b, a, power)
3398 3417 #print(power.shape)
3399 3418 #powerdB = 10* numpy.log10(power)
3400 3419
3401 th = p85
3420 th = p95 #* 1.1
3421 #th = mean_noise
3402 3422 index = numpy.where(fpower > th )
3403 #print("Noise ",mean_noise, p85)
3423 #print("Noise ",mean_noise, p95)
3404 3424 #print(index)
3405 3425
3406 if index[0].size < int(self.navg*profiles):
3426
3427 if index[0].size <= int(self.navg*profiles):
3407 3428 indexes = numpy.append(indexes, index[0])
3408 3429
3409 3430 #print("sdas ", noise_ref.mean())
3410 3431
3411 3432 if self.remYagi :
3412 3433 #print(self.minHJULIA_idx, self.maxHJULIA_idx)
3413 3434 powerJULIA = (data[c,:,self.minHJULIA_idx:self.maxHJULIA_idx] * numpy.conjugate(data[c,:,self.minHJULIA_idx:self.maxHJULIA_idx])).real
3414 3435 powerJULIA = powerJULIA.mean(axis=1)
3415 3436 th_JULIA = powerJULIA.mean()*0.85
3416 3437 indexJULIA = numpy.where(powerJULIA >= th_JULIA )
3417 3438
3418 3439 indexesJULIA= numpy.append(indexesJULIA, indexJULIA[0])
3419 3440
3420 3441 # fig, ax = plt.subplots()
3421 3442 # ax.plot(powerJULIA)
3422 3443 # ax.axhline(th_JULIA, color='r')
3423 3444 # plt.grid()
3424 3445 # plt.show()
3425 3446
3426 # fig, ax = plt.subplots()
3427 # ax.plot(fpower)
3428 # ax.axhline(th, color='r')
3429 # ax.axhline(std, color='b')
3430 # plt.grid()
3431 # plt.show()
3447 if self.debug:
3448 fig, ax = plt.subplots()
3449 ax.plot(fpower, label="power")
3450 #ax.plot(fnoise, label="noise ref")
3451 ax.axhline(th, color='g', label="th")
3452 #ax.axhline(std, color='b', label="mean")
3453 ax.legend()
3454 plt.grid()
3455 plt.show()
3432 3456
3433 3457 #print(indexes)
3434 3458
3435 3459 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
3436 3460 #outliers_IDs = numpy.unique(outliers_IDs)
3437 3461 # print(indexesJULIA)
3438 3462 if len(indexesJULIA > 1):
3439 3463 iJ = indexesJULIA
3440 3464 locs = [ (iJ[n]-iJ[n-1]) > 5 for n in range(len(iJ))]
3441 3465 locs_2 = numpy.where(locs)[0]
3442 3466 #print(locs_2, indexesJULIA[locs_2-1])
3443 3467 indexesYagi_up = numpy.append(indexesYagi_up, indexesJULIA[locs_2-1])
3444 3468 indexesYagi_down = numpy.append(indexesYagi_down, indexesJULIA[locs_2])
3445 3469
3446 3470
3447 3471 indexesYagi_up = numpy.append(indexesYagi_up,indexesJULIA[-1])
3448 3472 indexesYagi_down = numpy.append(indexesYagi_down,indexesJULIA[0])
3449 3473
3450 3474 indexesYagi_up = numpy.unique(indexesYagi_up)
3451 3475 indexesYagi_down = numpy.unique(indexesYagi_down)
3452 3476
3453 3477
3454 3478 aux_ind = [ numpy.arange( (self.offYagi + k)+1, (self.offYagi + k + self.nTxYagi)+1, 1, dtype=int) for k in indexesYagi_up]
3455 3479 indexesYagi_up = (numpy.asarray(aux_ind)).flatten()
3456 3480
3457 3481 aux_ind2 = [ numpy.arange( (k - self.nTxYagi)+1, k+1 , 1, dtype=int) for k in indexesYagi_down]
3458 3482 indexesYagi_down = (numpy.asarray(aux_ind2)).flatten()
3459 3483
3460 3484 indexesYagi = numpy.append(indexesYagi,indexesYagi_up)
3461 3485 indexesYagi = numpy.append(indexesYagi,indexesYagi_down)
3462 3486
3463 3487
3464 3488 indexesYagi = indexesYagi[ (indexesYagi >= 0) & (indexesYagi<profiles)]
3465 3489 indexesYagi = numpy.unique(indexesYagi)
3466 3490
3491 #print("indexes: " ,indexes)
3467 3492 outs_lines = numpy.unique(indexes)
3468
3493 #print(outs_lines)
3469 3494
3470 3495 #Agrupando el histograma de outliers,
3471 3496 my_bins = numpy.linspace(0,int(profiles), div, endpoint=True)
3472 3497
3473 3498
3474 3499 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
3475 hist_outliers_indexes = numpy.where(hist >= self.thHistOutlier) #es outlier
3476 hist_outliers_indexes = hist_outliers_indexes[0]
3500 #print("hist: ",hist)
3501 hist_outliers_indexes = numpy.where(hist >= self.thHistOutlier)[0] #es outlier
3502 #print(hist_outliers_indexes)
3477 3503 if len(hist_outliers_indexes>0):
3478 3504 hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1)
3479 3505
3480 3506 bins_outliers_indexes = [int(i) for i in (bins[hist_outliers_indexes])] #
3481 3507 outlier_loc_index = []
3482 3508 #print("out indexes ", bins_outliers_indexes)
3483 3509 if len(bins_outliers_indexes) <= 3:
3484 3510 extprof = 0
3485 3511 else:
3486 3512 extprof = self.profileMargin
3487 3513 outlier_loc_index = [e for n in range(len(bins_outliers_indexes)) for e in range(bins_outliers_indexes[n]-extprof,bins_outliers_indexes[n] + extprof) ]
3488 3514 outlier_loc_index = numpy.asarray(outlier_loc_index)
3489 3515 # if len(outlier_loc_index)>1:
3490 3516 # ipmax = numpy.where(fpower==fpower.max())[0]
3491 3517 # print("pmax: ",ipmax)
3492 3518
3493 3519
3494 3520
3495 3521
3496 3522 #print("outliers Ids: ", outlier_loc_index, outlier_loc_index.shape)
3497 3523 outlier_loc_index = outlier_loc_index[ (outlier_loc_index >= 0) & (outlier_loc_index<profiles)]
3498 3524 #print("outliers final: ", outlier_loc_index)
3499 3525
3500 3526
3501 3527 if self.debug:
3502 3528 x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
3503 3529 fig, ax = plt.subplots(nChannels,2,figsize=(8, 6))
3504 3530
3505 3531 for i in range(nChannels):
3506 3532 dat = data[i,:,:].real
3507 3533 dat = 10* numpy.log10((data[i,:,:] * numpy.conjugate(data[i,:,:])).real)
3508 3534 m = numpy.nanmean(dat)
3509 3535 o = numpy.nanstd(dat)
3510 3536 if nChannels>1:
3511 3537 c = ax[i][0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70)
3512 3538 ax[i][0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
3513 3539 #fig.colorbar(c)
3514 3540 ax[i][0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
3515 3541 ax[i][1].hist(outs_lines,bins=my_bins)
3516 3542 if self.remYagi :
3517 3543 ax[0].vlines(indexesYagi,750,850, linestyles='dashed', label = 'yagi', color='m')
3518 3544 else:
3519 c = ax[0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70)
3545 c = ax[0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = (70+2*self.cohFactor))
3520 3546 ax[0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
3521 3547 #fig.colorbar(c)
3522 3548 ax[0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
3523 3549
3524 3550 ax[1].hist(outs_lines,bins=my_bins)
3525 3551 if self.remYagi :
3526 3552 ax[0].vlines(indexesYagi,750,850, linestyles='dashed', label = 'yagi', color='m')
3527 3553 plt.show()
3528 3554
3529 3555
3530 3556
3531 3557
3532 3558 if self.remYagi and (self.currentTime < self.startTime and self.currentTime < self.endTime):
3533 3559 outlier_loc_index = numpy.append(outlier_loc_index,indexesYagi)
3534 3560
3535 3561 self.outliers_IDs_list = numpy.unique(outlier_loc_index)
3536 3562
3537 3563 #print("outs list: ", self.outliers_IDs_list)
3538 3564 return self.__buffer_data
3539 3565
3540 3566
3541 3567
3542 3568 def fillBuffer(self, data, datatime):
3543 3569
3544 3570 if self.__profIndex == 0:
3545 3571 self.__buffer_data = data.copy()
3546 3572
3547 3573 else:
3548 3574 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
3549 3575 self.__profIndex += 1
3550 3576 self.__buffer_times.append(datatime)
3551 3577
3552 3578 def getData(self, data, datatime=None):
3553 3579
3554 3580 if self.__profIndex == 0:
3555 3581 self.__initime = datatime
3556 3582
3557 3583
3558 3584 self.__dataReady = False
3559 3585
3560 3586 self.fillBuffer(data, datatime)
3561 3587 dataBlock = None
3562 3588
3563 3589 if self.__profIndex == self.n:
3564 3590 #print("apnd : ",data)
3565 3591 dataBlock = self.filterSatsProfiles()
3566 3592 self.__dataReady = True
3567 3593
3568 3594 return dataBlock
3569 3595
3570 3596 if dataBlock is None:
3571 3597 return None, None
3572 3598
3573 3599
3574 3600
3575 3601 return dataBlock
3576 3602
3577 3603 def releaseBlock(self):
3578 3604
3579 3605 if self.n % self.lenProfileOut != 0:
3580 3606 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
3581 3607 return None
3582 3608
3583 3609 data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
3584 3610
3585 3611 self.init_prof = self.end_prof
3586 3612 self.end_prof += self.lenProfileOut
3587 3613 #print("data release shape: ",dataOut.data.shape, self.end_prof)
3588 3614 self.n_prof_released += 1
3589 3615
3590 3616 return data
3591 3617
3592 3618 def run(self, dataOut, n=None, navg=0.9, nProfilesOut=1, profile_margin=50, th_hist_outlier=15,minHei=None,nBins=10,
3593 3619 maxHei=None, minRef=None, maxRef=None, debug=False, remYagi=False, nProfYagi = 0, offYagi=0, minHJULIA=None, maxHJULIA=None,
3594 3620 idate=None,startH=None,endH=None):
3595 3621
3596 3622 if not self.isConfig:
3597 3623 #print("init p idx: ", dataOut.profileIndex )
3598 3624 self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,thHistOutlier=th_hist_outlier,minHei=minHei,
3599 3625 nBins=10, maxHei=maxHei, minRef=minRef, maxRef=maxRef, debug=debug, remYagi=remYagi, nProfYagi = nProfYagi,
3600 3626 offYagi=offYagi, minHJULIA=minHJULIA,maxHJULIA=maxHJULIA,idate=idate,startH=startH,endH=endH)
3601 3627
3602 3628 self.isConfig = True
3603 3629
3604 3630 dataBlock = None
3605 3631 self.currentTime = datetime.datetime.fromtimestamp(dataOut.utctime)
3606 3632
3607 3633 if not dataOut.buffer_empty: #hay datos acumulados
3608 3634
3609 3635 if self.init_prof == 0:
3610 3636 self.n_prof_released = 0
3611 3637 self.lenProfileOut = nProfilesOut
3612 3638 dataOut.flagNoData = False
3613 3639 #print("tp 2 ",dataOut.data.shape)
3614 3640
3615 3641 self.init_prof = 0
3616 3642 self.end_prof = self.lenProfileOut
3617 3643
3618 3644 dataOut.nProfiles = self.lenProfileOut
3619 3645 if nProfilesOut == 1:
3620 3646 dataOut.flagDataAsBlock = False
3621 3647 else:
3622 3648 dataOut.flagDataAsBlock = True
3623 3649 #print("prof: ",self.init_prof)
3624 3650 dataOut.flagNoData = False
3625 3651 if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
3626 3652 #print("omitting: ", self.n_prof_released)
3627 3653 dataOut.flagNoData = True
3628 3654 dataOut.ippSeconds = self._ipp
3629 3655 dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
3630 3656 # print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
3631 3657 #dataOut.data = self.releaseBlock()
3632 3658 #########################################################3
3633 3659 if self.n % self.lenProfileOut != 0:
3634 3660 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
3635 3661 return None
3636 3662
3637 3663 dataOut.data = None
3638 3664
3639 3665 if nProfilesOut == 1:
3640 3666 dataOut.data = self.buffer[:,self.end_prof-1,:] #ch, prof, alt
3641 3667 else:
3642 3668 dataOut.data = self.buffer[:,self.init_prof:self.end_prof,:] #ch, prof, alt
3643 3669
3644 3670 self.init_prof = self.end_prof
3645 3671 self.end_prof += self.lenProfileOut
3646 3672 #print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)
3647 3673 self.n_prof_released += 1
3648 3674
3649 3675 if self.end_prof >= (self.n +self.lenProfileOut):
3650 3676
3651 3677 self.init_prof = 0
3652 3678 self.__profIndex = 0
3653 3679 self.buffer = None
3654 3680 dataOut.buffer_empty = True
3655 3681 self.outliers_IDs_list = []
3656 3682 self.n_prof_released = 0
3657 3683 dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
3658 3684 #print("cleaning...", dataOut.buffer_empty)
3659 3685 dataOut.profileIndex = self.__profIndex
3660 3686 ####################################################################
3661 3687 return dataOut
3662 3688
3663 3689
3664 3690 #print("tp 223 ",dataOut.data.shape)
3665 3691 dataOut.flagNoData = True
3666 3692
3667 3693
3668 3694
3669 3695 try:
3670 3696 #dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
3671 3697 dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
3672 3698 self.__count_exec +=1
3673 3699 except Exception as e:
3674 3700 print("Error getting profiles data",self.__count_exec )
3675 3701 print(e)
3676 3702 sys.exit()
3677 3703
3678 3704 if self.__dataReady:
3679 3705 #print("omitting: ", len(self.outliers_IDs_list))
3680 3706 self.__count_exec = 0
3681 3707 #dataOut.data =
3682 3708 #self.buffer = numpy.flip(dataBlock, axis=1)
3683 3709 self.buffer = dataBlock
3684 3710 self.first_utcBlock = self.__initime
3685 3711 dataOut.utctime = self.__initime
3686 3712 dataOut.nProfiles = self.__profIndex
3687 3713 #dataOut.flagNoData = False
3688 3714 self.init_prof = 0
3689 3715 self.__profIndex = 0
3690 3716 self.__initime = None
3691 3717 dataBlock = None
3692 3718 self.__buffer_times = []
3693 3719 dataOut.error = False
3694 3720 dataOut.useInputBuffer = True
3695 3721 dataOut.buffer_empty = False
3696 3722 #print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
3697 3723
3698 3724
3699 3725
3700 3726 #print(self.__count_exec)
3701 3727
3702 3728 return dataOut
3703 3729
3704 3730
3705 3731
3706 3732
3707 3733 class remHeightsIppInterf(Operation):
3708 3734
3709 3735 def __init__(self, **kwargs):
3710 3736
3711 3737
3712 3738 Operation.__init__(self, **kwargs)
3713 3739
3714 3740 self.isConfig = False
3715 3741
3716 3742 self.heights_indx = None
3717 3743 self.heightsList = []
3718 3744
3719 3745 self.ipp1 = None
3720 3746 self.ipp2 = None
3721 3747 self.tx1 = None
3722 3748 self.tx2 = None
3723 3749 self.dh1 = None
3724 3750
3725 3751
3726 3752 def setup(self, dataOut, ipp1=None, ipp2=None, tx1=None, tx2=None, dh1=None,
3727 3753 idate=None, startH=None, endH=None):
3728 3754
3729 3755
3730 3756 self.ipp1 = ipp1
3731 3757 self.ipp2 = ipp2
3732 3758 self.tx1 = tx1
3733 3759 self.tx2 = tx2
3734 3760 self.dh1 = dh1
3735 3761
3736 3762 _maxIpp1R = dataOut.heightList.max()
3737 3763
3738 3764 _n_repeats = int(_maxIpp1R / ipp2)
3739 3765 _init_hIntf = (tx1 + ipp2/2)+ dh1
3740 3766 _n_hIntf = int(tx2 / dh1)
3741 3767
3742 3768 self.heightsList = [_init_hIntf+n*ipp2 for n in range(_n_repeats) ]
3743 3769 heiList = dataOut.heightList
3744 3770 self.heights_indx = [getHei_index(h,h,heiList)[0] for h in self.heightsList]
3745 3771
3746 3772 self.heights_indx = [ numpy.asarray([k for k in range(_n_hIntf+2)])+(getHei_index(h,h,heiList)[0] -1) for h in self.heightsList]
3747 3773
3748 3774 self.heights_indx = numpy.asarray(self.heights_indx )
3749 3775 self.isConfig = True
3750 3776 self.startTime = datetime.datetime.combine(idate,startH)
3751 3777 self.endTime = datetime.datetime.combine(idate,endH)
3752 3778 #print(self.startTime, self.endTime)
3753 3779 #print("nrepeats: ", _n_repeats, " _nH: ",_n_hIntf )
3754 3780
3755 3781 log.warning("Heights set to zero (km): ", self.name)
3756 3782 log.warning(str((dataOut.heightList[self.heights_indx].flatten())), self.name)
3757 3783 log.warning("Be careful with the selection of heights for noise calculation!")
3758 3784
3759 3785 def run(self, dataOut, ipp1=None, ipp2=None, tx1=None, tx2=None, dh1=None, idate=None,
3760 3786 startH=None, endH=None):
3761 3787 #print(locals().values())
3762 3788 if None in locals().values():
3763 3789 log.warning('Missing kwargs, invalid values """None""" ', self.name)
3764 3790 return dataOut
3765 3791
3766 3792
3767 3793 if not self.isConfig:
3768 3794 self.setup(dataOut, ipp1=ipp1, ipp2=ipp2, tx1=tx1, tx2=tx2, dh1=dh1,
3769 3795 idate=idate, startH=startH, endH=endH)
3770 3796
3771 3797 dataOut.flagProfilesByRange = False
3772 3798 currentTime = datetime.datetime.fromtimestamp(dataOut.utctime)
3773 3799
3774 3800 if currentTime < self.startTime or currentTime > self.endTime:
3775 3801 return dataOut
3776 3802
3777 3803 for ch in range(dataOut.data.shape[0]):
3778 3804
3779 3805 for hk in self.heights_indx.flatten():
3780 3806 if dataOut.data.ndim < 3:
3781 3807 dataOut.data[ch,hk] = 0.0 + 0.0j
3782 3808 else:
3783 3809 dataOut.data[ch,:,hk] = 0.0 + 0.0j
3784 3810
3785 3811 dataOut.flagProfilesByRange = True
3786 3812
3787 3813 return dataOut No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now