##// END OF EJS Templates
schain
RSS
Bugs fixed: Eliminacion de numpy.abs() en el calculo de las cross-correlaciones y mejoras en los graficos de CrossSpectra
Miguel Valdez -
r227:3959fdcb8502
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,285 +1,202
1 1 import numpy
2 2 import datetime
3 3 import matplotlib
4 4 matplotlib.use("TKAgg")
5 5 import matplotlib.pyplot
6 6 import matplotlib.dates
7 7 #import scitools.numpyutils
8 8 from mpl_toolkits.axes_grid1 import make_axes_locatable
9 9
10 10 from matplotlib.dates import DayLocator, HourLocator, MinuteLocator, SecondLocator, DateFormatter
11 11 from matplotlib.ticker import FuncFormatter
12 12 from matplotlib.ticker import *
13 13
14 #def init(idfigure, wintitle, width, height, facecolor="w"):
15 #
16 # matplotlib.pyplot.ioff()
17 # fig = matplotlib.pyplot.matplotlib.pyplot.figure(num=idfigure, facecolor=facecolor)
18 # fig.canvas.manager.set_window_title(wintitle)
19 # fig.canvas.manager.resize(width, height)
20 # matplotlib.pyplot.ion()
21 #
22 # return fig
23 #
24 #def setWinTitle(fig, title):
25 #
26 # fig.canvas.manager.set_window_title(title)
27 #
28 #def setTitle(idfigure, title):
29 # fig = matplotlib.pyplot.figure(idfigure)
30 # fig.suptitle(title)
31 #
32 #def makeAxes(idfigure, nrow, ncol, xpos, ypos, colspan, rowspan):
33 # fig = matplotlib.pyplot.figure(idfigure)
34 # ax = matplotlib.pyplot.subplot2grid((nrow, ncol), (xpos, ypos), colspan=colspan, rowspan=rowspan)
35 # return ax
36 #
37 #def setTextFromAxes(idfigure, ax, title):
38 # fig = matplotlib.pyplot.figure(idfigure)
39 # ax.annotate(title, xy=(.1, .99),
40 # xycoords='figure fraction',
41 # horizontalalignment='left', verticalalignment='top',
42 # fontsize=10)
43 #
44 #def pline(ax, x, y, xmin, xmax, ymin, ymax, xlabel, ylabel, title, firsttime):
45 #
46 # if firsttime:
47 # ax.plot(x, y)
48 # ax.set_xlim([xmin,xmax])
49 # ax.set_ylim([ymin,ymax])
50 # ax.set_xlabel(xlabel, size=8)
51 # ax.set_ylabel(ylabel, size=8)
52 # ax.set_title(title, size=10)
53 # matplotlib.pyplot.tight_layout()
54 # else:
55 # ax.lines[0].set_data(x,y)
56 #
57 #def draw(idfigure):
58 #
59 # fig = matplotlib.pyplot.figure(idfigure)
60 # fig.canvas.draw()
61 #
62 #def pcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax, xlabel, ylabel, title, firsttime, mesh):
63 #
64 # if firsttime:
65 # divider = make_axes_locatable(ax)
66 # ax_cb = divider.new_horizontal(size="4%", pad=0.05)
67 # fig1 = ax.get_figure()
68 # fig1.add_axes(ax_cb)
69 #
70 # ax.set_xlim([xmin,xmax])
71 # ax.set_ylim([ymin,ymax])
72 # ax.set_xlabel(xlabel)
73 # ax.set_ylabel(ylabel)
74 # ax.set_title(title)
75 # print x
76 # imesh=ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax)
77 # matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
78 # ax_cb.yaxis.tick_right()
79 # for tl in ax_cb.get_yticklabels():
80 # tl.set_visible(True)
81 # ax_cb.yaxis.tick_right()
82 # matplotlib.pyplot.tight_layout()
83 # return imesh
84 # else:
85 ## ax.set_xlim([xmin,xmax])
86 ## ax.set_ylim([ymin,ymax])
87 # ax.set_xlabel(xlabel)
88 # ax.set_ylabel(ylabel)
89 # ax.set_title(title)
90 #
91 # z = z.T
92 ## z = z[0:-1,0:-1]
93 # mesh.set_array(z.ravel())
94 #
95 # return mesh
96
97 14 ###########################################
98 15 #Actualizacion de las funciones del driver
99 16 ###########################################
100 17
101 18 def createFigure(idfigure, wintitle, width, height, facecolor="w"):
102 19
103 20 matplotlib.pyplot.ioff()
104 fig = matplotlib.pyplot.matplotlib.pyplot.figure(num=idfigure, facecolor=facecolor)
21 fig = matplotlib.pyplot.figure(num=idfigure, facecolor=facecolor)
105 22 fig.canvas.manager.set_window_title(wintitle)
106 23 fig.canvas.manager.resize(width, height)
107 24 matplotlib.pyplot.ion()
108 25
109 26 return fig
110 27
111 28 def closeFigure():
112 29
113 30 matplotlib.pyplot.ioff()
114 31 matplotlib.pyplot.show()
115 32
116 33 return
117 34
118 35 def saveFigure(fig, filename):
119 36 fig.savefig(filename)
120 37
121 38 def setWinTitle(fig, title):
122 39
123 40 fig.canvas.manager.set_window_title(title)
124 41
125 42 def setTitle(fig, title):
126 43
127 44 fig.suptitle(title)
128 45
129 46 def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan):
130 47
131 48 matplotlib.pyplot.figure(fig.number)
132 49 axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
133 50 (xpos, ypos),
134 51 colspan=colspan,
135 52 rowspan=rowspan)
136 53 return axes
137 54
138 55 def setAxesText(ax, text):
139 56
140 57 ax.annotate(text,
141 58 xy = (.1, .99),
142 59 xycoords = 'figure fraction',
143 60 horizontalalignment = 'left',
144 61 verticalalignment = 'top',
145 62 fontsize = 10)
146 63
147 64 def printLabels(ax, xlabel, ylabel, title):
148 65
149 66 ax.set_xlabel(xlabel, size=11)
150 67 ax.set_ylabel(ylabel, size=11)
151 68 ax.set_title(title, size=12)
152 69
153 70 def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='',
154 71 ticksize=9, xtick_visible=True, ytick_visible=True,
155 72 nxticks=4, nyticks=10,
156 73 grid=None):
157 74
158 75 """
159 76
160 77 Input:
161 78 grid : None, 'both', 'x', 'y'
162 79 """
163 80
164 81 ax.plot(x, y)
165 82 ax.set_xlim([xmin,xmax])
166 83 ax.set_ylim([ymin,ymax])
167 84
168 85 printLabels(ax, xlabel, ylabel, title)
169 86
170 87 ######################################################
171 88 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/nxticks) + int(xmin)
172 89 ax.set_xticks(xtickspos)
173 90
174 91 for tick in ax.get_xticklabels():
175 92 tick.set_visible(xtick_visible)
176 93
177 94 for tick in ax.xaxis.get_major_ticks():
178 95 tick.label.set_fontsize(ticksize)
179 96
180 97 ######################################################
181 98 for tick in ax.get_yticklabels():
182 99 tick.set_visible(ytick_visible)
183 100
184 101 for tick in ax.yaxis.get_major_ticks():
185 102 tick.label.set_fontsize(ticksize)
186 103
187 104 iplot = ax.lines[-1]
188 105
189 106 ######################################################
190 107 if '0.' in matplotlib.__version__[0:2]:
191 108 print "The matplotlib version has to be updated to 1.1 or newer"
192 109 return iplot
193 110
194 111 if '1.0.' in matplotlib.__version__[0:4]:
195 112 print "The matplotlib version has to be updated to 1.1 or newer"
196 113 return iplot
197 114
198 115 if grid != None:
199 116 ax.grid(b=True, which='major', axis=grid)
200 117
201 118 matplotlib.pyplot.tight_layout()
202 119
203 120 return iplot
204 121
205 122 def pline(iplot, x, y, xlabel='', ylabel='', title=''):
206 123
207 124 ax = iplot.get_axes()
208 125
209 126 printLabels(ax, xlabel, ylabel, title)
210 127
211 128 iplot.set_data(x, y)
212 129
213 130 def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax,
214 131 xlabel='', ylabel='', title='', ticksize = 9,
215 cblabel='', cbsize="5%",
132 colormap='jet',cblabel='', cbsize="5%",
216 133 XAxisAsTime=False):
217 134
218 135 divider = make_axes_locatable(ax)
219 136 ax_cb = divider.new_horizontal(size=cbsize, pad=0.05)
220 137 fig = ax.get_figure()
221 138 fig.add_axes(ax_cb)
222 139
223 140 ax.set_xlim([xmin,xmax])
224 141 ax.set_ylim([ymin,ymax])
225 142
226 143 printLabels(ax, xlabel, ylabel, title)
227 144
228 imesh = ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax)
145 imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
229 146 cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
230 147 cb.set_label(cblabel)
231 148
232 149 # for tl in ax_cb.get_yticklabels():
233 150 # tl.set_visible(True)
234 151
235 152 for tick in ax.yaxis.get_major_ticks():
236 153 tick.label.set_fontsize(ticksize)
237 154
238 155 for tick in ax.xaxis.get_major_ticks():
239 156 tick.label.set_fontsize(ticksize)
240 157
241 158 for tick in cb.ax.get_yticklabels():
242 159 tick.set_fontsize(ticksize)
243 160
244 161 ax_cb.yaxis.tick_right()
245 162
246 163 if '0.' in matplotlib.__version__[0:2]:
247 164 print "The matplotlib version has to be updated to 1.1 or newer"
248 165 return imesh
249 166
250 167 if '1.0.' in matplotlib.__version__[0:4]:
251 168 print "The matplotlib version has to be updated to 1.1 or newer"
252 169 return imesh
253 170
254 171 matplotlib.pyplot.tight_layout()
255 172
256 173 if XAxisAsTime:
257 174
258 175 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
259 176 ax.xaxis.set_major_formatter(FuncFormatter(func))
260 177 ax.xaxis.set_major_locator(LinearLocator(7))
261 178
262 179 return imesh
263 180
264 181 def pcolor(imesh, z, xlabel='', ylabel='', title=''):
265 182
266 183 z = z.T
267 184
268 185 ax = imesh.get_axes()
269 186
270 187 printLabels(ax, xlabel, ylabel, title)
271 188
272 189 imesh.set_array(z.ravel())
273 190
274 191 def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title=''):
275 192
276 193 printLabels(ax, xlabel, ylabel, title)
277 194
278 195 imesh = ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax)
279 196
280 197 def draw(fig):
281 198
282 199 if type(fig) == 'int':
283 200 raise ValueError, "This parameter should be of tpye matplotlib figure"
284 201
285 202 fig.canvas.draw() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,512 +1,512
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import os, sys
8 8 import copy
9 9 import numpy
10 10 import datetime
11 11
12 12 from jroheaderIO import SystemHeader, RadarControllerHeader
13 13
14 14 def hildebrand_sekhon(data, navg):
15 15 """
16 16 This method is for the objective determination of de noise level in Doppler spectra. This
17 17 implementation technique is based on the fact that the standard deviation of the spectral
18 18 densities is equal to the mean spectral density for white Gaussian noise
19 19
20 20 Inputs:
21 21 Data : heights
22 22 navg : numbers of averages
23 23
24 24 Return:
25 25 -1 : any error
26 26 anoise : noise's level
27 27 """
28 28
29 29 dataflat = data.copy().reshape(-1)
30 30 dataflat.sort()
31 31 npts = dataflat.size #numbers of points of the data
32 32 npts_noise = 0.2*npts
33 33
34 34 if npts < 32:
35 35 print "error in noise - requires at least 32 points"
36 36 return -1.0
37 37
38 38 dataflat2 = numpy.power(dataflat,2)
39 39
40 40 cs = numpy.cumsum(dataflat)
41 41 cs2 = numpy.cumsum(dataflat2)
42 42
43 43 # data sorted in ascending order
44 44 nmin = int((npts + 7.)/8)
45 45
46 46 for i in range(nmin, npts):
47 47 s = cs[i]
48 48 s2 = cs2[i]
49 49 p = s / float(i);
50 50 p2 = p**2;
51 51 q = s2 / float(i) - p2;
52 52 leftc = p2;
53 53 rightc = q * float(navg);
54 54 R2 = leftc/rightc
55 55
56 56 # Signal detect: R2 < 1 (R2 = leftc/rightc)
57 57 if R2 < 1:
58 58 npts_noise = i
59 59 break
60 60
61 61
62 62 anoise = numpy.average(dataflat[0:npts_noise])
63 63
64 64 return anoise;
65 65
66 66 def sorting_bruce(data, navg):
67 67
68 68 data = data.copy()
69 69
70 70 sortdata = numpy.sort(data)
71 71 lenOfData = len(data)
72 72 nums_min = lenOfData/10
73 73
74 74 if (lenOfData/10) > 0:
75 75 nums_min = lenOfData/10
76 76 else:
77 77 nums_min = 0
78 78
79 79 rtest = 1.0 + 1.0/navg
80 80
81 81 sum = 0.
82 82
83 83 sumq = 0.
84 84
85 85 j = 0
86 86
87 87 cont = 1
88 88
89 89 while((cont==1)and(j<lenOfData)):
90 90
91 91 sum += sortdata[j]
92 92
93 93 sumq += sortdata[j]**2
94 94
95 95 j += 1
96 96
97 97 if j > nums_min:
98 98 if ((sumq*j) <= (rtest*sum**2)):
99 99 lnoise = sum / j
100 100 else:
101 101 j = j - 1
102 102 sum = sum - sordata[j]
103 103 sumq = sumq - sordata[j]**2
104 104 cont = 0
105 105
106 106 if j == nums_min:
107 107 lnoise = sum /j
108 108
109 109 return lnoise
110 110
111 111 class JROData:
112 112
113 113 # m_BasicHeader = BasicHeader()
114 114 # m_ProcessingHeader = ProcessingHeader()
115 115
116 116 systemHeaderObj = SystemHeader()
117 117
118 118 radarControllerHeaderObj = RadarControllerHeader()
119 119
120 120 # data = None
121 121
122 122 type = None
123 123
124 124 dtype = None
125 125
126 126 # nChannels = None
127 127
128 128 # nHeights = None
129 129
130 130 nProfiles = None
131 131
132 132 heightList = None
133 133
134 134 channelList = None
135 135
136 136 flagNoData = True
137 137
138 138 flagTimeBlock = False
139 139
140 140 utctime = None
141 141
142 142 blocksize = None
143 143
144 144 nCode = None
145 145
146 146 nBaud = None
147 147
148 148 code = None
149 149
150 150 flagDecodeData = True #asumo q la data esta decodificada
151 151
152 152 flagDeflipData = True #asumo q la data esta sin flip
153 153
154 154 flagShiftFFT = False
155 155
156 156 ippSeconds = None
157 157
158 158 timeInterval = None
159 159
160 160 nCohInt = None
161 161
162 162 noise = None
163 163
164 164 #Speed of ligth
165 165 C = 3e8
166 166
167 167 frequency = 49.92e6
168 168
169 169 def __init__(self):
170 170
171 171 raise ValueError, "This class has not been implemented"
172 172
173 173 def copy(self, inputObj=None):
174 174
175 175 if inputObj == None:
176 176 return copy.deepcopy(self)
177 177
178 178 for key in inputObj.__dict__.keys():
179 179 self.__dict__[key] = inputObj.__dict__[key]
180 180
181 181 def deepcopy(self):
182 182
183 183 return copy.deepcopy(self)
184 184
185 185 def isEmpty(self):
186 186
187 187 return self.flagNoData
188 188
189 189 def getNoise(self):
190 190
191 191 raise ValueError, "Not implemented"
192 192
193 193 def getNChannels(self):
194 194
195 195 return len(self.channelList)
196 196
197 197 def getChannelIndexList(self):
198 198
199 199 return range(self.nChannels)
200 200
201 201 def getNHeights(self):
202 202
203 203 return len(self.heightList)
204 204
205 205 def getHeiRange(self, extrapoints=0):
206 206
207 207 heis = self.heightList
208 208 # deltah = self.heightList[1] - self.heightList[0]
209 209 #
210 210 # heis.append(self.heightList[-1])
211 211
212 212 return heis
213 213
214 214 def getDatatime(self):
215 215
216 216 datatime = datetime.datetime.utcfromtimestamp(self.utctime)
217 217 return datatime
218 218
219 219 def getTimeRange(self):
220 220
221 221 datatime = []
222 222
223 223 datatime.append(self.utctime)
224 224 datatime.append(self.utctime + self.timeInterval)
225 225
226 226 datatime = numpy.array(datatime)
227 227
228 228 return datatime
229 229
230 230 def getFmax(self):
231 231
232 232 PRF = 1./(self.ippSeconds * self.nCohInt)
233 233
234 234 fmax = PRF/2.
235 235
236 236 return fmax
237 237
238 238 def getVmax(self):
239 239
240 240 _lambda = self.C/self.frequency
241 241
242 vmax = self.getFmax() * _lambda / 2.
242 vmax = self.getFmax() * _lambda
243 243
244 244 return vmax
245 245
246 246 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
247 247 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
248 248 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
249 249 noise = property(getNoise, "I'm the 'nHeights' property.")
250 250 datatime = property(getDatatime, "I'm the 'datatime' property")
251 251
252 252 class Voltage(JROData):
253 253
254 254 #data es un numpy array de 2 dmensiones (canales, alturas)
255 255 data = None
256 256
257 257 def __init__(self):
258 258 '''
259 259 Constructor
260 260 '''
261 261
262 262 self.radarControllerHeaderObj = RadarControllerHeader()
263 263
264 264 self.systemHeaderObj = SystemHeader()
265 265
266 266 self.type = "Voltage"
267 267
268 268 self.data = None
269 269
270 270 self.dtype = None
271 271
272 272 # self.nChannels = 0
273 273
274 274 # self.nHeights = 0
275 275
276 276 self.nProfiles = None
277 277
278 278 self.heightList = None
279 279
280 280 self.channelList = None
281 281
282 282 # self.channelIndexList = None
283 283
284 284 self.flagNoData = True
285 285
286 286 self.flagTimeBlock = False
287 287
288 288 self.utctime = None
289 289
290 290 self.nCohInt = None
291 291
292 292 self.blocksize = None
293 293
294 294 def getNoisebyHildebrand(self):
295 295 """
296 296 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
297 297
298 298 Return:
299 299 noiselevel
300 300 """
301 301
302 302 for channel in range(self.nChannels):
303 303 daux = self.data_spc[channel,:,:]
304 304 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
305 305
306 306 return self.noise
307 307
308 308 def getNoise(self, type = 1):
309 309
310 310 self.noise = numpy.zeros(self.nChannels)
311 311
312 312 if type == 1:
313 313 noise = self.getNoisebyHildebrand()
314 314
315 315 return 10*numpy.log10(noise)
316 316
317 317 class Spectra(JROData):
318 318
319 319 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
320 320 data_spc = None
321 321
322 322 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
323 323 data_cspc = None
324 324
325 325 #data es un numpy array de 2 dmensiones (canales, alturas)
326 326 data_dc = None
327 327
328 328 nFFTPoints = None
329 329
330 330 nPairs = None
331 331
332 332 pairsList = None
333 333
334 334 nIncohInt = None
335 335
336 336 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
337 337
338 338 nCohInt = None #se requiere para determinar el valor de timeInterval
339 339
340 340 def __init__(self):
341 341 '''
342 342 Constructor
343 343 '''
344 344
345 345 self.radarControllerHeaderObj = RadarControllerHeader()
346 346
347 347 self.systemHeaderObj = SystemHeader()
348 348
349 349 self.type = "Spectra"
350 350
351 351 # self.data = None
352 352
353 353 self.dtype = None
354 354
355 355 # self.nChannels = 0
356 356
357 357 # self.nHeights = 0
358 358
359 359 self.nProfiles = None
360 360
361 361 self.heightList = None
362 362
363 363 self.channelList = None
364 364
365 365 # self.channelIndexList = None
366 366
367 367 self.flagNoData = True
368 368
369 369 self.flagTimeBlock = False
370 370
371 371 self.utctime = None
372 372
373 373 self.nCohInt = None
374 374
375 375 self.nIncohInt = None
376 376
377 377 self.blocksize = None
378 378
379 379 self.nFFTPoints = None
380 380
381 381 self.wavelength = None
382 382
383 383 def getNoisebyHildebrand(self):
384 384 """
385 385 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
386 386
387 387 Return:
388 388 noiselevel
389 389 """
390 390
391 391 for channel in range(self.nChannels):
392 392 daux = self.data_spc[channel,:,:]
393 393 self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
394 394
395 395 return self.noise
396 396
397 397 def getNoisebyWindow(self, heiIndexMin=0, heiIndexMax=-1, freqIndexMin=0, freqIndexMax=-1):
398 398 """
399 399 Determina el ruido del canal utilizando la ventana indicada con las coordenadas:
400 400 (heiIndexMIn, freqIndexMin) hasta (heiIndexMax, freqIndexMAx)
401 401
402 402 Inputs:
403 403 heiIndexMin: Limite inferior del eje de alturas
404 404 heiIndexMax: Limite superior del eje de alturas
405 405 freqIndexMin: Limite inferior del eje de frecuencia
406 406 freqIndexMax: Limite supoerior del eje de frecuencia
407 407 """
408 408
409 409 data = self.data_spc[:, heiIndexMin:heiIndexMax, freqIndexMin:freqIndexMax]
410 410
411 411 for channel in range(self.nChannels):
412 412 daux = data[channel,:,:]
413 413 self.noise[channel] = numpy.average(daux)
414 414
415 415 return self.noise
416 416
417 417 def getNoisebySort(self):
418 418
419 419 for channel in range(self.nChannels):
420 420 daux = self.data_spc[channel,:,:]
421 421 self.noise[channel] = sorting_bruce(daux, self.nIncohInt)
422 422
423 423 return self.noise
424 424
425 425 def getNoise(self, type = 1):
426 426
427 427 self.noise = numpy.zeros(self.nChannels)
428 428
429 429 if type == 1:
430 430 noise = self.getNoisebyHildebrand()
431 431
432 432 if type == 2:
433 433 noise = self.getNoisebySort()
434 434
435 435 if type == 3:
436 436 noise = self.getNoisebyWindow()
437 437
438 438 return 10*numpy.log10(noise)
439 439
440 440
441 441 def getFreqRange(self, extrapoints=0):
442 442
443 delfreq = 2 * self.getFmax() / self.nFFTPoints
444 freqrange = deltafreqs*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
443 deltafreq = self.getFmax() / self.nFFTPoints
444 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
445 445
446 446 return freqrange
447 447
448 448 def getVelRange(self, extrapoints=0):
449 449
450 deltav = 2 * self.getVmax() / self.nFFTPoints
450 deltav = self.getVmax() / self.nFFTPoints
451 451 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
452 452
453 453 return velrange
454 454
455 455 def getNPairs(self):
456 456
457 457 return len(self.pairsList)
458 458
459 459 def getPairsIndexList(self):
460 460
461 461 return range(self.nPairs)
462 462
463 463 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
464 464 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
465 465
466 466 class SpectraHeis(JROData):
467 467
468 468 data_spc = None
469 469
470 470 data_cspc = None
471 471
472 472 data_dc = None
473 473
474 474 nFFTPoints = None
475 475
476 476 nPairs = None
477 477
478 478 pairsList = None
479 479
480 480 nIncohInt = None
481 481
482 482 def __init__(self):
483 483
484 484 self.radarControllerHeaderObj = RadarControllerHeader()
485 485
486 486 self.systemHeaderObj = SystemHeader()
487 487
488 488 self.type = "SpectraHeis"
489 489
490 490 self.dtype = None
491 491
492 492 # self.nChannels = 0
493 493
494 494 # self.nHeights = 0
495 495
496 496 self.nProfiles = None
497 497
498 498 self.heightList = None
499 499
500 500 self.channelList = None
501 501
502 502 # self.channelIndexList = None
503 503
504 504 self.flagNoData = True
505 505
506 506 self.flagTimeBlock = False
507 507
508 508 self.nPairs = 0
509 509
510 510 self.utctime = None
511 511
512 512 self.blocksize = None
Show file before | Show file after | Hide comments
@@ -1,589 +1,589
1 1 import numpy
2 2 import time, datetime
3 3 from graphics.figure import *
4 4
5 5 class CrossSpectraPlot(Figure):
6 6
7 7 __isConfig = None
8 8 __nsubplots = None
9 9
10 10 WIDTHPROF = None
11 11 HEIGHTPROF = None
12 12 PREFIX = 'spc'
13 13
14 14 def __init__(self):
15 15
16 16 self.__isConfig = False
17 17 self.__nsubplots = 4
18 18
19 19 self.WIDTH = 300
20 20 self.HEIGHT = 400
21 21 self.WIDTHPROF = 0
22 22 self.HEIGHTPROF = 0
23 23
24 24 def getSubplots(self):
25 25
26 26 ncol = 4
27 27 nrow = self.nplots
28 28
29 29 return nrow, ncol
30 30
31 31 def setup(self, idfigure, nplots, wintitle, showprofile=True):
32 32
33 33 self.__showprofile = showprofile
34 34 self.nplots = nplots
35 35
36 36 ncolspan = 1
37 37 colspan = 1
38 38
39 39 self.createFigure(idfigure = idfigure,
40 40 wintitle = wintitle,
41 41 widthplot = self.WIDTH + self.WIDTHPROF,
42 42 heightplot = self.HEIGHT + self.HEIGHTPROF)
43 43
44 44 nrow, ncol = self.getSubplots()
45 45
46 46 counter = 0
47 47 for y in range(nrow):
48 48 for x in range(ncol):
49 49 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
50 50
51 51 counter += 1
52 52
53 53 def run(self, dataOut, idfigure, wintitle="", pairsList=None, showprofile='True',
54 54 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
55 55 save=False, figpath='./', figfile=None):
56 56
57 57 """
58 58
59 59 Input:
60 60 dataOut :
61 61 idfigure :
62 62 wintitle :
63 63 channelList :
64 64 showProfile :
65 65 xmin : None,
66 66 xmax : None,
67 67 ymin : None,
68 68 ymax : None,
69 69 zmin : None,
70 70 zmax : None
71 71 """
72 72
73 73 if pairsList == None:
74 74 pairsIndexList = dataOut.pairsIndexList
75 75 else:
76 76 pairsIndexList = []
77 77 for pair in pairsList:
78 78 if pair not in dataOut.pairsList:
79 79 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
80 80 pairsIndexList.append(dataOut.pairsList.index(pair))
81 81
82 if pairIndexList == []:
82 if pairsIndexList == []:
83 83 return
84 84
85 x = dataOut.getVelRange(1)
85 x = dataOut.getFreqRange(1)
86 86 y = dataOut.getHeiRange()
87 87 z = 10.*numpy.log10(dataOut.data_spc[:,:,:])
88 88 avg = numpy.average(numpy.abs(z), axis=1)
89 89
90 90 noise = dataOut.getNoise()
91 91
92 92 if not self.__isConfig:
93 93
94 94 nplots = len(pairsIndexList)
95 95
96 96 self.setup(idfigure=idfigure,
97 97 nplots=nplots,
98 98 wintitle=wintitle,
99 99 showprofile=showprofile)
100 100
101 101 if xmin == None: xmin = numpy.nanmin(x)
102 102 if xmax == None: xmax = numpy.nanmax(x)
103 103 if ymin == None: ymin = numpy.nanmin(y)
104 104 if ymax == None: ymax = numpy.nanmax(y)
105 105 if zmin == None: zmin = numpy.nanmin(avg)*0.9
106 106 if zmax == None: zmax = numpy.nanmax(avg)*0.9
107 107
108 108 self.__isConfig = True
109 109
110 110 thisDatetime = dataOut.datatime
111 111 title = "Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
112 112 xlabel = "Velocity (m/s)"
113 113 ylabel = "Range (Km)"
114 114
115 115 self.setWinTitle(title)
116 116
117 117 for i in range(self.nplots):
118 118 pair = dataOut.pairsList[pairsIndexList[i]]
119 119
120 120 title = "Channel %d: %4.2fdB" %(pair[0], noise[pair[0]])
121 121 z = 10.*numpy.log10(dataOut.data_spc[pair[0],:,:])
122 122 axes0 = self.axesList[i*self.__nsubplots]
123 123 axes0.pcolor(x, y, z,
124 124 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
125 125 xlabel=xlabel, ylabel=ylabel, title=title,
126 126 ticksize=9, cblabel='')
127 127
128 128 title = "Channel %d: %4.2fdB" %(pair[1], noise[pair[1]])
129 129 z = 10.*numpy.log10(dataOut.data_spc[pair[1],:,:])
130 130 axes0 = self.axesList[i*self.__nsubplots+1]
131 131 axes0.pcolor(x, y, z,
132 132 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
133 133 xlabel=xlabel, ylabel=ylabel, title=title,
134 134 ticksize=9, cblabel='')
135 135
136 136 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:])
137 137 coherence = numpy.abs(coherenceComplex)
138 138 phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
139 139
140 140
141 141 title = "Coherence %d%d" %(pair[0], pair[1])
142 142 axes0 = self.axesList[i*self.__nsubplots+2]
143 143 axes0.pcolor(x, y, coherence,
144 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=-1, zmax=1,
144 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=1,
145 145 xlabel=xlabel, ylabel=ylabel, title=title,
146 146 ticksize=9, cblabel='')
147 147
148 148 title = "Phase %d%d" %(pair[0], pair[1])
149 149 axes0 = self.axesList[i*self.__nsubplots+3]
150 150 axes0.pcolor(x, y, phase,
151 151 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
152 152 xlabel=xlabel, ylabel=ylabel, title=title,
153 ticksize=9, cblabel='')
153 ticksize=9, cblabel='', colormap='RdBu')
154 154
155 155
156 156
157 157 self.draw()
158 158
159 159 if save:
160 160 date = thisDatetime.strftime("%Y%m%d")
161 161 if figfile == None:
162 162 figfile = self.getFilename(name = date)
163 163
164 164 self.saveFigure(figpath, figfile)
165 165
166 166
167 167 class RTIPlot(Figure):
168 168
169 169 __isConfig = None
170 170 __nsubplots = None
171 171
172 172 WIDTHPROF = None
173 173 HEIGHTPROF = None
174 174 PREFIX = 'rti'
175 175
176 176 def __init__(self):
177 177
178 178 self.__timerange = 24*60*60
179 179 self.__isConfig = False
180 180 self.__nsubplots = 1
181 181
182 182 self.WIDTH = 800
183 183 self.HEIGHT = 200
184 184 self.WIDTHPROF = 120
185 185 self.HEIGHTPROF = 0
186 186
187 187 def getSubplots(self):
188 188
189 189 ncol = 1
190 190 nrow = self.nplots
191 191
192 192 return nrow, ncol
193 193
194 194 def setup(self, idfigure, nplots, wintitle, showprofile=True):
195 195
196 196 self.__showprofile = showprofile
197 197 self.nplots = nplots
198 198
199 199 ncolspan = 1
200 200 colspan = 1
201 201 if showprofile:
202 202 ncolspan = 7
203 203 colspan = 6
204 204 self.__nsubplots = 2
205 205
206 206 self.createFigure(idfigure = idfigure,
207 207 wintitle = wintitle,
208 208 widthplot = self.WIDTH + self.WIDTHPROF,
209 209 heightplot = self.HEIGHT + self.HEIGHTPROF)
210 210
211 211 nrow, ncol = self.getSubplots()
212 212
213 213 counter = 0
214 214 for y in range(nrow):
215 215 for x in range(ncol):
216 216
217 217 if counter >= self.nplots:
218 218 break
219 219
220 220 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
221 221
222 222 if showprofile:
223 223 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
224 224
225 225 counter += 1
226 226
227 227 def __getTimeLim(self, x, xmin, xmax):
228 228
229 229 thisdatetime = datetime.datetime.utcfromtimestamp(numpy.min(x))
230 230 thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0))
231 231
232 232 ####################################################
233 233 #If the x is out of xrange
234 234 if xmax < (thisdatetime - thisdate).seconds/(60*60.):
235 235 xmin = None
236 236 xmax = None
237 237
238 238 if xmin == None:
239 239 td = thisdatetime - thisdate
240 240 xmin = td.seconds/(60*60.)
241 241
242 242 if xmax == None:
243 243 xmax = xmin + self.__timerange/(60*60.)
244 244
245 245 mindt = thisdate + datetime.timedelta(0,0,0,0,0, xmin)
246 246 tmin = time.mktime(mindt.timetuple())
247 247
248 248 maxdt = thisdate + datetime.timedelta(0,0,0,0,0, xmax)
249 249 tmax = time.mktime(maxdt.timetuple())
250 250
251 251 self.__timerange = tmax - tmin
252 252
253 253 return tmin, tmax
254 254
255 255 def run(self, dataOut, idfigure, wintitle="", channelList=None, showprofile='True',
256 256 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
257 257 timerange=None,
258 258 save=False, figpath='./', figfile=None):
259 259
260 260 """
261 261
262 262 Input:
263 263 dataOut :
264 264 idfigure :
265 265 wintitle :
266 266 channelList :
267 267 showProfile :
268 268 xmin : None,
269 269 xmax : None,
270 270 ymin : None,
271 271 ymax : None,
272 272 zmin : None,
273 273 zmax : None
274 274 """
275 275
276 276 if channelList == None:
277 277 channelIndexList = dataOut.channelIndexList
278 278 else:
279 279 channelIndexList = []
280 280 for channel in channelList:
281 281 if channel not in dataOut.channelList:
282 282 raise ValueError, "Channel %d is not in dataOut.channelList"
283 283 channelIndexList.append(dataOut.channelList.index(channel))
284 284
285 285 if timerange != None:
286 286 self.__timerange = timerange
287 287
288 288 tmin = None
289 289 tmax = None
290 290 x = dataOut.getTimeRange()
291 291 y = dataOut.getHeiRange()
292 292 z = 10.*numpy.log10(dataOut.data_spc[channelIndexList,:,:])
293 293 avg = numpy.average(z, axis=1)
294 294
295 295 noise = dataOut.getNoise()
296 296
297 297 if not self.__isConfig:
298 298
299 299 nplots = len(channelIndexList)
300 300
301 301 self.setup(idfigure=idfigure,
302 302 nplots=nplots,
303 303 wintitle=wintitle,
304 304 showprofile=showprofile)
305 305
306 306 tmin, tmax = self.__getTimeLim(x, xmin, xmax)
307 307 if ymin == None: ymin = numpy.nanmin(y)
308 308 if ymax == None: ymax = numpy.nanmax(y)
309 309 if zmin == None: zmin = numpy.nanmin(avg)*0.9
310 310 if zmax == None: zmax = numpy.nanmax(avg)*0.9
311 311
312 312 self.__isConfig = True
313 313
314 314 thisDatetime = dataOut.datatime
315 315 title = "RTI: %s" %(thisDatetime.strftime("%d-%b-%Y"))
316 316 xlabel = "Velocity (m/s)"
317 317 ylabel = "Range (Km)"
318 318
319 319 self.setWinTitle(title)
320 320
321 321 for i in range(self.nplots):
322 322 title = "Channel %d: %s" %(dataOut.channelList[i], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
323 323 axes = self.axesList[i*self.__nsubplots]
324 324 z = avg[i].reshape((1,-1))
325 325 axes.pcolor(x, y, z,
326 326 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
327 327 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
328 328 ticksize=9, cblabel='', cbsize="1%")
329 329
330 330 if self.__showprofile:
331 331 axes = self.axesList[i*self.__nsubplots +1]
332 332 axes.pline(avg[i], y,
333 333 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
334 334 xlabel='dB', ylabel='', title='',
335 335 ytick_visible=False,
336 336 grid='x')
337 337
338 338 self.draw()
339 339
340 340 if save:
341 341 date = thisDatetime.strftime("%Y%m%d")
342 342 if figfile == None:
343 343 figfile = self.getFilename(name = date)
344 344
345 345 self.saveFigure(figpath, figfile)
346 346
347 347 if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax:
348 348 self.__isConfig = False
349 349
350 350 class SpectraPlot(Figure):
351 351
352 352 __isConfig = None
353 353 __nsubplots = None
354 354
355 355 WIDTHPROF = None
356 356 HEIGHTPROF = None
357 357 PREFIX = 'spc'
358 358
359 359 def __init__(self):
360 360
361 361 self.__isConfig = False
362 362 self.__nsubplots = 1
363 363
364 364 self.WIDTH = 300
365 365 self.HEIGHT = 400
366 366 self.WIDTHPROF = 120
367 367 self.HEIGHTPROF = 0
368 368
369 369 def getSubplots(self):
370 370
371 371 ncol = int(numpy.sqrt(self.nplots)+0.9)
372 372 nrow = int(self.nplots*1./ncol + 0.9)
373 373
374 374 return nrow, ncol
375 375
376 376 def setup(self, idfigure, nplots, wintitle, showprofile=True):
377 377
378 378 self.__showprofile = showprofile
379 379 self.nplots = nplots
380 380
381 381 ncolspan = 1
382 382 colspan = 1
383 383 if showprofile:
384 384 ncolspan = 3
385 385 colspan = 2
386 386 self.__nsubplots = 2
387 387
388 388 self.createFigure(idfigure = idfigure,
389 389 wintitle = wintitle,
390 390 widthplot = self.WIDTH + self.WIDTHPROF,
391 391 heightplot = self.HEIGHT + self.HEIGHTPROF)
392 392
393 393 nrow, ncol = self.getSubplots()
394 394
395 395 counter = 0
396 396 for y in range(nrow):
397 397 for x in range(ncol):
398 398
399 399 if counter >= self.nplots:
400 400 break
401 401
402 402 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
403 403
404 404 if showprofile:
405 405 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
406 406
407 407 counter += 1
408 408
409 409 def run(self, dataOut, idfigure, wintitle="", channelList=None, showprofile='True',
410 410 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
411 411 save=False, figpath='./', figfile=None):
412 412
413 413 """
414 414
415 415 Input:
416 416 dataOut :
417 417 idfigure :
418 418 wintitle :
419 419 channelList :
420 420 showProfile :
421 421 xmin : None,
422 422 xmax : None,
423 423 ymin : None,
424 424 ymax : None,
425 425 zmin : None,
426 426 zmax : None
427 427 """
428 428
429 429 if channelList == None:
430 430 channelIndexList = dataOut.channelIndexList
431 431 else:
432 432 channelIndexList = []
433 433 for channel in channelList:
434 434 if channel not in dataOut.channelList:
435 435 raise ValueError, "Channel %d is not in dataOut.channelList"
436 436 channelIndexList.append(dataOut.channelList.index(channel))
437 437
438 438 x = dataOut.getVelRange(1)
439 439 y = dataOut.getHeiRange()
440 440 z = 10.*numpy.log10(dataOut.data_spc[channelIndexList,:,:])
441 441 avg = numpy.average(z, axis=1)
442 442
443 443 noise = dataOut.getNoise()
444 444
445 445 if not self.__isConfig:
446 446
447 447 nplots = len(channelIndexList)
448 448
449 449 self.setup(idfigure=idfigure,
450 450 nplots=nplots,
451 451 wintitle=wintitle,
452 452 showprofile=showprofile)
453 453
454 454 if xmin == None: xmin = numpy.nanmin(x)
455 455 if xmax == None: xmax = numpy.nanmax(x)
456 456 if ymin == None: ymin = numpy.nanmin(y)
457 457 if ymax == None: ymax = numpy.nanmax(y)
458 458 if zmin == None: zmin = numpy.nanmin(avg)*0.9
459 459 if zmax == None: zmax = numpy.nanmax(avg)*0.9
460 460
461 461 self.__isConfig = True
462 462
463 463 thisDatetime = dataOut.datatime
464 464 title = "Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
465 465 xlabel = "Velocity (m/s)"
466 466 ylabel = "Range (Km)"
467 467
468 468 self.setWinTitle(title)
469 469
470 470 for i in range(self.nplots):
471 471 title = "Channel %d: %4.2fdB" %(dataOut.channelList[i], noise[i])
472 472 axes = self.axesList[i*self.__nsubplots]
473 473 axes.pcolor(x, y, z[i,:,:],
474 474 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
475 475 xlabel=xlabel, ylabel=ylabel, title=title,
476 476 ticksize=9, cblabel='')
477 477
478 478 if self.__showprofile:
479 479 axes = self.axesList[i*self.__nsubplots +1]
480 480 axes.pline(avg[i], y,
481 481 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
482 482 xlabel='dB', ylabel='', title='',
483 483 ytick_visible=False,
484 484 grid='x')
485 485
486 486 self.draw()
487 487
488 488 if save:
489 489 date = thisDatetime.strftime("%Y%m%d")
490 490 if figfile == None:
491 491 figfile = self.getFilename(name = date)
492 492
493 493 self.saveFigure(figpath, figfile)
494 494
495 495 class Scope(Figure):
496 496
497 497 __isConfig = None
498 498
499 499 def __init__(self):
500 500
501 501 self.__isConfig = False
502 502 self.WIDTH = 600
503 503 self.HEIGHT = 200
504 504
505 505 def getSubplots(self):
506 506
507 507 nrow = self.nplots
508 508 ncol = 3
509 509 return nrow, ncol
510 510
511 511 def setup(self, idfigure, nplots, wintitle):
512 512
513 513 self.createFigure(idfigure, wintitle)
514 514
515 515 nrow,ncol = self.getSubplots()
516 516 colspan = 3
517 517 rowspan = 1
518 518
519 519 for i in range(nplots):
520 520 self.addAxes(nrow, ncol, i, 0, colspan, rowspan)
521 521
522 522 self.nplots = nplots
523 523
524 524 def run(self, dataOut, idfigure, wintitle="", channelList=None,
525 525 xmin=None, xmax=None, ymin=None, ymax=None, save=False, filename=None):
526 526
527 527 """
528 528
529 529 Input:
530 530 dataOut :
531 531 idfigure :
532 532 wintitle :
533 533 channelList :
534 534 xmin : None,
535 535 xmax : None,
536 536 ymin : None,
537 537 ymax : None,
538 538 """
539 539
540 540 if channelList == None:
541 541 channelIndexList = dataOut.channelIndexList
542 542 else:
543 543 channelIndexList = []
544 544 for channel in channelList:
545 545 if channel not in dataOut.channelList:
546 546 raise ValueError, "Channel %d is not in dataOut.channelList"
547 547 channelIndexList.append(dataOut.channelList.index(channel))
548 548
549 549 x = dataOut.heightList
550 550 y = dataOut.data[channelList,:] * numpy.conjugate(dataOut.data[channelList,:])
551 551 y = y.real
552 552
553 553 noise = dataOut.getNoise()
554 554
555 555 if not self.__isConfig:
556 556 nplots = len(channelList)
557 557
558 558 self.setup(idfigure=idfigure,
559 559 nplots=nplots,
560 560 wintitle=wintitle)
561 561
562 562 if xmin == None: xmin = numpy.nanmin(x)
563 563 if xmax == None: xmax = numpy.nanmax(x)
564 564 if ymin == None: ymin = numpy.nanmin(y)
565 565 if ymax == None: ymax = numpy.nanmax(y)
566 566
567 567 self.__isConfig = True
568 568
569 569
570 570 thisDatetime = dataOut.datatime
571 571 title = "Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
572 572 xlabel = "Range (Km)"
573 573 ylabel = "Intensity"
574 574
575 575 self.setWinTitle(title)
576 576
577 577 for i in range(len(self.axesList)):
578 578 title = "Channel %d: %4.2fdB" %(i, noise[i])
579 579 axes = self.axesList[i]
580 580 ychannel = y[i,:]
581 581 axes.pline(x, ychannel,
582 582 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
583 583 xlabel=xlabel, ylabel=ylabel, title=title)
584 584
585 585 self.draw()
586 586
587 587 if save:
588 588 self.saveFigure(filename)
589 589 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1145 +1,1145
1 1 '''
2 2
3 3 $Author: dsuarez $
4 4 $Id: Processor.py 1 2012-11-12 18:56:07Z dsuarez $
5 5 '''
6 6 import os
7 7 import numpy
8 8 import datetime
9 9 import time
10 10
11 11 from jrodata import *
12 12 from jrodataIO import *
13 13 from jroplot import *
14 14
15 15 class ProcessingUnit:
16 16
17 17 """
18 18 Esta es la clase base para el procesamiento de datos.
19 19
20 20 Contiene el metodo "call" para llamar operaciones. Las operaciones pueden ser:
21 21 - Metodos internos (callMethod)
22 22 - Objetos del tipo Operation (callObject). Antes de ser llamados, estos objetos
23 23 tienen que ser agreagados con el metodo "add".
24 24
25 25 """
26 26 # objeto de datos de entrada (Voltage, Spectra o Correlation)
27 27 dataIn = None
28 28
29 29 # objeto de datos de entrada (Voltage, Spectra o Correlation)
30 30 dataOut = None
31 31
32 32
33 33 objectDict = None
34 34
35 35 def __init__(self):
36 36
37 37 self.objectDict = {}
38 38
39 39 def init(self):
40 40
41 41 raise ValueError, "Not implemented"
42 42
43 43 def addOperation(self, object, objId):
44 44
45 45 """
46 46 Agrega el objeto "object" a la lista de objetos "self.objectList" y retorna el
47 47 identificador asociado a este objeto.
48 48
49 49 Input:
50 50
51 51 object : objeto de la clase "Operation"
52 52
53 53 Return:
54 54
55 55 objId : identificador del objeto, necesario para ejecutar la operacion
56 56 """
57 57
58 58 self.objectDict[objId] = object
59 59
60 60 return objId
61 61
62 62 def operation(self, **kwargs):
63 63
64 64 """
65 65 Operacion directa sobre la data (dataout.data). Es necesario actualizar los valores de los
66 66 atributos del objeto dataOut
67 67
68 68 Input:
69 69
70 70 **kwargs : Diccionario de argumentos de la funcion a ejecutar
71 71 """
72 72
73 73 raise ValueError, "ImplementedError"
74 74
75 75 def callMethod(self, name, **kwargs):
76 76
77 77 """
78 78 Ejecuta el metodo con el nombre "name" y con argumentos **kwargs de la propia clase.
79 79
80 80 Input:
81 81 name : nombre del metodo a ejecutar
82 82
83 83 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
84 84
85 85 """
86 86 if name != 'run':
87 87
88 88 if name == 'init' and self.dataIn.isEmpty():
89 89 self.dataOut.flagNoData = True
90 90 return False
91 91
92 92 if name != 'init' and self.dataOut.isEmpty():
93 93 return False
94 94
95 95 methodToCall = getattr(self, name)
96 96
97 97 methodToCall(**kwargs)
98 98
99 99 if name != 'run':
100 100 return True
101 101
102 102 if self.dataOut.isEmpty():
103 103 return False
104 104
105 105 return True
106 106
107 107 def callObject(self, objId, **kwargs):
108 108
109 109 """
110 110 Ejecuta la operacion asociada al identificador del objeto "objId"
111 111
112 112 Input:
113 113
114 114 objId : identificador del objeto a ejecutar
115 115
116 116 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
117 117
118 118 Return:
119 119
120 120 None
121 121 """
122 122
123 123 if self.dataOut.isEmpty():
124 124 return False
125 125
126 126 object = self.objectDict[objId]
127 127
128 128 object.run(self.dataOut, **kwargs)
129 129
130 130 return True
131 131
132 132 def call(self, operationConf, **kwargs):
133 133
134 134 """
135 135 Return True si ejecuta la operacion "operationConf.name" con los
136 136 argumentos "**kwargs". False si la operacion no se ha ejecutado.
137 137 La operacion puede ser de dos tipos:
138 138
139 139 1. Un metodo propio de esta clase:
140 140
141 141 operation.type = "self"
142 142
143 143 2. El metodo "run" de un objeto del tipo Operation o de un derivado de ella:
144 144 operation.type = "other".
145 145
146 146 Este objeto de tipo Operation debe de haber sido agregado antes con el metodo:
147 147 "addOperation" e identificado con el operation.id
148 148
149 149
150 150 con el id de la operacion.
151 151
152 152 Input:
153 153
154 154 Operation : Objeto del tipo operacion con los atributos: name, type y id.
155 155
156 156 """
157 157
158 158 if operationConf.type == 'self':
159 159 sts = self.callMethod(operationConf.name, **kwargs)
160 160
161 161 if operationConf.type == 'other':
162 162 sts = self.callObject(operationConf.id, **kwargs)
163 163
164 164 return sts
165 165
166 166 def setInput(self, dataIn):
167 167
168 168 self.dataIn = dataIn
169 169
170 170 def getOutput(self):
171 171
172 172 return self.dataOut
173 173
174 174 class Operation():
175 175
176 176 """
177 177 Clase base para definir las operaciones adicionales que se pueden agregar a la clase ProcessingUnit
178 178 y necesiten acumular informacion previa de los datos a procesar. De preferencia usar un buffer de
179 179 acumulacion dentro de esta clase
180 180
181 181 Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer)
182 182
183 183 """
184 184
185 185 __buffer = None
186 186 __isConfig = False
187 187
188 188 def __init__(self):
189 189
190 190 pass
191 191
192 192 def run(self, dataIn, **kwargs):
193 193
194 194 """
195 195 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los atributos del objeto dataIn.
196 196
197 197 Input:
198 198
199 199 dataIn : objeto del tipo JROData
200 200
201 201 Return:
202 202
203 203 None
204 204
205 205 Affected:
206 206 __buffer : buffer de recepcion de datos.
207 207
208 208 """
209 209
210 210 raise ValueError, "ImplementedError"
211 211
212 212 class VoltageProc(ProcessingUnit):
213 213
214 214
215 215 def __init__(self):
216 216
217 217 self.objectDict = {}
218 218 self.dataOut = Voltage()
219 219
220 220 def init(self):
221 221
222 222 self.dataOut.copy(self.dataIn)
223 223 # No necesita copiar en cada init() los atributos de dataIn
224 224 # la copia deberia hacerse por cada nuevo bloque de datos
225 225
226 226 def selectChannels(self, channelList):
227 227
228 228 channelIndexList = []
229 229
230 230 for channel in channelList:
231 231 index = self.dataOut.channelList.index(channel)
232 232 channelIndexList.append(index)
233 233
234 234 self.selectChannelsByIndex(channelIndexList)
235 235
236 236 def selectChannelsByIndex(self, channelIndexList):
237 237 """
238 238 Selecciona un bloque de datos en base a canales segun el channelIndexList
239 239
240 240 Input:
241 241 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
242 242
243 243 Affected:
244 244 self.dataOut.data
245 245 self.dataOut.channelIndexList
246 246 self.dataOut.nChannels
247 247 self.dataOut.m_ProcessingHeader.totalSpectra
248 248 self.dataOut.systemHeaderObj.numChannels
249 249 self.dataOut.m_ProcessingHeader.blockSize
250 250
251 251 Return:
252 252 None
253 253 """
254 254
255 255 for channelIndex in channelIndexList:
256 256 if channelIndex not in self.dataOut.channelIndexList:
257 257 print channelIndexList
258 258 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
259 259
260 260 nChannels = len(channelIndexList)
261 261
262 262 data = self.dataOut.data[channelIndexList,:]
263 263
264 264 self.dataOut.data = data
265 265 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
266 266 # self.dataOut.nChannels = nChannels
267 267
268 268 return 1
269 269
270 270 def selectHeights(self, minHei, maxHei):
271 271 """
272 272 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
273 273 minHei <= height <= maxHei
274 274
275 275 Input:
276 276 minHei : valor minimo de altura a considerar
277 277 maxHei : valor maximo de altura a considerar
278 278
279 279 Affected:
280 280 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
281 281
282 282 Return:
283 283 1 si el metodo se ejecuto con exito caso contrario devuelve 0
284 284 """
285 285 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
286 286 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
287 287
288 288 if (maxHei > self.dataOut.heightList[-1]):
289 289 maxHei = self.dataOut.heightList[-1]
290 290 # raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
291 291
292 292 minIndex = 0
293 293 maxIndex = 0
294 294 data = self.dataOut.heightList
295 295
296 296 for i,val in enumerate(data):
297 297 if val < minHei:
298 298 continue
299 299 else:
300 300 minIndex = i;
301 301 break
302 302
303 303 for i,val in enumerate(data):
304 304 if val <= maxHei:
305 305 maxIndex = i;
306 306 else:
307 307 break
308 308
309 309 self.selectHeightsByIndex(minIndex, maxIndex)
310 310
311 311 return 1
312 312
313 313
314 314 def selectHeightsByIndex(self, minIndex, maxIndex):
315 315 """
316 316 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
317 317 minIndex <= index <= maxIndex
318 318
319 319 Input:
320 320 minIndex : valor de indice minimo de altura a considerar
321 321 maxIndex : valor de indice maximo de altura a considerar
322 322
323 323 Affected:
324 324 self.dataOut.data
325 325 self.dataOut.heightList
326 326
327 327 Return:
328 328 1 si el metodo se ejecuto con exito caso contrario devuelve 0
329 329 """
330 330
331 331 if (minIndex < 0) or (minIndex > maxIndex):
332 332 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
333 333
334 334 if (maxIndex >= self.dataOut.nHeights):
335 335 maxIndex = self.dataOut.nHeights-1
336 336 # raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
337 337
338 338 nHeights = maxIndex - minIndex + 1
339 339
340 340 #voltage
341 341 data = self.dataOut.data[:,minIndex:maxIndex+1]
342 342
343 343 firstHeight = self.dataOut.heightList[minIndex]
344 344
345 345 self.dataOut.data = data
346 346 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
347 347
348 348 return 1
349 349
350 350
351 351 def filterByHeights(self, window):
352 352 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
353 353
354 354 if window == None:
355 355 window = self.dataOut.radarControllerHeaderObj.txA / deltaHeight
356 356
357 357 newdelta = deltaHeight * window
358 358 r = self.dataOut.data.shape[1] % window
359 359 buffer = self.dataOut.data[:,0:self.dataOut.data.shape[1]-r]
360 360 buffer = buffer.reshape(self.dataOut.data.shape[0],self.dataOut.data.shape[1]/window,window)
361 361 buffer = numpy.sum(buffer,2)
362 362 self.dataOut.data = buffer
363 363 self.dataOut.heightList = numpy.arange(self.dataOut.heightList[0],newdelta*self.dataOut.nHeights/window,newdelta)
364 364
365 365
366 366 class CohInt(Operation):
367 367
368 368 __profIndex = 0
369 369 __withOverapping = False
370 370
371 371 __byTime = False
372 372 __initime = None
373 373 __lastdatatime = None
374 374 __integrationtime = None
375 375
376 376 __buffer = None
377 377
378 378 __dataReady = False
379 379
380 380 n = None
381 381
382 382
383 383 def __init__(self):
384 384
385 385 self.__isConfig = False
386 386
387 387 def setup(self, n=None, timeInterval=None, overlapping=False):
388 388 """
389 389 Set the parameters of the integration class.
390 390
391 391 Inputs:
392 392
393 393 n : Number of coherent integrations
394 394 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
395 395 overlapping :
396 396
397 397 """
398 398
399 399 self.__initime = None
400 400 self.__lastdatatime = 0
401 401 self.__buffer = None
402 402 self.__dataReady = False
403 403
404 404
405 405 if n == None and timeInterval == None:
406 406 raise ValueError, "n or timeInterval should be specified ..."
407 407
408 408 if n != None:
409 409 self.n = n
410 410 self.__byTime = False
411 411 else:
412 412 self.__integrationtime = timeInterval * 60. #if (type(timeInterval)!=integer) -> change this line
413 413 self.n = 9999
414 414 self.__byTime = True
415 415
416 416 if overlapping:
417 417 self.__withOverapping = True
418 418 self.__buffer = None
419 419 else:
420 420 self.__withOverapping = False
421 421 self.__buffer = 0
422 422
423 423 self.__profIndex = 0
424 424
425 425 def putData(self, data):
426 426
427 427 """
428 428 Add a profile to the __buffer and increase in one the __profileIndex
429 429
430 430 """
431 431
432 432 if not self.__withOverapping:
433 433 self.__buffer += data.copy()
434 434 self.__profIndex += 1
435 435 return
436 436
437 437 #Overlapping data
438 438 nChannels, nHeis = data.shape
439 439 data = numpy.reshape(data, (1, nChannels, nHeis))
440 440
441 441 #If the buffer is empty then it takes the data value
442 442 if self.__buffer == None:
443 443 self.__buffer = data
444 444 self.__profIndex += 1
445 445 return
446 446
447 447 #If the buffer length is lower than n then stakcing the data value
448 448 if self.__profIndex < self.n:
449 449 self.__buffer = numpy.vstack((self.__buffer, data))
450 450 self.__profIndex += 1
451 451 return
452 452
453 453 #If the buffer length is equal to n then replacing the last buffer value with the data value
454 454 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
455 455 self.__buffer[self.n-1] = data
456 456 self.__profIndex = self.n
457 457 return
458 458
459 459
460 460 def pushData(self):
461 461 """
462 462 Return the sum of the last profiles and the profiles used in the sum.
463 463
464 464 Affected:
465 465
466 466 self.__profileIndex
467 467
468 468 """
469 469
470 470 if not self.__withOverapping:
471 471 data = self.__buffer
472 472 n = self.__profIndex
473 473
474 474 self.__buffer = 0
475 475 self.__profIndex = 0
476 476
477 477 return data, n
478 478
479 479 #Integration with Overlapping
480 480 data = numpy.sum(self.__buffer, axis=0)
481 481 n = self.__profIndex
482 482
483 483 return data, n
484 484
485 485 def byProfiles(self, data):
486 486
487 487 self.__dataReady = False
488 488 avgdata = None
489 489 n = None
490 490
491 491 self.putData(data)
492 492
493 493 if self.__profIndex == self.n:
494 494
495 495 avgdata, n = self.pushData()
496 496 self.__dataReady = True
497 497
498 498 return avgdata
499 499
500 500 def byTime(self, data, datatime):
501 501
502 502 self.__dataReady = False
503 503 avgdata = None
504 504 n = None
505 505
506 506 self.putData(data)
507 507
508 508 if (datatime - self.__initime) >= self.__integrationtime:
509 509 avgdata, n = self.pushData()
510 510 self.n = n
511 511 self.__dataReady = True
512 512
513 513 return avgdata
514 514
515 515 def integrate(self, data, datatime=None):
516 516
517 517 if self.__initime == None:
518 518 self.__initime = datatime
519 519
520 520 if self.__byTime:
521 521 avgdata = self.byTime(data, datatime)
522 522 else:
523 523 avgdata = self.byProfiles(data)
524 524
525 525
526 526 self.__lastdatatime = datatime
527 527
528 528 if avgdata == None:
529 529 return None, None
530 530
531 531 avgdatatime = self.__initime
532 532
533 533 deltatime = datatime -self.__lastdatatime
534 534
535 535 if not self.__withOverapping:
536 536 self.__initime = datatime
537 537 else:
538 538 self.__initime += deltatime
539 539
540 540 return avgdata, avgdatatime
541 541
542 542 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
543 543
544 544 if not self.__isConfig:
545 545 self.setup(n, timeInterval, overlapping)
546 546 self.__isConfig = True
547 547
548 548 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
549 549
550 550 # dataOut.timeInterval *= n
551 551 dataOut.flagNoData = True
552 552
553 553 if self.__dataReady:
554 554 dataOut.data = avgdata
555 555 dataOut.nCohInt *= self.n
556 556 dataOut.utctime = avgdatatime
557 557 dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
558 558 dataOut.flagNoData = False
559 559
560 560
561 561 class SpectraProc(ProcessingUnit):
562 562
563 563 def __init__(self):
564 564
565 565 self.objectDict = {}
566 566 self.buffer = None
567 567 self.firstdatatime = None
568 568 self.profIndex = 0
569 569 self.dataOut = Spectra()
570 570
571 571 def __updateObjFromInput(self):
572 572
573 573 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
574 574 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
575 575 self.dataOut.channelList = self.dataIn.channelList
576 576 self.dataOut.heightList = self.dataIn.heightList
577 577 self.dataOut.dtype = self.dataIn.dtype
578 578 # self.dataOut.nHeights = self.dataIn.nHeights
579 579 # self.dataOut.nChannels = self.dataIn.nChannels
580 580 self.dataOut.nBaud = self.dataIn.nBaud
581 581 self.dataOut.nCode = self.dataIn.nCode
582 582 self.dataOut.code = self.dataIn.code
583 583 self.dataOut.nProfiles = self.dataOut.nFFTPoints
584 584 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
585 585 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
586 586 self.dataOut.utctime = self.firstdatatime
587 587 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
588 588 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
589 589 self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
590 590 self.dataOut.nCohInt = self.dataIn.nCohInt
591 591 self.dataOut.nIncohInt = 1
592 592 self.dataOut.ippSeconds = self.dataIn.ippSeconds
593 593
594 594 self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt
595 595
596 596 def __getFft(self):
597 597 """
598 598 Convierte valores de Voltaje a Spectra
599 599
600 600 Affected:
601 601 self.dataOut.data_spc
602 602 self.dataOut.data_cspc
603 603 self.dataOut.data_dc
604 604 self.dataOut.heightList
605 605 self.profIndex
606 606 self.buffer
607 607 self.dataOut.flagNoData
608 608 """
609 fft_volt = numpy.fft.fft(self.buffer,axis=1)
609 fft_volt = numpy.fft.fft(self.buffer,axis=1)/numpy.sqrt(self.dataOut.nFFTPoints)
610 610 dc = fft_volt[:,0,:]
611 611
612 612 #calculo de self-spectra
613 613 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
614 614 spc = fft_volt * numpy.conjugate(fft_volt)
615 615 spc = spc.real
616 616
617 617 blocksize = 0
618 618 blocksize += dc.size
619 619 blocksize += spc.size
620 620
621 621 cspc = None
622 622 pairIndex = 0
623 623 if self.dataOut.pairsList != None:
624 624 #calculo de cross-spectra
625 625 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
626 626 for pair in self.dataOut.pairsList:
627 cspc[pairIndex,:,:] = numpy.abs(fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:]))
627 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
628 628 pairIndex += 1
629 629 blocksize += cspc.size
630 630
631 631 self.dataOut.data_spc = spc
632 632 self.dataOut.data_cspc = cspc
633 633 self.dataOut.data_dc = dc
634 634 self.dataOut.blockSize = blocksize
635 635
636 636 def init(self, nFFTPoints=None, pairsList=None):
637 637
638 638 if self.dataIn.type == "Spectra":
639 639 self.dataOut.copy(self.dataIn)
640 640 return
641 641
642 642 if self.dataIn.type == "Voltage":
643 643
644 644 if nFFTPoints == None:
645 645 raise ValueError, "This SpectraProc.init() need nFFTPoints input variable"
646 646
647 647 if pairsList == None:
648 648 nPairs = 0
649 649 else:
650 650 nPairs = len(pairsList)
651 651
652 652 self.dataOut.nFFTPoints = nFFTPoints
653 653 self.dataOut.pairsList = pairsList
654 654 self.dataOut.nPairs = nPairs
655 655
656 656 if self.buffer == None:
657 657 self.buffer = numpy.zeros((self.dataIn.nChannels,
658 658 self.dataOut.nFFTPoints,
659 659 self.dataIn.nHeights),
660 660 dtype='complex')
661 661
662 662
663 663 self.buffer[:,self.profIndex,:] = self.dataIn.data
664 664 self.profIndex += 1
665 665
666 666 if self.firstdatatime == None:
667 667 self.firstdatatime = self.dataIn.utctime
668 668
669 669 if self.profIndex == self.dataOut.nFFTPoints:
670 670 self.__updateObjFromInput()
671 671 self.__getFft()
672 672
673 673 self.dataOut.flagNoData = False
674 674
675 675 self.buffer = None
676 676 self.firstdatatime = None
677 677 self.profIndex = 0
678 678
679 679 return
680 680
681 681 raise ValuError, "The type object %s is not valid"%(self.dataIn.type)
682 682
683 683 def selectChannels(self, channelList):
684 684
685 685 channelIndexList = []
686 686
687 687 for channel in channelList:
688 688 index = self.dataOut.channelList.index(channel)
689 689 channelIndexList.append(index)
690 690
691 691 self.selectChannelsByIndex(channelIndexList)
692 692
693 693 def selectChannelsByIndex(self, channelIndexList):
694 694 """
695 695 Selecciona un bloque de datos en base a canales segun el channelIndexList
696 696
697 697 Input:
698 698 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
699 699
700 700 Affected:
701 701 self.dataOut.data_spc
702 702 self.dataOut.channelIndexList
703 703 self.dataOut.nChannels
704 704
705 705 Return:
706 706 None
707 707 """
708 708
709 709 for channelIndex in channelIndexList:
710 710 if channelIndex not in self.dataOut.channelIndexList:
711 711 print channelIndexList
712 712 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
713 713
714 714 nChannels = len(channelIndexList)
715 715
716 716 data_spc = self.dataOut.data_spc[channelIndexList,:]
717 717
718 718 self.dataOut.data_spc = data_spc
719 719 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
720 720 # self.dataOut.nChannels = nChannels
721 721
722 722 return 1
723 723
724 724
725 725 class IncohInt(Operation):
726 726
727 727
728 728 __profIndex = 0
729 729 __withOverapping = False
730 730
731 731 __byTime = False
732 732 __initime = None
733 733 __lastdatatime = None
734 734 __integrationtime = None
735 735
736 736 __buffer_spc = None
737 737 __buffer_cspc = None
738 738 __buffer_dc = None
739 739
740 740 __dataReady = False
741 741
742 742 n = None
743 743
744 744
745 745 def __init__(self):
746 746
747 747 self.__isConfig = False
748 748
749 749 def setup(self, n=None, timeInterval=None, overlapping=False):
750 750 """
751 751 Set the parameters of the integration class.
752 752
753 753 Inputs:
754 754
755 755 n : Number of coherent integrations
756 756 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
757 757 overlapping :
758 758
759 759 """
760 760
761 761 self.__initime = None
762 762 self.__lastdatatime = 0
763 763 self.__buffer_spc = None
764 764 self.__buffer_cspc = None
765 765 self.__buffer_dc = None
766 766 self.__dataReady = False
767 767
768 768
769 769 if n == None and timeInterval == None:
770 770 raise ValueError, "n or timeInterval should be specified ..."
771 771
772 772 if n != None:
773 773 self.n = n
774 774 self.__byTime = False
775 775 else:
776 776 self.__integrationtime = timeInterval * 60. #if (type(timeInterval)!=integer) -> change this line
777 777 self.n = 9999
778 778 self.__byTime = True
779 779
780 780 if overlapping:
781 781 self.__withOverapping = True
782 782 else:
783 783 self.__withOverapping = False
784 784 self.__buffer_spc = 0
785 785 self.__buffer_cspc = 0
786 786 self.__buffer_dc = 0
787 787
788 788 self.__profIndex = 0
789 789
790 790 def putData(self, data_spc, data_cspc, data_dc):
791 791
792 792 """
793 793 Add a profile to the __buffer_spc and increase in one the __profileIndex
794 794
795 795 """
796 796
797 797 if not self.__withOverapping:
798 798 self.__buffer_spc += data_spc
799 799
800 800 if data_cspc == None:
801 801 self.__buffer_cspc = None
802 802 else:
803 803 self.__buffer_cspc += data_cspc
804 804
805 805 if data_dc == None:
806 806 self.__buffer_dc = None
807 807 else:
808 808 self.__buffer_dc += data_dc
809 809
810 810 self.__profIndex += 1
811 811 return
812 812
813 813 #Overlapping data
814 814 nChannels, nFFTPoints, nHeis = data_spc.shape
815 815 data_spc = numpy.reshape(data_spc, (1, nChannels, nFFTPoints, nHeis))
816 816 if data_cspc != None:
817 817 data_cspc = numpy.reshape(data_cspc, (1, -1, nFFTPoints, nHeis))
818 818 if data_dc != None:
819 819 data_dc = numpy.reshape(data_dc, (1, -1, nHeis))
820 820
821 821 #If the buffer is empty then it takes the data value
822 822 if self.__buffer_spc == None:
823 823 self.__buffer_spc = data_spc
824 824
825 825 if data_cspc == None:
826 826 self.__buffer_cspc = None
827 827 else:
828 828 self.__buffer_cspc += data_cspc
829 829
830 830 if data_dc == None:
831 831 self.__buffer_dc = None
832 832 else:
833 833 self.__buffer_dc += data_dc
834 834
835 835 self.__profIndex += 1
836 836 return
837 837
838 838 #If the buffer length is lower than n then stakcing the data value
839 839 if self.__profIndex < self.n:
840 840 self.__buffer_spc = numpy.vstack((self.__buffer_spc, data_spc))
841 841
842 842 if data_cspc != None:
843 843 self.__buffer_cspc = numpy.vstack((self.__buffer_cspc, data_cspc))
844 844
845 845 if data_dc != None:
846 846 self.__buffer_dc = numpy.vstack((self.__buffer_dc, data_dc))
847 847
848 848 self.__profIndex += 1
849 849 return
850 850
851 851 #If the buffer length is equal to n then replacing the last buffer value with the data value
852 852 self.__buffer_spc = numpy.roll(self.__buffer_spc, -1, axis=0)
853 853 self.__buffer_spc[self.n-1] = data_spc
854 854
855 855 if data_cspc != None:
856 856 self.__buffer_cspc = numpy.roll(self.__buffer_cspc, -1, axis=0)
857 857 self.__buffer_cspc[self.n-1] = data_cspc
858 858
859 859 if data_dc != None:
860 860 self.__buffer_dc = numpy.roll(self.__buffer_dc, -1, axis=0)
861 861 self.__buffer_dc[self.n-1] = data_dc
862 862
863 863 self.__profIndex = self.n
864 864 return
865 865
866 866
867 867 def pushData(self):
868 868 """
869 869 Return the sum of the last profiles and the profiles used in the sum.
870 870
871 871 Affected:
872 872
873 873 self.__profileIndex
874 874
875 875 """
876 876 data_spc = None
877 877 data_cspc = None
878 878 data_dc = None
879 879
880 880 if not self.__withOverapping:
881 881 data_spc = self.__buffer_spc
882 882 data_cspc = self.__buffer_cspc
883 883 data_dc = self.__buffer_dc
884 884
885 885 n = self.__profIndex
886 886
887 887 self.__buffer_spc = 0
888 888 self.__buffer_cspc = 0
889 889 self.__buffer_dc = 0
890 890 self.__profIndex = 0
891 891
892 892 return data_spc, data_cspc, data_dc, n
893 893
894 894 #Integration with Overlapping
895 895 data_spc = numpy.sum(self.__buffer_spc, axis=0)
896 896
897 897 if self.__buffer_cspc != None:
898 898 data_cspc = numpy.sum(self.__buffer_cspc, axis=0)
899 899
900 900 if self.__buffer_dc != None:
901 901 data_dc = numpy.sum(self.__buffer_dc, axis=0)
902 902
903 903 n = self.__profIndex
904 904
905 905 return data_spc, data_cspc, data_dc, n
906 906
907 907 def byProfiles(self, *args):
908 908
909 909 self.__dataReady = False
910 910 avgdata_spc = None
911 911 avgdata_cspc = None
912 912 avgdata_dc = None
913 913 n = None
914 914
915 915 self.putData(*args)
916 916
917 917 if self.__profIndex == self.n:
918 918
919 919 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
920 920 self.__dataReady = True
921 921
922 922 return avgdata_spc, avgdata_cspc, avgdata_dc
923 923
924 924 def byTime(self, datatime, *args):
925 925
926 926 self.__dataReady = False
927 927 avgdata_spc = None
928 928 avgdata_cspc = None
929 929 avgdata_dc = None
930 930 n = None
931 931
932 932 self.putData(*args)
933 933
934 934 if (datatime - self.__initime) >= self.__integrationtime:
935 935 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
936 936 self.n = n
937 937 self.__dataReady = True
938 938
939 939 return avgdata_spc, avgdata_cspc, avgdata_dc
940 940
941 941 def integrate(self, datatime, *args):
942 942
943 943 if self.__initime == None:
944 944 self.__initime = datatime
945 945
946 946 if self.__byTime:
947 947 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
948 948 else:
949 949 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
950 950
951 951 self.__lastdatatime = datatime
952 952
953 953 if avgdata_spc == None:
954 954 return None, None, None, None
955 955
956 956 avgdatatime = self.__initime
957 957
958 958 deltatime = datatime -self.__lastdatatime
959 959
960 960 if not self.__withOverapping:
961 961 self.__initime = datatime
962 962 else:
963 963 self.__initime += deltatime
964 964
965 965 return avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc
966 966
967 967 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
968 968
969 969 if not self.__isConfig:
970 970 self.setup(n, timeInterval, overlapping)
971 971 self.__isConfig = True
972 972
973 973 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
974 974 dataOut.data_spc,
975 975 dataOut.data_cspc,
976 976 dataOut.data_dc)
977 977
978 978 # dataOut.timeInterval *= n
979 979 dataOut.flagNoData = True
980 980
981 981 if self.__dataReady:
982 982 dataOut.data_spc = avgdata_spc
983 983 dataOut.data_cspc = avgdata_cspc
984 984 dataOut.data_dc = avgdata_dc
985 985
986 986 dataOut.nIncohInt *= self.n
987 987 dataOut.utctime = avgdatatime
988 988 dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt * dataOut.nIncohInt * dataOut.nFFTPoints
989 989 dataOut.flagNoData = False
990 990
991 991
992 992 class ProfileSelector(Operation):
993 993
994 994 profileIndex = None
995 995 # Tamanho total de los perfiles
996 996 nProfiles = None
997 997
998 998 def __init__(self):
999 999
1000 1000 self.profileIndex = 0
1001 1001
1002 1002 def incIndex(self):
1003 1003 self.profileIndex += 1
1004 1004
1005 1005 if self.profileIndex >= self.nProfiles:
1006 1006 self.profileIndex = 0
1007 1007
1008 1008 def isProfileInRange(self, minIndex, maxIndex):
1009 1009
1010 1010 if self.profileIndex < minIndex:
1011 1011 return False
1012 1012
1013 1013 if self.profileIndex > maxIndex:
1014 1014 return False
1015 1015
1016 1016 return True
1017 1017
1018 1018 def isProfileInList(self, profileList):
1019 1019
1020 1020 if self.profileIndex not in profileList:
1021 1021 return False
1022 1022
1023 1023 return True
1024 1024
1025 1025 def run(self, dataOut, profileList=None, profileRangeList=None):
1026 1026
1027 1027 self.nProfiles = dataOut.nProfiles
1028 1028
1029 1029 if profileList != None:
1030 1030 if not(self.isProfileInList(profileList)):
1031 1031 dataOut.flagNoData = True
1032 1032 else:
1033 1033 dataOut.flagNoData = False
1034 1034 self.incIndex()
1035 1035 return 1
1036 1036
1037 1037
1038 1038 elif profileRangeList != None:
1039 1039 minIndex = profileRangeList[0]
1040 1040 maxIndex = profileRangeList[1]
1041 1041 if not(self.isProfileInRange(minIndex, maxIndex)):
1042 1042 dataOut.flagNoData = True
1043 1043 else:
1044 1044 dataOut.flagNoData = False
1045 1045 self.incIndex()
1046 1046 return 1
1047 1047 else:
1048 1048 raise ValueError, "ProfileSelector needs profileList or profileRangeList"
1049 1049
1050 1050 return 0
1051 1051
1052 1052 class Decoder:
1053 1053
1054 1054 data = None
1055 1055 profCounter = None
1056 1056 code = None
1057 1057 ncode = None
1058 1058 nbaud = None
1059 1059 codeIndex = None
1060 1060 flag = False
1061 1061
1062 1062 def __init__(self):
1063 1063
1064 1064 self.data = None
1065 1065 self.ndata = None
1066 1066 self.profCounter = 1
1067 1067 self.codeIndex = 0
1068 1068 self.flag = False
1069 1069 self.code = None
1070 1070 self.ncode = None
1071 1071 self.nbaud = None
1072 1072 self.__isConfig = False
1073 1073
1074 1074 def convolutionInFreq(self, data, ndata):
1075 1075
1076 1076 newcode = numpy.zeros(ndata)
1077 1077 newcode[0:self.nbaud] = self.code[self.codeIndex]
1078 1078
1079 1079 self.codeIndex += 1
1080 1080
1081 1081 fft_data = numpy.fft.fft(data, axis=1)
1082 1082 fft_code = numpy.conj(numpy.fft.fft(newcode))
1083 1083 fft_code = fft_code.reshape(1,len(fft_code))
1084 1084
1085 1085 conv = fft_data.copy()
1086 1086 conv.fill(0)
1087 1087
1088 1088 conv = fft_data*fft_code
1089 1089
1090 1090 data = numpy.fft.ifft(conv,axis=1)
1091 1091 self.data = data[:,:-self.nbaud+1]
1092 1092 self.flag = True
1093 1093
1094 1094 if self.profCounter == self.ncode:
1095 1095 self.profCounter = 0
1096 1096 self.codeIndex = 0
1097 1097
1098 1098 self.profCounter += 1
1099 1099
1100 1100 def convolutionInTime(self, data, ndata):
1101 1101
1102 1102 nchannel = data.shape[1]
1103 1103 newcode = self.code[self.codeIndex]
1104 1104 self.codeIndex += 1
1105 1105 conv = data.copy()
1106 1106 for i in range(nchannel):
1107 1107 conv[i,:] = numpy.correlate(data[i,:], newcode)
1108 1108
1109 1109 self.data = conv
1110 1110 self.flag = True
1111 1111
1112 1112 if self.profCounter == self.ncode:
1113 1113 self.profCounter = 0
1114 1114 self.codeIndex = 0
1115 1115
1116 1116 self.profCounter += 1
1117 1117
1118 1118 def run(self, dataOut, code=None, mode = 0):
1119 1119
1120 1120 if not(self.__isConfig):
1121 1121 if code == None:
1122 1122 code = dataOut.radarControllerHeaderObj.code
1123 1123 # code = dataOut.code
1124 1124
1125 1125 ncode, nbaud = code.shape
1126 1126 self.code = code
1127 1127 self.ncode = ncode
1128 1128 self.nbaud = nbaud
1129 1129 self.__isConfig = True
1130 1130
1131 1131 ndata = dataOut.data.shape[1]
1132 1132
1133 1133 if mode == 0:
1134 1134 self.convolutionInFreq(dataOut.data, ndata)
1135 1135
1136 1136 if mode == 1:
1137 1137 self.convolutionInTime(dataOut.data, ndata)
1138 1138
1139 1139 self.ndata = ndata - self.nbaud + 1
1140 1140
1141 1141 dataOut.data = self.data
1142 1142
1143 1143 dataOut.heightList = dataOut.heightList[:self.ndata]
1144 1144
1145 1145 dataOut.flagNoData = False No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now