##// END OF EJS Templates
schain-jc
RSS
Udating from v2.3
George Yong -
r1174:0a82420b353b merge
parent child
Show More
Side by Side Unified
Context file:
r1174:0a82420b353b Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1148 +1,1154
1 1
2 2 import os
3 3 import time
4 4 import glob
5 5 import datetime
6 6 from multiprocessing import Process
7 7
8 8 import zmq
9 9 import numpy
10 10 import matplotlib
11 11 import matplotlib.pyplot as plt
12 12 from matplotlib.patches import Polygon
13 13 from mpl_toolkits.axes_grid1 import make_axes_locatable
14 14 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
15 15
16 16 from schainpy.model.proc.jroproc_base import Operation
17 17 from schainpy.utils import log
18 18
19 19 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
20 20 blu_values = matplotlib.pyplot.get_cmap(
21 21 'seismic_r', 20)(numpy.arange(20))[10:15]
22 22 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
23 23 'jro', numpy.vstack((blu_values, jet_values)))
24 24 matplotlib.pyplot.register_cmap(cmap=ncmap)
25 25
26 26 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis', 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
27 27
28 28 EARTH_RADIUS = 6.3710e3
29 29
30 30 def ll2xy(lat1, lon1, lat2, lon2):
31 31
32 32 p = 0.017453292519943295
33 33 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
34 34 r = 12742 * numpy.arcsin(numpy.sqrt(a))
35 35 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)*numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
36 36 theta = -theta + numpy.pi/2
37 37 return r*numpy.cos(theta), r*numpy.sin(theta)
38 38
39 39 def km2deg(km):
40 40 '''
41 41 Convert distance in km to degrees
42 42 '''
43 43
44 44 return numpy.rad2deg(km/EARTH_RADIUS)
45 45
46 46 def figpause(interval):
47 47 backend = plt.rcParams['backend']
48 48 if backend in matplotlib.rcsetup.interactive_bk:
49 49 figManager = matplotlib._pylab_helpers.Gcf.get_active()
50 50 if figManager is not None:
51 51 canvas = figManager.canvas
52 52 if canvas.figure.stale:
53 53 canvas.draw()
54 54 try:
55 55 canvas.start_event_loop(interval)
56 56 except:
57 57 pass
58 58 return
59 59
60 60 def popup(message):
61 61 '''
62 62 '''
63 63
64 64 fig = plt.figure(figsize=(12, 8), facecolor='r')
65 65 text = '\n'.join([s.strip() for s in message.split(':')])
66 66 fig.text(0.01, 0.5, text, ha='left', va='center', size='20', weight='heavy', color='w')
67 67 fig.show()
68 68 figpause(1000)
69 69
70 70
71 71 class PlotData(Operation, Process):
72 72 '''
73 73 Base class for Schain plotting operations
74 74 '''
75 75
76 76 CODE = 'Figure'
77 77 colormap = 'jro'
78 78 bgcolor = 'white'
79 79 CONFLATE = False
80 80 __missing = 1E30
81 81
82 82 __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax',
83 83 'zlimits', 'xlabel', 'ylabel', 'xaxis','cb_label', 'title',
84 84 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure',
85 85 'showprofile', 'decimation', 'ftp']
86 86
87 87 def __init__(self, **kwargs):
88 88
89 89 Operation.__init__(self, plot=True, **kwargs)
90 90 Process.__init__(self)
91 91
92 92 self.kwargs['code'] = self.CODE
93 93 self.mp = False
94 94 self.data = None
95 95 self.isConfig = False
96 96 self.figures = []
97 97 self.axes = []
98 98 self.cb_axes = []
99 99 self.localtime = kwargs.pop('localtime', True)
100 100 self.show = kwargs.get('show', True)
101 101 self.save = kwargs.get('save', False)
102 102 self.ftp = kwargs.get('ftp', False)
103 103 self.colormap = kwargs.get('colormap', self.colormap)
104 104 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
105 105 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
106 106 self.colormaps = kwargs.get('colormaps', None)
107 107 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
108 108 self.showprofile = kwargs.get('showprofile', False)
109 109 self.title = kwargs.get('wintitle', self.CODE.upper())
110 110 self.cb_label = kwargs.get('cb_label', None)
111 111 self.cb_labels = kwargs.get('cb_labels', None)
112 112 self.labels = kwargs.get('labels', None)
113 113 self.xaxis = kwargs.get('xaxis', 'frequency')
114 114 self.zmin = kwargs.get('zmin', None)
115 115 self.zmax = kwargs.get('zmax', None)
116 116 self.zlimits = kwargs.get('zlimits', None)
117 117 self.xmin = kwargs.get('xmin', None)
118 118 self.xmax = kwargs.get('xmax', None)
119 119 self.xrange = kwargs.get('xrange', 24)
120 120 self.xscale = kwargs.get('xscale', None)
121 121 self.ymin = kwargs.get('ymin', None)
122 122 self.ymax = kwargs.get('ymax', None)
123 123 self.yscale = kwargs.get('yscale', None)
124 124 self.xlabel = kwargs.get('xlabel', None)
125 125 self.decimation = kwargs.get('decimation', None)
126 126 self.showSNR = kwargs.get('showSNR', False)
127 127 self.oneFigure = kwargs.get('oneFigure', True)
128 128 self.width = kwargs.get('width', None)
129 129 self.height = kwargs.get('height', None)
130 130 self.colorbar = kwargs.get('colorbar', True)
131 131 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
132 132 self.channels = kwargs.get('channels', None)
133 133 self.titles = kwargs.get('titles', [])
134 134 self.polar = False
135 135 self.grid = kwargs.get('grid', False)
136 136
137 137 def __fmtTime(self, x, pos):
138 138 '''
139 139 '''
140 140
141 141 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
142 142
143 143 def __setup(self):
144 144 '''
145 145 Common setup for all figures, here figures and axes are created
146 146 '''
147 147
148 148 if self.CODE not in self.data:
149 149 raise ValueError(log.error('Missing data for {}'.format(self.CODE),
150 150 self.name))
151 151
152 152 self.setup()
153 153
154 154 self.time_label = 'LT' if self.localtime else 'UTC'
155 155 if self.data.localtime:
156 156 self.getDateTime = datetime.datetime.fromtimestamp
157 157 else:
158 158 self.getDateTime = datetime.datetime.utcfromtimestamp
159 159
160 160 if self.width is None:
161 161 self.width = 8
162 162
163 163 self.figures = []
164 164 self.axes = []
165 165 self.cb_axes = []
166 166 self.pf_axes = []
167 167 self.cmaps = []
168 168
169 169 size = '15%' if self.ncols == 1 else '30%'
170 170 pad = '4%' if self.ncols == 1 else '8%'
171 171
172 172 if self.oneFigure:
173 173 if self.height is None:
174 174 self.height = 1.4 * self.nrows + 1
175 175 fig = plt.figure(figsize=(self.width, self.height),
176 176 edgecolor='k',
177 177 facecolor='w')
178 178 self.figures.append(fig)
179 179 for n in range(self.nplots):
180 180 ax = fig.add_subplot(self.nrows, self.ncols,
181 181 n + 1, polar=self.polar)
182 182 ax.tick_params(labelsize=8)
183 183 ax.firsttime = True
184 184 ax.index = 0
185 185 ax.press = None
186 186 self.axes.append(ax)
187 187 if self.showprofile:
188 188 cax = self.__add_axes(ax, size=size, pad=pad)
189 189 cax.tick_params(labelsize=8)
190 190 self.pf_axes.append(cax)
191 191 else:
192 192 if self.height is None:
193 193 self.height = 3
194 194 for n in range(self.nplots):
195 195 fig = plt.figure(figsize=(self.width, self.height),
196 196 edgecolor='k',
197 197 facecolor='w')
198 198 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
199 199 ax.tick_params(labelsize=8)
200 200 ax.firsttime = True
201 201 ax.index = 0
202 202 ax.press = None
203 203 self.figures.append(fig)
204 204 self.axes.append(ax)
205 205 if self.showprofile:
206 206 cax = self.__add_axes(ax, size=size, pad=pad)
207 207 cax.tick_params(labelsize=8)
208 208 self.pf_axes.append(cax)
209 209
210 210 for n in range(self.nrows):
211 211 if self.colormaps is not None:
212 212 cmap = plt.get_cmap(self.colormaps[n])
213 213 else:
214 214 cmap = plt.get_cmap(self.colormap)
215 215 cmap.set_bad(self.bgcolor, 1.)
216 216 self.cmaps.append(cmap)
217 217
218 218 for fig in self.figures:
219 219 fig.canvas.mpl_connect('key_press_event', self.OnKeyPress)
220 220 fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll)
221 221 fig.canvas.mpl_connect('button_press_event', self.onBtnPress)
222 222 fig.canvas.mpl_connect('motion_notify_event', self.onMotion)
223 223 fig.canvas.mpl_connect('button_release_event', self.onBtnRelease)
224 224 if self.show:
225 225 fig.show()
226 226
227 227 def OnKeyPress(self, event):
228 228 '''
229 229 Event for pressing keys (up, down) change colormap
230 230 '''
231 231 ax = event.inaxes
232 232 if ax in self.axes:
233 233 if event.key == 'down':
234 234 ax.index += 1
235 235 elif event.key == 'up':
236 236 ax.index -= 1
237 237 if ax.index < 0:
238 238 ax.index = len(CMAPS) - 1
239 239 elif ax.index == len(CMAPS):
240 240 ax.index = 0
241 241 cmap = CMAPS[ax.index]
242 242 ax.cbar.set_cmap(cmap)
243 243 ax.cbar.draw_all()
244 244 ax.plt.set_cmap(cmap)
245 245 ax.cbar.patch.figure.canvas.draw()
246 246 self.colormap = cmap.name
247 247
248 248 def OnBtnScroll(self, event):
249 249 '''
250 250 Event for scrolling, scale figure
251 251 '''
252 252 cb_ax = event.inaxes
253 253 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
254 254 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
255 255 pt = ax.cbar.ax.bbox.get_points()[:, 1]
256 256 nrm = ax.cbar.norm
257 257 vmin, vmax, p0, p1, pS = (
258 258 nrm.vmin, nrm.vmax, pt[0], pt[1], event.y)
259 259 scale = 2 if event.step == 1 else 0.5
260 260 point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0)
261 261 ax.cbar.norm.vmin = point - scale * (point - vmin)
262 262 ax.cbar.norm.vmax = point - scale * (point - vmax)
263 263 ax.plt.set_norm(ax.cbar.norm)
264 264 ax.cbar.draw_all()
265 265 ax.cbar.patch.figure.canvas.draw()
266 266
267 267 def onBtnPress(self, event):
268 268 '''
269 269 Event for mouse button press
270 270 '''
271 271 cb_ax = event.inaxes
272 272 if cb_ax is None:
273 273 return
274 274
275 275 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
276 276 cb_ax.press = event.x, event.y
277 277 else:
278 278 cb_ax.press = None
279 279
280 280 def onMotion(self, event):
281 281 '''
282 282 Event for move inside colorbar
283 283 '''
284 284 cb_ax = event.inaxes
285 285 if cb_ax is None:
286 286 return
287 287 if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]:
288 288 return
289 289 if cb_ax.press is None:
290 290 return
291 291
292 292 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
293 293 xprev, yprev = cb_ax.press
294 294 dx = event.x - xprev
295 295 dy = event.y - yprev
296 296 cb_ax.press = event.x, event.y
297 297 scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin
298 298 perc = 0.03
299 299
300 300 if event.button == 1:
301 301 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
302 302 ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy)
303 303 elif event.button == 3:
304 304 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
305 305 ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy)
306 306
307 307 ax.cbar.draw_all()
308 308 ax.plt.set_norm(ax.cbar.norm)
309 309 ax.cbar.patch.figure.canvas.draw()
310 310
311 311 def onBtnRelease(self, event):
312 312 '''
313 313 Event for mouse button release
314 314 '''
315 315 cb_ax = event.inaxes
316 316 if cb_ax is not None:
317 317 cb_ax.press = None
318 318
319 319 def __add_axes(self, ax, size='30%', pad='8%'):
320 320 '''
321 321 Add new axes to the given figure
322 322 '''
323 323 divider = make_axes_locatable(ax)
324 324 nax = divider.new_horizontal(size=size, pad=pad)
325 325 ax.figure.add_axes(nax)
326 326 return nax
327 327
328 328 self.setup()
329 329
330 330 def setup(self):
331 331 '''
332 332 This method should be implemented in the child class, the following
333 333 attributes should be set:
334 334
335 335 self.nrows: number of rows
336 336 self.ncols: number of cols
337 337 self.nplots: number of plots (channels or pairs)
338 338 self.ylabel: label for Y axes
339 339 self.titles: list of axes title
340 340
341 341 '''
342 342 raise NotImplementedError
343 343
344 344 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
345 345 '''
346 346 Create a masked array for missing data
347 347 '''
348 348 if x_buffer.shape[0] < 2:
349 349 return x_buffer, y_buffer, z_buffer
350 350
351 351 deltas = x_buffer[1:] - x_buffer[0:-1]
352 352 x_median = numpy.median(deltas)
353 353
354 354 index = numpy.where(deltas > 5 * x_median)
355 355
356 356 if len(index[0]) != 0:
357 357 z_buffer[::, index[0], ::] = self.__missing
358 358 z_buffer = numpy.ma.masked_inside(z_buffer,
359 359 0.99 * self.__missing,
360 360 1.01 * self.__missing)
361 361
362 362 return x_buffer, y_buffer, z_buffer
363 363
364 364 def decimate(self):
365 365
366 366 # dx = int(len(self.x)/self.__MAXNUMX) + 1
367 367 dy = int(len(self.y) / self.decimation) + 1
368 368
369 369 # x = self.x[::dx]
370 370 x = self.x
371 371 y = self.y[::dy]
372 372 z = self.z[::, ::, ::dy]
373 373
374 374 return x, y, z
375 375
376 376 def format(self):
377 377 '''
378 378 Set min and max values, labels, ticks and titles
379 379 '''
380 380
381 381 if self.xmin is None:
382 382 xmin = self.min_time
383 383 else:
384 384 if self.xaxis is 'time':
385 385 dt = self.getDateTime(self.min_time)
386 386 xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) -
387 387 datetime.datetime(1970, 1, 1)).total_seconds()
388 388 if self.data.localtime:
389 389 xmin += time.timezone
390 390 else:
391 391 xmin = self.xmin
392 392
393 393 if self.xmax is None:
394 394 xmax = xmin + self.xrange * 60 * 60
395 395 else:
396 396 if self.xaxis is 'time':
397 397 dt = self.getDateTime(self.max_time)
398 398 xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) -
399 399 datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds()
400 400 if self.data.localtime:
401 401 xmax += time.timezone
402 402 else:
403 403 xmax = self.xmax
404 404
405 405 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
406 406 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
407 407
408 408 Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])
409 409 i = 1 if numpy.where(abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0]
410 410 ystep = Y[i] / 10.
411 411
412 if self.xaxis is not 'time':
413 X = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])/2.
414 i = 1 if numpy.where(abs(xmax-xmin) <= X)[0][0] < 0 else numpy.where(abs(xmax-xmin) <= X)[0][0]
415 xstep = X[i] / 10.
416
412 417 for n, ax in enumerate(self.axes):
413 418 if ax.firsttime:
414 419 ax.set_facecolor(self.bgcolor)
415 420 ax.yaxis.set_major_locator(MultipleLocator(ystep))
416 ax.xaxis.set_major_locator(MultipleLocator(ystep))
417 421 if self.xscale:
418 422 ax.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{0:g}'.format(x*self.xscale)))
419 423 if self.xscale:
420 424 ax.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{0:g}'.format(x*self.yscale)))
421 425 if self.xaxis is 'time':
422 426 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
423 427 ax.xaxis.set_major_locator(LinearLocator(9))
428 else:
429 ax.xaxis.set_major_locator(MultipleLocator(xstep))
424 430 if self.xlabel is not None:
425 431 ax.set_xlabel(self.xlabel)
426 432 ax.set_ylabel(self.ylabel)
427 433 ax.firsttime = False
428 434 if self.showprofile:
429 435 self.pf_axes[n].set_ylim(ymin, ymax)
430 436 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
431 437 self.pf_axes[n].set_xlabel('dB')
432 438 self.pf_axes[n].grid(b=True, axis='x')
433 439 [tick.set_visible(False)
434 440 for tick in self.pf_axes[n].get_yticklabels()]
435 441 if self.colorbar:
436 442 ax.cbar = plt.colorbar(
437 443 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
438 444 ax.cbar.ax.tick_params(labelsize=8)
439 445 ax.cbar.ax.press = None
440 446 if self.cb_label:
441 447 ax.cbar.set_label(self.cb_label, size=8)
442 448 elif self.cb_labels:
443 449 ax.cbar.set_label(self.cb_labels[n], size=8)
444 450 else:
445 451 ax.cbar = None
446 452 if self.grid:
447 453 ax.grid(True)
448 454
449 455 if not self.polar:
450 456 ax.set_xlim(xmin, xmax)
451 457 ax.set_ylim(ymin, ymax)
452 458 ax.set_title('{} {} {}'.format(
453 459 self.titles[n],
454 460 self.getDateTime(self.max_time).strftime('%Y-%m-%dT%H:%M:%S'),
455 461 self.time_label),
456 462 size=8)
457 463 else:
458 464 ax.set_title('{}'.format(self.titles[n]), size=8)
459 465 ax.set_ylim(0, 90)
460 466 ax.set_yticks(numpy.arange(0, 90, 20))
461 467 ax.yaxis.labelpad = 40
462 468
463 469 def __plot(self):
464 470 '''
465 471 '''
466 472 log.log('Plotting', self.name)
467 473
468 474 try:
469 475 self.plot()
470 476 self.format()
471 477 except Exception as e:
472 478 log.warning('{} Plot could not be updated... check data'.format(self.CODE), self.name)
473 479 log.error(str(e), '')
474 480 return
475 481
476 482 for n, fig in enumerate(self.figures):
477 483 if self.nrows == 0 or self.nplots == 0:
478 484 log.warning('No data', self.name)
479 485 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
480 486 fig.canvas.manager.set_window_title(self.CODE)
481 487 continue
482 488
483 489 fig.tight_layout()
484 490 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
485 491 self.getDateTime(self.max_time).strftime('%Y/%m/%d')))
486 492 fig.canvas.draw()
487 493
488 494 if self.save and (self.data.ended or not self.data.buffering):
489 495
490 496 if self.save_labels:
491 497 labels = self.save_labels
492 498 else:
493 499 labels = list(range(self.nrows))
494 500
495 501 if self.oneFigure:
496 502 label = ''
497 503 else:
498 504 label = '-{}'.format(labels[n])
499 505 figname = os.path.join(
500 506 self.save,
501 507 self.CODE,
502 508 '{}{}_{}.png'.format(
503 509 self.CODE,
504 510 label,
505 511 self.getDateTime(self.saveTime).strftime(
506 512 '%Y%m%d_%H%M%S'),
507 513 )
508 514 )
509 515 log.log('Saving figure: {}'.format(figname), self.name)
510 516 if not os.path.isdir(os.path.dirname(figname)):
511 517 os.makedirs(os.path.dirname(figname))
512 518 fig.savefig(figname)
513 519
514 520 def plot(self):
515 521 '''
516 522 '''
517 523 raise NotImplementedError
518 524
519 525 def run(self):
520 526
521 527 log.log('Starting', self.name)
522 528
523 529 context = zmq.Context()
524 530 receiver = context.socket(zmq.SUB)
525 531 receiver.setsockopt(zmq.SUBSCRIBE, '')
526 532 receiver.setsockopt(zmq.CONFLATE, self.CONFLATE)
527 533
528 534 if 'server' in self.kwargs['parent']:
529 535 receiver.connect(
530 536 'ipc:///tmp/{}.plots'.format(self.kwargs['parent']['server']))
531 537 else:
532 538 receiver.connect("ipc:///tmp/zmq.plots")
533 539
534 540 while True:
535 541 try:
536 542 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)
537 543 if self.data.localtime and self.localtime:
538 544 self.times = self.data.times
539 545 elif self.data.localtime and not self.localtime:
540 546 self.times = self.data.times + time.timezone
541 547 elif not self.data.localtime and self.localtime:
542 548 self.times = self.data.times - time.timezone
543 549 else:
544 550 self.times = self.data.times
545 551
546 552 self.min_time = self.times[0]
547 553 self.max_time = self.times[-1]
548 554
549 555 if self.isConfig is False:
550 556 self.__setup()
551 557 self.isConfig = True
552 558
553 559 self.__plot()
554 560
555 561 except zmq.Again as e:
556 562 if self.data and self.data.ended:
557 563 break
558 564 log.log('Waiting for data...')
559 565 if self.data:
560 566 figpause(self.data.throttle)
561 567 else:
562 568 time.sleep(2)
563 569
564 570 def close(self):
565 571 if self.data:
566 572 self.__plot()
567 573
568 574
569 575 class PlotSpectraData(PlotData):
570 576 '''
571 577 Plot for Spectra data
572 578 '''
573 579
574 580 CODE = 'spc'
575 581 colormap = 'jro'
576 582
577 583 def setup(self):
578 584 self.nplots = len(self.data.channels)
579 585 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
580 586 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
581 587 self.width = 3.4 * self.ncols
582 588 self.height = 3 * self.nrows
583 589 self.cb_label = 'dB'
584 590 if self.showprofile:
585 591 self.width += 0.8 * self.ncols
586 592
587 593 self.ylabel = 'Range [km]'
588 594
589 595 def plot(self):
590 596 if self.xaxis == "frequency":
591 597 x = self.data.xrange[0]
592 598 self.xlabel = "Frequency (kHz)"
593 599 elif self.xaxis == "time":
594 600 x = self.data.xrange[1]
595 601 self.xlabel = "Time (ms)"
596 602 else:
597 603 x = self.data.xrange[2]
598 604 self.xlabel = "Velocity (m/s)"
599 605
600 606 if self.CODE == 'spc_mean':
601 607 x = self.data.xrange[2]
602 608 self.xlabel = "Velocity (m/s)"
603 609
604 610 self.titles = []
605 611
606 612 y = self.data.heights
607 613 self.y = y
608 614 z = self.data['spc']
609 615
610 616 for n, ax in enumerate(self.axes):
611 617 noise = self.data['noise'][n][-1]
612 618 if self.CODE == 'spc_mean':
613 619 mean = self.data['mean'][n][-1]
614 620 if ax.firsttime:
615 621 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
616 622 self.xmin = self.xmin if self.xmin else -self.xmax
617 623 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
618 624 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
619 625 ax.plt = ax.pcolormesh(x, y, z[n].T,
620 626 vmin=self.zmin,
621 627 vmax=self.zmax,
622 628 cmap=plt.get_cmap(self.colormap)
623 629 )
624 630
625 631 if self.showprofile:
626 632 ax.plt_profile = self.pf_axes[n].plot(
627 633 self.data['rti'][n][-1], y)[0]
628 634 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
629 635 color="k", linestyle="dashed", lw=1)[0]
630 636 if self.CODE == 'spc_mean':
631 637 ax.plt_mean = ax.plot(mean, y, color='k')[0]
632 638 else:
633 639 ax.plt.set_array(z[n].T.ravel())
634 640 if self.showprofile:
635 641 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
636 642 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
637 643 if self.CODE == 'spc_mean':
638 644 ax.plt_mean.set_data(mean, y)
639 645
640 646 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
641 647 self.saveTime = self.max_time
642 648
643 649
644 650 class PlotCrossSpectraData(PlotData):
645 651
646 652 CODE = 'cspc'
647 653 zmin_coh = None
648 654 zmax_coh = None
649 655 zmin_phase = None
650 656 zmax_phase = None
651 657
652 658 def setup(self):
653 659
654 660 self.ncols = 4
655 661 self.nrows = len(self.data.pairs)
656 662 self.nplots = self.nrows * 4
657 663 self.width = 3.4 * self.ncols
658 664 self.height = 3 * self.nrows
659 665 self.ylabel = 'Range [km]'
660 666 self.showprofile = False
661 667
662 668 def plot(self):
663 669
664 670 if self.xaxis == "frequency":
665 671 x = self.data.xrange[0]
666 672 self.xlabel = "Frequency (kHz)"
667 673 elif self.xaxis == "time":
668 674 x = self.data.xrange[1]
669 675 self.xlabel = "Time (ms)"
670 676 else:
671 677 x = self.data.xrange[2]
672 678 self.xlabel = "Velocity (m/s)"
673 679
674 680 self.titles = []
675 681
676 682 y = self.data.heights
677 683 self.y = y
678 684 spc = self.data['spc']
679 685 cspc = self.data['cspc']
680 686
681 687 for n in range(self.nrows):
682 688 noise = self.data['noise'][n][-1]
683 689 pair = self.data.pairs[n]
684 690 ax = self.axes[4 * n]
685 691 ax3 = self.axes[4 * n + 3]
686 692 if ax.firsttime:
687 693 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
688 694 self.xmin = self.xmin if self.xmin else -self.xmax
689 695 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
690 696 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
691 697 ax.plt = ax.pcolormesh(x, y, spc[pair[0]].T,
692 698 vmin=self.zmin,
693 699 vmax=self.zmax,
694 700 cmap=plt.get_cmap(self.colormap)
695 701 )
696 702 else:
697 703 ax.plt.set_array(spc[pair[0]].T.ravel())
698 704 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
699 705
700 706 ax = self.axes[4 * n + 1]
701 707 if ax.firsttime:
702 708 ax.plt = ax.pcolormesh(x, y, spc[pair[1]].T,
703 709 vmin=self.zmin,
704 710 vmax=self.zmax,
705 711 cmap=plt.get_cmap(self.colormap)
706 712 )
707 713 else:
708 714 ax.plt.set_array(spc[pair[1]].T.ravel())
709 715 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
710 716
711 717 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
712 718 coh = numpy.abs(out)
713 719 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
714 720
715 721 ax = self.axes[4 * n + 2]
716 722 if ax.firsttime:
717 723 ax.plt = ax.pcolormesh(x, y, coh.T,
718 724 vmin=0,
719 725 vmax=1,
720 726 cmap=plt.get_cmap(self.colormap_coh)
721 727 )
722 728 else:
723 729 ax.plt.set_array(coh.T.ravel())
724 730 self.titles.append(
725 731 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
726 732
727 733 ax = self.axes[4 * n + 3]
728 734 if ax.firsttime:
729 735 ax.plt = ax.pcolormesh(x, y, phase.T,
730 736 vmin=-180,
731 737 vmax=180,
732 738 cmap=plt.get_cmap(self.colormap_phase)
733 739 )
734 740 else:
735 741 ax.plt.set_array(phase.T.ravel())
736 742 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
737 743
738 744 self.saveTime = self.max_time
739 745
740 746
741 747 class PlotSpectraMeanData(PlotSpectraData):
742 748 '''
743 749 Plot for Spectra and Mean
744 750 '''
745 751 CODE = 'spc_mean'
746 752 colormap = 'jro'
747 753
748 754
749 755 class PlotRTIData(PlotData):
750 756 '''
751 757 Plot for RTI data
752 758 '''
753 759
754 760 CODE = 'rti'
755 761 colormap = 'jro'
756 762
757 763 def setup(self):
758 764 self.xaxis = 'time'
759 765 self.ncols = 1
760 766 self.nrows = len(self.data.channels)
761 767 self.nplots = len(self.data.channels)
762 768 self.ylabel = 'Range [km]'
763 769 self.cb_label = 'dB'
764 770 self.titles = ['{} Channel {}'.format(
765 771 self.CODE.upper(), x) for x in range(self.nrows)]
766 772
767 773 def plot(self):
768 774 self.x = self.times
769 775 self.y = self.data.heights
770 776 self.z = self.data[self.CODE]
771 777 self.z = numpy.ma.masked_invalid(self.z)
772 778
773 779 if self.decimation is None:
774 780 x, y, z = self.fill_gaps(self.x, self.y, self.z)
775 781 else:
776 782 x, y, z = self.fill_gaps(*self.decimate())
777 783
778 784 for n, ax in enumerate(self.axes):
779 785 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
780 786 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
781 787 if ax.firsttime:
782 788 ax.plt = ax.pcolormesh(x, y, z[n].T,
783 789 vmin=self.zmin,
784 790 vmax=self.zmax,
785 791 cmap=plt.get_cmap(self.colormap)
786 792 )
787 793 if self.showprofile:
788 794 ax.plot_profile = self.pf_axes[n].plot(
789 795 self.data['rti'][n][-1], self.y)[0]
790 796 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
791 797 color="k", linestyle="dashed", lw=1)[0]
792 798 else:
793 799 ax.collections.remove(ax.collections[0])
794 800 ax.plt = ax.pcolormesh(x, y, z[n].T,
795 801 vmin=self.zmin,
796 802 vmax=self.zmax,
797 803 cmap=plt.get_cmap(self.colormap)
798 804 )
799 805 if self.showprofile:
800 806 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
801 807 ax.plot_noise.set_data(numpy.repeat(
802 808 self.data['noise'][n][-1], len(self.y)), self.y)
803 809
804 810 self.saveTime = self.min_time
805 811
806 812
807 813 class PlotCOHData(PlotRTIData):
808 814 '''
809 815 Plot for Coherence data
810 816 '''
811 817
812 818 CODE = 'coh'
813 819
814 820 def setup(self):
815 821 self.xaxis = 'time'
816 822 self.ncols = 1
817 823 self.nrows = len(self.data.pairs)
818 824 self.nplots = len(self.data.pairs)
819 825 self.ylabel = 'Range [km]'
820 826 if self.CODE == 'coh':
821 827 self.cb_label = ''
822 828 self.titles = [
823 829 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
824 830 else:
825 831 self.cb_label = 'Degrees'
826 832 self.titles = [
827 833 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
828 834
829 835
830 836 class PlotPHASEData(PlotCOHData):
831 837 '''
832 838 Plot for Phase map data
833 839 '''
834 840
835 841 CODE = 'phase'
836 842 colormap = 'seismic'
837 843
838 844
839 845 class PlotNoiseData(PlotData):
840 846 '''
841 847 Plot for noise
842 848 '''
843 849
844 850 CODE = 'noise'
845 851
846 852 def setup(self):
847 853 self.xaxis = 'time'
848 854 self.ncols = 1
849 855 self.nrows = 1
850 856 self.nplots = 1
851 857 self.ylabel = 'Intensity [dB]'
852 858 self.titles = ['Noise']
853 859 self.colorbar = False
854 860
855 861 def plot(self):
856 862
857 863 x = self.times
858 864 xmin = self.min_time
859 865 xmax = xmin + self.xrange * 60 * 60
860 866 Y = self.data[self.CODE]
861 867
862 868 if self.axes[0].firsttime:
863 869 for ch in self.data.channels:
864 870 y = Y[ch]
865 871 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
866 872 plt.legend()
867 873 else:
868 874 for ch in self.data.channels:
869 875 y = Y[ch]
870 876 self.axes[0].lines[ch].set_data(x, y)
871 877
872 878 self.ymin = numpy.nanmin(Y) - 5
873 879 self.ymax = numpy.nanmax(Y) + 5
874 880 self.saveTime = self.min_time
875 881
876 882
877 883 class PlotSNRData(PlotRTIData):
878 884 '''
879 885 Plot for SNR Data
880 886 '''
881 887
882 888 CODE = 'snr'
883 889 colormap = 'jet'
884 890
885 891
886 892 class PlotDOPData(PlotRTIData):
887 893 '''
888 894 Plot for DOPPLER Data
889 895 '''
890 896
891 897 CODE = 'dop'
892 898 colormap = 'jet'
893 899
894 900
895 901 class PlotSkyMapData(PlotData):
896 902 '''
897 903 Plot for meteors detection data
898 904 '''
899 905
900 906 CODE = 'param'
901 907
902 908 def setup(self):
903 909
904 910 self.ncols = 1
905 911 self.nrows = 1
906 912 self.width = 7.2
907 913 self.height = 7.2
908 914 self.nplots = 1
909 915 self.xlabel = 'Zonal Zenith Angle (deg)'
910 916 self.ylabel = 'Meridional Zenith Angle (deg)'
911 917 self.polar = True
912 918 self.ymin = -180
913 919 self.ymax = 180
914 920 self.colorbar = False
915 921
916 922 def plot(self):
917 923
918 924 arrayParameters = numpy.concatenate(self.data['param'])
919 925 error = arrayParameters[:, -1]
920 926 indValid = numpy.where(error == 0)[0]
921 927 finalMeteor = arrayParameters[indValid, :]
922 928 finalAzimuth = finalMeteor[:, 3]
923 929 finalZenith = finalMeteor[:, 4]
924 930
925 931 x = finalAzimuth * numpy.pi / 180
926 932 y = finalZenith
927 933
928 934 ax = self.axes[0]
929 935
930 936 if ax.firsttime:
931 937 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
932 938 else:
933 939 ax.plot.set_data(x, y)
934 940
935 941 dt1 = self.getDateTime(self.min_time).strftime('%y/%m/%d %H:%M:%S')
936 942 dt2 = self.getDateTime(self.max_time).strftime('%y/%m/%d %H:%M:%S')
937 943 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
938 944 dt2,
939 945 len(x))
940 946 self.titles[0] = title
941 947 self.saveTime = self.max_time
942 948
943 949
944 950 class PlotParamData(PlotRTIData):
945 951 '''
946 952 Plot for data_param object
947 953 '''
948 954
949 955 CODE = 'param'
950 956 colormap = 'seismic'
951 957
952 958 def setup(self):
953 959 self.xaxis = 'time'
954 960 self.ncols = 1
955 961 self.nrows = self.data.shape(self.CODE)[0]
956 962 self.nplots = self.nrows
957 963 if self.showSNR:
958 964 self.nrows += 1
959 965 self.nplots += 1
960 966
961 967 self.ylabel = 'Height [km]'
962 968 if not self.titles:
963 969 self.titles = self.data.parameters \
964 970 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
965 971 if self.showSNR:
966 972 self.titles.append('SNR')
967 973
968 974 def plot(self):
969 975 self.data.normalize_heights()
970 976 self.x = self.times
971 977 self.y = self.data.heights
972 978 if self.showSNR:
973 979 self.z = numpy.concatenate(
974 980 (self.data[self.CODE], self.data['snr'])
975 981 )
976 982 else:
977 983 self.z = self.data[self.CODE]
978 984
979 985 self.z = numpy.ma.masked_invalid(self.z)
980 986
981 987 if self.decimation is None:
982 988 x, y, z = self.fill_gaps(self.x, self.y, self.z)
983 989 else:
984 990 x, y, z = self.fill_gaps(*self.decimate())
985 991
986 992 for n, ax in enumerate(self.axes):
987 993
988 994 self.zmax = self.zmax if self.zmax is not None else numpy.max(
989 995 self.z[n])
990 996 self.zmin = self.zmin if self.zmin is not None else numpy.min(
991 997 self.z[n])
992 998
993 999 if ax.firsttime:
994 1000 if self.zlimits is not None:
995 1001 self.zmin, self.zmax = self.zlimits[n]
996 1002
997 1003 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
998 1004 vmin=self.zmin,
999 1005 vmax=self.zmax,
1000 1006 cmap=self.cmaps[n]
1001 1007 )
1002 1008 else:
1003 1009 if self.zlimits is not None:
1004 1010 self.zmin, self.zmax = self.zlimits[n]
1005 1011 ax.collections.remove(ax.collections[0])
1006 1012 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
1007 1013 vmin=self.zmin,
1008 1014 vmax=self.zmax,
1009 1015 cmap=self.cmaps[n]
1010 1016 )
1011 1017
1012 1018 self.saveTime = self.min_time
1013 1019
1014 1020
1015 1021 class PlotOutputData(PlotParamData):
1016 1022 '''
1017 1023 Plot data_output object
1018 1024 '''
1019 1025
1020 1026 CODE = 'output'
1021 1027 colormap = 'seismic'
1022 1028
1023 1029
1024 1030 class PlotPolarMapData(PlotData):
1025 1031 '''
1026 1032 Plot for meteors detection data
1027 1033 '''
1028 1034
1029 1035 CODE = 'param'
1030 1036 colormap = 'seismic'
1031 1037
1032 1038 def setup(self):
1033 1039 self.ncols = 1
1034 1040 self.nrows = 1
1035 1041 self.width = 9
1036 1042 self.height = 8
1037 1043 self.mode = self.data.meta['mode']
1038 1044 if self.channels is not None:
1039 1045 self.nplots = len(self.channels)
1040 1046 self.nrows = len(self.channels)
1041 1047 else:
1042 1048 self.nplots = self.data.shape(self.CODE)[0]
1043 1049 self.nrows = self.nplots
1044 1050 self.channels = list(range(self.nplots))
1045 1051 if self.mode == 'E':
1046 1052 self.xlabel = 'Longitude'
1047 1053 self.ylabel = 'Latitude'
1048 1054 else:
1049 1055 self.xlabel = 'Range (km)'
1050 1056 self.ylabel = 'Height (km)'
1051 1057 self.bgcolor = 'white'
1052 1058 self.cb_labels = self.data.meta['units']
1053 1059 self.lat = self.data.meta['latitude']
1054 1060 self.lon = self.data.meta['longitude']
1055 1061 self.xmin, self.xmax = float(km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
1056 1062 self.ymin, self.ymax = float(km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
1057 1063 # self.polar = True
1058 1064
1059 1065 def plot(self):
1060 1066
1061 1067 for n, ax in enumerate(self.axes):
1062 1068 data = self.data['param'][self.channels[n]]
1063 1069
1064 1070 zeniths = numpy.linspace(0, self.data.meta['max_range'], data.shape[1])
1065 1071 if self.mode == 'E':
1066 1072 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
1067 1073 r, theta = numpy.meshgrid(zeniths, azimuths)
1068 1074 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
1069 1075 x = km2deg(x) + self.lon
1070 1076 y = km2deg(y) + self.lat
1071 1077 else:
1072 1078 azimuths = numpy.radians(self.data.heights)
1073 1079 r, theta = numpy.meshgrid(zeniths, azimuths)
1074 1080 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
1075 1081 self.y = zeniths
1076 1082
1077 1083 if ax.firsttime:
1078 1084 if self.zlimits is not None:
1079 1085 self.zmin, self.zmax = self.zlimits[n]
1080 1086 ax.plt = ax.pcolormesh(#r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
1081 1087 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
1082 1088 vmin=self.zmin,
1083 1089 vmax=self.zmax,
1084 1090 cmap=self.cmaps[n])
1085 1091 else:
1086 1092 if self.zlimits is not None:
1087 1093 self.zmin, self.zmax = self.zlimits[n]
1088 1094 ax.collections.remove(ax.collections[0])
1089 1095 ax.plt = ax.pcolormesh(# r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
1090 1096 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
1091 1097 vmin=self.zmin,
1092 1098 vmax=self.zmax,
1093 1099 cmap=self.cmaps[n])
1094 1100
1095 1101 if self.mode == 'A':
1096 1102 continue
1097 1103
1098 1104 # plot district names
1099 1105 f = open('/data/workspace/schain_scripts/distrito.csv')
1100 1106 for line in f:
1101 1107 label, lon, lat = [s.strip() for s in line.split(',') if s]
1102 1108 lat = float(lat)
1103 1109 lon = float(lon)
1104 1110 # ax.plot(lon, lat, '.b', ms=2)
1105 1111 ax.text(lon, lat, label.decode('utf8'), ha='center', va='bottom', size='8', color='black')
1106 1112
1107 1113 # plot limites
1108 1114 limites =[]
1109 1115 tmp = []
1110 1116 for line in open('/data/workspace/schain_scripts/lima.csv'):
1111 1117 if '#' in line:
1112 1118 if tmp:
1113 1119 limites.append(tmp)
1114 1120 tmp = []
1115 1121 continue
1116 1122 values = line.strip().split(',')
1117 1123 tmp.append((float(values[0]), float(values[1])))
1118 1124 for points in limites:
1119 1125 ax.add_patch(Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
1120 1126
1121 1127 # plot Cuencas
1122 1128 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
1123 1129 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
1124 1130 values = [line.strip().split(',') for line in f]
1125 1131 points = [(float(s[0]), float(s[1])) for s in values]
1126 1132 ax.add_patch(Polygon(points, ec='b', fc='none'))
1127 1133
1128 1134 # plot grid
1129 1135 for r in (15, 30, 45, 60):
1130 1136 ax.add_artist(plt.Circle((self.lon, self.lat), km2deg(r), color='0.6', fill=False, lw=0.2))
1131 1137 ax.text(
1132 1138 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
1133 1139 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
1134 1140 '{}km'.format(r),
1135 1141 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
1136 1142
1137 1143 if self.mode == 'E':
1138 1144 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
1139 1145 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
1140 1146 else:
1141 1147 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
1142 1148 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
1143 1149
1144 1150 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
1145 1151 self.titles = ['{} {}'.format(self.data.parameters[x], title) for x in self.channels]
1146 1152 self.saveTime = self.max_time
1147 1153
1148 1154 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,680 +1,680
1 1 '''
2 2 Created on Jul 2, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import numpy
7 7
8 8 from schainpy.model.io.jroIO_base import LOCALTIME, JRODataReader, JRODataWriter
9 9 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
10 10 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
11 11 from schainpy.model.data.jrodata import Spectra
12 12
13 13 @MPDecorator
14 14 class SpectraReader(JRODataReader, ProcessingUnit):
15 15 """
16 16 Esta clase permite leer datos de espectros desde archivos procesados (.pdata). La lectura
17 17 de los datos siempre se realiza por bloques. Los datos leidos (array de 3 dimensiones)
18 18 son almacenados en tres buffer's para el Self Spectra, el Cross Spectra y el DC Channel.
19 19
20 20 paresCanalesIguales * alturas * perfiles (Self Spectra)
21 21 paresCanalesDiferentes * alturas * perfiles (Cross Spectra)
22 22 canales * alturas (DC Channels)
23 23
24 24 Esta clase contiene instancias (objetos) de las clases BasicHeader, SystemHeader,
25 25 RadarControllerHeader y Spectra. Los tres primeros se usan para almacenar informacion de la
26 26 cabecera de datos (metadata), y el cuarto (Spectra) para obtener y almacenar un bloque de
27 27 datos desde el "buffer" cada vez que se ejecute el metodo "getData".
28 28
29 29 Example:
30 30 dpath = "/home/myuser/data"
31 31
32 32 startTime = datetime.datetime(2010,1,20,0,0,0,0,0,0)
33 33
34 34 endTime = datetime.datetime(2010,1,21,23,59,59,0,0,0)
35 35
36 36 readerObj = SpectraReader()
37 37
38 38 readerObj.setup(dpath, startTime, endTime)
39 39
40 40 while(True):
41 41
42 42 readerObj.getData()
43 43
44 44 print readerObj.data_spc
45 45
46 46 print readerObj.data_cspc
47 47
48 48 print readerObj.data_dc
49 49
50 50 if readerObj.flagNoMoreFiles:
51 51 break
52 52
53 53 """
54 54
55 55 pts2read_SelfSpectra = 0
56 56
57 57 pts2read_CrossSpectra = 0
58 58
59 59 pts2read_DCchannels = 0
60 60
61 61 ext = ".pdata"
62 62
63 63 optchar = "P"
64 64
65 65 dataOut = None
66 66
67 67 nRdChannels = None
68 68
69 69 nRdPairs = None
70 70
71 71 rdPairList = []
72 72
73 73 def __init__(self):#, **kwargs):
74 74 """
75 75 Inicializador de la clase SpectraReader para la lectura de datos de espectros.
76 76
77 77 Inputs:
78 78 dataOut : Objeto de la clase Spectra. Este objeto sera utilizado para
79 79 almacenar un perfil de datos cada vez que se haga un requerimiento
80 80 (getData). El perfil sera obtenido a partir del buffer de datos,
81 81 si el buffer esta vacio se hara un nuevo proceso de lectura de un
82 82 bloque de datos.
83 83 Si este parametro no es pasado se creara uno internamente.
84 84
85 85 Affected:
86 86 self.dataOut
87 87
88 88 Return : None
89 89 """
90 90
91 91 #Eliminar de la base la herencia
92 92 ProcessingUnit.__init__(self)#, **kwargs)
93 93
94 94 # self.isConfig = False
95 95
96 96 self.pts2read_SelfSpectra = 0
97 97
98 98 self.pts2read_CrossSpectra = 0
99 99
100 100 self.pts2read_DCchannels = 0
101 101
102 102 self.datablock = None
103 103
104 104 self.utc = None
105 105
106 106 self.ext = ".pdata"
107 107
108 108 self.optchar = "P"
109 109
110 110 self.basicHeaderObj = BasicHeader(LOCALTIME)
111 111
112 112 self.systemHeaderObj = SystemHeader()
113 113
114 114 self.radarControllerHeaderObj = RadarControllerHeader()
115 115
116 116 self.processingHeaderObj = ProcessingHeader()
117 117
118 118 self.online = 0
119 119
120 120 self.fp = None
121 121
122 122 self.idFile = None
123 123
124 124 self.dtype = None
125 125
126 126 self.fileSizeByHeader = None
127 127
128 128 self.filenameList = []
129 129
130 130 self.filename = None
131 131
132 132 self.fileSize = None
133 133
134 134 self.firstHeaderSize = 0
135 135
136 136 self.basicHeaderSize = 24
137 137
138 138 self.pathList = []
139 139
140 140 self.lastUTTime = 0
141 141
142 142 self.maxTimeStep = 30
143 143
144 144 self.flagNoMoreFiles = 0
145 145
146 146 self.set = 0
147 147
148 148 self.path = None
149 149
150 150 self.delay = 60 #seconds
151 151
152 152 self.nTries = 3 #quantity tries
153 153
154 154 self.nFiles = 3 #number of files for searching
155 155
156 156 self.nReadBlocks = 0
157 157
158 158 self.flagIsNewFile = 1
159 159
160 160 self.__isFirstTimeOnline = 1
161 161
162 162 # self.ippSeconds = 0
163 163
164 164 self.flagDiscontinuousBlock = 0
165 165
166 166 self.flagIsNewBlock = 0
167 167
168 168 self.nTotalBlocks = 0
169 169
170 170 self.blocksize = 0
171 171
172 172 self.dataOut = self.createObjByDefault()
173 173
174 174 self.profileIndex = 1 #Always
175 175
176 176
177 177 def createObjByDefault(self):
178 178
179 179 dataObj = Spectra()
180 180
181 181 return dataObj
182 182
183 183 def __hasNotDataInBuffer(self):
184 184 return 1
185 185
186 186
187 187 def getBlockDimension(self):
188 188 """
189 189 Obtiene la cantidad de puntos a leer por cada bloque de datos
190 190
191 191 Affected:
192 192 self.nRdChannels
193 193 self.nRdPairs
194 194 self.pts2read_SelfSpectra
195 195 self.pts2read_CrossSpectra
196 196 self.pts2read_DCchannels
197 197 self.blocksize
198 198 self.dataOut.nChannels
199 199 self.dataOut.nPairs
200 200
201 201 Return:
202 202 None
203 203 """
204 204 self.nRdChannels = 0
205 205 self.nRdPairs = 0
206 206 self.rdPairList = []
207 207
208 208 for i in range(0, self.processingHeaderObj.totalSpectra*2, 2):
209 209 if self.processingHeaderObj.spectraComb[i] == self.processingHeaderObj.spectraComb[i+1]:
210 210 self.nRdChannels = self.nRdChannels + 1 #par de canales iguales
211 211 else:
212 212 self.nRdPairs = self.nRdPairs + 1 #par de canales diferentes
213 213 self.rdPairList.append((self.processingHeaderObj.spectraComb[i], self.processingHeaderObj.spectraComb[i+1]))
214 214
215 215 pts2read = self.processingHeaderObj.nHeights * self.processingHeaderObj.profilesPerBlock
216 216
217 217 self.pts2read_SelfSpectra = int(self.nRdChannels * pts2read)
218 218 self.blocksize = self.pts2read_SelfSpectra
219 219
220 220 if self.processingHeaderObj.flag_cspc:
221 221 self.pts2read_CrossSpectra = int(self.nRdPairs * pts2read)
222 222 self.blocksize += self.pts2read_CrossSpectra
223 223
224 224 if self.processingHeaderObj.flag_dc:
225 225 self.pts2read_DCchannels = int(self.systemHeaderObj.nChannels * self.processingHeaderObj.nHeights)
226 226 self.blocksize += self.pts2read_DCchannels
227 227
228 228 # self.blocksize = self.pts2read_SelfSpectra + self.pts2read_CrossSpectra + self.pts2read_DCchannels
229 229
230 230
231 231 def readBlock(self):
232 232 """
233 233 Lee el bloque de datos desde la posicion actual del puntero del archivo
234 234 (self.fp) y actualiza todos los parametros relacionados al bloque de datos
235 235 (metadata + data). La data leida es almacenada en el buffer y el contador del buffer
236 236 es seteado a 0
237 237
238 238 Return: None
239 239
240 240 Variables afectadas:
241 241
242 242 self.flagIsNewFile
243 243 self.flagIsNewBlock
244 244 self.nTotalBlocks
245 245 self.data_spc
246 246 self.data_cspc
247 247 self.data_dc
248 248
249 249 Exceptions:
250 250 Si un bloque leido no es un bloque valido
251 251 """
252 252 blockOk_flag = False
253 253 fpointer = self.fp.tell()
254 254
255 255 spc = numpy.fromfile( self.fp, self.dtype[0], self.pts2read_SelfSpectra )
256 256 spc = spc.reshape( (self.nRdChannels, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D
257 257
258 258 if self.processingHeaderObj.flag_cspc:
259 259 cspc = numpy.fromfile( self.fp, self.dtype, self.pts2read_CrossSpectra )
260 260 cspc = cspc.reshape( (self.nRdPairs, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D
261 261
262 262 if self.processingHeaderObj.flag_dc:
263 263 dc = numpy.fromfile( self.fp, self.dtype, self.pts2read_DCchannels ) #int(self.processingHeaderObj.nHeights*self.systemHeaderObj.nChannels) )
264 264 dc = dc.reshape( (self.systemHeaderObj.nChannels, self.processingHeaderObj.nHeights) ) #transforma a un arreglo 2D
265 265
266 266
267 if self.processingHeaderObj.shif_fft:
267 if not self.processingHeaderObj.shif_fft:
268 268 #desplaza a la derecha en el eje 2 determinadas posiciones
269 269 shift = int(self.processingHeaderObj.profilesPerBlock/2)
270 270 spc = numpy.roll( spc, shift , axis=2 )
271 271
272 272 if self.processingHeaderObj.flag_cspc:
273 273 #desplaza a la derecha en el eje 2 determinadas posiciones
274 274 cspc = numpy.roll( cspc, shift, axis=2 )
275 275
276 276 #Dimensions : nChannels, nProfiles, nSamples
277 277 spc = numpy.transpose( spc, (0,2,1) )
278 278 self.data_spc = spc
279 279
280 280 if self.processingHeaderObj.flag_cspc:
281 281 cspc = numpy.transpose( cspc, (0,2,1) )
282 282 self.data_cspc = cspc['real'] + cspc['imag']*1j
283 283 else:
284 284 self.data_cspc = None
285 285
286 286 if self.processingHeaderObj.flag_dc:
287 287 self.data_dc = dc['real'] + dc['imag']*1j
288 288 else:
289 289 self.data_dc = None
290 290
291 291 self.flagIsNewFile = 0
292 292 self.flagIsNewBlock = 1
293 293
294 294 self.nTotalBlocks += 1
295 295 self.nReadBlocks += 1
296 296
297 297 return 1
298 298
299 299 def getFirstHeader(self):
300 300
301 301 self.getBasicHeader()
302 302
303 303 self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
304 304
305 305 self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
306 306
307 307 # self.dataOut.ippSeconds = self.ippSeconds
308 308
309 309 # self.dataOut.timeInterval = self.radarControllerHeaderObj.ippSeconds * self.processingHeaderObj.nCohInt * self.processingHeaderObj.nIncohInt * self.processingHeaderObj.profilesPerBlock
310 310
311 311 self.dataOut.dtype = self.dtype
312 312
313 313 # self.dataOut.nPairs = self.nPairs
314 314
315 315 self.dataOut.pairsList = self.rdPairList
316 316
317 317 self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
318 318
319 319 self.dataOut.nFFTPoints = self.processingHeaderObj.profilesPerBlock
320 320
321 321 self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
322 322
323 323 self.dataOut.nIncohInt = self.processingHeaderObj.nIncohInt
324 324
325 325 xf = self.processingHeaderObj.firstHeight + self.processingHeaderObj.nHeights*self.processingHeaderObj.deltaHeight
326 326
327 327 self.dataOut.heightList = numpy.arange(self.processingHeaderObj.firstHeight, xf, self.processingHeaderObj.deltaHeight)
328 328
329 329 self.dataOut.channelList = list(range(self.systemHeaderObj.nChannels))
330 330
331 331 self.dataOut.flagShiftFFT = True #Data is always shifted
332 332
333 333 self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode #asumo q la data no esta decodificada
334 334
335 335 self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip #asumo q la data esta sin flip
336 336
337 337 def getData(self):
338 338 """
339 339 First method to execute before "RUN" is called.
340 340
341 341 Copia el buffer de lectura a la clase "Spectra",
342 342 con todos los parametros asociados a este (metadata). cuando no hay datos en el buffer de
343 343 lectura es necesario hacer una nueva lectura de los bloques de datos usando "readNextBlock"
344 344
345 345 Return:
346 346 0 : Si no hay mas archivos disponibles
347 347 1 : Si hizo una buena copia del buffer
348 348
349 349 Affected:
350 350 self.dataOut
351 351
352 352 self.flagDiscontinuousBlock
353 353 self.flagIsNewBlock
354 354 """
355 355
356 356 if self.flagNoMoreFiles:
357 357 self.dataOut.flagNoData = True
358 358 print('Process finished')
359 359 return 0
360 360
361 361 self.flagDiscontinuousBlock = 0
362 362 self.flagIsNewBlock = 0
363 363
364 364 if self.__hasNotDataInBuffer():
365 365
366 366 if not( self.readNextBlock() ):
367 367 self.dataOut.flagNoData = True
368 368 return 0
369 369
370 370 #data es un numpy array de 3 dmensiones (perfiles, alturas y canales)
371 371
372 372 if self.data_spc is None:
373 373 self.dataOut.flagNoData = True
374 374 return 0
375 375
376 376 self.getBasicHeader()
377 377
378 378 self.getFirstHeader()
379 379
380 380 self.dataOut.data_spc = self.data_spc
381 381
382 382 self.dataOut.data_cspc = self.data_cspc
383 383
384 384 self.dataOut.data_dc = self.data_dc
385 385
386 386 self.dataOut.flagNoData = False
387 387
388 388 self.dataOut.realtime = self.online
389 389
390 390 return self.dataOut.data_spc
391 391
392 392 class SpectraWriter(JRODataWriter, Operation):
393 393
394 394 """
395 395 Esta clase permite escribir datos de espectros a archivos procesados (.pdata). La escritura
396 396 de los datos siempre se realiza por bloques.
397 397 """
398 398
399 399 ext = ".pdata"
400 400
401 401 optchar = "P"
402 402
403 403 shape_spc_Buffer = None
404 404
405 405 shape_cspc_Buffer = None
406 406
407 407 shape_dc_Buffer = None
408 408
409 409 data_spc = None
410 410
411 411 data_cspc = None
412 412
413 413 data_dc = None
414 414
415 415 # dataOut = None
416 416
417 417 def __init__(self, **kwargs):
418 418 """
419 419 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
420 420
421 421 Affected:
422 422 self.dataOut
423 423 self.basicHeaderObj
424 424 self.systemHeaderObj
425 425 self.radarControllerHeaderObj
426 426 self.processingHeaderObj
427 427
428 428 Return: None
429 429 """
430 430
431 431 Operation.__init__(self, **kwargs)
432 432
433 433 self.isConfig = False
434 434
435 435 self.nTotalBlocks = 0
436 436
437 437 self.data_spc = None
438 438
439 439 self.data_cspc = None
440 440
441 441 self.data_dc = None
442 442
443 443 self.fp = None
444 444
445 445 self.flagIsNewFile = 1
446 446
447 447 self.nTotalBlocks = 0
448 448
449 449 self.flagIsNewBlock = 0
450 450
451 451 self.setFile = None
452 452
453 453 self.dtype = None
454 454
455 455 self.path = None
456 456
457 457 self.noMoreFiles = 0
458 458
459 459 self.filename = None
460 460
461 461 self.basicHeaderObj = BasicHeader(LOCALTIME)
462 462
463 463 self.systemHeaderObj = SystemHeader()
464 464
465 465 self.radarControllerHeaderObj = RadarControllerHeader()
466 466
467 467 self.processingHeaderObj = ProcessingHeader()
468 468
469 469
470 470 def hasAllDataInBuffer(self):
471 471 return 1
472 472
473 473
474 474 def setBlockDimension(self):
475 475 """
476 476 Obtiene las formas dimensionales del los subbloques de datos que componen un bloque
477 477
478 478 Affected:
479 479 self.shape_spc_Buffer
480 480 self.shape_cspc_Buffer
481 481 self.shape_dc_Buffer
482 482
483 483 Return: None
484 484 """
485 485 self.shape_spc_Buffer = (self.dataOut.nChannels,
486 486 self.processingHeaderObj.nHeights,
487 487 self.processingHeaderObj.profilesPerBlock)
488 488
489 489 self.shape_cspc_Buffer = (self.dataOut.nPairs,
490 490 self.processingHeaderObj.nHeights,
491 491 self.processingHeaderObj.profilesPerBlock)
492 492
493 493 self.shape_dc_Buffer = (self.dataOut.nChannels,
494 494 self.processingHeaderObj.nHeights)
495 495
496 496
497 497 def writeBlock(self):
498 498 """
499 499 Escribe el buffer en el file designado
500 500
501 501 Affected:
502 502 self.data_spc
503 503 self.data_cspc
504 504 self.data_dc
505 505 self.flagIsNewFile
506 506 self.flagIsNewBlock
507 507 self.nTotalBlocks
508 508 self.nWriteBlocks
509 509
510 510 Return: None
511 511 """
512 512
513 513 spc = numpy.transpose( self.data_spc, (0,2,1) )
514 if not( self.processingHeaderObj.shif_fft ):
514 if not self.processingHeaderObj.shif_fft:
515 515 spc = numpy.roll( spc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
516 516 data = spc.reshape((-1))
517 517 data = data.astype(self.dtype[0])
518 518 data.tofile(self.fp)
519 519
520 520 if self.data_cspc is not None:
521 521 data = numpy.zeros( self.shape_cspc_Buffer, self.dtype )
522 522 cspc = numpy.transpose( self.data_cspc, (0,2,1) )
523 if not( self.processingHeaderObj.shif_fft ):
523 if not self.processingHeaderObj.shif_fft:
524 524 cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
525 525 data['real'] = cspc.real
526 526 data['imag'] = cspc.imag
527 527 data = data.reshape((-1))
528 528 data.tofile(self.fp)
529 529
530 530 if self.data_dc is not None:
531 531 data = numpy.zeros( self.shape_dc_Buffer, self.dtype )
532 532 dc = self.data_dc
533 533 data['real'] = dc.real
534 534 data['imag'] = dc.imag
535 535 data = data.reshape((-1))
536 536 data.tofile(self.fp)
537 537
538 538 # self.data_spc.fill(0)
539 539 #
540 540 # if self.data_dc is not None:
541 541 # self.data_dc.fill(0)
542 542 #
543 543 # if self.data_cspc is not None:
544 544 # self.data_cspc.fill(0)
545 545
546 546 self.flagIsNewFile = 0
547 547 self.flagIsNewBlock = 1
548 548 self.nTotalBlocks += 1
549 549 self.nWriteBlocks += 1
550 550 self.blockIndex += 1
551 551
552 552 # print "[Writing] Block = %d04" %self.blockIndex
553 553
554 554 def putData(self):
555 555 """
556 556 Setea un bloque de datos y luego los escribe en un file
557 557
558 558 Affected:
559 559 self.data_spc
560 560 self.data_cspc
561 561 self.data_dc
562 562
563 563 Return:
564 564 0 : Si no hay data o no hay mas files que puedan escribirse
565 565 1 : Si se escribio la data de un bloque en un file
566 566 """
567 567
568 568 if self.dataOut.flagNoData:
569 569 return 0
570 570
571 571 self.flagIsNewBlock = 0
572 572
573 573 if self.dataOut.flagDiscontinuousBlock:
574 574 self.data_spc.fill(0)
575 575 if self.dataOut.data_cspc is not None:
576 576 self.data_cspc.fill(0)
577 577 if self.dataOut.data_dc is not None:
578 578 self.data_dc.fill(0)
579 579 self.setNextFile()
580 580
581 581 if self.flagIsNewFile == 0:
582 582 self.setBasicHeader()
583 583
584 584 self.data_spc = self.dataOut.data_spc.copy()
585 585
586 586 if self.dataOut.data_cspc is not None:
587 587 self.data_cspc = self.dataOut.data_cspc.copy()
588 588
589 589 if self.dataOut.data_dc is not None:
590 590 self.data_dc = self.dataOut.data_dc.copy()
591 591
592 592 # #self.processingHeaderObj.dataBlocksPerFile)
593 593 if self.hasAllDataInBuffer():
594 594 # self.setFirstHeader()
595 595 self.writeNextBlock()
596 596
597 597 return 1
598 598
599 599 def __getBlockSize(self):
600 600 '''
601 601 Este metodos determina el cantidad de bytes para un bloque de datos de tipo Spectra
602 602 '''
603 603
604 604 dtype_width = self.getDtypeWidth()
605 605
606 606 pts2write = self.dataOut.nHeights * self.dataOut.nFFTPoints
607 607
608 608 pts2write_SelfSpectra = int(self.dataOut.nChannels * pts2write)
609 609 blocksize = (pts2write_SelfSpectra*dtype_width)
610 610
611 611 if self.dataOut.data_cspc is not None:
612 612 pts2write_CrossSpectra = int(self.dataOut.nPairs * pts2write)
613 613 blocksize += (pts2write_CrossSpectra*dtype_width*2)
614 614
615 615 if self.dataOut.data_dc is not None:
616 616 pts2write_DCchannels = int(self.dataOut.nChannels * self.dataOut.nHeights)
617 617 blocksize += (pts2write_DCchannels*dtype_width*2)
618 618
619 619 # blocksize = blocksize #* datatypeValue * 2 #CORREGIR ESTO
620 620
621 621 return blocksize
622 622
623 623 def setFirstHeader(self):
624 624
625 625 """
626 626 Obtiene una copia del First Header
627 627
628 628 Affected:
629 629 self.systemHeaderObj
630 630 self.radarControllerHeaderObj
631 631 self.dtype
632 632
633 633 Return:
634 634 None
635 635 """
636 636
637 637 self.systemHeaderObj = self.dataOut.systemHeaderObj.copy()
638 638 self.systemHeaderObj.nChannels = self.dataOut.nChannels
639 639 self.radarControllerHeaderObj = self.dataOut.radarControllerHeaderObj.copy()
640 640
641 641 self.processingHeaderObj.dtype = 1 # Spectra
642 642 self.processingHeaderObj.blockSize = self.__getBlockSize()
643 643 self.processingHeaderObj.profilesPerBlock = self.dataOut.nFFTPoints
644 644 self.processingHeaderObj.dataBlocksPerFile = self.blocksPerFile
645 645 self.processingHeaderObj.nWindows = 1 #podria ser 1 o self.dataOut.processingHeaderObj.nWindows
646 646 self.processingHeaderObj.nCohInt = self.dataOut.nCohInt# Se requiere para determinar el valor de timeInterval
647 647 self.processingHeaderObj.nIncohInt = self.dataOut.nIncohInt
648 648 self.processingHeaderObj.totalSpectra = self.dataOut.nPairs + self.dataOut.nChannels
649 649 self.processingHeaderObj.shif_fft = self.dataOut.flagShiftFFT
650 650
651 651 if self.processingHeaderObj.totalSpectra > 0:
652 652 channelList = []
653 653 for channel in range(self.dataOut.nChannels):
654 654 channelList.append(channel)
655 655 channelList.append(channel)
656 656
657 657 pairsList = []
658 658 if self.dataOut.nPairs > 0:
659 659 for pair in self.dataOut.pairsList:
660 660 pairsList.append(pair[0])
661 661 pairsList.append(pair[1])
662 662
663 663 spectraComb = channelList + pairsList
664 664 spectraComb = numpy.array(spectraComb, dtype="u1")
665 665 self.processingHeaderObj.spectraComb = spectraComb
666 666
667 667 if self.dataOut.code is not None:
668 668 self.processingHeaderObj.code = self.dataOut.code
669 669 self.processingHeaderObj.nCode = self.dataOut.nCode
670 670 self.processingHeaderObj.nBaud = self.dataOut.nBaud
671 671
672 672 if self.processingHeaderObj.nWindows != 0:
673 673 self.processingHeaderObj.firstHeight = self.dataOut.heightList[0]
674 674 self.processingHeaderObj.deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
675 675 self.processingHeaderObj.nHeights = self.dataOut.nHeights
676 676 self.processingHeaderObj.samplesWin = self.dataOut.nHeights
677 677
678 678 self.processingHeaderObj.processFlags = self.getProcessFlags()
679 679
680 680 self.setBasicHeader() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,139 +1,139
1 1 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
2 2 #define NUM_CPY_THREADS 8
3 3 #include <Python.h>
4 4 #include <numpy/arrayobject.h>
5 5 #include <math.h>
6 6 #include <complex.h>
7 7 #include <time.h>
8 8
9 9 // void printArr(int *array);
10 10 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args);
11 11 static PyObject *correlateByBlock(PyObject *self, PyObject *args);
12 12 #ifndef PyMODINIT_FUNC /* declarations for DLL import/export */
13 13 #define PyMODINIT_FUNC void
14 14 #endif
15 15
16 16 static PyMethodDef extensionsMethods[] = {
17 17 { "correlateByBlock", (PyCFunction)correlateByBlock, METH_VARARGS, "get correlation by block" },
18 18 { "hildebrand_sekhon", (PyCFunction)hildebrand_sekhon, METH_VARARGS, "get noise with hildebrand_sekhon" },
19 19 { NULL, NULL, 0, NULL }
20 20 };
21 21
22 22 PyMODINIT_FUNC initcSchain() {
23 23 Py_InitModule("cSchain", extensionsMethods);
24 24 import_array();
25 25 }
26 26
27 27 static PyObject *correlateByBlock(PyObject *self, PyObject *args) {
28 28
29 29 // int *x = (int*) malloc(4000000 * 216 * sizeof(int));;
30 30 // int a = 5;
31 31 // x = &a;
32 32 // int b = 6;
33 33 // x = &b;
34 34 // printf("Antes de imprimir x \n");
35 35 // printf("%d \n", x[0]);
36 36
37 37 PyObject *data_obj1, *data_obj2;
38 38 PyArrayObject *data_array1, *data_array2, *correlateRow, *out, *dataRow, *codeRow; //, ,
39 39 int mode;
40 40
41 41 if (!PyArg_ParseTuple(args, "OOi", &data_obj1, &data_obj2, &mode)) return NULL;
42 42
43 data_array1 = (PyArrayObject *) PyArray_FROM_OTF(data_obj1, NPY_COMPLEX128, NPY_ARRAY_DEFAULT);
44 data_array2 = (PyArrayObject *) PyArray_FROM_OTF(data_obj2, NPY_FLOAT64, NPY_ARRAY_DEFAULT);
43 data_array1 = (PyArrayObject *) PyArray_FROM_OTF(data_obj1, NPY_COMPLEX128, NPY_ARRAY_IN_ARRAY);
44 data_array2 = (PyArrayObject *) PyArray_FROM_OTF(data_obj2, NPY_FLOAT64, NPY_ARRAY_IN_ARRAY);
45 45
46 46 npy_intp dims[1];
47 47 dims[0] = 200;
48 48 npy_intp dims_code[1];
49 49 dims_code[0] = 16;
50 50
51 51 double complex * dataRaw;
52 52 double * codeRaw;
53 53 dataRaw = (double complex*) PyArray_DATA(data_array1);
54 54 codeRaw = (double *) PyArray_DATA(data_array2);
55 55 double complex ** outC = malloc(40000*200*sizeof(double complex));
56 56 int i;
57 57
58 58 clock_t start = clock();
59 59 for(i=0; i<40000; i++){
60 60 // codeRow = PyArray_SimpleNewFromData(1, dims_code, NPY_FLOAT64, codeRaw + 16 * i);
61 61 // dataRow = PyArray_SimpleNewFromData(1, dims, NPY_COMPLEX128, dataRaw + 200 * i);
62 62 // Py_INCREF(codeRow);
63 63 // Py_INCREF(dataRow);
64 64 // PyArray_ENABLEFLAGS(codeRow, NPY_ARRAY_OWNDATA);
65 65 // PyArray_ENABLEFLAGS(dataRow, NPY_ARRAY_OWNDATA);
66 66 correlateRow = (PyArrayObject *) PyArray_Correlate2(PyArray_SimpleNewFromData(1, dims_code, NPY_FLOAT64, codeRaw + 16 * i), PyArray_SimpleNewFromData(1, dims, NPY_COMPLEX128, dataRaw + 200 * i), (npy_intp) 2);
67 67 //Py_INCREF(correlateRow);
68 68 // PyArray_ENABLEFLAGS(correlateRow, NPY_ARRAY_OWNDATA);
69 69 memcpy(outC + 200*i, (double complex*) PyArray_DATA(correlateRow), 200 * sizeof(double complex));
70 70
71 71 Py_DECREF(correlateRow);
72 72 // Py_DECREF(codeRow);
73 73 // Py_DECREF(dataRow);
74 74 }
75 75 clock_t end = clock();
76 76 float seconds = (float)(end - start) / CLOCKS_PER_SEC;
77 77 printf("%f", seconds);
78 78 //
79 79 npy_intp dimsret[2];
80 80 dimsret[0] = 40000;
81 81 dimsret[1] = 200;
82 82 out = PyArray_SimpleNewFromData(2, dimsret, NPY_COMPLEX128, outC);
83 83 PyArray_ENABLEFLAGS(out, NPY_ARRAY_OWNDATA);
84 84 //Py_INCREF(out);
85 85 Py_DECREF(data_array1);
86 86 Py_DECREF(data_array2);
87 87 // PyArray_DebugPrint(out);
88 88 // Py_DECREF(data_obj2);
89 89 // Py_DECREF(data_obj1);
90 90 // Py_DECREF(codeRow);
91 91 // Py_DECREF(dataRow);
92 92 // free(dataRaw);
93 93 // free(codeRaw);
94 94
95 95 return PyArray_Return(out);
96 96 }
97 97
98 98 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
99 99 double navg;
100 100 PyObject *data_obj, *data_array;
101 101
102 102 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) return NULL;
103 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_ARRAY_DEFAULT);
103 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_ARRAY_IN_ARRAY);
104 104 if (data_array == NULL) {
105 105 Py_XDECREF(data_array);
106 106 Py_XDECREF(data_obj);
107 107 return NULL;
108 108 }
109 109 double *sortdata = (double*)PyArray_DATA(data_array);
110 110 int lenOfData = (int)PyArray_SIZE(data_array) ;
111 111 double nums_min = lenOfData*0.2;
112 112 if (nums_min <= 5) nums_min = 5;
113 113 double sump = 0;
114 114 double sumq = 0;
115 115 int j = 0;
116 116 int cont = 1;
117 117 double rtest = 0;
118 118 while ((cont == 1) && (j < lenOfData)) {
119 119 sump = sump + sortdata[j];
120 120 sumq = sumq + pow(sortdata[j], 2);
121 121 if (j > nums_min) {
122 122 rtest = (double)j/(j-1) + 1/navg;
123 123 if ((sumq*j) > (rtest*pow(sump, 2))) {
124 124 j = j - 1;
125 125 sump = sump - sortdata[j];
126 126 sumq = sumq - pow(sortdata[j],2);
127 127 cont = 0;
128 128 }
129 129 }
130 130 j = j + 1;
131 131 }
132 132
133 133 double lnoise = sump / j;
134 134
135 135 Py_DECREF(data_array);
136 136
137 137 return Py_BuildValue("d", lnoise);
138 138 }
139 139
General Comments 0
You need to be logged in to leave comments. Login now