##// END OF EJS Templates
schain
RSS
sin matplotlib de los modulos de claire
José Chávez -
r1075:eb67833bff79
parent child
Show More
Side by Side Unified
Context file:
r1075:eb67833bff79 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,12 +1,12
1 1 #from schainpy.model.data.jrodata import *
2 2 # from schainpy.model.io.jrodataIO import *
3 3 # from schainpy.model.proc.jroprocessing import *
4 4 # from schainpy.model.graphics.jroplot import *
5 5 # from schainpy.model.utils.jroutils import *
6 6 # from schainpy.serializer import *
7 7
8 from graphics import *
8 9 from data import *
9 10 from io import *
10 11 from proc import *
11 from graphics import *
12 from utils import * No newline at end of file
12 from utils import *
Show file before | Show file after | Hide comments
@@ -1,468 +1,495
1 1 import numpy
2 2 import datetime
3 3 import sys
4 4 import matplotlib
5 5
6 6 if 'linux' in sys.platform:
7 matplotlib.use("TKAgg")
7 matplotlib.use("TkAgg")
8 8
9 9 if 'darwin' in sys.platform:
10 10 matplotlib.use('TKAgg')
11 #Qt4Agg', 'GTK', 'GTKAgg', 'ps', 'agg', 'cairo', 'MacOSX', 'GTKCairo', 'WXAgg', 'template', 'TkAgg', 'GTK3Cairo', 'GTK3Agg', 'svg', 'WebAgg', 'CocoaAgg', 'emf', 'gdk', 'WX'
11 # Qt4Agg', 'GTK', 'GTKAgg', 'ps', 'agg', 'cairo', 'MacOSX', 'GTKCairo', 'WXAgg', 'template', 'TkAgg', 'GTK3Cairo', 'GTK3Agg', 'svg', 'WebAgg', 'CocoaAgg', 'emf', 'gdk', 'WX'
12 12 import matplotlib.pyplot
13 13
14 14 from mpl_toolkits.axes_grid1 import make_axes_locatable
15 15 from matplotlib.ticker import FuncFormatter, LinearLocator
16 16
17 17 ###########################################
18 #Actualizacion de las funciones del driver
18 # Actualizacion de las funciones del driver
19 19 ###########################################
20 20
21 21 # create jro colormap
22 22
23 23 jet_values = matplotlib.pyplot.get_cmap("jet", 100)(numpy.arange(100))[10:90]
24 blu_values = matplotlib.pyplot.get_cmap("seismic_r", 20)(numpy.arange(20))[10:15]
25 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list("jro", numpy.vstack((blu_values, jet_values)))
24 blu_values = matplotlib.pyplot.get_cmap(
25 "seismic_r", 20)(numpy.arange(20))[10:15]
26 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
27 "jro", numpy.vstack((blu_values, jet_values)))
26 28 matplotlib.pyplot.register_cmap(cmap=ncmap)
27 29
28 def createFigure(id, wintitle, width, height, facecolor="w", show=True, dpi = 80):
30
31 def createFigure(id, wintitle, width, height, facecolor="w", show=True, dpi=80):
29 32
30 33 matplotlib.pyplot.ioff()
31 34
32 fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor, figsize=(1.0*width/dpi, 1.0*height/dpi))
35 fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor, figsize=(
36 1.0 * width / dpi, 1.0 * height / dpi))
33 37 fig.canvas.manager.set_window_title(wintitle)
34 38 # fig.canvas.manager.resize(width, height)
35 39 matplotlib.pyplot.ion()
36 40
37 41 if show:
38 42 matplotlib.pyplot.show()
39 43
40 44 return fig
41 45
46
42 47 def closeFigure(show=False, fig=None):
43 48
44 # matplotlib.pyplot.ioff()
45 # matplotlib.pyplot.pause(0)
49 # matplotlib.pyplot.ioff()
50 # matplotlib.pyplot.pause(0)
46 51
47 52 if show:
48 53 matplotlib.pyplot.show()
49 54
50 55 if fig != None:
51 56 matplotlib.pyplot.close(fig)
52 57 # matplotlib.pyplot.pause(0)
53 58 # matplotlib.pyplot.ion()
54 59
55 60 return
56 61
57 62 matplotlib.pyplot.close("all")
58 63 # matplotlib.pyplot.pause(0)
59 64 # matplotlib.pyplot.ion()
60 65
61 66 return
62 67
68
63 69 def saveFigure(fig, filename):
64 70
65 # matplotlib.pyplot.ioff()
71 # matplotlib.pyplot.ioff()
66 72 fig.savefig(filename, dpi=matplotlib.pyplot.gcf().dpi)
67 73 # matplotlib.pyplot.ion()
68 74
75
69 76 def clearFigure(fig):
70 77
71 78 fig.clf()
72 79
80
73 81 def setWinTitle(fig, title):
74 82
75 83 fig.canvas.manager.set_window_title(title)
76 84
85
77 86 def setTitle(fig, title):
78 87
79 88 fig.suptitle(title)
80 89
90
81 91 def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
82 92
83 93 matplotlib.pyplot.ioff()
84 94 matplotlib.pyplot.figure(fig.number)
85 95 axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
86 (xpos, ypos),
87 colspan=colspan,
88 rowspan=rowspan,
89 polar=polar)
96 (xpos, ypos),
97 colspan=colspan,
98 rowspan=rowspan,
99 polar=polar)
90 100
91 101 matplotlib.pyplot.ion()
92 102 return axes
93 103
104
94 105 def setAxesText(ax, text):
95 106
96 107 ax.annotate(text,
97 xy = (.1, .99),
98 xycoords = 'figure fraction',
99 horizontalalignment = 'left',
100 verticalalignment = 'top',
101 fontsize = 10)
108 xy=(.1, .99),
109 xycoords='figure fraction',
110 horizontalalignment='left',
111 verticalalignment='top',
112 fontsize=10)
113
102 114
103 115 def printLabels(ax, xlabel, ylabel, title):
104 116
105 117 ax.set_xlabel(xlabel, size=11)
106 118 ax.set_ylabel(ylabel, size=11)
107 119 ax.set_title(title, size=8)
108 120
121
109 122 def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='',
110 123 ticksize=9, xtick_visible=True, ytick_visible=True,
111 124 nxticks=4, nyticks=10,
112 grid=None,color='blue'):
113
125 grid=None, color='blue'):
114 126 """
115 127
116 128 Input:
117 129 grid : None, 'both', 'x', 'y'
118 130 """
119 131
120 132 matplotlib.pyplot.ioff()
121 133
122 ax.set_xlim([xmin,xmax])
123 ax.set_ylim([ymin,ymax])
134 ax.set_xlim([xmin, xmax])
135 ax.set_ylim([ymin, ymax])
124 136
125 137 printLabels(ax, xlabel, ylabel, title)
126 138
127 139 ######################################################
128 if (xmax-xmin)<=1:
129 xtickspos = numpy.linspace(xmin,xmax,nxticks)
130 xtickspos = numpy.array([float("%.1f"%i) for i in xtickspos])
140 if (xmax - xmin) <= 1:
141 xtickspos = numpy.linspace(xmin, xmax, nxticks)
142 xtickspos = numpy.array([float("%.1f" % i) for i in xtickspos])
131 143 ax.set_xticks(xtickspos)
132 144 else:
133 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
145 xtickspos = numpy.arange(nxticks) * \
146 int((xmax - xmin) / (nxticks)) + int(xmin)
134 147 # xtickspos = numpy.arange(nxticks)*float(xmax-xmin)/float(nxticks) + int(xmin)
135 148 ax.set_xticks(xtickspos)
136 149
137 150 for tick in ax.get_xticklabels():
138 151 tick.set_visible(xtick_visible)
139 152
140 153 for tick in ax.xaxis.get_major_ticks():
141 154 tick.label.set_fontsize(ticksize)
142 155
143 156 ######################################################
144 157 for tick in ax.get_yticklabels():
145 158 tick.set_visible(ytick_visible)
146 159
147 160 for tick in ax.yaxis.get_major_ticks():
148 161 tick.label.set_fontsize(ticksize)
149 162
150 163 ax.plot(x, y, color=color)
151 164 iplot = ax.lines[-1]
152 165
153 166 ######################################################
154 167 if '0.' in matplotlib.__version__[0:2]:
155 168 print "The matplotlib version has to be updated to 1.1 or newer"
156 169 return iplot
157 170
158 171 if '1.0.' in matplotlib.__version__[0:4]:
159 172 print "The matplotlib version has to be updated to 1.1 or newer"
160 173 return iplot
161 174
162 175 if grid != None:
163 176 ax.grid(b=True, which='major', axis=grid)
164 177
165 178 matplotlib.pyplot.tight_layout()
166 179
167 180 matplotlib.pyplot.ion()
168 181
169 182 return iplot
170 183
184
171 185 def set_linedata(ax, x, y, idline):
172 186
173 ax.lines[idline].set_data(x,y)
187 ax.lines[idline].set_data(x, y)
188
174 189
175 190 def pline(iplot, x, y, xlabel='', ylabel='', title=''):
176 191
177 192 ax = iplot.axes
178 193
179 194 printLabels(ax, xlabel, ylabel, title)
180 195
181 196 set_linedata(ax, x, y, idline=0)
182 197
198
183 199 def addpline(ax, x, y, color, linestyle, lw):
184 200
185 ax.plot(x,y,color=color,linestyle=linestyle,lw=lw)
201 ax.plot(x, y, color=color, linestyle=linestyle, lw=lw)
186 202
187 203
188 204 def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax,
189 xlabel='', ylabel='', title='', ticksize = 9,
190 colormap='jet',cblabel='', cbsize="5%",
205 xlabel='', ylabel='', title='', ticksize=9,
206 colormap='jet', cblabel='', cbsize="5%",
191 207 XAxisAsTime=False):
192 208
193 209 matplotlib.pyplot.ioff()
194 210
195 211 divider = make_axes_locatable(ax)
196 212 ax_cb = divider.new_horizontal(size=cbsize, pad=0.05)
197 213 fig = ax.get_figure()
198 214 fig.add_axes(ax_cb)
199 215
200 ax.set_xlim([xmin,xmax])
201 ax.set_ylim([ymin,ymax])
216 ax.set_xlim([xmin, xmax])
217 ax.set_ylim([ymin, ymax])
202 218
203 219 printLabels(ax, xlabel, ylabel, title)
204 220
205 221 z = numpy.ma.masked_invalid(z)
206 cmap=matplotlib.pyplot.get_cmap(colormap)
222 cmap = matplotlib.pyplot.get_cmap(colormap)
207 223 cmap.set_bad('black', 1.)
208 imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=cmap)
209 cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
224 imesh = ax.pcolormesh(x, y, z.T, vmin=zmin, vmax=zmax, cmap=cmap)
225 cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
210 226 cb.set_label(cblabel)
211 227
212 228 # for tl in ax_cb.get_yticklabels():
213 229 # tl.set_visible(True)
214 230
215 231 for tick in ax.yaxis.get_major_ticks():
216 232 tick.label.set_fontsize(ticksize)
217 233
218 234 for tick in ax.xaxis.get_major_ticks():
219 235 tick.label.set_fontsize(ticksize)
220 236
221 237 for tick in cb.ax.get_yticklabels():
222 238 tick.set_fontsize(ticksize)
223 239
224 240 ax_cb.yaxis.tick_right()
225 241
226 242 if '0.' in matplotlib.__version__[0:2]:
227 243 print "The matplotlib version has to be updated to 1.1 or newer"
228 244 return imesh
229 245
230 246 if '1.0.' in matplotlib.__version__[0:4]:
231 247 print "The matplotlib version has to be updated to 1.1 or newer"
232 248 return imesh
233 249
234 250 matplotlib.pyplot.tight_layout()
235 251
236 252 if XAxisAsTime:
237 253
238 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
254 def func(x, pos): return ('%s') % (
255 datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
239 256 ax.xaxis.set_major_formatter(FuncFormatter(func))
240 257 ax.xaxis.set_major_locator(LinearLocator(7))
241 258
242 259 matplotlib.pyplot.ion()
243 260 return imesh
244 261
262
245 263 def pcolor(imesh, z, xlabel='', ylabel='', title=''):
246 264
247 265 z = z.T
248 266 ax = imesh.axes
249 267 printLabels(ax, xlabel, ylabel, title)
250 268 imesh.set_array(z.ravel())
251 269
270
252 271 def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
253 272
254 273 printLabels(ax, xlabel, ylabel, title)
255 274
256 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
275 ax.pcolormesh(x, y, z.T, vmin=zmin, vmax=zmax,
276 cmap=matplotlib.pyplot.get_cmap(colormap))
277
257 278
258 279 def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
259 280
260 281 printLabels(ax, xlabel, ylabel, title)
261 282
262 283 ax.collections.remove(ax.collections[0])
263 284
264 285 z = numpy.ma.masked_invalid(z)
265 286
266 cmap=matplotlib.pyplot.get_cmap(colormap)
287 cmap = matplotlib.pyplot.get_cmap(colormap)
267 288 cmap.set_bad('black', 1.)
268 289
290 ax.pcolormesh(x, y, z.T, vmin=zmin, vmax=zmax, cmap=cmap)
269 291
270 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=cmap)
271 292
272 293 def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
273 ticksize=9, xtick_visible=True, ytick_visible=True,
274 nxticks=4, nyticks=10,
275 grid=None):
276
294 ticksize=9, xtick_visible=True, ytick_visible=True,
295 nxticks=4, nyticks=10,
296 grid=None):
277 297 """
278 298
279 299 Input:
280 300 grid : None, 'both', 'x', 'y'
281 301 """
282 302
283 303 matplotlib.pyplot.ioff()
284 304
285 305 lines = ax.plot(x.T, y)
286 306 leg = ax.legend(lines, legendlabels, loc='upper right')
287 307 leg.get_frame().set_alpha(0.5)
288 ax.set_xlim([xmin,xmax])
289 ax.set_ylim([ymin,ymax])
308 ax.set_xlim([xmin, xmax])
309 ax.set_ylim([ymin, ymax])
290 310 printLabels(ax, xlabel, ylabel, title)
291 311
292 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
312 xtickspos = numpy.arange(nxticks) * \
313 int((xmax - xmin) / (nxticks)) + int(xmin)
293 314 ax.set_xticks(xtickspos)
294 315
295 316 for tick in ax.get_xticklabels():
296 317 tick.set_visible(xtick_visible)
297 318
298 319 for tick in ax.xaxis.get_major_ticks():
299 320 tick.label.set_fontsize(ticksize)
300 321
301 322 for tick in ax.get_yticklabels():
302 323 tick.set_visible(ytick_visible)
303 324
304 325 for tick in ax.yaxis.get_major_ticks():
305 326 tick.label.set_fontsize(ticksize)
306 327
307 328 iplot = ax.lines[-1]
308 329
309 330 if '0.' in matplotlib.__version__[0:2]:
310 331 print "The matplotlib version has to be updated to 1.1 or newer"
311 332 return iplot
312 333
313 334 if '1.0.' in matplotlib.__version__[0:4]:
314 335 print "The matplotlib version has to be updated to 1.1 or newer"
315 336 return iplot
316 337
317 338 if grid != None:
318 339 ax.grid(b=True, which='major', axis=grid)
319 340
320 341 matplotlib.pyplot.tight_layout()
321 342
322 343 matplotlib.pyplot.ion()
323 344
324 345 return iplot
325 346
326 347
327 348 def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''):
328 349
329 350 ax = iplot.axes
330 351
331 352 printLabels(ax, xlabel, ylabel, title)
332 353
333 354 for i in range(len(ax.lines)):
334 355 line = ax.lines[i]
335 line.set_data(x[i,:],y)
356 line.set_data(x[i, :], y)
336 357
337 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
338 ticksize=9, xtick_visible=True, ytick_visible=True,
339 nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
340 grid=None, XAxisAsTime=False):
341 358
359 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
360 ticksize=9, xtick_visible=True, ytick_visible=True,
361 nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
362 grid=None, XAxisAsTime=False):
342 363 """
343 364
344 365 Input:
345 366 grid : None, 'both', 'x', 'y'
346 367 """
347 368
348 369 matplotlib.pyplot.ioff()
349 370
350 371 # lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
351 372 lines = ax.plot(x, y.T)
352 373 # leg = ax.legend(lines, legendlabels, loc=2, bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
353 374 # handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)
354 375
355 376 leg = ax.legend(lines, legendlabels,
356 377 loc='upper right', bbox_to_anchor=(1.16, 1), borderaxespad=0)
357 378
358 for label in leg.get_texts(): label.set_fontsize(9)
379 for label in leg.get_texts():
380 label.set_fontsize(9)
359 381
360 ax.set_xlim([xmin,xmax])
361 ax.set_ylim([ymin,ymax])
382 ax.set_xlim([xmin, xmax])
383 ax.set_ylim([ymin, ymax])
362 384 printLabels(ax, xlabel, ylabel, title)
363 385
364 386 # xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
365 387 # ax.set_xticks(xtickspos)
366 388
367 389 for tick in ax.get_xticklabels():
368 390 tick.set_visible(xtick_visible)
369 391
370 392 for tick in ax.xaxis.get_major_ticks():
371 393 tick.label.set_fontsize(ticksize)
372 394
373 395 for tick in ax.get_yticklabels():
374 396 tick.set_visible(ytick_visible)
375 397
376 398 for tick in ax.yaxis.get_major_ticks():
377 399 tick.label.set_fontsize(ticksize)
378 400
379 401 iplot = ax.lines[-1]
380 402
381 403 if '0.' in matplotlib.__version__[0:2]:
382 404 print "The matplotlib version has to be updated to 1.1 or newer"
383 405 return iplot
384 406
385 407 if '1.0.' in matplotlib.__version__[0:4]:
386 408 print "The matplotlib version has to be updated to 1.1 or newer"
387 409 return iplot
388 410
389 411 if grid != None:
390 412 ax.grid(b=True, which='major', axis=grid)
391 413
392 414 matplotlib.pyplot.tight_layout()
393 415
394 416 if XAxisAsTime:
395 417
396 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
418 def func(x, pos): return ('%s') % (
419 datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
397 420 ax.xaxis.set_major_formatter(FuncFormatter(func))
398 421 ax.xaxis.set_major_locator(LinearLocator(7))
399 422
400 423 matplotlib.pyplot.ion()
401 424
402 425 return iplot
403 426
427
404 428 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
405 429
406 430 ax = iplot.axes
407 431
408 432 printLabels(ax, xlabel, ylabel, title)
409 433
410 434 for i in range(len(ax.lines)):
411 435 line = ax.lines[i]
412 line.set_data(x,y[i,:])
436 line.set_data(x, y[i, :])
437
413 438
414 439 def createPolar(ax, x, y,
415 xlabel='', ylabel='', title='', ticksize = 9,
416 colormap='jet',cblabel='', cbsize="5%",
417 XAxisAsTime=False):
440 xlabel='', ylabel='', title='', ticksize=9,
441 colormap='jet', cblabel='', cbsize="5%",
442 XAxisAsTime=False):
418 443
419 444 matplotlib.pyplot.ioff()
420 445
421 ax.plot(x,y,'bo', markersize=5)
446 ax.plot(x, y, 'bo', markersize=5)
422 447 # ax.set_rmax(90)
423 ax.set_ylim(0,90)
424 ax.set_yticks(numpy.arange(0,90,20))
448 ax.set_ylim(0, 90)
449 ax.set_yticks(numpy.arange(0, 90, 20))
425 450 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
426 451 # ax.text(0, 50, ylabel, rotation='vertical', va ='center', ha = 'left' ,size='11')
427 452 # ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
428 453 ax.yaxis.labelpad = 40
429 454 printLabels(ax, xlabel, ylabel, title)
430 455 iplot = ax.lines[-1]
431 456
432 457 if '0.' in matplotlib.__version__[0:2]:
433 458 print "The matplotlib version has to be updated to 1.1 or newer"
434 459 return iplot
435 460
436 461 if '1.0.' in matplotlib.__version__[0:4]:
437 462 print "The matplotlib version has to be updated to 1.1 or newer"
438 463 return iplot
439 464
440 465 # if grid != None:
441 466 # ax.grid(b=True, which='major', axis=grid)
442 467
443 468 matplotlib.pyplot.tight_layout()
444 469
445 470 matplotlib.pyplot.ion()
446 471
447
448 472 return iplot
449 473
474
450 475 def polar(iplot, x, y, xlabel='', ylabel='', title=''):
451 476
452 477 ax = iplot.axes
453 478
454 479 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
455 480 printLabels(ax, xlabel, ylabel, title)
456 481
457 482 set_linedata(ax, x, y, idline=0)
458 483
484
459 485 def draw(fig):
460 486
461 487 if type(fig) == 'int':
462 488 raise ValueError, "Error drawing: Fig parameter should be a matplotlib figure object figure"
463 489
464 490 fig.canvas.draw()
465 491
492
466 493 def pause(interval=0.000001):
467 494
468 495 matplotlib.pyplot.pause(interval)
Show file before | Show file after | Hide comments
@@ -1,321 +1,331
1 import os, sys
1 import os
2 import sys
2 3 import glob
3 4 import fnmatch
4 5 import datetime
5 6 import time
6 7 import re
7 8 import h5py
8 9 import numpy
9 import matplotlib.pyplot as plt
10 10
11 import pylab as plb
12 11 from scipy.optimize import curve_fit
13 from scipy import asarray as ar,exp
12 from scipy import asarray as ar, exp
14 13 from scipy import stats
15 14
16 15 from duplicity.path import Path
17 16 from numpy.ma.core import getdata
18 17
19 18 SPEED_OF_LIGHT = 299792458
20 19 SPEED_OF_LIGHT = 3e8
21 20
22 21 try:
23 22 from gevent import sleep
24 23 except:
25 24 from time import sleep
26 25
27 26 from schainpy.model.data.jrodata import Spectra
28 27 #from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader
29 28 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
30 29 #from schainpy.model.io.jroIO_bltr import BLTRReader
31 30 from numpy import imag, shape, NaN
32 31
33 32
34 startFp = open('/home/erick/Documents/MIRA35C/20160117/20160117_0000.zspc',"rb")
35
36
37 FILE_HEADER = numpy.dtype([ #HEADER 1024bytes
38 ('Hname',numpy.str_,32), #Original file name
39 ('Htime',numpy.str_,32), #Date and time when the file was created
40 ('Hoper',numpy.str_,64), #Name of operator who created the file
41 ('Hplace',numpy.str_,128), #Place where the measurements was carried out
42 ('Hdescr',numpy.str_,256), #Description of measurements
43 ('Hdummy',numpy.str_,512), #Reserved space
44 #Main chunk
45 ('Msign','<i4'), #Main chunk signature FZKF or NUIG
46 ('MsizeData','<i4'), #Size of data block main chunk
47 #Processing DSP parameters
48 ('PPARsign','<i4'), #PPAR signature
49 ('PPARsize','<i4'), #PPAR size of block
50 ('PPARprf','<i4'), #Pulse repetition frequency
51 ('PPARpdr','<i4'), #Pulse duration
52 ('PPARsft','<i4'), #FFT length
53 ('PPARavc','<i4'), #Number of spectral (in-coherent) averages
54 ('PPARihp','<i4'), #Number of lowest range gate for moment estimation
55 ('PPARchg','<i4'), #Count for gates for moment estimation
56 ('PPARpol','<i4'), #switch on/off polarimetric measurements. Should be 1.
57 #Service DSP parameters
58 ('SPARatt','<i4'), #STC attenuation on the lowest ranges on/off
59 ('SPARtx','<i4'), #OBSOLETE
60 ('SPARaddGain0','<f4'), #OBSOLETE
61 ('SPARaddGain1','<f4'), #OBSOLETE
62 ('SPARwnd','<i4'), #Debug only. It normal mode it is 0.
63 ('SPARpos','<i4'), #Delay between sync pulse and tx pulse for phase corr, ns
64 ('SPARadd','<i4'), #"add to pulse" to compensate for delay between the leading edge of driver pulse and envelope of the RF signal.
65 ('SPARlen','<i4'), #Time for measuring txn pulse phase. OBSOLETE
66 ('SPARcal','<i4'), #OBSOLETE
67 ('SPARnos','<i4'), #OBSOLETE
68 ('SPARof0','<i4'), #detection threshold
69 ('SPARof1','<i4'), #OBSOLETE
70 ('SPARswt','<i4'), #2nd moment estimation threshold
71 ('SPARsum','<i4'), #OBSOLETE
72 ('SPARosc','<i4'), #flag Oscillosgram mode
73 ('SPARtst','<i4'), #OBSOLETE
74 ('SPARcor','<i4'), #OBSOLETE
75 ('SPARofs','<i4'), #OBSOLETE
76 ('SPARhsn','<i4'), #Hildebrand div noise detection on noise gate
77 ('SPARhsa','<f4'), #Hildebrand div noise detection on all gates
78 ('SPARcalibPow_M','<f4'), #OBSOLETE
79 ('SPARcalibSNR_M','<f4'), #OBSOLETE
80 ('SPARcalibPow_S','<f4'), #OBSOLETE
81 ('SPARcalibSNR_S','<f4'), #OBSOLETE
82 ('SPARrawGate1','<i4'), #Lowest range gate for spectra saving Raw_Gate1 >=5
83 ('SPARrawGate2','<i4'), #Number of range gates with atmospheric signal
84 ('SPARraw','<i4'), #flag - IQ or spectra saving on/off
85 ('SPARprc','<i4'),]) #flag - Moment estimation switched on/off
86
87
88
89 self.Hname= None
90 self.Htime= None
91 self.Hoper= None
92 self.Hplace= None
93 self.Hdescr= None
94 self.Hdummy= None
95
96 self.Msign=None
97 self.MsizeData=None
98
99 self.PPARsign=None
100 self.PPARsize=None
101 self.PPARprf=None
102 self.PPARpdr=None
103 self.PPARsft=None
104 self.PPARavc=None
105 self.PPARihp=None
106 self.PPARchg=None
107 self.PPARpol=None
108 #Service DSP parameters
109 self.SPARatt=None
110 self.SPARtx=None
111 self.SPARaddGain0=None
112 self.SPARaddGain1=None
113 self.SPARwnd=None
114 self.SPARpos=None
115 self.SPARadd=None
116 self.SPARlen=None
117 self.SPARcal=None
118 self.SPARnos=None
119 self.SPARof0=None
120 self.SPARof1=None
121 self.SPARswt=None
122 self.SPARsum=None
123 self.SPARosc=None
124 self.SPARtst=None
125 self.SPARcor=None
126 self.SPARofs=None
127 self.SPARhsn=None
128 self.SPARhsa=None
129 self.SPARcalibPow_M=None
130 self.SPARcalibSNR_M=None
131 self.SPARcalibPow_S=None
132 self.SPARcalibSNR_S=None
133 self.SPARrawGate1=None
134 self.SPARrawGate2=None
135 self.SPARraw=None
136 self.SPARprc=None
137
138
139
140 header = numpy.fromfile(fp, FILE_HEADER,1)
33 startFp = open(
34 '/home/erick/Documents/MIRA35C/20160117/20160117_0000.zspc', "rb")
35
36
37 FILE_HEADER = numpy.dtype([ # HEADER 1024bytes
38 ('Hname', numpy.str_, 32), # Original file name
39 # Date and time when the file was created
40 ('Htime', numpy.str_, 32),
41 # Name of operator who created the file
42 ('Hoper', numpy.str_, 64),
43 # Place where the measurements was carried out
44 ('Hplace', numpy.str_, 128),
45 # Description of measurements
46 ('Hdescr', numpy.str_, 256),
47 ('Hdummy', numpy.str_, 512), # Reserved space
48 # Main chunk
49 ('Msign', '<i4'), # Main chunk signature FZKF or NUIG
50 ('MsizeData', '<i4'), # Size of data block main chunk
51 # Processing DSP parameters
52 ('PPARsign', '<i4'), # PPAR signature
53 ('PPARsize', '<i4'), # PPAR size of block
54 ('PPARprf', '<i4'), # Pulse repetition frequency
55 ('PPARpdr', '<i4'), # Pulse duration
56 ('PPARsft', '<i4'), # FFT length
57 # Number of spectral (in-coherent) averages
58 ('PPARavc', '<i4'),
59 # Number of lowest range gate for moment estimation
60 ('PPARihp', '<i4'),
61 # Count for gates for moment estimation
62 ('PPARchg', '<i4'),
63 # switch on/off polarimetric measurements. Should be 1.
64 ('PPARpol', '<i4'),
65 # Service DSP parameters
66 # STC attenuation on the lowest ranges on/off
67 ('SPARatt', '<i4'),
68 ('SPARtx', '<i4'), # OBSOLETE
69 ('SPARaddGain0', '<f4'), # OBSOLETE
70 ('SPARaddGain1', '<f4'), # OBSOLETE
71 # Debug only. It normal mode it is 0.
72 ('SPARwnd', '<i4'),
73 # Delay between sync pulse and tx pulse for phase corr, ns
74 ('SPARpos', '<i4'),
75 # "add to pulse" to compensate for delay between the leading edge of driver pulse and envelope of the RF signal.
76 ('SPARadd', '<i4'),
77 # Time for measuring txn pulse phase. OBSOLETE
78 ('SPARlen', '<i4'),
79 ('SPARcal', '<i4'), # OBSOLETE
80 ('SPARnos', '<i4'), # OBSOLETE
81 ('SPARof0', '<i4'), # detection threshold
82 ('SPARof1', '<i4'), # OBSOLETE
83 ('SPARswt', '<i4'), # 2nd moment estimation threshold
84 ('SPARsum', '<i4'), # OBSOLETE
85 ('SPARosc', '<i4'), # flag Oscillosgram mode
86 ('SPARtst', '<i4'), # OBSOLETE
87 ('SPARcor', '<i4'), # OBSOLETE
88 ('SPARofs', '<i4'), # OBSOLETE
89 # Hildebrand div noise detection on noise gate
90 ('SPARhsn', '<i4'),
91 # Hildebrand div noise detection on all gates
92 ('SPARhsa', '<f4'),
93 ('SPARcalibPow_M', '<f4'), # OBSOLETE
94 ('SPARcalibSNR_M', '<f4'), # OBSOLETE
95 ('SPARcalibPow_S', '<f4'), # OBSOLETE
96 ('SPARcalibSNR_S', '<f4'), # OBSOLETE
97 # Lowest range gate for spectra saving Raw_Gate1 >=5
98 ('SPARrawGate1', '<i4'),
99 # Number of range gates with atmospheric signal
100 ('SPARrawGate2', '<i4'),
101 # flag - IQ or spectra saving on/off
102 ('SPARraw', '<i4'),
103 ('SPARprc', '<i4'), ]) # flag - Moment estimation switched on/off
104
105
106 self.Hname = None
107 self.Htime = None
108 self.Hoper = None
109 self.Hplace = None
110 self.Hdescr = None
111 self.Hdummy = None
112
113 self.Msign = None
114 self.MsizeData = None
115
116 self.PPARsign = None
117 self.PPARsize = None
118 self.PPARprf = None
119 self.PPARpdr = None
120 self.PPARsft = None
121 self.PPARavc = None
122 self.PPARihp = None
123 self.PPARchg = None
124 self.PPARpol = None
125 # Service DSP parameters
126 self.SPARatt = None
127 self.SPARtx = None
128 self.SPARaddGain0 = None
129 self.SPARaddGain1 = None
130 self.SPARwnd = None
131 self.SPARpos = None
132 self.SPARadd = None
133 self.SPARlen = None
134 self.SPARcal = None
135 self.SPARnos = None
136 self.SPARof0 = None
137 self.SPARof1 = None
138 self.SPARswt = None
139 self.SPARsum = None
140 self.SPARosc = None
141 self.SPARtst = None
142 self.SPARcor = None
143 self.SPARofs = None
144 self.SPARhsn = None
145 self.SPARhsa = None
146 self.SPARcalibPow_M = None
147 self.SPARcalibSNR_M = None
148 self.SPARcalibPow_S = None
149 self.SPARcalibSNR_S = None
150 self.SPARrawGate1 = None
151 self.SPARrawGate2 = None
152 self.SPARraw = None
153 self.SPARprc = None
154
155
156 header = numpy.fromfile(fp, FILE_HEADER, 1)
141 157 ''' numpy.fromfile(file, dtype, count, sep='')
142 158 file : file or str
143 159 Open file object or filename.
144 160
145 161 dtype : data-type
146 162 Data type of the returned array. For binary files, it is used to determine
147 163 the size and byte-order of the items in the file.
148 164
149 165 count : int
150 166 Number of items to read. -1 means all items (i.e., the complete file).
151 167
152 168 sep : str
153 169 Separator between items if file is a text file. Empty ("") separator means
154 170 the file should be treated as binary. Spaces (" ") in the separator match zero
155 171 or more whitespace characters. A separator consisting only of spaces must match
156 172 at least one whitespace.
157 173
158 174 '''
159 175
160 Hname= str(header['Hname'][0])
161 Htime= str(header['Htime'][0])
162 Hoper= str(header['Hoper'][0])
163 Hplace= str(header['Hplace'][0])
164 Hdescr= str(header['Hdescr'][0])
165 Hdummy= str(header['Hdummy'][0])
166
167 Msign=header['Msign'][0]
168 MsizeData=header['MsizeData'][0]
169
170 PPARsign=header['PPARsign'][0]
171 PPARsize=header['PPARsize'][0]
172 PPARprf=header['PPARprf'][0]
173 PPARpdr=header['PPARpdr'][0]
174 PPARsft=header['PPARsft'][0]
175 PPARavc=header['PPARavc'][0]
176 PPARihp=header['PPARihp'][0]
177 PPARchg=header['PPARchg'][0]
178 PPARpol=header['PPARpol'][0]
179 #Service DSP parameters
180 SPARatt=header['SPARatt'][0]
181 SPARtx=header['SPARtx'][0]
182 SPARaddGain0=header['SPARaddGain0'][0]
183 SPARaddGain1=header['SPARaddGain1'][0]
184 SPARwnd=header['SPARwnd'][0]
185 SPARpos=header['SPARpos'][0]
186 SPARadd=header['SPARadd'][0]
187 SPARlen=header['SPARlen'][0]
188 SPARcal=header['SPARcal'][0]
189 SPARnos=header['SPARnos'][0]
190 SPARof0=header['SPARof0'][0]
191 SPARof1=header['SPARof1'][0]
192 SPARswt=header['SPARswt'][0]
193 SPARsum=header['SPARsum'][0]
194 SPARosc=header['SPARosc'][0]
195 SPARtst=header['SPARtst'][0]
196 SPARcor=header['SPARcor'][0]
197 SPARofs=header['SPARofs'][0]
198 SPARhsn=header['SPARhsn'][0]
199 SPARhsa=header['SPARhsa'][0]
200 SPARcalibPow_M=header['SPARcalibPow_M'][0]
201 SPARcalibSNR_M=header['SPARcalibSNR_M'][0]
202 SPARcalibPow_S=header['SPARcalibPow_S'][0]
203 SPARcalibSNR_S=header['SPARcalibSNR_S'][0]
204 SPARrawGate1=header['SPARrawGate1'][0]
205 SPARrawGate2=header['SPARrawGate2'][0]
206 SPARraw=header['SPARraw'][0]
207 SPARprc=header['SPARprc'][0]
208
209
210
211 SRVI_STRUCTURE = numpy.dtype([
212 ('frame_cnt','<u4'),#
213 ('time_t','<u4'), #
214 ('tpow','<f4'), #
215 ('npw1','<f4'), #
216 ('npw2','<f4'), #
217 ('cpw1','<f4'), #
218 ('pcw2','<f4'), #
219 ('ps_err','<u4'), #
220 ('te_err','<u4'), #
221 ('rc_err','<u4'), #
222 ('grs1','<u4'), #
223 ('grs2','<u4'), #
224 ('azipos','<f4'), #
225 ('azivel','<f4'), #
226 ('elvpos','<f4'), #
227 ('elvvel','<f4'), #
228 ('northAngle','<f4'), #
229 ('microsec','<u4'), #
230 ('azisetvel','<f4'), #
231 ('elvsetpos','<f4'), #
232 ('RadarConst','<f4'),]) #
233
234 JUMP_STRUCTURE = numpy.dtype([
235 ('jump','<u140'),#
236 ('SizeOfDataBlock1',numpy.str_,32),#
237 ('jump','<i4'),#
238 ('DataBlockTitleSRVI1',numpy.str_,32),#
239 ('SizeOfSRVI1','<i4'),])#
240
241
242
243 #frame_cnt=0, time_t= 0, tpow=0, npw1=0, npw2=0,
244 #cpw1=0, pcw2=0, ps_err=0, te_err=0, rc_err=0, grs1=0,
245 #grs2=0, azipos=0, azivel=0, elvpos=0, elvvel=0, northangle=0,
246 #microsec=0, azisetvel=0, elvsetpos=0, RadarConst=0
247
248
176 Hname = str(header['Hname'][0])
177 Htime = str(header['Htime'][0])
178 Hoper = str(header['Hoper'][0])
179 Hplace = str(header['Hplace'][0])
180 Hdescr = str(header['Hdescr'][0])
181 Hdummy = str(header['Hdummy'][0])
182
183 Msign = header['Msign'][0]
184 MsizeData = header['MsizeData'][0]
185
186 PPARsign = header['PPARsign'][0]
187 PPARsize = header['PPARsize'][0]
188 PPARprf = header['PPARprf'][0]
189 PPARpdr = header['PPARpdr'][0]
190 PPARsft = header['PPARsft'][0]
191 PPARavc = header['PPARavc'][0]
192 PPARihp = header['PPARihp'][0]
193 PPARchg = header['PPARchg'][0]
194 PPARpol = header['PPARpol'][0]
195 # Service DSP parameters
196 SPARatt = header['SPARatt'][0]
197 SPARtx = header['SPARtx'][0]
198 SPARaddGain0 = header['SPARaddGain0'][0]
199 SPARaddGain1 = header['SPARaddGain1'][0]
200 SPARwnd = header['SPARwnd'][0]
201 SPARpos = header['SPARpos'][0]
202 SPARadd = header['SPARadd'][0]
203 SPARlen = header['SPARlen'][0]
204 SPARcal = header['SPARcal'][0]
205 SPARnos = header['SPARnos'][0]
206 SPARof0 = header['SPARof0'][0]
207 SPARof1 = header['SPARof1'][0]
208 SPARswt = header['SPARswt'][0]
209 SPARsum = header['SPARsum'][0]
210 SPARosc = header['SPARosc'][0]
211 SPARtst = header['SPARtst'][0]
212 SPARcor = header['SPARcor'][0]
213 SPARofs = header['SPARofs'][0]
214 SPARhsn = header['SPARhsn'][0]
215 SPARhsa = header['SPARhsa'][0]
216 SPARcalibPow_M = header['SPARcalibPow_M'][0]
217 SPARcalibSNR_M = header['SPARcalibSNR_M'][0]
218 SPARcalibPow_S = header['SPARcalibPow_S'][0]
219 SPARcalibSNR_S = header['SPARcalibSNR_S'][0]
220 SPARrawGate1 = header['SPARrawGate1'][0]
221 SPARrawGate2 = header['SPARrawGate2'][0]
222 SPARraw = header['SPARraw'][0]
223 SPARprc = header['SPARprc'][0]
224
225
226 SRVI_STRUCTURE = numpy.dtype([
227 ('frame_cnt', '<u4'),
228 ('time_t', '<u4'), #
229 ('tpow', '<f4'), #
230 ('npw1', '<f4'), #
231 ('npw2', '<f4'), #
232 ('cpw1', '<f4'), #
233 ('pcw2', '<f4'), #
234 ('ps_err', '<u4'), #
235 ('te_err', '<u4'), #
236 ('rc_err', '<u4'), #
237 ('grs1', '<u4'), #
238 ('grs2', '<u4'), #
239 ('azipos', '<f4'), #
240 ('azivel', '<f4'), #
241 ('elvpos', '<f4'), #
242 ('elvvel', '<f4'), #
243 ('northAngle', '<f4'),
244 ('microsec', '<u4'), #
245 ('azisetvel', '<f4'), #
246 ('elvsetpos', '<f4'), #
247 ('RadarConst', '<f4'), ]) #
248
249 JUMP_STRUCTURE = numpy.dtype([
250 ('jump', '<u140'),
251 ('SizeOfDataBlock1', numpy.str_, 32),
252 ('jump', '<i4'),
253 ('DataBlockTitleSRVI1', numpy.str_, 32),
254 ('SizeOfSRVI1', '<i4'), ])
255
256
257 # frame_cnt=0, time_t= 0, tpow=0, npw1=0, npw2=0,
258 # cpw1=0, pcw2=0, ps_err=0, te_err=0, rc_err=0, grs1=0,
259 # grs2=0, azipos=0, azivel=0, elvpos=0, elvvel=0, northangle=0,
260 # microsec=0, azisetvel=0, elvsetpos=0, RadarConst=0
261
262
249 263 frame_cnt = frame_cnt
250 264 dwell = time_t
251 265 tpow = tpow
252 266 npw1 = npw1
253 267 npw2 = npw2
254 268 cpw1 = cpw1
255 269 pcw2 = pcw2
256 270 ps_err = ps_err
257 271 te_err = te_err
258 272 rc_err = rc_err
259 273 grs1 = grs1
260 274 grs2 = grs2
261 275 azipos = azipos
262 276 azivel = azivel
263 277 elvpos = elvpos
264 278 elvvel = elvvel
265 279 northAngle = northAngle
266 280 microsec = microsec
267 281 azisetvel = azisetvel
268 282 elvsetpos = elvsetpos
269 283 RadarConst5 = RadarConst
270
271 284
272 285
273 #print fp
274 #startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
275 #startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
276 #RecCounter=0
277 #Off2StartNxtRec=811248
278 #print 'OffsetStartHeader ',self.OffsetStartHeader,'RecCounter ', self.RecCounter, 'Off2StartNxtRec ' , self.Off2StartNxtRec
286 # print fp
287 # startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
288 # startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
289 # RecCounter=0
290 # Off2StartNxtRec=811248
291 # print 'OffsetStartHeader ',self.OffsetStartHeader,'RecCounter ', self.RecCounter, 'Off2StartNxtRec ' , self.Off2StartNxtRec
279 292 #OffRHeader= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
280 293 #startFp.seek(OffRHeader, os.SEEK_SET)
281 print 'debe ser 48, RecCounter*811248', self.OffsetStartHeader,self.RecCounter,self.Off2StartNxtRec
282 print 'Posicion del bloque: ',OffRHeader
294 print 'debe ser 48, RecCounter*811248', self.OffsetStartHeader, self.RecCounter, self.Off2StartNxtRec
295 print 'Posicion del bloque: ', OffRHeader
283 296
284 header = numpy.fromfile(startFp,SRVI_STRUCTURE,1)
297 header = numpy.fromfile(startFp, SRVI_STRUCTURE, 1)
285 298
286 self.frame_cnt = header['frame_cnt'][0]#
299 self.frame_cnt = header['frame_cnt'][0]
287 300 self.time_t = header['frame_cnt'][0] #
288 301 self.tpow = header['frame_cnt'][0] #
289 302 self.npw1 = header['frame_cnt'][0] #
290 303 self.npw2 = header['frame_cnt'][0] #
291 304 self.cpw1 = header['frame_cnt'][0] #
292 305 self.pcw2 = header['frame_cnt'][0] #
293 306 self.ps_err = header['frame_cnt'][0] #
294 307 self.te_err = header['frame_cnt'][0] #
295 308 self.rc_err = header['frame_cnt'][0] #
296 309 self.grs1 = header['frame_cnt'][0] #
297 310 self.grs2 = header['frame_cnt'][0] #
298 311 self.azipos = header['frame_cnt'][0] #
299 312 self.azivel = header['frame_cnt'][0] #
300 313 self.elvpos = header['frame_cnt'][0] #
301 314 self.elvvel = header['frame_cnt'][0] #
302 315 self.northAngle = header['frame_cnt'][0] #
303 316 self.microsec = header['frame_cnt'][0] #
304 317 self.azisetvel = header['frame_cnt'][0] #
305 318 self.elvsetpos = header['frame_cnt'][0] #
306 self.RadarConst = header['frame_cnt'][0] #
319 self.RadarConst = header['frame_cnt'][0] #
307 320
308 321
309 self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
322 self.ipp = 0.5 * (SPEED_OF_LIGHT / self.PRFhz)
310 323
311 self.RHsize = 180+20*self.nChannels
312 self.Datasize= self.nProfiles*self.nChannels*self.nHeights*2*4
313 #print 'Datasize',self.Datasize
314 endFp = self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
324 self.RHsize = 180 + 20 * self.nChannels
325 self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
326 # print 'Datasize',self.Datasize
327 endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
315 328
316 329 print '=============================================='
317 330
318 331 print '=============================================='
319
320
321 No newline at end of file
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1270 lines changed) Show them Hide them
@@ -1,1154 +1,1182
1 import os, sys
1 import os
2 import sys
2 3 import glob
3 4 import fnmatch
4 5 import datetime
5 6 import time
6 7 import re
7 8 import h5py
8 9 import numpy
9 import matplotlib.pyplot as plt
10 10
11 11 import pylab as plb
12 12 from scipy.optimize import curve_fit
13 13 from scipy import asarray as ar, exp
14 14 from scipy import stats
15 15
16 16 from numpy.ma.core import getdata
17 17
18 18 SPEED_OF_LIGHT = 299792458
19 19 SPEED_OF_LIGHT = 3e8
20 20
21 21 try:
22 22 from gevent import sleep
23 23 except:
24 24 from time import sleep
25 25
26 26 from schainpy.model.data.jrodata import Spectra
27 27 #from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader
28 28 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
29 29 #from schainpy.model.io.jroIO_bltr import BLTRReader
30 30 from numpy import imag, shape, NaN
31 31
32 32 from jroIO_base import JRODataReader
33 33
34 34
35 35 class Header(object):
36
36
37 37 def __init__(self):
38 38 raise NotImplementedError
39
40
39
41 40 def read(self):
42
41
43 42 raise NotImplementedError
44
43
45 44 def write(self):
46
45
47 46 raise NotImplementedError
48
47
49 48 def printInfo(self):
50
51 message = "#"*50 + "\n"
49
50 message = "#" * 50 + "\n"
52 51 message += self.__class__.__name__.upper() + "\n"
53 message += "#"*50 + "\n"
54
52 message += "#" * 50 + "\n"
53
55 54 keyList = self.__dict__.keys()
56 55 keyList.sort()
57
56
58 57 for key in keyList:
59 message += "%s = %s" %(key, self.__dict__[key]) + "\n"
60
58 message += "%s = %s" % (key, self.__dict__[key]) + "\n"
59
61 60 if "size" not in keyList:
62 61 attr = getattr(self, "size")
63
64 if attr:
65 message += "%s = %s" %("size", attr) + "\n"
66
67 #print message
68 62
63 if attr:
64 message += "%s = %s" % ("size", attr) + "\n"
69 65
66 # print message
70 67
71 68
69 FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes
70 ('FileMgcNumber', '<u4'), # 0x23020100
71 # No Of FDT data records in this file (0 or more)
72 ('nFDTdataRecors', '<u4'),
73 ('OffsetStartHeader', '<u4'),
74 ('RadarUnitId', '<u4'),
75 ('SiteName', numpy.str_, 32), # Null terminated
76 ])
72 77
73 FILE_STRUCTURE = numpy.dtype([ #HEADER 48bytes
74 ('FileMgcNumber','<u4'), #0x23020100
75 ('nFDTdataRecors','<u4'), #No Of FDT data records in this file (0 or more)
76 ('OffsetStartHeader','<u4'),
77 ('RadarUnitId','<u4'),
78 ('SiteName',numpy.str_,32), #Null terminated
79 ])
80 78
81 79 class FileHeaderBLTR(Header):
82
80
83 81 def __init__(self):
84
85 self.FileMgcNumber= 0 #0x23020100
86 self.nFDTdataRecors=0 #No Of FDT data records in this file (0 or more)
87 self.RadarUnitId= 0
88 self.OffsetStartHeader=0
89 self.SiteName= ""
82
83 self.FileMgcNumber = 0 # 0x23020100
84 # No Of FDT data records in this file (0 or more)
85 self.nFDTdataRecors = 0
86 self.RadarUnitId = 0
87 self.OffsetStartHeader = 0
88 self.SiteName = ""
90 89 self.size = 48
91
90
92 91 def FHread(self, fp):
93 #try:
94 startFp = open(fp,"rb")
95
96 header = numpy.fromfile(startFp, FILE_STRUCTURE,1)
97
92 # try:
93 startFp = open(fp, "rb")
94
95 header = numpy.fromfile(startFp, FILE_STRUCTURE, 1)
96
98 97 print ' '
99 98 print 'puntero file header', startFp.tell()
100 99 print ' '
101
102
100
103 101 ''' numpy.fromfile(file, dtype, count, sep='')
104 102 file : file or str
105 103 Open file object or filename.
106 104
107 105 dtype : data-type
108 106 Data type of the returned array. For binary files, it is used to determine
109 107 the size and byte-order of the items in the file.
110 108
111 109 count : int
112 110 Number of items to read. -1 means all items (i.e., the complete file).
113 111
114 112 sep : str
115 113 Separator between items if file is a text file. Empty ("") separator means
116 114 the file should be treated as binary. Spaces (" ") in the separator match zero
117 115 or more whitespace characters. A separator consisting only of spaces must match
118 116 at least one whitespace.
119 117
120 118 '''
121 119
120 self.FileMgcNumber = hex(header['FileMgcNumber'][0])
121 # No Of FDT data records in this file (0 or more)
122 self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
123 self.RadarUnitId = int(header['RadarUnitId'][0])
124 self.OffsetStartHeader = int(header['OffsetStartHeader'][0])
125 self.SiteName = str(header['SiteName'][0])
126
127 # print 'Numero de bloques', self.nFDTdataRecors
122 128
123
124 self.FileMgcNumber= hex(header['FileMgcNumber'][0])
125 self.nFDTdataRecors=int(header['nFDTdataRecors'][0]) #No Of FDT data records in this file (0 or more)
126 self.RadarUnitId= int(header['RadarUnitId'][0])
127 self.OffsetStartHeader= int(header['OffsetStartHeader'][0])
128 self.SiteName= str(header['SiteName'][0])
129
130 #print 'Numero de bloques', self.nFDTdataRecors
131
132
133 if self.size <48:
129 if self.size < 48:
134 130 return 0
135
131
136 132 return 1
137
138
133
139 134 def write(self, fp):
140
135
141 136 headerTuple = (self.FileMgcNumber,
142 137 self.nFDTdataRecors,
143 138 self.RadarUnitId,
144 139 self.SiteName,
145 140 self.size)
146
147
141
148 142 header = numpy.array(headerTuple, FILE_STRUCTURE)
149 143 # numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0)
150 144 header.tofile(fp)
151 145 ''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
152 146
153 147 fid : file or str
154 148 An open file object, or a string containing a filename.
155 149
156 150 sep : str
157 151 Separator between array items for text output. If "" (empty), a binary file is written,
158 152 equivalent to file.write(a.tobytes()).
159 153
160 154 format : str
161 155 Format string for text file output. Each entry in the array is formatted to text by
162 156 first converting it to the closest Python type, and then using "format" % item.
163 157
164 158 '''
165
166 return 1
167 159
160 return 1
168 161
169 162
163 RECORD_STRUCTURE = numpy.dtype([ # RECORD HEADER 180+20N bytes
164 ('RecMgcNumber', '<u4'), # 0x23030001
165 ('RecCounter', '<u4'), # Record counter(0,1, ...)
166 # Offset to start of next record form start of this record
167 ('Off2StartNxtRec', '<u4'),
168 # Offset to start of data from start of this record
169 ('Off2StartData', '<u4'),
170 # Epoch time stamp of start of acquisition (seconds)
171 ('nUtime', '<i4'),
172 # Millisecond component of time stamp (0,...,999)
173 ('nMilisec', '<u4'),
174 # Experiment tag name (null terminated)
175 ('ExpTagName', numpy.str_, 32),
176 # Experiment comment (null terminated)
177 ('ExpComment', numpy.str_, 32),
178 # Site latitude (from GPS) in degrees (positive implies North)
179 ('SiteLatDegrees', '<f4'),
180 # Site longitude (from GPS) in degrees (positive implies East)
181 ('SiteLongDegrees', '<f4'),
182 # RTC GPS engine status (0=SEEK, 1=LOCK, 2=NOT FITTED, 3=UNAVAILABLE)
183 ('RTCgpsStatus', '<u4'),
184 ('TransmitFrec', '<u4'), # Transmit frequency (Hz)
185 ('ReceiveFrec', '<u4'), # Receive frequency
186 # First local oscillator frequency (Hz)
187 ('FirstOsciFrec', '<u4'),
188 # (0="O", 1="E", 2="linear 1", 3="linear2")
189 ('Polarisation', '<u4'),
190 # Receiver filter settings (0,1,2,3)
191 ('ReceiverFiltSett', '<u4'),
192 # Number of modes in use (1 or 2)
193 ('nModesInUse', '<u4'),
194 # Dual Mode index number for these data (0 or 1)
195 ('DualModeIndex', '<u4'),
196 # Dual Mode range correction for these data (m)
197 ('DualModeRange', '<u4'),
198 # Number of digital channels acquired (2*N)
199 ('nDigChannels', '<u4'),
200 # Sampling resolution (meters)
201 ('SampResolution', '<u4'),
202 # Number of range gates sampled
203 ('nHeights', '<u4'),
204 # Start range of sampling (meters)
205 ('StartRangeSamp', '<u4'),
206 ('PRFhz', '<u4'), # PRF (Hz)
207 ('nCohInt', '<u4'), # Integrations
208 # Number of data points transformed
209 ('nProfiles', '<u4'),
210 # Number of receive beams stored in file (1 or N)
211 ('nChannels', '<u4'),
212 ('nIncohInt', '<u4'), # Number of spectral averages
213 # FFT windowing index (0 = no window)
214 ('FFTwindowingInd', '<u4'),
215 # Beam steer angle (azimuth) in degrees (clockwise from true North)
216 ('BeamAngleAzim', '<f4'),
217 # Beam steer angle (zenith) in degrees (0=> vertical)
218 ('BeamAngleZen', '<f4'),
219 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
220 ('AntennaCoord0', '<f4'),
221 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
222 ('AntennaAngl0', '<f4'),
223 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
224 ('AntennaCoord1', '<f4'),
225 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
226 ('AntennaAngl1', '<f4'),
227 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
228 ('AntennaCoord2', '<f4'),
229 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
230 ('AntennaAngl2', '<f4'),
231 # Receiver phase calibration (degrees) - N values
232 ('RecPhaseCalibr0', '<f4'),
233 # Receiver phase calibration (degrees) - N values
234 ('RecPhaseCalibr1', '<f4'),
235 # Receiver phase calibration (degrees) - N values
236 ('RecPhaseCalibr2', '<f4'),
237 # Receiver amplitude calibration (ratio relative to receiver one) - N values
238 ('RecAmpCalibr0', '<f4'),
239 # Receiver amplitude calibration (ratio relative to receiver one) - N values
240 ('RecAmpCalibr1', '<f4'),
241 # Receiver amplitude calibration (ratio relative to receiver one) - N values
242 ('RecAmpCalibr2', '<f4'),
243 # Receiver gains in dB - N values
244 ('ReceiverGaindB0', '<i4'),
245 # Receiver gains in dB - N values
246 ('ReceiverGaindB1', '<i4'),
247 # Receiver gains in dB - N values
248 ('ReceiverGaindB2', '<i4'),
249 ])
170 250
171 251
172 RECORD_STRUCTURE = numpy.dtype([ #RECORD HEADER 180+20N bytes
173 ('RecMgcNumber','<u4'), #0x23030001
174 ('RecCounter','<u4'), #Record counter(0,1, ...)
175 ('Off2StartNxtRec','<u4'), #Offset to start of next record form start of this record
176 ('Off2StartData','<u4'), #Offset to start of data from start of this record
177 ('nUtime','<i4'), #Epoch time stamp of start of acquisition (seconds)
178 ('nMilisec','<u4'), #Millisecond component of time stamp (0,...,999)
179 ('ExpTagName',numpy.str_,32), #Experiment tag name (null terminated)
180 ('ExpComment',numpy.str_,32), #Experiment comment (null terminated)
181 ('SiteLatDegrees','<f4'), #Site latitude (from GPS) in degrees (positive implies North)
182 ('SiteLongDegrees','<f4'), #Site longitude (from GPS) in degrees (positive implies East)
183 ('RTCgpsStatus','<u4'), #RTC GPS engine status (0=SEEK, 1=LOCK, 2=NOT FITTED, 3=UNAVAILABLE)
184 ('TransmitFrec','<u4'), #Transmit frequency (Hz)
185 ('ReceiveFrec','<u4'), #Receive frequency
186 ('FirstOsciFrec','<u4'), #First local oscillator frequency (Hz)
187 ('Polarisation','<u4'), #(0="O", 1="E", 2="linear 1", 3="linear2")
188 ('ReceiverFiltSett','<u4'), #Receiver filter settings (0,1,2,3)
189 ('nModesInUse','<u4'), #Number of modes in use (1 or 2)
190 ('DualModeIndex','<u4'), #Dual Mode index number for these data (0 or 1)
191 ('DualModeRange','<u4'), #Dual Mode range correction for these data (m)
192 ('nDigChannels','<u4'), #Number of digital channels acquired (2*N)
193 ('SampResolution','<u4'), #Sampling resolution (meters)
194 ('nHeights','<u4'), #Number of range gates sampled
195 ('StartRangeSamp','<u4'), #Start range of sampling (meters)
196 ('PRFhz','<u4'), #PRF (Hz)
197 ('nCohInt','<u4'), #Integrations
198 ('nProfiles','<u4'), #Number of data points transformed
199 ('nChannels','<u4'), #Number of receive beams stored in file (1 or N)
200 ('nIncohInt','<u4'), #Number of spectral averages
201 ('FFTwindowingInd','<u4'), #FFT windowing index (0 = no window)
202 ('BeamAngleAzim','<f4'), #Beam steer angle (azimuth) in degrees (clockwise from true North)
203 ('BeamAngleZen','<f4'), #Beam steer angle (zenith) in degrees (0=> vertical)
204 ('AntennaCoord0','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
205 ('AntennaAngl0','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
206 ('AntennaCoord1','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
207 ('AntennaAngl1','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
208 ('AntennaCoord2','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
209 ('AntennaAngl2','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
210 ('RecPhaseCalibr0','<f4'), #Receiver phase calibration (degrees) - N values
211 ('RecPhaseCalibr1','<f4'), #Receiver phase calibration (degrees) - N values
212 ('RecPhaseCalibr2','<f4'), #Receiver phase calibration (degrees) - N values
213 ('RecAmpCalibr0','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
214 ('RecAmpCalibr1','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
215 ('RecAmpCalibr2','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
216 ('ReceiverGaindB0','<i4'), #Receiver gains in dB - N values
217 ('ReceiverGaindB1','<i4'), #Receiver gains in dB - N values
218 ('ReceiverGaindB2','<i4'), #Receiver gains in dB - N values
219 ])
252 class RecordHeaderBLTR(Header):
220 253
254 def __init__(self, RecMgcNumber=None, RecCounter=0, Off2StartNxtRec=811248,
255 nUtime=0, nMilisec=0, ExpTagName=None,
256 ExpComment=None, SiteLatDegrees=0, SiteLongDegrees=0,
257 RTCgpsStatus=0, TransmitFrec=0, ReceiveFrec=0,
258 FirstOsciFrec=0, Polarisation=0, ReceiverFiltSett=0,
259 nModesInUse=0, DualModeIndex=0, DualModeRange=0,
260 nDigChannels=0, SampResolution=0, nHeights=0,
261 StartRangeSamp=0, PRFhz=0, nCohInt=0,
262 nProfiles=0, nChannels=0, nIncohInt=0,
263 FFTwindowingInd=0, BeamAngleAzim=0, BeamAngleZen=0,
264 AntennaCoord0=0, AntennaCoord1=0, AntennaCoord2=0,
265 RecPhaseCalibr0=0, RecPhaseCalibr1=0, RecPhaseCalibr2=0,
266 RecAmpCalibr0=0, RecAmpCalibr1=0, RecAmpCalibr2=0,
267 AntennaAngl0=0, AntennaAngl1=0, AntennaAngl2=0,
268 ReceiverGaindB0=0, ReceiverGaindB1=0, ReceiverGaindB2=0, Off2StartData=0, OffsetStartHeader=0):
221 269
222 class RecordHeaderBLTR(Header):
223
224 def __init__(self, RecMgcNumber=None, RecCounter= 0, Off2StartNxtRec= 811248,
225 nUtime= 0, nMilisec= 0, ExpTagName= None,
226 ExpComment=None, SiteLatDegrees=0, SiteLongDegrees= 0,
227 RTCgpsStatus= 0, TransmitFrec= 0, ReceiveFrec= 0,
228 FirstOsciFrec= 0, Polarisation= 0, ReceiverFiltSett= 0,
229 nModesInUse= 0, DualModeIndex= 0, DualModeRange= 0,
230 nDigChannels= 0, SampResolution= 0, nHeights= 0,
231 StartRangeSamp= 0, PRFhz= 0, nCohInt= 0,
232 nProfiles= 0, nChannels= 0, nIncohInt= 0,
233 FFTwindowingInd= 0, BeamAngleAzim= 0, BeamAngleZen= 0,
234 AntennaCoord0= 0, AntennaCoord1= 0, AntennaCoord2= 0,
235 RecPhaseCalibr0= 0, RecPhaseCalibr1= 0, RecPhaseCalibr2= 0,
236 RecAmpCalibr0= 0, RecAmpCalibr1= 0, RecAmpCalibr2= 0,
237 AntennaAngl0=0, AntennaAngl1=0, AntennaAngl2=0,
238 ReceiverGaindB0= 0, ReceiverGaindB1= 0, ReceiverGaindB2= 0, Off2StartData=0, OffsetStartHeader=0):
239
240 self.RecMgcNumber = RecMgcNumber #0x23030001
241 self.RecCounter = RecCounter
242 self.Off2StartNxtRec = Off2StartNxtRec
270 self.RecMgcNumber = RecMgcNumber # 0x23030001
271 self.RecCounter = RecCounter
272 self.Off2StartNxtRec = Off2StartNxtRec
243 273 self.Off2StartData = Off2StartData
244 self.nUtime = nUtime
245 self.nMilisec = nMilisec
246 self.ExpTagName = ExpTagName
247 self.ExpComment = ExpComment
248 self.SiteLatDegrees = SiteLatDegrees
249 self.SiteLongDegrees = SiteLongDegrees
250 self.RTCgpsStatus = RTCgpsStatus
251 self.TransmitFrec = TransmitFrec
252 self.ReceiveFrec = ReceiveFrec
253 self.FirstOsciFrec = FirstOsciFrec
254 self.Polarisation = Polarisation
255 self.ReceiverFiltSett = ReceiverFiltSett
256 self.nModesInUse = nModesInUse
257 self.DualModeIndex = DualModeIndex
258 self.DualModeRange = DualModeRange
274 self.nUtime = nUtime
275 self.nMilisec = nMilisec
276 self.ExpTagName = ExpTagName
277 self.ExpComment = ExpComment
278 self.SiteLatDegrees = SiteLatDegrees
279 self.SiteLongDegrees = SiteLongDegrees
280 self.RTCgpsStatus = RTCgpsStatus
281 self.TransmitFrec = TransmitFrec
282 self.ReceiveFrec = ReceiveFrec
283 self.FirstOsciFrec = FirstOsciFrec
284 self.Polarisation = Polarisation
285 self.ReceiverFiltSett = ReceiverFiltSett
286 self.nModesInUse = nModesInUse
287 self.DualModeIndex = DualModeIndex
288 self.DualModeRange = DualModeRange
259 289 self.nDigChannels = nDigChannels
260 self.SampResolution = SampResolution
261 self.nHeights = nHeights
262 self.StartRangeSamp = StartRangeSamp
263 self.PRFhz = PRFhz
264 self.nCohInt = nCohInt
265 self.nProfiles = nProfiles
266 self.nChannels = nChannels
267 self.nIncohInt = nIncohInt
268 self.FFTwindowingInd = FFTwindowingInd
269 self.BeamAngleAzim = BeamAngleAzim
270 self.BeamAngleZen = BeamAngleZen
290 self.SampResolution = SampResolution
291 self.nHeights = nHeights
292 self.StartRangeSamp = StartRangeSamp
293 self.PRFhz = PRFhz
294 self.nCohInt = nCohInt
295 self.nProfiles = nProfiles
296 self.nChannels = nChannels
297 self.nIncohInt = nIncohInt
298 self.FFTwindowingInd = FFTwindowingInd
299 self.BeamAngleAzim = BeamAngleAzim
300 self.BeamAngleZen = BeamAngleZen
271 301 self.AntennaCoord0 = AntennaCoord0
272 302 self.AntennaAngl0 = AntennaAngl0
273 303 self.AntennaAngl1 = AntennaAngl1
274 304 self.AntennaAngl2 = AntennaAngl2
275 305 self.AntennaCoord1 = AntennaCoord1
276 self.AntennaCoord2 = AntennaCoord2
306 self.AntennaCoord2 = AntennaCoord2
277 307 self.RecPhaseCalibr0 = RecPhaseCalibr0
278 308 self.RecPhaseCalibr1 = RecPhaseCalibr1
279 self.RecPhaseCalibr2 = RecPhaseCalibr2
309 self.RecPhaseCalibr2 = RecPhaseCalibr2
280 310 self.RecAmpCalibr0 = RecAmpCalibr0
281 311 self.RecAmpCalibr1 = RecAmpCalibr1
282 self.RecAmpCalibr2 = RecAmpCalibr2
312 self.RecAmpCalibr2 = RecAmpCalibr2
283 313 self.ReceiverGaindB0 = ReceiverGaindB0
284 314 self.ReceiverGaindB1 = ReceiverGaindB1
285 self.ReceiverGaindB2 = ReceiverGaindB2
286 self.OffsetStartHeader = 48
287
288
289
315 self.ReceiverGaindB2 = ReceiverGaindB2
316 self.OffsetStartHeader = 48
317
290 318 def RHread(self, fp):
291 #print fp
292 #startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
293 startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
294 #RecCounter=0
295 #Off2StartNxtRec=811248
296 OffRHeader= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
319 # print fp
320 # startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
321 # The method tell() returns the current position of the file read/write pointer within the file.
322 startFp = open(fp, "rb")
323 # RecCounter=0
324 # Off2StartNxtRec=811248
325 OffRHeader = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
297 326 print ' '
298 327 print 'puntero Record Header', startFp.tell()
299 328 print ' '
300
301
329
302 330 startFp.seek(OffRHeader, os.SEEK_SET)
303
331
304 332 print ' '
305 333 print 'puntero Record Header con seek', startFp.tell()
306 334 print ' '
307
308 #print 'Posicion del bloque: ',OffRHeader
309
310 header = numpy.fromfile(startFp,RECORD_STRUCTURE,1)
311
335
336 # print 'Posicion del bloque: ',OffRHeader
337
338 header = numpy.fromfile(startFp, RECORD_STRUCTURE, 1)
339
312 340 print ' '
313 341 print 'puntero Record Header con seek', startFp.tell()
314 342 print ' '
315
343
316 344 print ' '
317 345 #
318 #print 'puntero Record Header despues de seek', header.tell()
346 # print 'puntero Record Header despues de seek', header.tell()
319 347 print ' '
320
321 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) #0x23030001
322 self.RecCounter = int(header['RecCounter'][0])
323 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
324 self.Off2StartData = int(header['Off2StartData'][0])
325 self.nUtime = header['nUtime'][0]
326 self.nMilisec = header['nMilisec'][0]
327 self.ExpTagName = str(header['ExpTagName'][0])
328 self.ExpComment = str(header['ExpComment'][0])
329 self.SiteLatDegrees = header['SiteLatDegrees'][0]
330 self.SiteLongDegrees = header['SiteLongDegrees'][0]
331 self.RTCgpsStatus = header['RTCgpsStatus'][0]
332 self.TransmitFrec = header['TransmitFrec'][0]
333 self.ReceiveFrec = header['ReceiveFrec'][0]
334 self.FirstOsciFrec = header['FirstOsciFrec'][0]
335 self.Polarisation = header['Polarisation'][0]
336 self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
337 self.nModesInUse = header['nModesInUse'][0]
338 self.DualModeIndex = header['DualModeIndex'][0]
339 self.DualModeRange = header['DualModeRange'][0]
348
349 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) # 0x23030001
350 self.RecCounter = int(header['RecCounter'][0])
351 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
352 self.Off2StartData = int(header['Off2StartData'][0])
353 self.nUtime = header['nUtime'][0]
354 self.nMilisec = header['nMilisec'][0]
355 self.ExpTagName = str(header['ExpTagName'][0])
356 self.ExpComment = str(header['ExpComment'][0])
357 self.SiteLatDegrees = header['SiteLatDegrees'][0]
358 self.SiteLongDegrees = header['SiteLongDegrees'][0]
359 self.RTCgpsStatus = header['RTCgpsStatus'][0]
360 self.TransmitFrec = header['TransmitFrec'][0]
361 self.ReceiveFrec = header['ReceiveFrec'][0]
362 self.FirstOsciFrec = header['FirstOsciFrec'][0]
363 self.Polarisation = header['Polarisation'][0]
364 self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
365 self.nModesInUse = header['nModesInUse'][0]
366 self.DualModeIndex = header['DualModeIndex'][0]
367 self.DualModeRange = header['DualModeRange'][0]
340 368 self.nDigChannels = header['nDigChannels'][0]
341 self.SampResolution = header['SampResolution'][0]
342 self.nHeights = header['nHeights'][0]
343 self.StartRangeSamp = header['StartRangeSamp'][0]
344 self.PRFhz = header['PRFhz'][0]
345 self.nCohInt = header['nCohInt'][0]
346 self.nProfiles = header['nProfiles'][0]
347 self.nChannels = header['nChannels'][0]
348 self.nIncohInt = header['nIncohInt'][0]
349 self.FFTwindowingInd = header['FFTwindowingInd'][0]
350 self.BeamAngleAzim = header['BeamAngleAzim'][0]
351 self.BeamAngleZen = header['BeamAngleZen'][0]
369 self.SampResolution = header['SampResolution'][0]
370 self.nHeights = header['nHeights'][0]
371 self.StartRangeSamp = header['StartRangeSamp'][0]
372 self.PRFhz = header['PRFhz'][0]
373 self.nCohInt = header['nCohInt'][0]
374 self.nProfiles = header['nProfiles'][0]
375 self.nChannels = header['nChannels'][0]
376 self.nIncohInt = header['nIncohInt'][0]
377 self.FFTwindowingInd = header['FFTwindowingInd'][0]
378 self.BeamAngleAzim = header['BeamAngleAzim'][0]
379 self.BeamAngleZen = header['BeamAngleZen'][0]
352 380 self.AntennaCoord0 = header['AntennaCoord0'][0]
353 381 self.AntennaAngl0 = header['AntennaAngl0'][0]
354 382 self.AntennaCoord1 = header['AntennaCoord1'][0]
355 self.AntennaAngl1 = header['AntennaAngl1'][0]
383 self.AntennaAngl1 = header['AntennaAngl1'][0]
356 384 self.AntennaCoord2 = header['AntennaCoord2'][0]
357 self.AntennaAngl2 = header['AntennaAngl2'][0]
385 self.AntennaAngl2 = header['AntennaAngl2'][0]
358 386 self.RecPhaseCalibr0 = header['RecPhaseCalibr0'][0]
359 387 self.RecPhaseCalibr1 = header['RecPhaseCalibr1'][0]
360 self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0]
388 self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0]
361 389 self.RecAmpCalibr0 = header['RecAmpCalibr0'][0]
362 390 self.RecAmpCalibr1 = header['RecAmpCalibr1'][0]
363 self.RecAmpCalibr2 = header['RecAmpCalibr2'][0]
391 self.RecAmpCalibr2 = header['RecAmpCalibr2'][0]
364 392 self.ReceiverGaindB0 = header['ReceiverGaindB0'][0]
365 393 self.ReceiverGaindB1 = header['ReceiverGaindB1'][0]
366 394 self.ReceiverGaindB2 = header['ReceiverGaindB2'][0]
367 395
368 self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
369
370 self.RHsize = 180+20*self.nChannels
371 self.Datasize= self.nProfiles*self.nChannels*self.nHeights*2*4
372 #print 'Datasize',self.Datasize
373 endFp = self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
374
396 self.ipp = 0.5 * (SPEED_OF_LIGHT / self.PRFhz)
397
398 self.RHsize = 180 + 20 * self.nChannels
399 self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
400 # print 'Datasize',self.Datasize
401 endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
402
375 403 print '=============================================='
376 print 'RecMgcNumber ',self.RecMgcNumber
377 print 'RecCounter ',self.RecCounter
378 print 'Off2StartNxtRec ',self.Off2StartNxtRec
379 print 'Off2StartData ',self.Off2StartData
380 print 'Range Resolution ',self.SampResolution
381 print 'First Height ',self.StartRangeSamp
382 print 'PRF (Hz) ',self.PRFhz
383 print 'Heights (K) ',self.nHeights
384 print 'Channels (N) ',self.nChannels
385 print 'Profiles (J) ',self.nProfiles
386 print 'iCoh ',self.nCohInt
387 print 'iInCoh ',self.nIncohInt
388 print 'BeamAngleAzim ',self.BeamAngleAzim
389 print 'BeamAngleZen ',self.BeamAngleZen
390
391 #print 'ModoEnUso ',self.DualModeIndex
392 #print 'UtcTime ',self.nUtime
393 #print 'MiliSec ',self.nMilisec
394 #print 'Exp TagName ',self.ExpTagName
395 #print 'Exp Comment ',self.ExpComment
396 #print 'FFT Window Index ',self.FFTwindowingInd
397 #print 'N Dig. Channels ',self.nDigChannels
398 print 'Size de bloque ',self.RHsize
399 print 'DataSize ',self.Datasize
400 print 'BeamAngleAzim ',self.BeamAngleAzim
401 #print 'AntennaCoord0 ',self.AntennaCoord0
402 #print 'AntennaAngl0 ',self.AntennaAngl0
403 #print 'AntennaCoord1 ',self.AntennaCoord1
404 #print 'AntennaAngl1 ',self.AntennaAngl1
405 #print 'AntennaCoord2 ',self.AntennaCoord2
406 #print 'AntennaAngl2 ',self.AntennaAngl2
407 print 'RecPhaseCalibr0 ',self.RecPhaseCalibr0
408 print 'RecPhaseCalibr1 ',self.RecPhaseCalibr1
409 print 'RecPhaseCalibr2 ',self.RecPhaseCalibr2
410 print 'RecAmpCalibr0 ',self.RecAmpCalibr0
411 print 'RecAmpCalibr1 ',self.RecAmpCalibr1
412 print 'RecAmpCalibr2 ',self.RecAmpCalibr2
413 print 'ReceiverGaindB0 ',self.ReceiverGaindB0
414 print 'ReceiverGaindB1 ',self.ReceiverGaindB1
415 print 'ReceiverGaindB2 ',self.ReceiverGaindB2
404 print 'RecMgcNumber ', self.RecMgcNumber
405 print 'RecCounter ', self.RecCounter
406 print 'Off2StartNxtRec ', self.Off2StartNxtRec
407 print 'Off2StartData ', self.Off2StartData
408 print 'Range Resolution ', self.SampResolution
409 print 'First Height ', self.StartRangeSamp
410 print 'PRF (Hz) ', self.PRFhz
411 print 'Heights (K) ', self.nHeights
412 print 'Channels (N) ', self.nChannels
413 print 'Profiles (J) ', self.nProfiles
414 print 'iCoh ', self.nCohInt
415 print 'iInCoh ', self.nIncohInt
416 print 'BeamAngleAzim ', self.BeamAngleAzim
417 print 'BeamAngleZen ', self.BeamAngleZen
418
419 # print 'ModoEnUso ',self.DualModeIndex
420 # print 'UtcTime ',self.nUtime
421 # print 'MiliSec ',self.nMilisec
422 # print 'Exp TagName ',self.ExpTagName
423 # print 'Exp Comment ',self.ExpComment
424 # print 'FFT Window Index ',self.FFTwindowingInd
425 # print 'N Dig. Channels ',self.nDigChannels
426 print 'Size de bloque ', self.RHsize
427 print 'DataSize ', self.Datasize
428 print 'BeamAngleAzim ', self.BeamAngleAzim
429 # print 'AntennaCoord0 ',self.AntennaCoord0
430 # print 'AntennaAngl0 ',self.AntennaAngl0
431 # print 'AntennaCoord1 ',self.AntennaCoord1
432 # print 'AntennaAngl1 ',self.AntennaAngl1
433 # print 'AntennaCoord2 ',self.AntennaCoord2
434 # print 'AntennaAngl2 ',self.AntennaAngl2
435 print 'RecPhaseCalibr0 ', self.RecPhaseCalibr0
436 print 'RecPhaseCalibr1 ', self.RecPhaseCalibr1
437 print 'RecPhaseCalibr2 ', self.RecPhaseCalibr2
438 print 'RecAmpCalibr0 ', self.RecAmpCalibr0
439 print 'RecAmpCalibr1 ', self.RecAmpCalibr1
440 print 'RecAmpCalibr2 ', self.RecAmpCalibr2
441 print 'ReceiverGaindB0 ', self.ReceiverGaindB0
442 print 'ReceiverGaindB1 ', self.ReceiverGaindB1
443 print 'ReceiverGaindB2 ', self.ReceiverGaindB2
416 444 print '=============================================='
417
445
418 446 if OffRHeader > endFp:
419 sys.stderr.write("Warning %s: Size value read from System Header is lower than it has to be\n" %fp)
447 sys.stderr.write(
448 "Warning %s: Size value read from System Header is lower than it has to be\n" % fp)
420 449 return 0
421
450
422 451 if OffRHeader < endFp:
423 sys.stderr.write("Warning %s: Size value read from System Header size is greater than it has to be\n" %fp)
452 sys.stderr.write(
453 "Warning %s: Size value read from System Header size is greater than it has to be\n" % fp)
424 454 return 0
425
455
426 456 return 1
427 457
428
458
429 459 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODataReader):
430
460
431 461 path = None
432 462 startDate = None
433 463 endDate = None
434 464 startTime = None
435 465 endTime = None
436 466 walk = None
437 467 isConfig = False
438
439
440 fileList= None
441
442 #metadata
443 TimeZone= None
444 Interval= None
445 heightList= None
446
447 #data
448 data= None
449 utctime= None
450
451
452
468
469 fileList = None
470
471 # metadata
472 TimeZone = None
473 Interval = None
474 heightList = None
475
476 # data
477 data = None
478 utctime = None
479
453 480 def __init__(self, **kwargs):
454
455 #Eliminar de la base la herencia
481
482 # Eliminar de la base la herencia
456 483 ProcessingUnit.__init__(self, **kwargs)
457
484
458 485 #self.isConfig = False
459
486
460 487 #self.pts2read_SelfSpectra = 0
461 488 #self.pts2read_CrossSpectra = 0
462 489 #self.pts2read_DCchannels = 0
463 490 #self.datablock = None
464 491 self.utc = None
465 492 self.ext = ".fdt"
466 493 self.optchar = "P"
467 self.fpFile=None
494 self.fpFile = None
468 495 self.fp = None
469 self.BlockCounter=0
496 self.BlockCounter = 0
470 497 self.dtype = None
471 498 self.fileSizeByHeader = None
472 499 self.filenameList = []
473 500 self.fileSelector = 0
474 self.Off2StartNxtRec=0
475 self.RecCounter=0
501 self.Off2StartNxtRec = 0
502 self.RecCounter = 0
476 503 self.flagNoMoreFiles = 0
477 self.data_spc=None
478 self.data_cspc=None
479 self.data_output=None
504 self.data_spc = None
505 self.data_cspc = None
506 self.data_output = None
480 507 self.path = None
481 self.OffsetStartHeader=0
482 self.Off2StartData=0
508 self.OffsetStartHeader = 0
509 self.Off2StartData = 0
483 510 self.ipp = 0
484 self.nFDTdataRecors=0
511 self.nFDTdataRecors = 0
485 512 self.blocksize = 0
486 513 self.dataOut = Spectra()
487 self.profileIndex = 1 #Always
488 self.dataOut.flagNoData=False
514 self.profileIndex = 1 # Always
515 self.dataOut.flagNoData = False
489 516 self.dataOut.nRdPairs = 0
490 517 self.dataOut.pairsList = []
491 self.dataOut.data_spc=None
492 self.dataOut.noise=[]
493 self.dataOut.velocityX=[]
494 self.dataOut.velocityY=[]
495 self.dataOut.velocityV=[]
496
497
518 self.dataOut.data_spc = None
519 self.dataOut.noise = []
520 self.dataOut.velocityX = []
521 self.dataOut.velocityY = []
522 self.dataOut.velocityV = []
498 523
499 524 def Files2Read(self, fp):
500 525 '''
501 526 Function that indicates the number of .fdt files that exist in the folder to be read.
502 527 It also creates an organized list with the names of the files to read.
503 528 '''
504 #self.__checkPath()
505
506 ListaData=os.listdir(fp) #Gets the list of files within the fp address
507 ListaData=sorted(ListaData) #Sort the list of files from least to largest by names
508 nFiles=0 #File Counter
509 FileList=[] #A list is created that will contain the .fdt files
510 for IndexFile in ListaData :
511 if '.fdt' in IndexFile:
529 # self.__checkPath()
530
531 # Gets the list of files within the fp address
532 ListaData = os.listdir(fp)
533 # Sort the list of files from least to largest by names
534 ListaData = sorted(ListaData)
535 nFiles = 0 # File Counter
536 FileList = [] # A list is created that will contain the .fdt files
537 for IndexFile in ListaData:
538 if '.fdt' in IndexFile:
512 539 FileList.append(IndexFile)
513 nFiles+=1
514
515 #print 'Files2Read'
516 #print 'Existen '+str(nFiles)+' archivos .fdt'
517
518 self.filenameList=FileList #List of files from least to largest by names
519
520
540 nFiles += 1
541
542 # print 'Files2Read'
543 # print 'Existen '+str(nFiles)+' archivos .fdt'
544
545 self.filenameList = FileList # List of files from least to largest by names
546
521 547 def run(self, **kwargs):
522 548 '''
523 549 This method will be the one that will initiate the data entry, will be called constantly.
524 550 You should first verify that your Setup () is set up and then continue to acquire
525 551 the data to be processed with getData ().
526 552 '''
527 553 if not self.isConfig:
528 554 self.setup(**kwargs)
529 555 self.isConfig = True
530
556
531 557 self.getData()
532 #print 'running'
533
534
558 # print 'running'
559
535 560 def setup(self, path=None,
536 startDate=None,
537 endDate=None,
538 startTime=None,
539 endTime=None,
540 walk=True,
541 timezone='utc',
542 code = None,
543 online=False,
544 ReadMode=None,
545 **kwargs):
546
561 startDate=None,
562 endDate=None,
563 startTime=None,
564 endTime=None,
565 walk=True,
566 timezone='utc',
567 code=None,
568 online=False,
569 ReadMode=None,
570 **kwargs):
571
547 572 self.isConfig = True
548
549 self.path=path
550 self.startDate=startDate
551 self.endDate=endDate
552 self.startTime=startTime
553 self.endTime=endTime
554 self.walk=walk
555 self.ReadMode=int(ReadMode)
556
573
574 self.path = path
575 self.startDate = startDate
576 self.endDate = endDate
577 self.startTime = startTime
578 self.endTime = endTime
579 self.walk = walk
580 self.ReadMode = int(ReadMode)
581
557 582 pass
558
559
583
560 584 def getData(self):
561 585 '''
562 586 Before starting this function, you should check that there is still an unread file,
563 587 If there are still blocks to read or if the data block is empty.
564
588
565 589 You should call the file "read".
566
590
567 591 '''
568
592
569 593 if self.flagNoMoreFiles:
570 594 self.dataOut.flagNoData = True
571 595 print 'NoData se vuelve true'
572 596 return 0
573 597
574 self.fp=self.path
598 self.fp = self.path
575 599 self.Files2Read(self.fp)
576 600 self.readFile(self.fp)
577 601 self.dataOut.data_spc = self.data_spc
578 self.dataOut.data_cspc =self.data_cspc
579 self.dataOut.data_output=self.data_output
580
602 self.dataOut.data_cspc = self.data_cspc
603 self.dataOut.data_output = self.data_output
604
581 605 print 'self.dataOut.data_output', shape(self.dataOut.data_output)
582
583 #self.removeDC()
584 return self.dataOut.data_spc
585
586
587 def readFile(self,fp):
606
607 # self.removeDC()
608 return self.dataOut.data_spc
609
610 def readFile(self, fp):
588 611 '''
589 612 You must indicate if you are reading in Online or Offline mode and load the
590 613 The parameters for this file reading mode.
591
614
592 615 Then you must do 2 actions:
593
616
594 617 1. Get the BLTR FileHeader.
595 618 2. Start reading the first block.
596 619 '''
597
598 #The address of the folder is generated the name of the .fdt file that will be read
599 print "File: ",self.fileSelector+1
600
620
621 # The address of the folder is generated the name of the .fdt file that will be read
622 print "File: ", self.fileSelector + 1
623
601 624 if self.fileSelector < len(self.filenameList):
602
603 self.fpFile=str(fp)+'/'+str(self.filenameList[self.fileSelector])
604 #print self.fpFile
625
626 self.fpFile = str(fp) + '/' + \
627 str(self.filenameList[self.fileSelector])
628 # print self.fpFile
605 629 fheader = FileHeaderBLTR()
606 fheader.FHread(self.fpFile) #Bltr FileHeader Reading
607 self.nFDTdataRecors=fheader.nFDTdataRecors
608
609 self.readBlock() #Block reading
630 fheader.FHread(self.fpFile) # Bltr FileHeader Reading
631 self.nFDTdataRecors = fheader.nFDTdataRecors
632
633 self.readBlock() # Block reading
610 634 else:
611 635 print 'readFile FlagNoData becomes true'
612 self.flagNoMoreFiles=True
636 self.flagNoMoreFiles = True
613 637 self.dataOut.flagNoData = True
614 return 0
615
638 return 0
639
616 640 def getVelRange(self, extrapoints=0):
617 Lambda= SPEED_OF_LIGHT/50000000
618 PRF = self.dataOut.PRF#1./(self.dataOut.ippSeconds * self.dataOut.nCohInt)
619 Vmax=-Lambda/(4.*(1./PRF)*self.dataOut.nCohInt*2.)
620 deltafreq = PRF / (self.nProfiles)
621 deltavel = (Vmax*2) / (self.nProfiles)
622 freqrange = deltafreq*(numpy.arange(self.nProfiles)-self.nProfiles/2.) - deltafreq/2
623 velrange = deltavel*(numpy.arange(self.nProfiles)-self.nProfiles/2.)
624 return velrange
625
626 def readBlock(self):
641 Lambda = SPEED_OF_LIGHT / 50000000
642 # 1./(self.dataOut.ippSeconds * self.dataOut.nCohInt)
643 PRF = self.dataOut.PRF
644 Vmax = -Lambda / (4. * (1. / PRF) * self.dataOut.nCohInt * 2.)
645 deltafreq = PRF / (self.nProfiles)
646 deltavel = (Vmax * 2) / (self.nProfiles)
647 freqrange = deltafreq * \
648 (numpy.arange(self.nProfiles) - self.nProfiles / 2.) - deltafreq / 2
649 velrange = deltavel * \
650 (numpy.arange(self.nProfiles) - self.nProfiles / 2.)
651 return velrange
652
653 def readBlock(self):
627 654 '''
628 655 It should be checked if the block has data, if it is not passed to the next file.
629
656
630 657 Then the following is done:
631
658
632 659 1. Read the RecordHeader
633 660 2. Fill the buffer with the current block number.
634
661
635 662 '''
636
637 if self.BlockCounter < self.nFDTdataRecors-2:
663
664 if self.BlockCounter < self.nFDTdataRecors - 2:
638 665 print self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
639 if self.ReadMode==1:
640 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1)
641 elif self.ReadMode==0:
666 if self.ReadMode == 1:
667 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter + 1)
668 elif self.ReadMode == 0:
642 669 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter)
643
644 rheader.RHread(self.fpFile) #Bltr FileHeader Reading
645
646 self.OffsetStartHeader=rheader.OffsetStartHeader
647 self.RecCounter=rheader.RecCounter
648 self.Off2StartNxtRec=rheader.Off2StartNxtRec
649 self.Off2StartData=rheader.Off2StartData
650 self.nProfiles=rheader.nProfiles
651 self.nChannels=rheader.nChannels
652 self.nHeights=rheader.nHeights
653 self.frequency=rheader.TransmitFrec
654 self.DualModeIndex=rheader.DualModeIndex
655
656 self.pairsList =[(0,1),(0,2),(1,2)]
670
671 rheader.RHread(self.fpFile) # Bltr FileHeader Reading
672
673 self.OffsetStartHeader = rheader.OffsetStartHeader
674 self.RecCounter = rheader.RecCounter
675 self.Off2StartNxtRec = rheader.Off2StartNxtRec
676 self.Off2StartData = rheader.Off2StartData
677 self.nProfiles = rheader.nProfiles
678 self.nChannels = rheader.nChannels
679 self.nHeights = rheader.nHeights
680 self.frequency = rheader.TransmitFrec
681 self.DualModeIndex = rheader.DualModeIndex
682
683 self.pairsList = [(0, 1), (0, 2), (1, 2)]
657 684 self.dataOut.pairsList = self.pairsList
658
659 self.nRdPairs=len(self.dataOut.pairsList)
685
686 self.nRdPairs = len(self.dataOut.pairsList)
660 687 self.dataOut.nRdPairs = self.nRdPairs
661
662 self.__firstHeigth=rheader.StartRangeSamp
663 self.__deltaHeigth=rheader.SampResolution
664 self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth
688
689 self.__firstHeigth = rheader.StartRangeSamp
690 self.__deltaHeigth = rheader.SampResolution
691 self.dataOut.heightList = self.__firstHeigth + \
692 numpy.array(range(self.nHeights)) * self.__deltaHeigth
665 693 self.dataOut.channelList = range(self.nChannels)
666 self.dataOut.nProfiles=rheader.nProfiles
667 self.dataOut.nIncohInt=rheader.nIncohInt
668 self.dataOut.nCohInt=rheader.nCohInt
669 self.dataOut.ippSeconds= 1/float(rheader.PRFhz)
670 self.dataOut.PRF=rheader.PRFhz
671 self.dataOut.nFFTPoints=rheader.nProfiles
672 self.dataOut.utctime=rheader.nUtime
673 self.dataOut.timeZone=0
674 self.dataOut.normFactor= self.dataOut.nProfiles*self.dataOut.nIncohInt*self.dataOut.nCohInt
675 self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
676
677 self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN
694 self.dataOut.nProfiles = rheader.nProfiles
695 self.dataOut.nIncohInt = rheader.nIncohInt
696 self.dataOut.nCohInt = rheader.nCohInt
697 self.dataOut.ippSeconds = 1 / float(rheader.PRFhz)
698 self.dataOut.PRF = rheader.PRFhz
699 self.dataOut.nFFTPoints = rheader.nProfiles
700 self.dataOut.utctime = rheader.nUtime
701 self.dataOut.timeZone = 0
702 self.dataOut.normFactor = self.dataOut.nProfiles * \
703 self.dataOut.nIncohInt * self.dataOut.nCohInt
704 self.dataOut.outputInterval = self.dataOut.ippSeconds * \
705 self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
706
707 self.data_output = numpy.ones([3, rheader.nHeights]) * numpy.NaN
678 708 print 'self.data_output', shape(self.data_output)
679 self.dataOut.velocityX=[]
680 self.dataOut.velocityY=[]
681 self.dataOut.velocityV=[]
682
709 self.dataOut.velocityX = []
710 self.dataOut.velocityY = []
711 self.dataOut.velocityV = []
712
683 713 '''Block Reading, the Block Data is received and Reshape is used to give it
684 714 shape.
685 '''
686
687 #Procedure to take the pointer to where the date block starts
688 startDATA = open(self.fpFile,"rb")
689 OffDATA= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec+self.Off2StartData
715 '''
716
717 # Procedure to take the pointer to where the date block starts
718 startDATA = open(self.fpFile, "rb")
719 OffDATA = self.OffsetStartHeader + self.RecCounter * \
720 self.Off2StartNxtRec + self.Off2StartData
690 721 startDATA.seek(OffDATA, os.SEEK_SET)
691
722
692 723 def moving_average(x, N=2):
693 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
694
695 def gaus(xSamples,a,x0,sigma):
696 return a*exp(-(xSamples-x0)**2/(2*sigma**2))
697
698 def Find(x,value):
724 return numpy.convolve(x, numpy.ones((N,)) / N)[(N - 1):]
725
726 def gaus(xSamples, a, x0, sigma):
727 return a * exp(-(xSamples - x0)**2 / (2 * sigma**2))
728
729 def Find(x, value):
699 730 for index in range(len(x)):
700 if x[index]==value:
701 return index
702
731 if x[index] == value:
732 return index
733
703 734 def pol2cart(rho, phi):
704 735 x = rho * numpy.cos(phi)
705 736 y = rho * numpy.sin(phi)
706 737 return(x, y)
707
708
709
710
711 if self.DualModeIndex==self.ReadMode:
712
713 self.data_fft = numpy.fromfile( startDATA, [('complex','<c8')],self.nProfiles*self.nChannels*self.nHeights )
714
715 self.data_fft=self.data_fft.astype(numpy.dtype('complex'))
716
717 self.data_block=numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles ))
718
719 self.data_block = numpy.transpose(self.data_block, (1,2,0))
720
738
739 if self.DualModeIndex == self.ReadMode:
740
741 self.data_fft = numpy.fromfile(
742 startDATA, [('complex', '<c8')], self.nProfiles * self.nChannels * self.nHeights)
743
744 self.data_fft = self.data_fft.astype(numpy.dtype('complex'))
745
746 self.data_block = numpy.reshape(
747 self.data_fft, (self.nHeights, self.nChannels, self.nProfiles))
748
749 self.data_block = numpy.transpose(self.data_block, (1, 2, 0))
750
721 751 copy = self.data_block.copy()
722 spc = copy * numpy.conjugate(copy)
723
724 self.data_spc = numpy.absolute(spc) # valor absoluto o magnitud
725
752 spc = copy * numpy.conjugate(copy)
753
754 self.data_spc = numpy.absolute(
755 spc) # valor absoluto o magnitud
756
726 757 factor = self.dataOut.normFactor
727
728
729 z = self.data_spc.copy()#/factor
758
759 z = self.data_spc.copy() # /factor
730 760 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
731 761 #zdB = 10*numpy.log10(z)
732 762 print ' '
733 763 print 'Z: '
734 print shape(z)
764 print shape(z)
735 765 print ' '
736 766 print ' '
737
738 self.dataOut.data_spc=self.data_spc
739
740 self.noise = self.dataOut.getNoise(ymin_index=80, ymax_index=132)#/factor
767
768 self.dataOut.data_spc = self.data_spc
769
770 self.noise = self.dataOut.getNoise(
771 ymin_index=80, ymax_index=132) # /factor
741 772 #noisedB = 10*numpy.log10(self.noise)
742
743
744 ySamples=numpy.ones([3,self.nProfiles])
745 phase=numpy.ones([3,self.nProfiles])
746 CSPCSamples=numpy.ones([3,self.nProfiles],dtype=numpy.complex_)
747 coherence=numpy.ones([3,self.nProfiles])
748 PhaseSlope=numpy.ones(3)
749 PhaseInter=numpy.ones(3)
750
773
774 ySamples = numpy.ones([3, self.nProfiles])
775 phase = numpy.ones([3, self.nProfiles])
776 CSPCSamples = numpy.ones(
777 [3, self.nProfiles], dtype=numpy.complex_)
778 coherence = numpy.ones([3, self.nProfiles])
779 PhaseSlope = numpy.ones(3)
780 PhaseInter = numpy.ones(3)
781
751 782 '''****** Getting CrossSpectra ******'''
752 cspc=self.data_block.copy()
753 self.data_cspc=self.data_block.copy()
754
755 xFrec=self.getVelRange(1)
756 VelRange=self.getVelRange(1)
757 self.dataOut.VelRange=VelRange
758 #print ' '
759 #print ' '
760 #print 'xFrec',xFrec
761 #print ' '
762 #print ' '
763 #Height=35
764 for i in range(self.nRdPairs):
765
783 cspc = self.data_block.copy()
784 self.data_cspc = self.data_block.copy()
785
786 xFrec = self.getVelRange(1)
787 VelRange = self.getVelRange(1)
788 self.dataOut.VelRange = VelRange
789 # print ' '
790 # print ' '
791 # print 'xFrec',xFrec
792 # print ' '
793 # print ' '
794 # Height=35
795 for i in range(self.nRdPairs):
796
766 797 chan_index0 = self.dataOut.pairsList[i][0]
767 798 chan_index1 = self.dataOut.pairsList[i][1]
768
769 self.data_cspc[i,:,:]=cspc[chan_index0,:,:] * numpy.conjugate(cspc[chan_index1,:,:])
770
771
799
800 self.data_cspc[i, :, :] = cspc[chan_index0, :,
801 :] * numpy.conjugate(cspc[chan_index1, :, :])
802
772 803 '''Getting Eij and Nij'''
773 (AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180)
774 (AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180)
775 (AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180)
776
777 E01=AntennaX0-AntennaX1
778 N01=AntennaY0-AntennaY1
779
780 E02=AntennaX0-AntennaX2
781 N02=AntennaY0-AntennaY2
782
783 E12=AntennaX1-AntennaX2
784 N12=AntennaY1-AntennaY2
785
786 self.ChanDist= numpy.array([[E01, N01],[E02,N02],[E12,N12]])
787
804 (AntennaX0, AntennaY0) = pol2cart(
805 rheader.AntennaCoord0, rheader.AntennaAngl0 * numpy.pi / 180)
806 (AntennaX1, AntennaY1) = pol2cart(
807 rheader.AntennaCoord1, rheader.AntennaAngl1 * numpy.pi / 180)
808 (AntennaX2, AntennaY2) = pol2cart(
809 rheader.AntennaCoord2, rheader.AntennaAngl2 * numpy.pi / 180)
810
811 E01 = AntennaX0 - AntennaX1
812 N01 = AntennaY0 - AntennaY1
813
814 E02 = AntennaX0 - AntennaX2
815 N02 = AntennaY0 - AntennaY2
816
817 E12 = AntennaX1 - AntennaX2
818 N12 = AntennaY1 - AntennaY2
819
820 self.ChanDist = numpy.array(
821 [[E01, N01], [E02, N02], [E12, N12]])
822
788 823 self.dataOut.ChanDist = self.ChanDist
789
790
824
825
791 826 # for Height in range(self.nHeights):
792 #
827 #
793 828 # for i in range(self.nRdPairs):
794 #
829 #
795 830 # '''****** Line of Data SPC ******'''
796 831 # zline=z[i,:,Height]
797 #
832 #
798 833 # '''****** DC is removed ******'''
799 834 # DC=Find(zline,numpy.amax(zline))
800 835 # zline[DC]=(zline[DC-1]+zline[DC+1])/2
801 #
802 #
836 #
837 #
803 838 # '''****** SPC is normalized ******'''
804 839 # FactNorm= zline.copy() / numpy.sum(zline.copy())
805 840 # FactNorm= FactNorm/numpy.sum(FactNorm)
806 #
841 #
807 842 # SmoothSPC=moving_average(FactNorm,N=3)
808 #
843 #
809 844 # xSamples = ar(range(len(SmoothSPC)))
810 845 # ySamples[i] = SmoothSPC-self.noise[i]
811 #
846 #
812 847 # for i in range(self.nRdPairs):
813 #
848 #
814 849 # '''****** Line of Data CSPC ******'''
815 850 # cspcLine=self.data_cspc[i,:,Height].copy()
816 #
817 #
818 #
851 #
852 #
853 #
819 854 # '''****** CSPC is normalized ******'''
820 855 # chan_index0 = self.dataOut.pairsList[i][0]
821 # chan_index1 = self.dataOut.pairsList[i][1]
856 # chan_index1 = self.dataOut.pairsList[i][1]
822 857 # CSPCFactor= numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])
823 #
824 #
858 #
859 #
825 860 # CSPCNorm= cspcLine.copy() / numpy.sqrt(CSPCFactor)
826 #
827 #
861 #
862 #
828 863 # CSPCSamples[i] = CSPCNorm-self.noise[i]
829 864 # coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
830 #
865 #
831 866 # '''****** DC is removed ******'''
832 867 # DC=Find(coherence[i],numpy.amax(coherence[i]))
833 868 # coherence[i][DC]=(coherence[i][DC-1]+coherence[i][DC+1])/2
834 869 # coherence[i]= moving_average(coherence[i],N=2)
835 #
870 #
836 871 # phase[i] = moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
837 #
838 #
872 #
873 #
839 874 # '''****** Getting fij width ******'''
840 #
875 #
841 876 # yMean=[]
842 # yMean2=[]
843 #
877 # yMean2=[]
878 #
844 879 # for j in range(len(ySamples[1])):
845 880 # yMean=numpy.append(yMean,numpy.average([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
846 #
881 #
847 882 # '''******* Getting fitting Gaussian ******'''
848 883 # meanGauss=sum(xSamples*yMean) / len(xSamples)
849 884 # sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
850 885 # #print 'Height',Height,'SNR', meanGauss/sigma**2
851 #
886 #
852 887 # if (abs(meanGauss/sigma**2) > 0.0001) :
853 #
854 # try:
888 #
889 # try:
855 890 # popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
856 #
891 #
857 892 # if numpy.amax(popt)>numpy.amax(yMean)*0.3:
858 893 # FitGauss=gaus(xSamples,*popt)
859 #
860 # else:
894 #
895 # else:
861 896 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
862 897 # print 'Verificador: Dentro', Height
863 898 # except RuntimeError:
864 #
899 #
865 900 # try:
866 901 # for j in range(len(ySamples[1])):
867 902 # yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]]))
868 903 # popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma])
869 904 # FitGauss=gaus(xSamples,*popt)
870 905 # print 'Verificador: Exepcion1', Height
871 906 # except RuntimeError:
872 #
907 #
873 908 # try:
874 909 # popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma])
875 910 # FitGauss=gaus(xSamples,*popt)
876 911 # print 'Verificador: Exepcion2', Height
877 912 # except RuntimeError:
878 913 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
879 914 # print 'Verificador: Exepcion3', Height
880 915 # else:
881 916 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
882 917 # #print 'Verificador: Fuera', Height
883 #
884 #
885 #
918 #
919 #
920 #
886 921 # Maximun=numpy.amax(yMean)
887 922 # eMinus1=Maximun*numpy.exp(-1)
888 #
923 #
889 924 # HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
890 925 # HalfWidth= xFrec[HWpos]
891 926 # GCpos=Find(FitGauss, numpy.amax(FitGauss))
892 927 # Vpos=Find(FactNorm, numpy.amax(FactNorm))
893 928 # #Vpos=numpy.sum(FactNorm)/len(FactNorm)
894 929 # #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) )))
895 930 # #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos
896 931 # '''****** Getting Fij ******'''
897 #
932 #
898 933 # GaussCenter=xFrec[GCpos]
899 934 # if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
900 935 # Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
901 936 # else:
902 937 # Fij=abs(GaussCenter-HalfWidth)+0.0000001
903 #
938 #
904 939 # '''****** Getting Frecuency range of significant data ******'''
905 #
940 #
906 941 # Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
907 #
942 #
908 943 # if Rangpos<GCpos:
909 944 # Range=numpy.array([Rangpos,2*GCpos-Rangpos])
910 945 # else:
911 946 # Range=numpy.array([2*GCpos-Rangpos,Rangpos])
912 #
947 #
913 948 # FrecRange=xFrec[Range[0]:Range[1]]
914 #
949 #
915 950 # #print 'FrecRange', FrecRange
916 951 # '''****** Getting SCPC Slope ******'''
917 #
952 #
918 953 # for i in range(self.nRdPairs):
919 #
954 #
920 955 # if len(FrecRange)>5 and len(FrecRange)<self.nProfiles*0.5:
921 # PhaseRange=moving_average(phase[i,Range[0]:Range[1]],N=3)
922 #
956 # PhaseRange=moving_average(phase[i,Range[0]:Range[1]],N=3)
957 #
923 958 # slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
924 959 # PhaseSlope[i]=slope
925 960 # PhaseInter[i]=intercept
926 961 # else:
927 962 # PhaseSlope[i]=0
928 963 # PhaseInter[i]=0
929 #
964 #
930 965 # # plt.figure(i+15)
931 966 # # plt.title('FASE ( CH%s*CH%s )' %(self.dataOut.pairsList[i][0],self.dataOut.pairsList[i][1]))
932 967 # # plt.xlabel('Frecuencia (KHz)')
933 968 # # plt.ylabel('Magnitud')
934 969 # # #plt.subplot(311+i)
935 970 # # plt.plot(FrecRange,PhaseRange,'b')
936 971 # # plt.plot(FrecRange,FrecRange*PhaseSlope[i]+PhaseInter[i],'r')
937 #
972 #
938 973 # #plt.axis([-0.6, 0.2, -3.2, 3.2])
939 #
940 #
974 #
975 #
941 976 # '''Getting constant C'''
942 977 # cC=(Fij*numpy.pi)**2
943 #
978 #
944 979 # # '''Getting Eij and Nij'''
945 980 # # (AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180)
946 981 # # (AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180)
947 982 # # (AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180)
948 # #
983 # #
949 984 # # E01=AntennaX0-AntennaX1
950 985 # # N01=AntennaY0-AntennaY1
951 # #
986 # #
952 987 # # E02=AntennaX0-AntennaX2
953 988 # # N02=AntennaY0-AntennaY2
954 # #
989 # #
955 990 # # E12=AntennaX1-AntennaX2
956 991 # # N12=AntennaY1-AntennaY2
957 #
992 #
958 993 # '''****** Getting constants F and G ******'''
959 994 # MijEijNij=numpy.array([[E02,N02], [E12,N12]])
960 995 # MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
961 # MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
996 # MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
962 997 # MijResults=numpy.array([MijResult0,MijResult1])
963 998 # (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
964 #
999 #
965 1000 # '''****** Getting constants A, B and H ******'''
966 1001 # W01=numpy.amax(coherence[0])
967 1002 # W02=numpy.amax(coherence[1])
968 1003 # W12=numpy.amax(coherence[2])
969 #
1004 #
970 1005 # WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
971 1006 # WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
972 1007 # WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
973 #
1008 #
974 1009 # WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
975 #
976 # WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
1010 #
1011 # WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
977 1012 # (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
978 #
1013 #
979 1014 # VxVy=numpy.array([[cA,cH],[cH,cB]])
980 #
1015 #
981 1016 # VxVyResults=numpy.array([-cF,-cG])
982 1017 # (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
983 1018 # Vzon = Vy
984 1019 # Vmer = Vx
985 1020 # Vmag=numpy.sqrt(Vzon**2+Vmer**2)
986 1021 # Vang=numpy.arctan2(Vmer,Vzon)
987 #
1022 #
988 1023 # if abs(Vy)<100 and abs(Vy)> 0.:
989 1024 # self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag
990 1025 # #print 'Vmag',Vmag
991 1026 # else:
992 1027 # self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN)
993 #
1028 #
994 1029 # if abs(Vx)<100 and abs(Vx) > 0.:
995 1030 # self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang
996 # #print 'Vang',Vang
1031 # #print 'Vang',Vang
997 1032 # else:
998 1033 # self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN)
999 #
1034 #
1000 1035 # if abs(GaussCenter)<2:
1001 1036 # self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos])
1002 #
1037 #
1003 1038 # else:
1004 1039 # self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN)
1005 #
1006 #
1040 #
1041 #
1007 1042 # # print '********************************************'
1008 1043 # # print 'HalfWidth ', HalfWidth
1009 1044 # # print 'Maximun ', Maximun
1010 1045 # # print 'eMinus1 ', eMinus1
1011 1046 # # print 'Rangpos ', Rangpos
1012 1047 # # print 'GaussCenter ',GaussCenter
1013 1048 # # print 'E01 ',E01
1014 1049 # # print 'N01 ',N01
1015 1050 # # print 'E02 ',E02
1016 1051 # # print 'N02 ',N02
1017 1052 # # print 'E12 ',E12
1018 1053 # # print 'N12 ',N12
1019 1054 # #print 'self.dataOut.velocityX ', self.dataOut.velocityX
1020 1055 # # print 'Fij ', Fij
1021 1056 # # print 'cC ', cC
1022 1057 # # print 'cF ', cF
1023 1058 # # print 'cG ', cG
1024 1059 # # print 'cA ', cA
1025 1060 # # print 'cB ', cB
1026 1061 # # print 'cH ', cH
1027 1062 # # print 'Vx ', Vx
1028 1063 # # print 'Vy ', Vy
1029 1064 # # print 'Vmag ', Vmag
1030 1065 # # print 'Vang ', Vang*180/numpy.pi
1031 1066 # # print 'PhaseSlope ',PhaseSlope[0]
1032 1067 # # print 'PhaseSlope ',PhaseSlope[1]
1033 1068 # # print 'PhaseSlope ',PhaseSlope[2]
1034 1069 # # print '********************************************'
1035 1070 # #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY)
1036 #
1071 #
1037 1072 # #print 'self.dataOut.velocityX', len(self.dataOut.velocityX)
1038 1073 # #print 'self.dataOut.velocityY', len(self.dataOut.velocityY)
1039 1074 # #print 'self.dataOut.velocityV', self.dataOut.velocityV
1040 #
1075 #
1041 1076 # self.data_output[0]=numpy.array(self.dataOut.velocityX)
1042 1077 # self.data_output[1]=numpy.array(self.dataOut.velocityY)
1043 1078 # self.data_output[2]=numpy.array(self.dataOut.velocityV)
1044 #
1079 #
1045 1080 # prin= self.data_output[0][~numpy.isnan(self.data_output[0])]
1046 1081 # print ' '
1047 # print 'VmagAverage',numpy.mean(prin)
1082 # print 'VmagAverage',numpy.mean(prin)
1048 1083 # print ' '
1049 1084 # # plt.figure(5)
1050 1085 # # plt.subplot(211)
1051 1086 # # plt.plot(self.dataOut.velocityX,'yo:')
1052 1087 # # plt.subplot(212)
1053 1088 # # plt.plot(self.dataOut.velocityY,'yo:')
1054 #
1089 #
1055 1090 # # plt.figure(1)
1056 1091 # # # plt.subplot(121)
1057 1092 # # # plt.plot(xFrec,ySamples[0],'k',label='Ch0')
1058 1093 # # # plt.plot(xFrec,ySamples[1],'g',label='Ch1')
1059 1094 # # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
1060 1095 # # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
1061 1096 # # # plt.legend()
1062 1097 # # plt.title('DATOS A ALTURA DE 2850 METROS')
1063 # #
1098 # #
1064 1099 # # plt.xlabel('Frecuencia (KHz)')
1065 1100 # # plt.ylabel('Magnitud')
1066 1101 # # # plt.subplot(122)
1067 1102 # # # plt.title('Fit for Time Constant')
1068 1103 # # #plt.plot(xFrec,zline)
1069 1104 # # #plt.plot(xFrec,SmoothSPC,'g')
1070 1105 # # plt.plot(xFrec,FactNorm)
1071 1106 # # plt.axis([-4, 4, 0, 0.15])
1072 1107 # # # plt.xlabel('SelfSpectra KHz')
1073 # #
1108 # #
1074 1109 # # plt.figure(10)
1075 1110 # # # plt.subplot(121)
1076 1111 # # plt.plot(xFrec,ySamples[0],'b',label='Ch0')
1077 1112 # # plt.plot(xFrec,ySamples[1],'y',label='Ch1')
1078 1113 # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
1079 1114 # # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
1080 1115 # # plt.legend()
1081 1116 # # plt.title('SELFSPECTRA EN CANALES')
1082 # #
1117 # #
1083 1118 # # plt.xlabel('Frecuencia (KHz)')
1084 1119 # # plt.ylabel('Magnitud')
1085 1120 # # # plt.subplot(122)
1086 1121 # # # plt.title('Fit for Time Constant')
1087 1122 # # #plt.plot(xFrec,zline)
1088 1123 # # #plt.plot(xFrec,SmoothSPC,'g')
1089 1124 # # # plt.plot(xFrec,FactNorm)
1090 1125 # # # plt.axis([-4, 4, 0, 0.15])
1091 1126 # # # plt.xlabel('SelfSpectra KHz')
1092 # #
1127 # #
1093 1128 # # plt.figure(9)
1094 # #
1095 # #
1129 # #
1130 # #
1096 1131 # # plt.title('DATOS SUAVIZADOS')
1097 1132 # # plt.xlabel('Frecuencia (KHz)')
1098 1133 # # plt.ylabel('Magnitud')
1099 1134 # # plt.plot(xFrec,SmoothSPC,'g')
1100 # #
1135 # #
1101 1136 # # #plt.plot(xFrec,FactNorm)
1102 1137 # # plt.axis([-4, 4, 0, 0.15])
1103 1138 # # # plt.xlabel('SelfSpectra KHz')
1104 # # #
1139 # # #
1105 1140 # # plt.figure(2)
1106 1141 # # # #plt.subplot(121)
1107 1142 # # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra')
1108 1143 # # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano')
1109 1144 # # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo')
1110 1145 # # # #plt.plot(xFrec,phase)
1111 1146 # # # plt.xlabel('Suavizado, promediado KHz')
1112 1147 # # plt.title('SELFSPECTRA PROMEDIADO')
1113 1148 # # # #plt.subplot(122)
1114 1149 # # # #plt.plot(xSamples,zline)
1115 1150 # # plt.xlabel('Frecuencia (KHz)')
1116 1151 # # plt.ylabel('Magnitud')
1117 1152 # # plt.legend()
1118 # # #
1153 # # #
1119 1154 # # # plt.figure(3)
1120 1155 # # # plt.subplot(311)
1121 1156 # # # #plt.plot(xFrec,phase[0])
1122 1157 # # # plt.plot(xFrec,phase[0],'g')
1123 1158 # # # plt.subplot(312)
1124 1159 # # # plt.plot(xFrec,phase[1],'g')
1125 1160 # # # plt.subplot(313)
1126 1161 # # # plt.plot(xFrec,phase[2],'g')
1127 1162 # # # #plt.plot(xFrec,phase[2])
1128 # # #
1163 # # #
1129 1164 # # # plt.figure(4)
1130 # # #
1165 # # #
1131 1166 # # # plt.plot(xSamples,coherence[0],'b')
1132 1167 # # # plt.plot(xSamples,coherence[1],'r')
1133 1168 # # # plt.plot(xSamples,coherence[2],'g')
1134 1169 # # plt.show()
1135 # # #
1170 # # #
1136 1171 # # # plt.clf()
1137 1172 # # # plt.cla()
1138 # # # plt.close()
1139 #
1140 # print ' '
1141
1142
1143
1144 self.BlockCounter+=2
1145
1173 # # # plt.close()
1174 #
1175 # print ' '
1176
1177 self.BlockCounter += 2
1178
1146 1179 else:
1147 self.fileSelector+=1
1148 self.BlockCounter=0
1180 self.fileSelector += 1
1181 self.BlockCounter = 0
1149 1182 print "Next File"
1150
1151
1152
1153
1154
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1027 lines changed) Show them Hide them
@@ -1,803 +1,802
1 import os, sys
1 import os
2 import sys
2 3 import glob
3 4 import fnmatch
4 5 import datetime
5 6 import time
6 7 import re
7 8 import h5py
8 9 import numpy
9 import matplotlib.pyplot as plt
10 10
11 import pylab as plb
12 11 from scipy.optimize import curve_fit
13 from scipy import asarray as ar,exp
12 from scipy import asarray as ar, exp
14 13 from scipy import stats
15 14
16 15 from numpy.ma.core import getdata
17 16
18 17 SPEED_OF_LIGHT = 299792458
19 18 SPEED_OF_LIGHT = 3e8
20 19
21 20 try:
22 21 from gevent import sleep
23 22 except:
24 23 from time import sleep
25 24
26 25 from schainpy.model.data.jrodata import Spectra
27 26 #from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader
28 27 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
29 28 #from schainpy.model.io.jroIO_bltr import BLTRReader
30 29 from numpy import imag, shape, NaN, empty
31 30
32 31
33
34 32 class Header(object):
35
33
36 34 def __init__(self):
37 35 raise NotImplementedError
38
39
36
40 37 def read(self):
41
38
42 39 raise NotImplementedError
43
40
44 41 def write(self):
45
42
46 43 raise NotImplementedError
47
44
48 45 def printInfo(self):
49
50 message = "#"*50 + "\n"
46
47 message = "#" * 50 + "\n"
51 48 message += self.__class__.__name__.upper() + "\n"
52 message += "#"*50 + "\n"
53
49 message += "#" * 50 + "\n"
50
54 51 keyList = self.__dict__.keys()
55 52 keyList.sort()
56
53
57 54 for key in keyList:
58 message += "%s = %s" %(key, self.__dict__[key]) + "\n"
59
55 message += "%s = %s" % (key, self.__dict__[key]) + "\n"
56
60 57 if "size" not in keyList:
61 58 attr = getattr(self, "size")
62
59
63 60 if attr:
64 message += "%s = %s" %("size", attr) + "\n"
65
66 #print message
67
68
69 FILE_HEADER = numpy.dtype([ #HEADER 1024bytes
70 ('Hname','a32'), #Original file name
71 ('Htime',numpy.str_,32), #Date and time when the file was created
72 ('Hoper',numpy.str_,64), #Name of operator who created the file
73 ('Hplace',numpy.str_,128), #Place where the measurements was carried out
74 ('Hdescr',numpy.str_,256), #Description of measurements
75 ('Hdummy',numpy.str_,512), #Reserved space
76 #Main chunk 8bytes
77 ('Msign',numpy.str_,4), #Main chunk signature FZKF or NUIG
78 ('MsizeData','<i4'), #Size of data block main chunk
79 #Processing DSP parameters 36bytes
80 ('PPARsign',numpy.str_,4), #PPAR signature
81 ('PPARsize','<i4'), #PPAR size of block
82 ('PPARprf','<i4'), #Pulse repetition frequency
83 ('PPARpdr','<i4'), #Pulse duration
84 ('PPARsft','<i4'), #FFT length
85 ('PPARavc','<i4'), #Number of spectral (in-coherent) averages
86 ('PPARihp','<i4'), #Number of lowest range gate for moment estimation
87 ('PPARchg','<i4'), #Count for gates for moment estimation
88 ('PPARpol','<i4'), #switch on/off polarimetric measurements. Should be 1.
89 #Service DSP parameters 112bytes
90 ('SPARatt','<i4'), #STC attenuation on the lowest ranges on/off
91 ('SPARtx','<i4'), #OBSOLETE
92 ('SPARaddGain0','<f4'), #OBSOLETE
93 ('SPARaddGain1','<f4'), #OBSOLETE
94 ('SPARwnd','<i4'), #Debug only. It normal mode it is 0.
95 ('SPARpos','<i4'), #Delay between sync pulse and tx pulse for phase corr, ns
96 ('SPARadd','<i4'), #"add to pulse" to compensate for delay between the leading edge of driver pulse and envelope of the RF signal.
97 ('SPARlen','<i4'), #Time for measuring txn pulse phase. OBSOLETE
98 ('SPARcal','<i4'), #OBSOLETE
99 ('SPARnos','<i4'), #OBSOLETE
100 ('SPARof0','<i4'), #detection threshold
101 ('SPARof1','<i4'), #OBSOLETE
102 ('SPARswt','<i4'), #2nd moment estimation threshold
103 ('SPARsum','<i4'), #OBSOLETE
104 ('SPARosc','<i4'), #flag Oscillosgram mode
105 ('SPARtst','<i4'), #OBSOLETE
106 ('SPARcor','<i4'), #OBSOLETE
107 ('SPARofs','<i4'), #OBSOLETE
108 ('SPARhsn','<i4'), #Hildebrand div noise detection on noise gate
109 ('SPARhsa','<f4'), #Hildebrand div noise detection on all gates
110 ('SPARcalibPow_M','<f4'), #OBSOLETE
111 ('SPARcalibSNR_M','<f4'), #OBSOLETE
112 ('SPARcalibPow_S','<f4'), #OBSOLETE
113 ('SPARcalibSNR_S','<f4'), #OBSOLETE
114 ('SPARrawGate1','<i4'), #Lowest range gate for spectra saving Raw_Gate1 >=5
115 ('SPARrawGate2','<i4'), #Number of range gates with atmospheric signal
116 ('SPARraw','<i4'), #flag - IQ or spectra saving on/off
117 ('SPARprc','<i4'),]) #flag - Moment estimation switched on/off
118
119
120
61 message += "%s = %s" % ("size", attr) + "\n"
62
63 # print message
64
65
66 FILE_HEADER = numpy.dtype([ # HEADER 1024bytes
67 ('Hname', 'a32'), # Original file name
68 # Date and time when the file was created
69 ('Htime', numpy.str_, 32),
70 # Name of operator who created the file
71 ('Hoper', numpy.str_, 64),
72 # Place where the measurements was carried out
73 ('Hplace', numpy.str_, 128),
74 # Description of measurements
75 ('Hdescr', numpy.str_, 256),
76 ('Hdummy', numpy.str_, 512), # Reserved space
77 # Main chunk 8bytes
78 # Main chunk signature FZKF or NUIG
79 ('Msign', numpy.str_, 4),
80 ('MsizeData', '<i4'), # Size of data block main chunk
81 # Processing DSP parameters 36bytes
82 ('PPARsign', numpy.str_, 4), # PPAR signature
83 ('PPARsize', '<i4'), # PPAR size of block
84 ('PPARprf', '<i4'), # Pulse repetition frequency
85 ('PPARpdr', '<i4'), # Pulse duration
86 ('PPARsft', '<i4'), # FFT length
87 # Number of spectral (in-coherent) averages
88 ('PPARavc', '<i4'),
89 # Number of lowest range gate for moment estimation
90 ('PPARihp', '<i4'),
91 # Count for gates for moment estimation
92 ('PPARchg', '<i4'),
93 # switch on/off polarimetric measurements. Should be 1.
94 ('PPARpol', '<i4'),
95 # Service DSP parameters 112bytes
96 # STC attenuation on the lowest ranges on/off
97 ('SPARatt', '<i4'),
98 ('SPARtx', '<i4'), # OBSOLETE
99 ('SPARaddGain0', '<f4'), # OBSOLETE
100 ('SPARaddGain1', '<f4'), # OBSOLETE
101 # Debug only. It normal mode it is 0.
102 ('SPARwnd', '<i4'),
103 # Delay between sync pulse and tx pulse for phase corr, ns
104 ('SPARpos', '<i4'),
105 # "add to pulse" to compensate for delay between the leading edge of driver pulse and envelope of the RF signal.
106 ('SPARadd', '<i4'),
107 # Time for measuring txn pulse phase. OBSOLETE
108 ('SPARlen', '<i4'),
109 ('SPARcal', '<i4'), # OBSOLETE
110 ('SPARnos', '<i4'), # OBSOLETE
111 ('SPARof0', '<i4'), # detection threshold
112 ('SPARof1', '<i4'), # OBSOLETE
113 ('SPARswt', '<i4'), # 2nd moment estimation threshold
114 ('SPARsum', '<i4'), # OBSOLETE
115 ('SPARosc', '<i4'), # flag Oscillosgram mode
116 ('SPARtst', '<i4'), # OBSOLETE
117 ('SPARcor', '<i4'), # OBSOLETE
118 ('SPARofs', '<i4'), # OBSOLETE
119 # Hildebrand div noise detection on noise gate
120 ('SPARhsn', '<i4'),
121 # Hildebrand div noise detection on all gates
122 ('SPARhsa', '<f4'),
123 ('SPARcalibPow_M', '<f4'), # OBSOLETE
124 ('SPARcalibSNR_M', '<f4'), # OBSOLETE
125 ('SPARcalibPow_S', '<f4'), # OBSOLETE
126 ('SPARcalibSNR_S', '<f4'), # OBSOLETE
127 # Lowest range gate for spectra saving Raw_Gate1 >=5
128 ('SPARrawGate1', '<i4'),
129 # Number of range gates with atmospheric signal
130 ('SPARrawGate2', '<i4'),
131 # flag - IQ or spectra saving on/off
132 ('SPARraw', '<i4'),
133 ('SPARprc', '<i4'), ]) # flag - Moment estimation switched on/off
134
135
121 136 class FileHeaderMIRA35c(Header):
122
137
123 138 def __init__(self):
124
125 self.Hname= None
126 self.Htime= None
127 self.Hoper= None
128 self.Hplace= None
129 self.Hdescr= None
130 self.Hdummy= None
131
132 self.Msign=None
133 self.MsizeData=None
134
135 self.PPARsign=None
136 self.PPARsize=None
137 self.PPARprf=None
138 self.PPARpdr=None
139 self.PPARsft=None
140 self.PPARavc=None
141 self.PPARihp=None
142 self.PPARchg=None
143 self.PPARpol=None
144 #Service DSP parameters
145 self.SPARatt=None
146 self.SPARtx=None
147 self.SPARaddGain0=None
148 self.SPARaddGain1=None
149 self.SPARwnd=None
150 self.SPARpos=None
151 self.SPARadd=None
152 self.SPARlen=None
153 self.SPARcal=None
154 self.SPARnos=None
155 self.SPARof0=None
156 self.SPARof1=None
157 self.SPARswt=None
158 self.SPARsum=None
159 self.SPARosc=None
160 self.SPARtst=None
161 self.SPARcor=None
162 self.SPARofs=None
163 self.SPARhsn=None
164 self.SPARhsa=None
165 self.SPARcalibPow_M=None
166 self.SPARcalibSNR_M=None
167 self.SPARcalibPow_S=None
168 self.SPARcalibSNR_S=None
169 self.SPARrawGate1=None
170 self.SPARrawGate2=None
171 self.SPARraw=None
172 self.SPARprc=None
173
174 self.FHsize=1180
175
139
140 self.Hname = None
141 self.Htime = None
142 self.Hoper = None
143 self.Hplace = None
144 self.Hdescr = None
145 self.Hdummy = None
146
147 self.Msign = None
148 self.MsizeData = None
149
150 self.PPARsign = None
151 self.PPARsize = None
152 self.PPARprf = None
153 self.PPARpdr = None
154 self.PPARsft = None
155 self.PPARavc = None
156 self.PPARihp = None
157 self.PPARchg = None
158 self.PPARpol = None
159 # Service DSP parameters
160 self.SPARatt = None
161 self.SPARtx = None
162 self.SPARaddGain0 = None
163 self.SPARaddGain1 = None
164 self.SPARwnd = None
165 self.SPARpos = None
166 self.SPARadd = None
167 self.SPARlen = None
168 self.SPARcal = None
169 self.SPARnos = None
170 self.SPARof0 = None
171 self.SPARof1 = None
172 self.SPARswt = None
173 self.SPARsum = None
174 self.SPARosc = None
175 self.SPARtst = None
176 self.SPARcor = None
177 self.SPARofs = None
178 self.SPARhsn = None
179 self.SPARhsa = None
180 self.SPARcalibPow_M = None
181 self.SPARcalibSNR_M = None
182 self.SPARcalibPow_S = None
183 self.SPARcalibSNR_S = None
184 self.SPARrawGate1 = None
185 self.SPARrawGate2 = None
186 self.SPARraw = None
187 self.SPARprc = None
188
189 self.FHsize = 1180
190
176 191 def FHread(self, fp):
177
178 header = numpy.fromfile(fp, FILE_HEADER,1)
192
193 header = numpy.fromfile(fp, FILE_HEADER, 1)
179 194 ''' numpy.fromfile(file, dtype, count, sep='')
180 195 file : file or str
181 196 Open file object or filename.
182 197
183 198 dtype : data-type
184 199 Data type of the returned array. For binary files, it is used to determine
185 200 the size and byte-order of the items in the file.
186 201
187 202 count : int
188 203 Number of items to read. -1 means all items (i.e., the complete file).
189 204
190 205 sep : str
191 206 Separator between items if file is a text file. Empty ("") separator means
192 207 the file should be treated as binary. Spaces (" ") in the separator match zero
193 208 or more whitespace characters. A separator consisting only of spaces must match
194 209 at least one whitespace.
195 210
196 211 '''
197
198
199 self.Hname= str(header['Hname'][0])
200 self.Htime= str(header['Htime'][0])
201 self.Hoper= str(header['Hoper'][0])
202 self.Hplace= str(header['Hplace'][0])
203 self.Hdescr= str(header['Hdescr'][0])
204 self.Hdummy= str(header['Hdummy'][0])
205 #1024
206
207 self.Msign=str(header['Msign'][0])
208 self.MsizeData=header['MsizeData'][0]
209 #8
210
211 self.PPARsign=str(header['PPARsign'][0])
212 self.PPARsize=header['PPARsize'][0]
213 self.PPARprf=header['PPARprf'][0]
214 self.PPARpdr=header['PPARpdr'][0]
215 self.PPARsft=header['PPARsft'][0]
216 self.PPARavc=header['PPARavc'][0]
217 self.PPARihp=header['PPARihp'][0]
218 self.PPARchg=header['PPARchg'][0]
219 self.PPARpol=header['PPARpol'][0]
220 #Service DSP parameters
221 #36
222
223 self.SPARatt=header['SPARatt'][0]
224 self.SPARtx=header['SPARtx'][0]
225 self.SPARaddGain0=header['SPARaddGain0'][0]
226 self.SPARaddGain1=header['SPARaddGain1'][0]
227 self.SPARwnd=header['SPARwnd'][0]
228 self.SPARpos=header['SPARpos'][0]
229 self.SPARadd=header['SPARadd'][0]
230 self.SPARlen=header['SPARlen'][0]
231 self.SPARcal=header['SPARcal'][0]
232 self.SPARnos=header['SPARnos'][0]
233 self.SPARof0=header['SPARof0'][0]
234 self.SPARof1=header['SPARof1'][0]
235 self.SPARswt=header['SPARswt'][0]
236 self.SPARsum=header['SPARsum'][0]
237 self.SPARosc=header['SPARosc'][0]
238 self.SPARtst=header['SPARtst'][0]
239 self.SPARcor=header['SPARcor'][0]
240 self.SPARofs=header['SPARofs'][0]
241 self.SPARhsn=header['SPARhsn'][0]
242 self.SPARhsa=header['SPARhsa'][0]
243 self.SPARcalibPow_M=header['SPARcalibPow_M'][0]
244 self.SPARcalibSNR_M=header['SPARcalibSNR_M'][0]
245 self.SPARcalibPow_S=header['SPARcalibPow_S'][0]
246 self.SPARcalibSNR_S=header['SPARcalibSNR_S'][0]
247 self.SPARrawGate1=header['SPARrawGate1'][0]
248 self.SPARrawGate2=header['SPARrawGate2'][0]
249 self.SPARraw=header['SPARraw'][0]
250 self.SPARprc=header['SPARprc'][0]
251 #112
252 #1180
253 #print 'Pointer fp header', fp.tell()
254 #print ' '
255 #print 'SPARrawGate'
256 #print self.SPARrawGate2 - self.SPARrawGate1
257
258 #print ' '
259 #print 'Hname'
260 #print self.Hname
261
262 #print ' '
263 #print 'Msign'
264 #print self.Msign
265
212
213 self.Hname = str(header['Hname'][0])
214 self.Htime = str(header['Htime'][0])
215 self.Hoper = str(header['Hoper'][0])
216 self.Hplace = str(header['Hplace'][0])
217 self.Hdescr = str(header['Hdescr'][0])
218 self.Hdummy = str(header['Hdummy'][0])
219 # 1024
220
221 self.Msign = str(header['Msign'][0])
222 self.MsizeData = header['MsizeData'][0]
223 # 8
224
225 self.PPARsign = str(header['PPARsign'][0])
226 self.PPARsize = header['PPARsize'][0]
227 self.PPARprf = header['PPARprf'][0]
228 self.PPARpdr = header['PPARpdr'][0]
229 self.PPARsft = header['PPARsft'][0]
230 self.PPARavc = header['PPARavc'][0]
231 self.PPARihp = header['PPARihp'][0]
232 self.PPARchg = header['PPARchg'][0]
233 self.PPARpol = header['PPARpol'][0]
234 # Service DSP parameters
235 # 36
236
237 self.SPARatt = header['SPARatt'][0]
238 self.SPARtx = header['SPARtx'][0]
239 self.SPARaddGain0 = header['SPARaddGain0'][0]
240 self.SPARaddGain1 = header['SPARaddGain1'][0]
241 self.SPARwnd = header['SPARwnd'][0]
242 self.SPARpos = header['SPARpos'][0]
243 self.SPARadd = header['SPARadd'][0]
244 self.SPARlen = header['SPARlen'][0]
245 self.SPARcal = header['SPARcal'][0]
246 self.SPARnos = header['SPARnos'][0]
247 self.SPARof0 = header['SPARof0'][0]
248 self.SPARof1 = header['SPARof1'][0]
249 self.SPARswt = header['SPARswt'][0]
250 self.SPARsum = header['SPARsum'][0]
251 self.SPARosc = header['SPARosc'][0]
252 self.SPARtst = header['SPARtst'][0]
253 self.SPARcor = header['SPARcor'][0]
254 self.SPARofs = header['SPARofs'][0]
255 self.SPARhsn = header['SPARhsn'][0]
256 self.SPARhsa = header['SPARhsa'][0]
257 self.SPARcalibPow_M = header['SPARcalibPow_M'][0]
258 self.SPARcalibSNR_M = header['SPARcalibSNR_M'][0]
259 self.SPARcalibPow_S = header['SPARcalibPow_S'][0]
260 self.SPARcalibSNR_S = header['SPARcalibSNR_S'][0]
261 self.SPARrawGate1 = header['SPARrawGate1'][0]
262 self.SPARrawGate2 = header['SPARrawGate2'][0]
263 self.SPARraw = header['SPARraw'][0]
264 self.SPARprc = header['SPARprc'][0]
265 # 112
266 # 1180
267 # print 'Pointer fp header', fp.tell()
268 # print ' '
269 # print 'SPARrawGate'
270 # print self.SPARrawGate2 - self.SPARrawGate1
271
272 # print ' '
273 # print 'Hname'
274 # print self.Hname
275
276 # print ' '
277 # print 'Msign'
278 # print self.Msign
279
266 280 def write(self, fp):
267
281
268 282 headerTuple = (self.Hname,
269 283 self.Htime,
270 284 self.Hoper,
271 285 self.Hplace,
272 286 self.Hdescr,
273 287 self.Hdummy)
274
275
288
276 289 header = numpy.array(headerTuple, FILE_HEADER)
277 290 # numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0)
278 291 header.tofile(fp)
279 292 ''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
280 293
281 294 fid : file or str
282 295 An open file object, or a string containing a filename.
283 296
284 297 sep : str
285 298 Separator between array items for text output. If "" (empty), a binary file is written,
286 299 equivalent to file.write(a.tobytes()).
287 300
288 301 format : str
289 302 Format string for text file output. Each entry in the array is formatted to text by
290 303 first converting it to the closest Python type, and then using "format" % item.
291 304
292 305 '''
293
306
294 307 return 1
295 308
296 SRVI_HEADER = numpy.dtype([
297 ('SignatureSRVI1',numpy.str_,4),#
298 ('SizeOfDataBlock1','<i4'),#
299 ('DataBlockTitleSRVI1',numpy.str_,4),#
300 ('SizeOfSRVI1','<i4'),])#
301
302 class SRVIHeader(Header):
309
310 SRVI_HEADER = numpy.dtype([
311 ('SignatureSRVI1', numpy.str_, 4),
312 ('SizeOfDataBlock1', '<i4'),
313 ('DataBlockTitleSRVI1', numpy.str_, 4),
314 ('SizeOfSRVI1', '<i4'), ])
315
316
317 class SRVIHeader(Header):
303 318 def __init__(self, SignatureSRVI1=0, SizeOfDataBlock1=0, DataBlockTitleSRVI1=0, SizeOfSRVI1=0):
304
319
305 320 self.SignatureSRVI1 = SignatureSRVI1
306 321 self.SizeOfDataBlock1 = SizeOfDataBlock1
307 322 self.DataBlockTitleSRVI1 = DataBlockTitleSRVI1
308 self.SizeOfSRVI1 = SizeOfSRVI1
309
310 self.SRVIHsize=16
311
312 def SRVIread(self, fp):
313
314 header = numpy.fromfile(fp, SRVI_HEADER,1)
323 self.SizeOfSRVI1 = SizeOfSRVI1
324
325 self.SRVIHsize = 16
326
327 def SRVIread(self, fp):
328
329 header = numpy.fromfile(fp, SRVI_HEADER, 1)
315 330
316 331 self.SignatureSRVI1 = str(header['SignatureSRVI1'][0])
317 332 self.SizeOfDataBlock1 = header['SizeOfDataBlock1'][0]
318 333 self.DataBlockTitleSRVI1 = str(header['DataBlockTitleSRVI1'][0])
319 334 self.SizeOfSRVI1 = header['SizeOfSRVI1'][0]
320 #16
335 # 16
321 336 print 'Pointer fp SRVIheader', fp.tell()
322
323
324 SRVI_STRUCTURE = numpy.dtype([
325 ('frame_cnt','<u4'),#
326 ('time_t','<u4'), #
327 ('tpow','<f4'), #
328 ('npw1','<f4'), #
329 ('npw2','<f4'), #
330 ('cpw1','<f4'), #
331 ('pcw2','<f4'), #
332 ('ps_err','<u4'), #
333 ('te_err','<u4'), #
334 ('rc_err','<u4'), #
335 ('grs1','<u4'), #
336 ('grs2','<u4'), #
337 ('azipos','<f4'), #
338 ('azivel','<f4'), #
339 ('elvpos','<f4'), #
340 ('elvvel','<f4'), #
341 ('northAngle','<f4'), #
342 ('microsec','<u4'), #
343 ('azisetvel','<f4'), #
344 ('elvsetpos','<f4'), #
345 ('RadarConst','<f4'),]) #
346 337
347 338
339 SRVI_STRUCTURE = numpy.dtype([
340 ('frame_cnt', '<u4'),
341 ('time_t', '<u4'), #
342 ('tpow', '<f4'), #
343 ('npw1', '<f4'), #
344 ('npw2', '<f4'), #
345 ('cpw1', '<f4'), #
346 ('pcw2', '<f4'), #
347 ('ps_err', '<u4'), #
348 ('te_err', '<u4'), #
349 ('rc_err', '<u4'), #
350 ('grs1', '<u4'), #
351 ('grs2', '<u4'), #
352 ('azipos', '<f4'), #
353 ('azivel', '<f4'), #
354 ('elvpos', '<f4'), #
355 ('elvvel', '<f4'), #
356 ('northAngle', '<f4'),
357 ('microsec', '<u4'), #
358 ('azisetvel', '<f4'), #
359 ('elvsetpos', '<f4'), #
360 ('RadarConst', '<f4'), ]) #
348 361
349 362
350 363 class RecordHeader(Header):
351
352
353 def __init__(self, frame_cnt=0, time_t= 0, tpow=0, npw1=0, npw2=0,
354 cpw1=0, pcw2=0, ps_err=0, te_err=0, rc_err=0, grs1=0,
355 grs2=0, azipos=0, azivel=0, elvpos=0, elvvel=0, northangle=0,
356 microsec=0, azisetvel=0, elvsetpos=0, RadarConst=0 , RecCounter=0, Off2StartNxtRec=0):
357
358
364
365 def __init__(self, frame_cnt=0, time_t=0, tpow=0, npw1=0, npw2=0,
366 cpw1=0, pcw2=0, ps_err=0, te_err=0, rc_err=0, grs1=0,
367 grs2=0, azipos=0, azivel=0, elvpos=0, elvvel=0, northangle=0,
368 microsec=0, azisetvel=0, elvsetpos=0, RadarConst=0, RecCounter=0, Off2StartNxtRec=0):
369
359 370 self.frame_cnt = frame_cnt
360 371 self.dwell = time_t
361 372 self.tpow = tpow
362 373 self.npw1 = npw1
363 374 self.npw2 = npw2
364 375 self.cpw1 = cpw1
365 376 self.pcw2 = pcw2
366 377 self.ps_err = ps_err
367 378 self.te_err = te_err
368 379 self.rc_err = rc_err
369 380 self.grs1 = grs1
370 381 self.grs2 = grs2
371 382 self.azipos = azipos
372 383 self.azivel = azivel
373 384 self.elvpos = elvpos
374 385 self.elvvel = elvvel
375 386 self.northAngle = northangle
376 387 self.microsec = microsec
377 388 self.azisetvel = azisetvel
378 389 self.elvsetpos = elvsetpos
379 390 self.RadarConst = RadarConst
380 self.RHsize=84
391 self.RHsize = 84
381 392 self.RecCounter = RecCounter
382 self.Off2StartNxtRec=Off2StartNxtRec
383
393 self.Off2StartNxtRec = Off2StartNxtRec
394
384 395 def RHread(self, fp):
385
386 #startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
387
396
397 # startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
398
388 399 #OffRHeader= 1180 + self.RecCounter*(self.Off2StartNxtRec)
389 400 #startFp.seek(OffRHeader, os.SEEK_SET)
390
391 #print 'Posicion del bloque: ',OffRHeader
392
393 header = numpy.fromfile(fp,SRVI_STRUCTURE,1)
394
395 self.frame_cnt = header['frame_cnt'][0]#
401
402 # print 'Posicion del bloque: ',OffRHeader
403
404 header = numpy.fromfile(fp, SRVI_STRUCTURE, 1)
405
406 self.frame_cnt = header['frame_cnt'][0]
396 407 self.time_t = header['time_t'][0] #
397 408 self.tpow = header['tpow'][0] #
398 409 self.npw1 = header['npw1'][0] #
399 410 self.npw2 = header['npw2'][0] #
400 411 self.cpw1 = header['cpw1'][0] #
401 412 self.pcw2 = header['pcw2'][0] #
402 413 self.ps_err = header['ps_err'][0] #
403 414 self.te_err = header['te_err'][0] #
404 415 self.rc_err = header['rc_err'][0] #
405 416 self.grs1 = header['grs1'][0] #
406 417 self.grs2 = header['grs2'][0] #
407 418 self.azipos = header['azipos'][0] #
408 419 self.azivel = header['azivel'][0] #
409 420 self.elvpos = header['elvpos'][0] #
410 421 self.elvvel = header['elvvel'][0] #
411 422 self.northAngle = header['northAngle'][0] #
412 423 self.microsec = header['microsec'][0] #
413 424 self.azisetvel = header['azisetvel'][0] #
414 425 self.elvsetpos = header['elvsetpos'][0] #
415 426 self.RadarConst = header['RadarConst'][0] #
416 #84
417
418 #print 'Pointer fp RECheader', fp.tell()
419
427 # 84
428
429 # print 'Pointer fp RECheader', fp.tell()
430
420 431 #self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
421
432
422 433 #self.RHsize = 180+20*self.nChannels
423 434 #self.Datasize= self.nProfiles*self.nChannels*self.nHeights*2*4
424 #print 'Datasize',self.Datasize
435 # print 'Datasize',self.Datasize
425 436 #endFp = self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
426
437
427 438 print '=============================================='
428
439
429 440 print '=============================================='
430 441
431
432 442 return 1
433
434 class MIRA35CReader (ProcessingUnit,FileHeaderMIRA35c,SRVIHeader,RecordHeader):
435
443
444
445 class MIRA35CReader (ProcessingUnit, FileHeaderMIRA35c, SRVIHeader, RecordHeader):
446
436 447 path = None
437 448 startDate = None
438 449 endDate = None
439 450 startTime = None
440 451 endTime = None
441 452 walk = None
442 453 isConfig = False
443
444
445 fileList= None
446
447 #metadata
448 TimeZone= None
449 Interval= None
450 heightList= None
451
452 #data
453 data= None
454 utctime= None
455
456
457
458 def __init__(self, **kwargs):
459
460 #Eliminar de la base la herencia
454
455 fileList = None
456
457 # metadata
458 TimeZone = None
459 Interval = None
460 heightList = None
461
462 # data
463 data = None
464 utctime = None
465
466 def __init__(self, **kwargs):
467
468 # Eliminar de la base la herencia
461 469 ProcessingUnit.__init__(self, **kwargs)
462 470 self.PointerReader = 0
463 471 self.FileHeaderFlag = False
464 472 self.utc = None
465 473 self.ext = ".zspca"
466 474 self.optchar = "P"
467 self.fpFile=None
475 self.fpFile = None
468 476 self.fp = None
469 self.BlockCounter=0
477 self.BlockCounter = 0
470 478 self.dtype = None
471 479 self.fileSizeByHeader = None
472 480 self.filenameList = []
473 481 self.fileSelector = 0
474 self.Off2StartNxtRec=0
475 self.RecCounter=0
482 self.Off2StartNxtRec = 0
483 self.RecCounter = 0
476 484 self.flagNoMoreFiles = 0
477 self.data_spc=None
478 #self.data_cspc=None
479 self.data_output=None
485 self.data_spc = None
486 # self.data_cspc=None
487 self.data_output = None
480 488 self.path = None
481 self.OffsetStartHeader=0
482 self.Off2StartData=0
489 self.OffsetStartHeader = 0
490 self.Off2StartData = 0
483 491 self.ipp = 0
484 self.nFDTdataRecors=0
492 self.nFDTdataRecors = 0
485 493 self.blocksize = 0
486 494 self.dataOut = Spectra()
487 self.profileIndex = 1 #Always
488 self.dataOut.flagNoData=False
495 self.profileIndex = 1 # Always
496 self.dataOut.flagNoData = False
489 497 self.dataOut.nRdPairs = 0
490 498 self.dataOut.pairsList = []
491 self.dataOut.data_spc=None
492
493 self.dataOut.normFactor=1
499 self.dataOut.data_spc = None
500
501 self.dataOut.normFactor = 1
494 502 self.nextfileflag = True
495 503 self.dataOut.RadarConst = 0
496 504 self.dataOut.HSDV = []
497 505 self.dataOut.NPW = []
498 506 self.dataOut.COFA = []
499 507 self.dataOut.noise = 0
500
501
508
502 509 def Files2Read(self, fp):
503 510 '''
504 511 Function that indicates the number of .fdt files that exist in the folder to be read.
505 512 It also creates an organized list with the names of the files to read.
506 513 '''
507 #self.__checkPath()
508
509 ListaData=os.listdir(fp) #Gets the list of files within the fp address
510 ListaData=sorted(ListaData) #Sort the list of files from least to largest by names
511 nFiles=0 #File Counter
512 FileList=[] #A list is created that will contain the .fdt files
513 for IndexFile in ListaData :
514 if '.zspca' in IndexFile and '.gz' not in IndexFile:
514 # self.__checkPath()
515
516 # Gets the list of files within the fp address
517 ListaData = os.listdir(fp)
518 # Sort the list of files from least to largest by names
519 ListaData = sorted(ListaData)
520 nFiles = 0 # File Counter
521 FileList = [] # A list is created that will contain the .fdt files
522 for IndexFile in ListaData:
523 if '.zspca' in IndexFile and '.gz' not in IndexFile:
515 524 FileList.append(IndexFile)
516 nFiles+=1
517
518 #print 'Files2Read'
519 #print 'Existen '+str(nFiles)+' archivos .fdt'
520
521 self.filenameList=FileList #List of files from least to largest by names
522
523
525 nFiles += 1
526
527 # print 'Files2Read'
528 # print 'Existen '+str(nFiles)+' archivos .fdt'
529
530 self.filenameList = FileList # List of files from least to largest by names
531
524 532 def run(self, **kwargs):
525 533 '''
526 534 This method will be the one that will initiate the data entry, will be called constantly.
527 535 You should first verify that your Setup () is set up and then continue to acquire
528 536 the data to be processed with getData ().
529 537 '''
530 538 if not self.isConfig:
531 539 self.setup(**kwargs)
532 540 self.isConfig = True
533
541
534 542 self.getData()
535
536
543
537 544 def setup(self, path=None,
538 startDate=None,
539 endDate=None,
540 startTime=None,
541 endTime=None,
542 walk=True,
543 timezone='utc',
544 code = None,
545 online=False,
546 ReadMode=None, **kwargs):
547
545 startDate=None,
546 endDate=None,
547 startTime=None,
548 endTime=None,
549 walk=True,
550 timezone='utc',
551 code=None,
552 online=False,
553 ReadMode=None, **kwargs):
554
548 555 self.isConfig = True
549
550 self.path=path
551 self.startDate=startDate
552 self.endDate=endDate
553 self.startTime=startTime
554 self.endTime=endTime
555 self.walk=walk
556 #self.ReadMode=int(ReadMode)
557
556
557 self.path = path
558 self.startDate = startDate
559 self.endDate = endDate
560 self.startTime = startTime
561 self.endTime = endTime
562 self.walk = walk
563 # self.ReadMode=int(ReadMode)
564
558 565 pass
559
560
566
561 567 def getData(self):
562 568 '''
563 569 Before starting this function, you should check that there is still an unread file,
564 570 If there are still blocks to read or if the data block is empty.
565
571
566 572 You should call the file "read".
567
573
568 574 '''
569
575
570 576 if self.flagNoMoreFiles:
571 577 self.dataOut.flagNoData = True
572 578 print 'NoData se vuelve true'
573 579 return 0
574
575 self.fp=self.path
580
581 self.fp = self.path
576 582 self.Files2Read(self.fp)
577 583 self.readFile(self.fp)
578
579 self.dataOut.data_spc = self.dataOut_spc#self.data_spc.copy()
584
585 self.dataOut.data_spc = self.dataOut_spc # self.data_spc.copy()
580 586 self.dataOut.RadarConst = self.RadarConst
581 self.dataOut.data_output=self.data_output
587 self.dataOut.data_output = self.data_output
582 588 self.dataOut.noise = self.dataOut.getNoise()
583 #print 'ACAAAAAA', self.dataOut.noise
584 self.dataOut.data_spc = self.dataOut.data_spc+self.dataOut.noise
585 #print 'self.dataOut.noise',self.dataOut.noise
586
587
589 # print 'ACAAAAAA', self.dataOut.noise
590 self.dataOut.data_spc = self.dataOut.data_spc + self.dataOut.noise
591 # print 'self.dataOut.noise',self.dataOut.noise
592
588 593 return self.dataOut.data_spc
589
590
591 def readFile(self,fp):
594
595 def readFile(self, fp):
592 596 '''
593 597 You must indicate if you are reading in Online or Offline mode and load the
594 598 The parameters for this file reading mode.
595
599
596 600 Then you must do 2 actions:
597
601
598 602 1. Get the BLTR FileHeader.
599 603 2. Start reading the first block.
600 604 '''
601
602 #The address of the folder is generated the name of the .fdt file that will be read
603 print "File: ",self.fileSelector+1
604
605
606 # The address of the folder is generated the name of the .fdt file that will be read
607 print "File: ", self.fileSelector + 1
608
605 609 if self.fileSelector < len(self.filenameList):
606
607 self.fpFile=str(fp)+'/'+str(self.filenameList[self.fileSelector])
608
609 if self.nextfileflag==True:
610 self.fp = open(self.fpFile,"rb")
611 self.nextfileflag==False
612
610
611 self.fpFile = str(fp) + '/' + \
612 str(self.filenameList[self.fileSelector])
613
614 if self.nextfileflag == True:
615 self.fp = open(self.fpFile, "rb")
616 self.nextfileflag == False
617
613 618 '''HERE STARTING THE FILE READING'''
614
615
619
616 620 self.fheader = FileHeaderMIRA35c()
617 self.fheader.FHread(self.fp) #Bltr FileHeader Reading
618
621 self.fheader.FHread(self.fp) # Bltr FileHeader Reading
619 622
620 623 self.SPARrawGate1 = self.fheader.SPARrawGate1
621 624 self.SPARrawGate2 = self.fheader.SPARrawGate2
622 625 self.Num_Hei = self.SPARrawGate2 - self.SPARrawGate1
623 626 self.Num_Bins = self.fheader.PPARsft
624 627 self.dataOut.nFFTPoints = self.fheader.PPARsft
625
626
628
627 629 self.Num_inCoh = self.fheader.PPARavc
628 630 self.dataOut.PRF = self.fheader.PPARprf
629 self.dataOut.frequency = 34.85*10**9
630 self.Lambda = SPEED_OF_LIGHT/self.dataOut.frequency
631 self.dataOut.ippSeconds= 1./float(self.dataOut.PRF)
632
633 pulse_width = self.fheader.PPARpdr * 10**-9
631 self.dataOut.frequency = 34.85 * 10**9
632 self.Lambda = SPEED_OF_LIGHT / self.dataOut.frequency
633 self.dataOut.ippSeconds = 1. / float(self.dataOut.PRF)
634
635 pulse_width = self.fheader.PPARpdr * 10**-9
634 636 self.__deltaHeigth = 0.5 * SPEED_OF_LIGHT * pulse_width
635
636 self.data_spc = numpy.zeros((self.Num_Hei, self.Num_Bins,2))#
637
638 self.data_spc = numpy.zeros((self.Num_Hei, self.Num_Bins, 2))
637 639 self.dataOut.HSDV = numpy.zeros((self.Num_Hei, 2))
638
640
639 641 self.Ze = numpy.zeros(self.Num_Hei)
640 self.ETA = numpy.zeros(([2,self.Num_Hei]))
641
642
643
644 self.readBlock() #Block reading
645
642 self.ETA = numpy.zeros(([2, self.Num_Hei]))
643
644 self.readBlock() # Block reading
645
646 646 else:
647 647 print 'readFile FlagNoData becomes true'
648 self.flagNoMoreFiles=True
648 self.flagNoMoreFiles = True
649 649 self.dataOut.flagNoData = True
650 650 self.FileHeaderFlag == True
651 return 0
652
653
654
655 def readBlock(self):
651 return 0
652
653 def readBlock(self):
656 654 '''
657 655 It should be checked if the block has data, if it is not passed to the next file.
658
656
659 657 Then the following is done:
660
658
661 659 1. Read the RecordHeader
662 660 2. Fill the buffer with the current block number.
663
661
664 662 '''
665
663
666 664 if self.PointerReader > 1180:
667 self.fp.seek(self.PointerReader , os.SEEK_SET)
665 self.fp.seek(self.PointerReader, os.SEEK_SET)
668 666 self.FirstPoint = self.PointerReader
669
670 else :
667
668 else:
671 669 self.FirstPoint = 1180
672
673
674
670
675 671 self.srviHeader = SRVIHeader()
676
677 self.srviHeader.SRVIread(self.fp) #Se obtiene la cabecera del SRVI
678
679 self.blocksize = self.srviHeader.SizeOfDataBlock1 # Se obtiene el tamao del bloque
680
672
673 self.srviHeader.SRVIread(self.fp) # Se obtiene la cabecera del SRVI
674
675 self.blocksize = self.srviHeader.SizeOfDataBlock1 # Se obtiene el tamao del bloque
676
681 677 if self.blocksize == 148:
682 678 print 'blocksize == 148 bug'
683 jump = numpy.fromfile(self.fp,[('jump',numpy.str_,140)] ,1)
684
685 self.srviHeader.SRVIread(self.fp) #Se obtiene la cabecera del SRVI
686
679 jump = numpy.fromfile(self.fp, [('jump', numpy.str_, 140)], 1)
680
681 # Se obtiene la cabecera del SRVI
682 self.srviHeader.SRVIread(self.fp)
683
687 684 if not self.srviHeader.SizeOfSRVI1:
688 self.fileSelector+=1
689 self.nextfileflag==True
685 self.fileSelector += 1
686 self.nextfileflag == True
690 687 self.FileHeaderFlag == True
691
688
692 689 self.recordheader = RecordHeader()
693 690 self.recordheader.RHread(self.fp)
694 691 self.RadarConst = self.recordheader.RadarConst
695 692 dwell = self.recordheader.time_t
696 693 npw1 = self.recordheader.npw1
697 694 npw2 = self.recordheader.npw2
698
699
695
700 696 self.dataOut.channelList = range(1)
701 self.dataOut.nIncohInt = self.Num_inCoh
697 self.dataOut.nIncohInt = self.Num_inCoh
702 698 self.dataOut.nProfiles = self.Num_Bins
703 699 self.dataOut.nCohInt = 1
704 700 self.dataOut.windowOfFilter = 1
705 701 self.dataOut.utctime = dwell
706 self.dataOut.timeZone=0
707
702 self.dataOut.timeZone = 0
703
708 704 self.dataOut.outputInterval = self.dataOut.getTimeInterval()
709 self.dataOut.heightList = self.SPARrawGate1*self.__deltaHeigth + numpy.array(range(self.Num_Hei))*self.__deltaHeigth
710
711
712
713 self.HSDVsign = numpy.fromfile( self.fp, [('HSDV',numpy.str_,4)],1)
714 self.SizeHSDV = numpy.fromfile( self.fp, [('SizeHSDV','<i4')],1)
715 self.HSDV_Co = numpy.fromfile( self.fp, [('HSDV_Co','<f4')],self.Num_Hei)
716 self.HSDV_Cx = numpy.fromfile( self.fp, [('HSDV_Cx','<f4')],self.Num_Hei)
717
718 self.COFAsign = numpy.fromfile( self.fp, [('COFA',numpy.str_,4)],1)
719 self.SizeCOFA = numpy.fromfile( self.fp, [('SizeCOFA','<i4')],1)
720 self.COFA_Co = numpy.fromfile( self.fp, [('COFA_Co','<f4')],self.Num_Hei)
721 self.COFA_Cx = numpy.fromfile( self.fp, [('COFA_Cx','<f4')],self.Num_Hei)
722
723 self.ZSPCsign = numpy.fromfile(self.fp, [('ZSPCsign',numpy.str_,4)],1)
724 self.SizeZSPC = numpy.fromfile(self.fp, [('SizeZSPC','<i4')],1)
725
726 self.dataOut.HSDV[0]=self.HSDV_Co[:][0]
727 self.dataOut.HSDV[1]=self.HSDV_Cx[:][0]
728
705 self.dataOut.heightList = self.SPARrawGate1 * self.__deltaHeigth + \
706 numpy.array(range(self.Num_Hei)) * self.__deltaHeigth
707
708 self.HSDVsign = numpy.fromfile(self.fp, [('HSDV', numpy.str_, 4)], 1)
709 self.SizeHSDV = numpy.fromfile(self.fp, [('SizeHSDV', '<i4')], 1)
710 self.HSDV_Co = numpy.fromfile(
711 self.fp, [('HSDV_Co', '<f4')], self.Num_Hei)
712 self.HSDV_Cx = numpy.fromfile(
713 self.fp, [('HSDV_Cx', '<f4')], self.Num_Hei)
714
715 self.COFAsign = numpy.fromfile(self.fp, [('COFA', numpy.str_, 4)], 1)
716 self.SizeCOFA = numpy.fromfile(self.fp, [('SizeCOFA', '<i4')], 1)
717 self.COFA_Co = numpy.fromfile(
718 self.fp, [('COFA_Co', '<f4')], self.Num_Hei)
719 self.COFA_Cx = numpy.fromfile(
720 self.fp, [('COFA_Cx', '<f4')], self.Num_Hei)
721
722 self.ZSPCsign = numpy.fromfile(
723 self.fp, [('ZSPCsign', numpy.str_, 4)], 1)
724 self.SizeZSPC = numpy.fromfile(self.fp, [('SizeZSPC', '<i4')], 1)
725
726 self.dataOut.HSDV[0] = self.HSDV_Co[:][0]
727 self.dataOut.HSDV[1] = self.HSDV_Cx[:][0]
728
729 729 for irg in range(self.Num_Hei):
730 nspc = numpy.fromfile(self.fp, [('nspc','int16')],1)[0][0] # Number of spectral sub pieces containing significant power
731
730 # Number of spectral sub pieces containing significant power
731 nspc = numpy.fromfile(self.fp, [('nspc', 'int16')], 1)[0][0]
732
732 733 for k in range(nspc):
733 binIndex = numpy.fromfile(self.fp, [('binIndex','int16')],1)[0][0] # Index of the spectral bin where the piece is beginning
734 nbins = numpy.fromfile(self.fp, [('nbins','int16')],1)[0][0] # Number of bins of the piece
735
736 #Co_Channel
737 jbin = numpy.fromfile(self.fp, [('jbin','uint16')],nbins)[0][0] # Spectrum piece to be normaliced
738 jmax = numpy.fromfile(self.fp, [('jmax','float32')],1)[0][0] # Maximun piece to be normaliced
739
740
741 self.data_spc[irg,binIndex:binIndex+nbins,0] = self.data_spc[irg,binIndex:binIndex+nbins,0]+jbin/65530.*jmax
742
743 #Cx_Channel
744 jbin = numpy.fromfile(self.fp, [('jbin','uint16')],nbins)[0][0]
745 jmax = numpy.fromfile(self.fp, [('jmax','float32')],1)[0][0]
746
747
748 self.data_spc[irg,binIndex:binIndex+nbins,1] = self.data_spc[irg,binIndex:binIndex+nbins,1]+jbin/65530.*jmax
749
750 for bin in range(self.Num_Bins):
751
752 self.data_spc[:,bin,0] = self.data_spc[:,bin,0] - self.dataOut.HSDV[:,0]
753
754 self.data_spc[:,bin,1] = self.data_spc[:,bin,1] - self.dataOut.HSDV[:,1]
755
756
734 # Index of the spectral bin where the piece is beginning
735 binIndex = numpy.fromfile(
736 self.fp, [('binIndex', 'int16')], 1)[0][0]
737 nbins = numpy.fromfile(self.fp, [('nbins', 'int16')], 1)[
738 0][0] # Number of bins of the piece
739
740 # Co_Channel
741 jbin = numpy.fromfile(self.fp, [('jbin', 'uint16')], nbins)[
742 0][0] # Spectrum piece to be normaliced
743 jmax = numpy.fromfile(self.fp, [('jmax', 'float32')], 1)[
744 0][0] # Maximun piece to be normaliced
745
746 self.data_spc[irg, binIndex:binIndex + nbins, 0] = self.data_spc[irg,
747 binIndex:binIndex + nbins, 0] + jbin / 65530. * jmax
748
749 # Cx_Channel
750 jbin = numpy.fromfile(
751 self.fp, [('jbin', 'uint16')], nbins)[0][0]
752 jmax = numpy.fromfile(self.fp, [('jmax', 'float32')], 1)[0][0]
753
754 self.data_spc[irg, binIndex:binIndex + nbins, 1] = self.data_spc[irg,
755 binIndex:binIndex + nbins, 1] + jbin / 65530. * jmax
756
757 for bin in range(self.Num_Bins):
758
759 self.data_spc[:, bin, 0] = self.data_spc[:,
760 bin, 0] - self.dataOut.HSDV[:, 0]
761
762 self.data_spc[:, bin, 1] = self.data_spc[:,
763 bin, 1] - self.dataOut.HSDV[:, 1]
764
757 765 numpy.set_printoptions(threshold='nan')
758
759 self.data_spc = numpy.where(self.data_spc > 0. , self.data_spc, 0)
760
761 self.dataOut.COFA = numpy.array([self.COFA_Co , self.COFA_Cx])
762
766
767 self.data_spc = numpy.where(self.data_spc > 0., self.data_spc, 0)
768
769 self.dataOut.COFA = numpy.array([self.COFA_Co, self.COFA_Cx])
770
763 771 print ' '
764 print 'SPC',numpy.shape(self.dataOut.data_spc)
765 #print 'SPC',self.dataOut.data_spc
766
772 print 'SPC', numpy.shape(self.dataOut.data_spc)
773 # print 'SPC',self.dataOut.data_spc
774
767 775 noinor1 = 713031680
768 776 noinor2 = 30
769
770 npw1 = 1#0**(npw1/10) * noinor1 * noinor2
771 npw2 = 1#0**(npw2/10) * noinor1 * noinor2
777
778 npw1 = 1 # 0**(npw1/10) * noinor1 * noinor2
779 npw2 = 1 # 0**(npw2/10) * noinor1 * noinor2
772 780 self.dataOut.NPW = numpy.array([npw1, npw2])
773
781
774 782 print ' '
775
776 self.data_spc = numpy.transpose(self.data_spc, (2,1,0))
777 self.data_spc = numpy.fft.fftshift(self.data_spc, axes = 1)
778
783
784 self.data_spc = numpy.transpose(self.data_spc, (2, 1, 0))
785 self.data_spc = numpy.fft.fftshift(self.data_spc, axes=1)
786
779 787 self.data_spc = numpy.fliplr(self.data_spc)
780
781 self.data_spc = numpy.where(self.data_spc > 0. , self.data_spc, 0)
782 self.dataOut_spc= numpy.ones([1, self.Num_Bins , self.Num_Hei])
783 self.dataOut_spc[0,:,:] = self.data_spc[0,:,:]
784 #print 'SHAPE', self.dataOut_spc.shape
785 #For nyquist correction:
786 #fix = 20 # ~3m/s
788
789 self.data_spc = numpy.where(self.data_spc > 0., self.data_spc, 0)
790 self.dataOut_spc = numpy.ones([1, self.Num_Bins, self.Num_Hei])
791 self.dataOut_spc[0, :, :] = self.data_spc[0, :, :]
792 # print 'SHAPE', self.dataOut_spc.shape
793 # For nyquist correction:
794 # fix = 20 # ~3m/s
787 795 #shift = self.Num_Bins/2 + fix
788 796 #self.data_spc = numpy.array([ self.data_spc[: , self.Num_Bins-shift+1: , :] , self.data_spc[: , 0:self.Num_Bins-shift , :]])
789
790
791
797
792 798 '''Block Reading, the Block Data is received and Reshape is used to give it
793 799 shape.
794 800 '''
795
801
796 802 self.PointerReader = self.fp.tell()
797
798
799
800
801
802
803 No newline at end of file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now