##// 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 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 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 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
30
28 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 49 # matplotlib.pyplot.ioff()
45 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 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 96 (xpos, ypos),
87 97 colspan=colspan,
88 98 rowspan=rowspan,
89 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 108 xy = (.1, .99),
98 109 xycoords = 'figure fraction',
99 110 horizontalalignment = 'left',
100 111 verticalalignment = 'top',
101 112 fontsize = 10)
102 113
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 125 grid=None,color='blue'):
113
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 134 ax.set_xlim([xmin,xmax])
123 135 ax.set_ylim([ymin,ymax])
124 136
125 137 printLabels(ax, xlabel, ylabel, title)
126 138
127 139 ######################################################
128 140 if (xmax-xmin)<=1:
129 141 xtickspos = numpy.linspace(xmin,xmax,nxticks)
130 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 187 ax.lines[idline].set_data(x,y)
174 188
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 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 205 xlabel='', ylabel='', title='', ticksize = 9,
190 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 216 ax.set_xlim([xmin,xmax])
201 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 222 cmap=matplotlib.pyplot.get_cmap(colormap)
207 223 cmap.set_bad('black', 1.)
208 224 imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=cmap)
209 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 287 cmap=matplotlib.pyplot.get_cmap(colormap)
267 288 cmap.set_bad('black', 1.)
268 289
269
270 290 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=cmap)
271 291
292
272 293 def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
273 294 ticksize=9, xtick_visible=True, ytick_visible=True,
274 295 nxticks=4, nyticks=10,
275 296 grid=None):
276
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 308 ax.set_xlim([xmin,xmax])
289 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 356 line.set_data(x[i,:],y)
336 357
358
337 359 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
338 360 ticksize=9, xtick_visible=True, ytick_visible=True,
339 361 nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
340 362 grid=None, XAxisAsTime=False):
341
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 382 ax.set_xlim([xmin,xmax])
361 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 436 line.set_data(x,y[i,:])
413 437
438
414 439 def createPolar(ax, x, y,
415 440 xlabel='', ylabel='', title='', ticksize = 9,
416 441 colormap='jet',cblabel='', cbsize="5%",
417 442 XAxisAsTime=False):
418 443
419 444 matplotlib.pyplot.ioff()
420 445
421 446 ax.plot(x,y,'bo', markersize=5)
422 447 # ax.set_rmax(90)
423 448 ax.set_ylim(0,90)
424 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 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")
33 startFp = open(
34 '/home/erick/Documents/MIRA35C/20160117/20160117_0000.zspc', "rb")
35 35
36 36
37 37 FILE_HEADER = numpy.dtype([ #HEADER 1024bytes
38 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
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),
43 47 ('Hdummy',numpy.str_,512), #Reserved space
44 48 #Main chunk
45 49 ('Msign','<i4'), #Main chunk signature FZKF or NUIG
46 50 ('MsizeData','<i4'), #Size of data block main chunk
47 51 #Processing DSP parameters
48 52 ('PPARsign','<i4'), #PPAR signature
49 53 ('PPARsize','<i4'), #PPAR size of block
50 54 ('PPARprf','<i4'), #Pulse repetition frequency
51 55 ('PPARpdr','<i4'), #Pulse duration
52 56 ('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 # 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'),
57 65 #Service DSP parameters
58 ('SPARatt','<i4'), #STC attenuation on the lowest ranges on/off
66 # STC attenuation on the lowest ranges on/off
67 ('SPARatt', '<i4'),
59 68 ('SPARtx','<i4'), #OBSOLETE
60 69 ('SPARaddGain0','<f4'), #OBSOLETE
61 70 ('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
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'),
66 79 ('SPARcal','<i4'), #OBSOLETE
67 80 ('SPARnos','<i4'), #OBSOLETE
68 81 ('SPARof0','<i4'), #detection threshold
69 82 ('SPARof1','<i4'), #OBSOLETE
70 83 ('SPARswt','<i4'), #2nd moment estimation threshold
71 84 ('SPARsum','<i4'), #OBSOLETE
72 85 ('SPARosc','<i4'), #flag Oscillosgram mode
73 86 ('SPARtst','<i4'), #OBSOLETE
74 87 ('SPARcor','<i4'), #OBSOLETE
75 88 ('SPARofs','<i4'), #OBSOLETE
76 ('SPARhsn','<i4'), #Hildebrand div noise detection on noise gate
77 ('SPARhsa','<f4'), #Hildebrand div noise detection on all gates
89 # Hildebrand div noise detection on noise gate
90 ('SPARhsn', '<i4'),
91 # Hildebrand div noise detection on all gates
92 ('SPARhsa', '<f4'),
78 93 ('SPARcalibPow_M','<f4'), #OBSOLETE
79 94 ('SPARcalibSNR_M','<f4'), #OBSOLETE
80 95 ('SPARcalibPow_S','<f4'), #OBSOLETE
81 96 ('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
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'),
85 103 ('SPARprc','<i4'),]) #flag - Moment estimation switched on/off
86 104
87 105
88
89 106 self.Hname= None
90 107 self.Htime= None
91 108 self.Hoper= None
92 109 self.Hplace= None
93 110 self.Hdescr= None
94 111 self.Hdummy= None
95 112
96 113 self.Msign=None
97 114 self.MsizeData=None
98 115
99 116 self.PPARsign=None
100 117 self.PPARsize=None
101 118 self.PPARprf=None
102 119 self.PPARpdr=None
103 120 self.PPARsft=None
104 121 self.PPARavc=None
105 122 self.PPARihp=None
106 123 self.PPARchg=None
107 124 self.PPARpol=None
108 125 #Service DSP parameters
109 126 self.SPARatt=None
110 127 self.SPARtx=None
111 128 self.SPARaddGain0=None
112 129 self.SPARaddGain1=None
113 130 self.SPARwnd=None
114 131 self.SPARpos=None
115 132 self.SPARadd=None
116 133 self.SPARlen=None
117 134 self.SPARcal=None
118 135 self.SPARnos=None
119 136 self.SPARof0=None
120 137 self.SPARof1=None
121 138 self.SPARswt=None
122 139 self.SPARsum=None
123 140 self.SPARosc=None
124 141 self.SPARtst=None
125 142 self.SPARcor=None
126 143 self.SPARofs=None
127 144 self.SPARhsn=None
128 145 self.SPARhsa=None
129 146 self.SPARcalibPow_M=None
130 147 self.SPARcalibSNR_M=None
131 148 self.SPARcalibPow_S=None
132 149 self.SPARcalibSNR_S=None
133 150 self.SPARrawGate1=None
134 151 self.SPARrawGate2=None
135 152 self.SPARraw=None
136 153 self.SPARprc=None
137 154
138 155
139
140 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 176 Hname= str(header['Hname'][0])
161 177 Htime= str(header['Htime'][0])
162 178 Hoper= str(header['Hoper'][0])
163 179 Hplace= str(header['Hplace'][0])
164 180 Hdescr= str(header['Hdescr'][0])
165 181 Hdummy= str(header['Hdummy'][0])
166 182
167 183 Msign=header['Msign'][0]
168 184 MsizeData=header['MsizeData'][0]
169 185
170 186 PPARsign=header['PPARsign'][0]
171 187 PPARsize=header['PPARsize'][0]
172 188 PPARprf=header['PPARprf'][0]
173 189 PPARpdr=header['PPARpdr'][0]
174 190 PPARsft=header['PPARsft'][0]
175 191 PPARavc=header['PPARavc'][0]
176 192 PPARihp=header['PPARihp'][0]
177 193 PPARchg=header['PPARchg'][0]
178 194 PPARpol=header['PPARpol'][0]
179 195 #Service DSP parameters
180 196 SPARatt=header['SPARatt'][0]
181 197 SPARtx=header['SPARtx'][0]
182 198 SPARaddGain0=header['SPARaddGain0'][0]
183 199 SPARaddGain1=header['SPARaddGain1'][0]
184 200 SPARwnd=header['SPARwnd'][0]
185 201 SPARpos=header['SPARpos'][0]
186 202 SPARadd=header['SPARadd'][0]
187 203 SPARlen=header['SPARlen'][0]
188 204 SPARcal=header['SPARcal'][0]
189 205 SPARnos=header['SPARnos'][0]
190 206 SPARof0=header['SPARof0'][0]
191 207 SPARof1=header['SPARof1'][0]
192 208 SPARswt=header['SPARswt'][0]
193 209 SPARsum=header['SPARsum'][0]
194 210 SPARosc=header['SPARosc'][0]
195 211 SPARtst=header['SPARtst'][0]
196 212 SPARcor=header['SPARcor'][0]
197 213 SPARofs=header['SPARofs'][0]
198 214 SPARhsn=header['SPARhsn'][0]
199 215 SPARhsa=header['SPARhsa'][0]
200 216 SPARcalibPow_M=header['SPARcalibPow_M'][0]
201 217 SPARcalibSNR_M=header['SPARcalibSNR_M'][0]
202 218 SPARcalibPow_S=header['SPARcalibPow_S'][0]
203 219 SPARcalibSNR_S=header['SPARcalibSNR_S'][0]
204 220 SPARrawGate1=header['SPARrawGate1'][0]
205 221 SPARrawGate2=header['SPARrawGate2'][0]
206 222 SPARraw=header['SPARraw'][0]
207 223 SPARprc=header['SPARprc'][0]
208 224
209 225
210
211 226 SRVI_STRUCTURE = numpy.dtype([
212 ('frame_cnt','<u4'),#
227 ('frame_cnt', '<u4'),
213 228 ('time_t','<u4'), #
214 229 ('tpow','<f4'), #
215 230 ('npw1','<f4'), #
216 231 ('npw2','<f4'), #
217 232 ('cpw1','<f4'), #
218 233 ('pcw2','<f4'), #
219 234 ('ps_err','<u4'), #
220 235 ('te_err','<u4'), #
221 236 ('rc_err','<u4'), #
222 237 ('grs1','<u4'), #
223 238 ('grs2','<u4'), #
224 239 ('azipos','<f4'), #
225 240 ('azivel','<f4'), #
226 241 ('elvpos','<f4'), #
227 242 ('elvvel','<f4'), #
228 ('northAngle','<f4'), #
243 ('northAngle', '<f4'),
229 244 ('microsec','<u4'), #
230 245 ('azisetvel','<f4'), #
231 246 ('elvsetpos','<f4'), #
232 247 ('RadarConst','<f4'),]) #
233 248
234 249 JUMP_STRUCTURE = numpy.dtype([
235 ('jump','<u140'),#
236 ('SizeOfDataBlock1',numpy.str_,32),#
237 ('jump','<i4'),#
238 ('DataBlockTitleSRVI1',numpy.str_,32),#
239 ('SizeOfSRVI1','<i4'),])#
240
250 ('jump', '<u140'),
251 ('SizeOfDataBlock1', numpy.str_, 32),
252 ('jump', '<i4'),
253 ('DataBlockTitleSRVI1', numpy.str_, 32),
254 ('SizeOfSRVI1', '<i4'), ])
241 255
242 256
243 257 #frame_cnt=0, time_t= 0, tpow=0, npw1=0, npw2=0,
244 258 #cpw1=0, pcw2=0, ps_err=0, te_err=0, rc_err=0, grs1=0,
245 259 #grs2=0, azipos=0, azivel=0, elvpos=0, elvvel=0, northangle=0,
246 260 #microsec=0, azisetvel=0, elvsetpos=0, RadarConst=0
247 261
248 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 284
271 285
272
273 286 #print fp
274 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.
275 288 #startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
276 289 #RecCounter=0
277 290 #Off2StartNxtRec=811248
278 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 294 print 'debe ser 48, RecCounter*811248', self.OffsetStartHeader,self.RecCounter,self.Off2StartNxtRec
282 295 print 'Posicion del bloque: ',OffRHeader
283 296
284 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 319 self.RadarConst = header['frame_cnt'][0] #
307 320
308 321
309 322 self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
310 323
311 324 self.RHsize = 180+20*self.nChannels
312 325 self.Datasize= self.nProfiles*self.nChannels*self.nHeights*2*4
313 326 #print 'Datasize',self.Datasize
314 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
@@ -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 39
40
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 49
51 50 message = "#"*50 + "\n"
52 51 message += self.__class__.__name__.upper() + "\n"
53 52 message += "#"*50 + "\n"
54 53
55 54 keyList = self.__dict__.keys()
56 55 keyList.sort()
57 56
58 57 for key in keyList:
59 58 message += "%s = %s" %(key, self.__dict__[key]) + "\n"
60 59
61 60 if "size" not in keyList:
62 61 attr = getattr(self, "size")
63 62
64 63 if attr:
65 64 message += "%s = %s" %("size", attr) + "\n"
66 65
67 66 #print message
68 67
69 68
70
71
72
73 69 FILE_STRUCTURE = numpy.dtype([ #HEADER 48bytes
74 70 ('FileMgcNumber','<u4'), #0x23020100
75 ('nFDTdataRecors','<u4'), #No Of FDT data records in this file (0 or more)
71 # No Of FDT data records in this file (0 or more)
72 ('nFDTdataRecors', '<u4'),
76 73 ('OffsetStartHeader','<u4'),
77 74 ('RadarUnitId','<u4'),
78 75 ('SiteName',numpy.str_,32), #Null terminated
79 76 ])
80 77
78
81 79 class FileHeaderBLTR(Header):
82 80
83 81 def __init__(self):
84 82
85 83 self.FileMgcNumber= 0 #0x23020100
86 self.nFDTdataRecors=0 #No Of FDT data records in this file (0 or more)
84 # No Of FDT data records in this file (0 or more)
85 self.nFDTdataRecors = 0
87 86 self.RadarUnitId= 0
88 87 self.OffsetStartHeader=0
89 88 self.SiteName= ""
90 89 self.size = 48
91 90
92 91 def FHread(self, fp):
93 92 #try:
94 93 startFp = open(fp,"rb")
95 94
96 95 header = numpy.fromfile(startFp, FILE_STRUCTURE,1)
97 96
98 97 print ' '
99 98 print 'puntero file header', startFp.tell()
100 99 print ' '
101 100
102
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
122
123
124 120 self.FileMgcNumber= hex(header['FileMgcNumber'][0])
125 self.nFDTdataRecors=int(header['nFDTdataRecors'][0]) #No Of FDT data records in this file (0 or more)
121 # No Of FDT data records in this file (0 or more)
122 self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
126 123 self.RadarUnitId= int(header['RadarUnitId'][0])
127 124 self.OffsetStartHeader= int(header['OffsetStartHeader'][0])
128 125 self.SiteName= str(header['SiteName'][0])
129 126
130 127 #print 'Numero de bloques', self.nFDTdataRecors
131 128
132
133 129 if self.size <48:
134 130 return 0
135 131
136 132 return 1
137 133
138
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 141
147
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 159
166 160 return 1
167 161
168 162
169
170
171
172 163 RECORD_STRUCTURE = numpy.dtype([ #RECORD HEADER 180+20N bytes
173 164 ('RecMgcNumber','<u4'), #0x23030001
174 165 ('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)
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 184 ('TransmitFrec','<u4'), #Transmit frequency (Hz)
185 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)
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'),
196 206 ('PRFhz','<u4'), #PRF (Hz)
197 207 ('nCohInt','<u4'), #Integrations
198 ('nProfiles','<u4'), #Number of data points transformed
199 ('nChannels','<u4'), #Number of receive beams stored in file (1 or N)
208 # Number of data points transformed
209 ('nProfiles', '<u4'),
210 # Number of receive beams stored in file (1 or N)
211 ('nChannels', '<u4'),
200 212 ('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
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'),
219 249 ])
220 250
221 251
222 252 class RecordHeaderBLTR(Header):
223 253
224 254 def __init__(self, RecMgcNumber=None, RecCounter= 0, Off2StartNxtRec= 811248,
225 255 nUtime= 0, nMilisec= 0, ExpTagName= None,
226 256 ExpComment=None, SiteLatDegrees=0, SiteLongDegrees= 0,
227 257 RTCgpsStatus= 0, TransmitFrec= 0, ReceiveFrec= 0,
228 258 FirstOsciFrec= 0, Polarisation= 0, ReceiverFiltSett= 0,
229 259 nModesInUse= 0, DualModeIndex= 0, DualModeRange= 0,
230 260 nDigChannels= 0, SampResolution= 0, nHeights= 0,
231 261 StartRangeSamp= 0, PRFhz= 0, nCohInt= 0,
232 262 nProfiles= 0, nChannels= 0, nIncohInt= 0,
233 263 FFTwindowingInd= 0, BeamAngleAzim= 0, BeamAngleZen= 0,
234 264 AntennaCoord0= 0, AntennaCoord1= 0, AntennaCoord2= 0,
235 265 RecPhaseCalibr0= 0, RecPhaseCalibr1= 0, RecPhaseCalibr2= 0,
236 266 RecAmpCalibr0= 0, RecAmpCalibr1= 0, RecAmpCalibr2= 0,
237 267 AntennaAngl0=0, AntennaAngl1=0, AntennaAngl2=0,
238 268 ReceiverGaindB0= 0, ReceiverGaindB1= 0, ReceiverGaindB2= 0, Off2StartData=0, OffsetStartHeader=0):
239 269
240 270 self.RecMgcNumber = RecMgcNumber #0x23030001
241 271 self.RecCounter = RecCounter
242 272 self.Off2StartNxtRec = Off2StartNxtRec
243 273 self.Off2StartData = Off2StartData
244 274 self.nUtime = nUtime
245 275 self.nMilisec = nMilisec
246 276 self.ExpTagName = ExpTagName
247 277 self.ExpComment = ExpComment
248 278 self.SiteLatDegrees = SiteLatDegrees
249 279 self.SiteLongDegrees = SiteLongDegrees
250 280 self.RTCgpsStatus = RTCgpsStatus
251 281 self.TransmitFrec = TransmitFrec
252 282 self.ReceiveFrec = ReceiveFrec
253 283 self.FirstOsciFrec = FirstOsciFrec
254 284 self.Polarisation = Polarisation
255 285 self.ReceiverFiltSett = ReceiverFiltSett
256 286 self.nModesInUse = nModesInUse
257 287 self.DualModeIndex = DualModeIndex
258 288 self.DualModeRange = DualModeRange
259 289 self.nDigChannels = nDigChannels
260 290 self.SampResolution = SampResolution
261 291 self.nHeights = nHeights
262 292 self.StartRangeSamp = StartRangeSamp
263 293 self.PRFhz = PRFhz
264 294 self.nCohInt = nCohInt
265 295 self.nProfiles = nProfiles
266 296 self.nChannels = nChannels
267 297 self.nIncohInt = nIncohInt
268 298 self.FFTwindowingInd = FFTwindowingInd
269 299 self.BeamAngleAzim = BeamAngleAzim
270 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 306 self.AntennaCoord2 = AntennaCoord2
277 307 self.RecPhaseCalibr0 = RecPhaseCalibr0
278 308 self.RecPhaseCalibr1 = RecPhaseCalibr1
279 309 self.RecPhaseCalibr2 = RecPhaseCalibr2
280 310 self.RecAmpCalibr0 = RecAmpCalibr0
281 311 self.RecAmpCalibr1 = RecAmpCalibr1
282 312 self.RecAmpCalibr2 = RecAmpCalibr2
283 313 self.ReceiverGaindB0 = ReceiverGaindB0
284 314 self.ReceiverGaindB1 = ReceiverGaindB1
285 315 self.ReceiverGaindB2 = ReceiverGaindB2
286 316 self.OffsetStartHeader = 48
287 317
288
289
290 318 def RHread(self, fp):
291 319 #print fp
292 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.
293 startFp = open(fp,"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")
294 323 #RecCounter=0
295 324 #Off2StartNxtRec=811248
296 325 OffRHeader= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
297 326 print ' '
298 327 print 'puntero Record Header', startFp.tell()
299 328 print ' '
300 329
301
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 335
308 336 #print 'Posicion del bloque: ',OffRHeader
309 337
310 338 header = numpy.fromfile(startFp,RECORD_STRUCTURE,1)
311 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 346 #print 'puntero Record Header despues de seek', header.tell()
319 347 print ' '
320 348
321 349 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) #0x23030001
322 350 self.RecCounter = int(header['RecCounter'][0])
323 351 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
324 352 self.Off2StartData = int(header['Off2StartData'][0])
325 353 self.nUtime = header['nUtime'][0]
326 354 self.nMilisec = header['nMilisec'][0]
327 355 self.ExpTagName = str(header['ExpTagName'][0])
328 356 self.ExpComment = str(header['ExpComment'][0])
329 357 self.SiteLatDegrees = header['SiteLatDegrees'][0]
330 358 self.SiteLongDegrees = header['SiteLongDegrees'][0]
331 359 self.RTCgpsStatus = header['RTCgpsStatus'][0]
332 360 self.TransmitFrec = header['TransmitFrec'][0]
333 361 self.ReceiveFrec = header['ReceiveFrec'][0]
334 362 self.FirstOsciFrec = header['FirstOsciFrec'][0]
335 363 self.Polarisation = header['Polarisation'][0]
336 364 self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
337 365 self.nModesInUse = header['nModesInUse'][0]
338 366 self.DualModeIndex = header['DualModeIndex'][0]
339 367 self.DualModeRange = header['DualModeRange'][0]
340 368 self.nDigChannels = header['nDigChannels'][0]
341 369 self.SampResolution = header['SampResolution'][0]
342 370 self.nHeights = header['nHeights'][0]
343 371 self.StartRangeSamp = header['StartRangeSamp'][0]
344 372 self.PRFhz = header['PRFhz'][0]
345 373 self.nCohInt = header['nCohInt'][0]
346 374 self.nProfiles = header['nProfiles'][0]
347 375 self.nChannels = header['nChannels'][0]
348 376 self.nIncohInt = header['nIncohInt'][0]
349 377 self.FFTwindowingInd = header['FFTwindowingInd'][0]
350 378 self.BeamAngleAzim = header['BeamAngleAzim'][0]
351 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 383 self.AntennaAngl1 = header['AntennaAngl1'][0]
356 384 self.AntennaCoord2 = header['AntennaCoord2'][0]
357 385 self.AntennaAngl2 = header['AntennaAngl2'][0]
358 386 self.RecPhaseCalibr0 = header['RecPhaseCalibr0'][0]
359 387 self.RecPhaseCalibr1 = header['RecPhaseCalibr1'][0]
360 388 self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0]
361 389 self.RecAmpCalibr0 = header['RecAmpCalibr0'][0]
362 390 self.RecAmpCalibr1 = header['RecAmpCalibr1'][0]
363 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 396 self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
369 397
370 398 self.RHsize = 180+20*self.nChannels
371 399 self.Datasize= self.nProfiles*self.nChannels*self.nHeights*2*4
372 400 #print 'Datasize',self.Datasize
373 401 endFp = self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
374 402
375 403 print '=============================================='
376 404 print 'RecMgcNumber ',self.RecMgcNumber
377 405 print 'RecCounter ',self.RecCounter
378 406 print 'Off2StartNxtRec ',self.Off2StartNxtRec
379 407 print 'Off2StartData ',self.Off2StartData
380 408 print 'Range Resolution ',self.SampResolution
381 409 print 'First Height ',self.StartRangeSamp
382 410 print 'PRF (Hz) ',self.PRFhz
383 411 print 'Heights (K) ',self.nHeights
384 412 print 'Channels (N) ',self.nChannels
385 413 print 'Profiles (J) ',self.nProfiles
386 414 print 'iCoh ',self.nCohInt
387 415 print 'iInCoh ',self.nIncohInt
388 416 print 'BeamAngleAzim ',self.BeamAngleAzim
389 417 print 'BeamAngleZen ',self.BeamAngleZen
390 418
391 419 #print 'ModoEnUso ',self.DualModeIndex
392 420 #print 'UtcTime ',self.nUtime
393 421 #print 'MiliSec ',self.nMilisec
394 422 #print 'Exp TagName ',self.ExpTagName
395 423 #print 'Exp Comment ',self.ExpComment
396 424 #print 'FFT Window Index ',self.FFTwindowingInd
397 425 #print 'N Dig. Channels ',self.nDigChannels
398 426 print 'Size de bloque ',self.RHsize
399 427 print 'DataSize ',self.Datasize
400 428 print 'BeamAngleAzim ',self.BeamAngleAzim
401 429 #print 'AntennaCoord0 ',self.AntennaCoord0
402 430 #print 'AntennaAngl0 ',self.AntennaAngl0
403 431 #print 'AntennaCoord1 ',self.AntennaCoord1
404 432 #print 'AntennaAngl1 ',self.AntennaAngl1
405 433 #print 'AntennaCoord2 ',self.AntennaCoord2
406 434 #print 'AntennaAngl2 ',self.AntennaAngl2
407 435 print 'RecPhaseCalibr0 ',self.RecPhaseCalibr0
408 436 print 'RecPhaseCalibr1 ',self.RecPhaseCalibr1
409 437 print 'RecPhaseCalibr2 ',self.RecPhaseCalibr2
410 438 print 'RecAmpCalibr0 ',self.RecAmpCalibr0
411 439 print 'RecAmpCalibr1 ',self.RecAmpCalibr1
412 440 print 'RecAmpCalibr2 ',self.RecAmpCalibr2
413 441 print 'ReceiverGaindB0 ',self.ReceiverGaindB0
414 442 print 'ReceiverGaindB1 ',self.ReceiverGaindB1
415 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 468
439
440 469 fileList= None
441 470
442 471 #metadata
443 472 TimeZone= None
444 473 Interval= None
445 474 heightList= None
446 475
447 476 #data
448 477 data= None
449 478 utctime= None
450 479
451
452
453 480 def __init__(self, **kwargs):
454 481
455 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 494 self.fpFile=None
468 495 self.fp = None
469 496 self.BlockCounter=0
470 497 self.dtype = None
471 498 self.fileSizeByHeader = None
472 499 self.filenameList = []
473 500 self.fileSelector = 0
474 501 self.Off2StartNxtRec=0
475 502 self.RecCounter=0
476 503 self.flagNoMoreFiles = 0
477 504 self.data_spc=None
478 505 self.data_cspc=None
479 506 self.data_output=None
480 507 self.path = None
481 508 self.OffsetStartHeader=0
482 509 self.Off2StartData=0
483 510 self.ipp = 0
484 511 self.nFDTdataRecors=0
485 512 self.blocksize = 0
486 513 self.dataOut = Spectra()
487 514 self.profileIndex = 1 #Always
488 515 self.dataOut.flagNoData=False
489 516 self.dataOut.nRdPairs = 0
490 517 self.dataOut.pairsList = []
491 518 self.dataOut.data_spc=None
492 519 self.dataOut.noise=[]
493 520 self.dataOut.velocityX=[]
494 521 self.dataOut.velocityY=[]
495 522 self.dataOut.velocityV=[]
496 523
497
498
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 529 #self.__checkPath()
505 530
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
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)
508 535 nFiles=0 #File Counter
509 536 FileList=[] #A list is created that will contain the .fdt files
510 537 for IndexFile in ListaData :
511 538 if '.fdt' in IndexFile:
512 539 FileList.append(IndexFile)
513 540 nFiles+=1
514 541
515 542 #print 'Files2Read'
516 543 #print 'Existen '+str(nFiles)+' archivos .fdt'
517 544
518 545 self.filenameList=FileList #List of files from least to largest by names
519 546
520
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 558 #print 'running'
533 559
534
535 560 def setup(self, path=None,
536 561 startDate=None,
537 562 endDate=None,
538 563 startTime=None,
539 564 endTime=None,
540 565 walk=True,
541 566 timezone='utc',
542 567 code = None,
543 568 online=False,
544 569 ReadMode=None,
545 570 **kwargs):
546 571
547 572 self.isConfig = True
548 573
549 574 self.path=path
550 575 self.startDate=startDate
551 576 self.endDate=endDate
552 577 self.startTime=startTime
553 578 self.endTime=endTime
554 579 self.walk=walk
555 580 self.ReadMode=int(ReadMode)
556 581
557 582 pass
558 583
559
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 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 602 self.dataOut.data_cspc =self.data_cspc
579 603 self.dataOut.data_output=self.data_output
580 604
581 605 print 'self.dataOut.data_output', shape(self.dataOut.data_output)
582 606
583 607 #self.removeDC()
584 608 return self.dataOut.data_spc
585 609
586
587 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 620
598 621 #The address of the folder is generated the name of the .fdt file that will be read
599 622 print "File: ",self.fileSelector+1
600 623
601 624 if self.fileSelector < len(self.filenameList):
602 625
603 self.fpFile=str(fp)+'/'+str(self.filenameList[self.fileSelector])
626 self.fpFile = str(fp) + '/' + \
627 str(self.filenameList[self.fileSelector])
604 628 #print self.fpFile
605 629 fheader = FileHeaderBLTR()
606 630 fheader.FHread(self.fpFile) #Bltr FileHeader Reading
607 631 self.nFDTdataRecors=fheader.nFDTdataRecors
608 632
609 633 self.readBlock() #Block reading
610 634 else:
611 635 print 'readFile FlagNoData becomes true'
612 636 self.flagNoMoreFiles=True
613 637 self.dataOut.flagNoData = True
614 638 return 0
615 639
616 640 def getVelRange(self, extrapoints=0):
617 641 Lambda= SPEED_OF_LIGHT/50000000
618 PRF = self.dataOut.PRF#1./(self.dataOut.ippSeconds * self.dataOut.nCohInt)
642 # 1./(self.dataOut.ippSeconds * self.dataOut.nCohInt)
643 PRF = self.dataOut.PRF
619 644 Vmax=-Lambda/(4.*(1./PRF)*self.dataOut.nCohInt*2.)
620 645 deltafreq = PRF / (self.nProfiles)
621 646 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.)
647 freqrange = deltafreq * \
648 (numpy.arange(self.nProfiles) - self.nProfiles / 2.) - deltafreq / 2
649 velrange = deltavel * \
650 (numpy.arange(self.nProfiles) - self.nProfiles / 2.)
624 651 return velrange
625 652
626 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 663
637 664 if self.BlockCounter < self.nFDTdataRecors-2:
638 665 print self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
639 666 if self.ReadMode==1:
640 667 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1)
641 668 elif self.ReadMode==0:
642 669 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter)
643 670
644 671 rheader.RHread(self.fpFile) #Bltr FileHeader Reading
645 672
646 673 self.OffsetStartHeader=rheader.OffsetStartHeader
647 674 self.RecCounter=rheader.RecCounter
648 675 self.Off2StartNxtRec=rheader.Off2StartNxtRec
649 676 self.Off2StartData=rheader.Off2StartData
650 677 self.nProfiles=rheader.nProfiles
651 678 self.nChannels=rheader.nChannels
652 679 self.nHeights=rheader.nHeights
653 680 self.frequency=rheader.TransmitFrec
654 681 self.DualModeIndex=rheader.DualModeIndex
655 682
656 683 self.pairsList =[(0,1),(0,2),(1,2)]
657 684 self.dataOut.pairsList = self.pairsList
658 685
659 686 self.nRdPairs=len(self.dataOut.pairsList)
660 687 self.dataOut.nRdPairs = self.nRdPairs
661 688
662 689 self.__firstHeigth=rheader.StartRangeSamp
663 690 self.__deltaHeigth=rheader.SampResolution
664 self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth
691 self.dataOut.heightList = self.__firstHeigth + \
692 numpy.array(range(self.nHeights)) * self.__deltaHeigth
665 693 self.dataOut.channelList = range(self.nChannels)
666 694 self.dataOut.nProfiles=rheader.nProfiles
667 695 self.dataOut.nIncohInt=rheader.nIncohInt
668 696 self.dataOut.nCohInt=rheader.nCohInt
669 697 self.dataOut.ippSeconds= 1/float(rheader.PRFhz)
670 698 self.dataOut.PRF=rheader.PRFhz
671 699 self.dataOut.nFFTPoints=rheader.nProfiles
672 700 self.dataOut.utctime=rheader.nUtime
673 701 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
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
676 706
677 707 self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN
678 708 print 'self.data_output', shape(self.data_output)
679 709 self.dataOut.velocityX=[]
680 710 self.dataOut.velocityY=[]
681 711 self.dataOut.velocityV=[]
682 712
683 713 '''Block Reading, the Block Data is received and Reshape is used to give it
684 714 shape.
685 715 '''
686 716
687 717 #Procedure to take the pointer to where the date block starts
688 718 startDATA = open(self.fpFile,"rb")
689 OffDATA= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec+self.Off2StartData
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 724 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
694 725
695 726 def gaus(xSamples,a,x0,sigma):
696 727 return a*exp(-(xSamples-x0)**2/(2*sigma**2))
697 728
698 729 def Find(x,value):
699 730 for index in range(len(x)):
700 731 if x[index]==value:
701 732 return index
702 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 738
708
709
710
711 739 if self.DualModeIndex==self.ReadMode:
712 740
713 self.data_fft = numpy.fromfile( startDATA, [('complex','<c8')],self.nProfiles*self.nChannels*self.nHeights )
741 self.data_fft = numpy.fromfile(
742 startDATA, [('complex', '<c8')], self.nProfiles * self.nChannels * self.nHeights)
714 743
715 744 self.data_fft=self.data_fft.astype(numpy.dtype('complex'))
716 745
717 self.data_block=numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles ))
746 self.data_block = numpy.reshape(
747 self.data_fft, (self.nHeights, self.nChannels, self.nProfiles))
718 748
719 749 self.data_block = numpy.transpose(self.data_block, (1,2,0))
720 750
721 751 copy = self.data_block.copy()
722 752 spc = copy * numpy.conjugate(copy)
723 753
724 self.data_spc = numpy.absolute(spc) # valor absoluto o magnitud
754 self.data_spc = numpy.absolute(
755 spc) # valor absoluto o magnitud
725 756
726 757 factor = self.dataOut.normFactor
727 758
728
729 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 764 print shape(z)
735 765 print ' '
736 766 print ' '
737 767
738 768 self.dataOut.data_spc=self.data_spc
739 769
740 self.noise = self.dataOut.getNoise(ymin_index=80, ymax_index=132)#/factor
770 self.noise = self.dataOut.getNoise(
771 ymin_index=80, ymax_index=132) # /factor
741 772 #noisedB = 10*numpy.log10(self.noise)
742 773
743
744 774 ySamples=numpy.ones([3,self.nProfiles])
745 775 phase=numpy.ones([3,self.nProfiles])
746 CSPCSamples=numpy.ones([3,self.nProfiles],dtype=numpy.complex_)
776 CSPCSamples = numpy.ones(
777 [3, self.nProfiles], dtype=numpy.complex_)
747 778 coherence=numpy.ones([3,self.nProfiles])
748 779 PhaseSlope=numpy.ones(3)
749 780 PhaseInter=numpy.ones(3)
750 781
751 782 '''****** Getting CrossSpectra ******'''
752 783 cspc=self.data_block.copy()
753 784 self.data_cspc=self.data_block.copy()
754 785
755 786 xFrec=self.getVelRange(1)
756 787 VelRange=self.getVelRange(1)
757 788 self.dataOut.VelRange=VelRange
758 789 #print ' '
759 790 #print ' '
760 791 #print 'xFrec',xFrec
761 792 #print ' '
762 793 #print ' '
763 794 #Height=35
764 795 for i in range(self.nRdPairs):
765 796
766 797 chan_index0 = self.dataOut.pairsList[i][0]
767 798 chan_index1 = self.dataOut.pairsList[i][1]
768 799
769 self.data_cspc[i,:,:]=cspc[chan_index0,:,:] * numpy.conjugate(cspc[chan_index1,:,:])
770
800 self.data_cspc[i, :, :] = cspc[chan_index0, :,
801 :] * numpy.conjugate(cspc[chan_index1, :, :])
771 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)
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)
776 810
777 811 E01=AntennaX0-AntennaX1
778 812 N01=AntennaY0-AntennaY1
779 813
780 814 E02=AntennaX0-AntennaX2
781 815 N02=AntennaY0-AntennaY2
782 816
783 817 E12=AntennaX1-AntennaX2
784 818 N12=AntennaY1-AntennaY2
785 819
786 self.ChanDist= numpy.array([[E01, N01],[E02,N02],[E12,N12]])
820 self.ChanDist = numpy.array(
821 [[E01, N01], [E02, N02], [E12, N12]])
787 822
788 823 self.dataOut.ChanDist = self.ChanDist
789 824
790 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 836 #
802 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 851 #
817 852 #
818 853 #
819 854 # '''****** CSPC is normalized ******'''
820 855 # chan_index0 = self.dataOut.pairsList[i][0]
821 856 # chan_index1 = self.dataOut.pairsList[i][1]
822 857 # CSPCFactor= numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])
823 858 #
824 859 #
825 860 # CSPCNorm= cspcLine.copy() / numpy.sqrt(CSPCFactor)
826 861 #
827 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 872 #
838 873 #
839 874 # '''****** Getting fij width ******'''
840 875 #
841 876 # yMean=[]
842 877 # yMean2=[]
843 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 888 #
854 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 894 #
860 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 918 #
884 919 #
885 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 956 # PhaseRange=moving_average(phase[i,Range[0]:Range[1]],N=3)
922 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 974 #
940 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 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 1010 #
976 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 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 1040 #
1006 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 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 1129 # #
1095 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 1173 # # # plt.close()
1139 1174 #
1140 1175 # print ' '
1141 1176
1142
1143
1144 1177 self.BlockCounter+=2
1145 1178
1146 1179 else:
1147 1180 self.fileSelector+=1
1148 1181 self.BlockCounter=0
1149 1182 print "Next File"
1150
1151
1152
1153
1154
Show file before | Show file after | Hide comments
@@ -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 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 36
39
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 46
50 47 message = "#"*50 + "\n"
51 48 message += self.__class__.__name__.upper() + "\n"
52 49 message += "#"*50 + "\n"
53 50
54 51 keyList = self.__dict__.keys()
55 52 keyList.sort()
56 53
57 54 for key in keyList:
58 55 message += "%s = %s" %(key, self.__dict__[key]) + "\n"
59 56
60 57 if "size" not in keyList:
61 58 attr = getattr(self, "size")
62 59
63 60 if attr:
64 61 message += "%s = %s" %("size", attr) + "\n"
65 62
66 63 #print message
67 64
68 65
69 66 FILE_HEADER = numpy.dtype([ #HEADER 1024bytes
70 67 ('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
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),
75 76 ('Hdummy',numpy.str_,512), #Reserved space
76 77 #Main chunk 8bytes
77 ('Msign',numpy.str_,4), #Main chunk signature FZKF or NUIG
78 # Main chunk signature FZKF or NUIG
79 ('Msign', numpy.str_, 4),
78 80 ('MsizeData','<i4'), #Size of data block main chunk
79 81 #Processing DSP parameters 36bytes
80 82 ('PPARsign',numpy.str_,4), #PPAR signature
81 83 ('PPARsize','<i4'), #PPAR size of block
82 84 ('PPARprf','<i4'), #Pulse repetition frequency
83 85 ('PPARpdr','<i4'), #Pulse duration
84 86 ('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.
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'),
89 95 #Service DSP parameters 112bytes
90 ('SPARatt','<i4'), #STC attenuation on the lowest ranges on/off
96 # STC attenuation on the lowest ranges on/off
97 ('SPARatt', '<i4'),
91 98 ('SPARtx','<i4'), #OBSOLETE
92 99 ('SPARaddGain0','<f4'), #OBSOLETE
93 100 ('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
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'),
98 109 ('SPARcal','<i4'), #OBSOLETE
99 110 ('SPARnos','<i4'), #OBSOLETE
100 111 ('SPARof0','<i4'), #detection threshold
101 112 ('SPARof1','<i4'), #OBSOLETE
102 113 ('SPARswt','<i4'), #2nd moment estimation threshold
103 114 ('SPARsum','<i4'), #OBSOLETE
104 115 ('SPARosc','<i4'), #flag Oscillosgram mode
105 116 ('SPARtst','<i4'), #OBSOLETE
106 117 ('SPARcor','<i4'), #OBSOLETE
107 118 ('SPARofs','<i4'), #OBSOLETE
108 ('SPARhsn','<i4'), #Hildebrand div noise detection on noise gate
109 ('SPARhsa','<f4'), #Hildebrand div noise detection on all gates
119 # Hildebrand div noise detection on noise gate
120 ('SPARhsn', '<i4'),
121 # Hildebrand div noise detection on all gates
122 ('SPARhsa', '<f4'),
110 123 ('SPARcalibPow_M','<f4'), #OBSOLETE
111 124 ('SPARcalibSNR_M','<f4'), #OBSOLETE
112 125 ('SPARcalibPow_S','<f4'), #OBSOLETE
113 126 ('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
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'),
117 133 ('SPARprc','<i4'),]) #flag - Moment estimation switched on/off
118 134
119 135
120
121 136 class FileHeaderMIRA35c(Header):
122 137
123 138 def __init__(self):
124 139
125 140 self.Hname= None
126 141 self.Htime= None
127 142 self.Hoper= None
128 143 self.Hplace= None
129 144 self.Hdescr= None
130 145 self.Hdummy= None
131 146
132 147 self.Msign=None
133 148 self.MsizeData=None
134 149
135 150 self.PPARsign=None
136 151 self.PPARsize=None
137 152 self.PPARprf=None
138 153 self.PPARpdr=None
139 154 self.PPARsft=None
140 155 self.PPARavc=None
141 156 self.PPARihp=None
142 157 self.PPARchg=None
143 158 self.PPARpol=None
144 159 #Service DSP parameters
145 160 self.SPARatt=None
146 161 self.SPARtx=None
147 162 self.SPARaddGain0=None
148 163 self.SPARaddGain1=None
149 164 self.SPARwnd=None
150 165 self.SPARpos=None
151 166 self.SPARadd=None
152 167 self.SPARlen=None
153 168 self.SPARcal=None
154 169 self.SPARnos=None
155 170 self.SPARof0=None
156 171 self.SPARof1=None
157 172 self.SPARswt=None
158 173 self.SPARsum=None
159 174 self.SPARosc=None
160 175 self.SPARtst=None
161 176 self.SPARcor=None
162 177 self.SPARofs=None
163 178 self.SPARhsn=None
164 179 self.SPARhsa=None
165 180 self.SPARcalibPow_M=None
166 181 self.SPARcalibSNR_M=None
167 182 self.SPARcalibPow_S=None
168 183 self.SPARcalibSNR_S=None
169 184 self.SPARrawGate1=None
170 185 self.SPARrawGate2=None
171 186 self.SPARraw=None
172 187 self.SPARprc=None
173 188
174 189 self.FHsize=1180
175 190
176 191 def FHread(self, fp):
177 192
178 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 212
198
199 213 self.Hname= str(header['Hname'][0])
200 214 self.Htime= str(header['Htime'][0])
201 215 self.Hoper= str(header['Hoper'][0])
202 216 self.Hplace= str(header['Hplace'][0])
203 217 self.Hdescr= str(header['Hdescr'][0])
204 218 self.Hdummy= str(header['Hdummy'][0])
205 219 #1024
206 220
207 221 self.Msign=str(header['Msign'][0])
208 222 self.MsizeData=header['MsizeData'][0]
209 223 #8
210 224
211 225 self.PPARsign=str(header['PPARsign'][0])
212 226 self.PPARsize=header['PPARsize'][0]
213 227 self.PPARprf=header['PPARprf'][0]
214 228 self.PPARpdr=header['PPARpdr'][0]
215 229 self.PPARsft=header['PPARsft'][0]
216 230 self.PPARavc=header['PPARavc'][0]
217 231 self.PPARihp=header['PPARihp'][0]
218 232 self.PPARchg=header['PPARchg'][0]
219 233 self.PPARpol=header['PPARpol'][0]
220 234 #Service DSP parameters
221 235 #36
222 236
223 237 self.SPARatt=header['SPARatt'][0]
224 238 self.SPARtx=header['SPARtx'][0]
225 239 self.SPARaddGain0=header['SPARaddGain0'][0]
226 240 self.SPARaddGain1=header['SPARaddGain1'][0]
227 241 self.SPARwnd=header['SPARwnd'][0]
228 242 self.SPARpos=header['SPARpos'][0]
229 243 self.SPARadd=header['SPARadd'][0]
230 244 self.SPARlen=header['SPARlen'][0]
231 245 self.SPARcal=header['SPARcal'][0]
232 246 self.SPARnos=header['SPARnos'][0]
233 247 self.SPARof0=header['SPARof0'][0]
234 248 self.SPARof1=header['SPARof1'][0]
235 249 self.SPARswt=header['SPARswt'][0]
236 250 self.SPARsum=header['SPARsum'][0]
237 251 self.SPARosc=header['SPARosc'][0]
238 252 self.SPARtst=header['SPARtst'][0]
239 253 self.SPARcor=header['SPARcor'][0]
240 254 self.SPARofs=header['SPARofs'][0]
241 255 self.SPARhsn=header['SPARhsn'][0]
242 256 self.SPARhsa=header['SPARhsa'][0]
243 257 self.SPARcalibPow_M=header['SPARcalibPow_M'][0]
244 258 self.SPARcalibSNR_M=header['SPARcalibSNR_M'][0]
245 259 self.SPARcalibPow_S=header['SPARcalibPow_S'][0]
246 260 self.SPARcalibSNR_S=header['SPARcalibSNR_S'][0]
247 261 self.SPARrawGate1=header['SPARrawGate1'][0]
248 262 self.SPARrawGate2=header['SPARrawGate2'][0]
249 263 self.SPARraw=header['SPARraw'][0]
250 264 self.SPARprc=header['SPARprc'][0]
251 265 #112
252 266 #1180
253 267 #print 'Pointer fp header', fp.tell()
254 268 #print ' '
255 269 #print 'SPARrawGate'
256 270 #print self.SPARrawGate2 - self.SPARrawGate1
257 271
258 272 #print ' '
259 273 #print 'Hname'
260 274 #print self.Hname
261 275
262 276 #print ' '
263 277 #print 'Msign'
264 278 #print self.Msign
265 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 288
275
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
309
296 310 SRVI_HEADER = numpy.dtype([
297 ('SignatureSRVI1',numpy.str_,4),#
298 ('SizeOfDataBlock1','<i4'),#
299 ('DataBlockTitleSRVI1',numpy.str_,4),#
300 ('SizeOfSRVI1','<i4'),])#
311 ('SignatureSRVI1', numpy.str_, 4),
312 ('SizeOfDataBlock1', '<i4'),
313 ('DataBlockTitleSRVI1', numpy.str_, 4),
314 ('SizeOfSRVI1', '<i4'), ])
315
301 316
302 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 323 self.SizeOfSRVI1 = SizeOfSRVI1
309 324
310 325 self.SRVIHsize=16
311 326
312 327 def SRVIread(self, fp):
313 328
314 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 335 #16
321 336 print 'Pointer fp SRVIheader', fp.tell()
322 337
323 338
324 339 SRVI_STRUCTURE = numpy.dtype([
325 ('frame_cnt','<u4'),#
340 ('frame_cnt', '<u4'),
326 341 ('time_t','<u4'), #
327 342 ('tpow','<f4'), #
328 343 ('npw1','<f4'), #
329 344 ('npw2','<f4'), #
330 345 ('cpw1','<f4'), #
331 346 ('pcw2','<f4'), #
332 347 ('ps_err','<u4'), #
333 348 ('te_err','<u4'), #
334 349 ('rc_err','<u4'), #
335 350 ('grs1','<u4'), #
336 351 ('grs2','<u4'), #
337 352 ('azipos','<f4'), #
338 353 ('azivel','<f4'), #
339 354 ('elvpos','<f4'), #
340 355 ('elvvel','<f4'), #
341 ('northAngle','<f4'), #
356 ('northAngle', '<f4'),
342 357 ('microsec','<u4'), #
343 358 ('azisetvel','<f4'), #
344 359 ('elvsetpos','<f4'), #
345 360 ('RadarConst','<f4'),]) #
346 361
347 362
348
349
350 363 class RecordHeader(Header):
351 364
352
353 365 def __init__(self, frame_cnt=0, time_t= 0, tpow=0, npw1=0, npw2=0,
354 366 cpw1=0, pcw2=0, ps_err=0, te_err=0, rc_err=0, grs1=0,
355 367 grs2=0, azipos=0, azivel=0, elvpos=0, elvvel=0, northangle=0,
356 368 microsec=0, azisetvel=0, elvsetpos=0, RadarConst=0 , RecCounter=0, Off2StartNxtRec=0):
357 369
358
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 391 self.RHsize=84
381 392 self.RecCounter = RecCounter
382 393 self.Off2StartNxtRec=Off2StartNxtRec
383 394
384 395 def RHread(self, fp):
385 396
386 397 #startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
387 398
388 399 #OffRHeader= 1180 + self.RecCounter*(self.Off2StartNxtRec)
389 400 #startFp.seek(OffRHeader, os.SEEK_SET)
390 401
391 402 #print 'Posicion del bloque: ',OffRHeader
392 403
393 404 header = numpy.fromfile(fp,SRVI_STRUCTURE,1)
394 405
395 self.frame_cnt = header['frame_cnt'][0]#
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 427 #84
417 428
418 429 #print 'Pointer fp RECheader', fp.tell()
419 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 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 443
444
434 445 class MIRA35CReader (ProcessingUnit,FileHeaderMIRA35c,SRVIHeader,RecordHeader):
435 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 454
444
445 455 fileList= None
446 456
447 457 #metadata
448 458 TimeZone= None
449 459 Interval= None
450 460 heightList= None
451 461
452 462 #data
453 463 data= None
454 464 utctime= None
455 465
456
457
458 466 def __init__(self, **kwargs):
459 467
460 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 475 self.fpFile=None
468 476 self.fp = None
469 477 self.BlockCounter=0
470 478 self.dtype = None
471 479 self.fileSizeByHeader = None
472 480 self.filenameList = []
473 481 self.fileSelector = 0
474 482 self.Off2StartNxtRec=0
475 483 self.RecCounter=0
476 484 self.flagNoMoreFiles = 0
477 485 self.data_spc=None
478 486 #self.data_cspc=None
479 487 self.data_output=None
480 488 self.path = None
481 489 self.OffsetStartHeader=0
482 490 self.Off2StartData=0
483 491 self.ipp = 0
484 492 self.nFDTdataRecors=0
485 493 self.blocksize = 0
486 494 self.dataOut = Spectra()
487 495 self.profileIndex = 1 #Always
488 496 self.dataOut.flagNoData=False
489 497 self.dataOut.nRdPairs = 0
490 498 self.dataOut.pairsList = []
491 499 self.dataOut.data_spc=None
492 500
493 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 508
501
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 514 #self.__checkPath()
508 515
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
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)
511 520 nFiles=0 #File Counter
512 521 FileList=[] #A list is created that will contain the .fdt files
513 522 for IndexFile in ListaData :
514 523 if '.zspca' in IndexFile and '.gz' not in IndexFile:
515 524 FileList.append(IndexFile)
516 525 nFiles+=1
517 526
518 527 #print 'Files2Read'
519 528 #print 'Existen '+str(nFiles)+' archivos .fdt'
520 529
521 530 self.filenameList=FileList #List of files from least to largest by names
522 531
523
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 543
536
537 544 def setup(self, path=None,
538 545 startDate=None,
539 546 endDate=None,
540 547 startTime=None,
541 548 endTime=None,
542 549 walk=True,
543 550 timezone='utc',
544 551 code = None,
545 552 online=False,
546 553 ReadMode=None, **kwargs):
547 554
548 555 self.isConfig = True
549 556
550 557 self.path=path
551 558 self.startDate=startDate
552 559 self.endDate=endDate
553 560 self.startTime=startTime
554 561 self.endTime=endTime
555 562 self.walk=walk
556 563 #self.ReadMode=int(ReadMode)
557 564
558 565 pass
559 566
560
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 580
575 581 self.fp=self.path
576 582 self.Files2Read(self.fp)
577 583 self.readFile(self.fp)
578 584
579 585 self.dataOut.data_spc = self.dataOut_spc#self.data_spc.copy()
580 586 self.dataOut.RadarConst = self.RadarConst
581 587 self.dataOut.data_output=self.data_output
582 588 self.dataOut.noise = self.dataOut.getNoise()
583 589 #print 'ACAAAAAA', self.dataOut.noise
584 590 self.dataOut.data_spc = self.dataOut.data_spc+self.dataOut.noise
585 591 #print 'self.dataOut.noise',self.dataOut.noise
586 592
587
588 593 return self.dataOut.data_spc
589 594
590
591 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 605
602 606 #The address of the folder is generated the name of the .fdt file that will be read
603 607 print "File: ",self.fileSelector+1
604 608
605 609 if self.fileSelector < len(self.filenameList):
606 610
607 self.fpFile=str(fp)+'/'+str(self.filenameList[self.fileSelector])
611 self.fpFile = str(fp) + '/' + \
612 str(self.filenameList[self.fileSelector])
608 613
609 614 if self.nextfileflag==True:
610 615 self.fp = open(self.fpFile,"rb")
611 616 self.nextfileflag==False
612 617
613 618 '''HERE STARTING THE FILE READING'''
614 619
615
616 620 self.fheader = FileHeaderMIRA35c()
617 621 self.fheader.FHread(self.fp) #Bltr FileHeader Reading
618 622
619
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 628
626
627 629 self.Num_inCoh = self.fheader.PPARavc
628 630 self.dataOut.PRF = self.fheader.PPARprf
629 631 self.dataOut.frequency = 34.85*10**9
630 632 self.Lambda = SPEED_OF_LIGHT/self.dataOut.frequency
631 633 self.dataOut.ippSeconds= 1./float(self.dataOut.PRF)
632 634
633 635 pulse_width = self.fheader.PPARpdr * 10**-9
634 636 self.__deltaHeigth = 0.5 * SPEED_OF_LIGHT * pulse_width
635 637
636 self.data_spc = numpy.zeros((self.Num_Hei, self.Num_Bins,2))#
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 642 self.ETA = numpy.zeros(([2,self.Num_Hei]))
641 643
642
643
644 644 self.readBlock() #Block reading
645 645
646 646 else:
647 647 print 'readFile FlagNoData becomes true'
648 648 self.flagNoMoreFiles=True
649 649 self.dataOut.flagNoData = True
650 650 self.FileHeaderFlag == True
651 651 return 0
652 652
653
654
655 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 665 self.fp.seek(self.PointerReader , os.SEEK_SET)
668 666 self.FirstPoint = self.PointerReader
669 667
670 668 else :
671 669 self.FirstPoint = 1180
672 670
673
674
675 671 self.srviHeader = SRVIHeader()
676 672
677 673 self.srviHeader.SRVIread(self.fp) #Se obtiene la cabecera del SRVI
678 674
679 675 self.blocksize = self.srviHeader.SizeOfDataBlock1 # Se obtiene el tamao del bloque
680 676
681 677 if self.blocksize == 148:
682 678 print 'blocksize == 148 bug'
683 679 jump = numpy.fromfile(self.fp,[('jump',numpy.str_,140)] ,1)
684 680
685 self.srviHeader.SRVIread(self.fp) #Se obtiene la cabecera del SRVI
681 # Se obtiene la cabecera del SRVI
682 self.srviHeader.SRVIread(self.fp)
686 683
687 684 if not self.srviHeader.SizeOfSRVI1:
688 685 self.fileSelector+=1
689 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 695
699
700 696 self.dataOut.channelList = range(1)
701 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 702 self.dataOut.timeZone=0
707 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
705 self.dataOut.heightList = self.SPARrawGate1 * self.__deltaHeigth + \
706 numpy.array(range(self.Num_Hei)) * self.__deltaHeigth
712 707
713 708 self.HSDVsign = numpy.fromfile( self.fp, [('HSDV',numpy.str_,4)],1)
714 709 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)
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)
717 714
718 715 self.COFAsign = numpy.fromfile( self.fp, [('COFA',numpy.str_,4)],1)
719 716 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)
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)
722 721
723 self.ZSPCsign = numpy.fromfile(self.fp, [('ZSPCsign',numpy.str_,4)],1)
722 self.ZSPCsign = numpy.fromfile(
723 self.fp, [('ZSPCsign', numpy.str_, 4)], 1)
724 724 self.SizeZSPC = numpy.fromfile(self.fp, [('SizeZSPC','<i4')],1)
725 725
726 726 self.dataOut.HSDV[0]=self.HSDV_Co[:][0]
727 727 self.dataOut.HSDV[1]=self.HSDV_Cx[:][0]
728 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
730 # Number of spectral sub pieces containing significant power
731 nspc = numpy.fromfile(self.fp, [('nspc', 'int16')], 1)[0][0]
731 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
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
735 739
736 740 #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
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
740 745
741 self.data_spc[irg,binIndex:binIndex+nbins,0] = self.data_spc[irg,binIndex:binIndex+nbins,0]+jbin/65530.*jmax
746 self.data_spc[irg, binIndex:binIndex + nbins, 0] = self.data_spc[irg,
747 binIndex:binIndex + nbins, 0] + jbin / 65530. * jmax
742 748
743 749 #Cx_Channel
744 jbin = numpy.fromfile(self.fp, [('jbin','uint16')],nbins)[0][0]
750 jbin = numpy.fromfile(
751 self.fp, [('jbin', 'uint16')], nbins)[0][0]
745 752 jmax = numpy.fromfile(self.fp, [('jmax','float32')],1)[0][0]
746 753
747
748 self.data_spc[irg,binIndex:binIndex+nbins,1] = self.data_spc[irg,binIndex:binIndex+nbins,1]+jbin/65530.*jmax
754 self.data_spc[irg, binIndex:binIndex + nbins, 1] = self.data_spc[irg,
755 binIndex:binIndex + nbins, 1] + jbin / 65530. * jmax
749 756
750 757 for bin in range(self.Num_Bins):
751 758
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]
759 self.data_spc[:, bin, 0] = self.data_spc[:,
760 bin, 0] - self.dataOut.HSDV[:, 0]
755 761
762 self.data_spc[:, bin, 1] = self.data_spc[:,
763 bin, 1] - self.dataOut.HSDV[:, 1]
756 764
757 765 numpy.set_printoptions(threshold='nan')
758 766
759 767 self.data_spc = numpy.where(self.data_spc > 0. , self.data_spc, 0)
760 768
761 769 self.dataOut.COFA = numpy.array([self.COFA_Co , self.COFA_Cx])
762 770
763 771 print ' '
764 772 print 'SPC',numpy.shape(self.dataOut.data_spc)
765 773 #print 'SPC',self.dataOut.data_spc
766 774
767 775 noinor1 = 713031680
768 776 noinor2 = 30
769 777
770 778 npw1 = 1#0**(npw1/10) * noinor1 * noinor2
771 779 npw2 = 1#0**(npw2/10) * noinor1 * noinor2
772 780 self.dataOut.NPW = numpy.array([npw1, npw2])
773 781
774 782 print ' '
775 783
776 784 self.data_spc = numpy.transpose(self.data_spc, (2,1,0))
777 785 self.data_spc = numpy.fft.fftshift(self.data_spc, axes = 1)
778 786
779 787 self.data_spc = numpy.fliplr(self.data_spc)
780 788
781 789 self.data_spc = numpy.where(self.data_spc > 0. , self.data_spc, 0)
782 790 self.dataOut_spc= numpy.ones([1, self.Num_Bins , self.Num_Hei])
783 791 self.dataOut_spc[0,:,:] = self.data_spc[0,:,:]
784 792 #print 'SHAPE', self.dataOut_spc.shape
785 793 #For nyquist correction:
786 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 797
790
791
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