##// END OF EJS Templates
schain
RSS
last_update primeras correciones rhi
avaldezp -
r1418:2ccc6450afd4
parent child
Show More
Side by Side Unified
Context file:
r1418:2ccc6450afd4 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -0,0 +1,64
1 import numpy as np
2 import matplotlib.pyplot as pl
3 import warnings
4 #export WRADLIB_DATA=/path/to/wradlib-data
5 warnings.filterwarnings('ignore')
6
7 from wradlib.io import read_generic_netcdf
8 from wradlib.util import get_wradlib_data_file
9 import os
10
11
12
13 # A little helper function for repeated tasks
14 def read_and_overview(filename):
15 """Read NetCDF using read_generic_netcdf and print upper level dictionary keys
16 """
17 test = read_generic_netcdf(filename)
18 print("\nPrint keys for file %s" % os.path.basename(filename))
19 for key in test.keys():
20 print("\t%s" % key)
21 return test
22
23
24 filename = '/home/soporte/Downloads/PX1000/PX-20180220-174014-E0.0-Z.nc'
25 filename = get_wradlib_data_file(filename)
26 test= read_and_overview(filename)
27 print("Height",test['Height'])
28 print("Azimuth",test['Azimuth'])
29 print("Elevation",test['Elevation'])
30 print("CalibH-value",test['CalibH-value'])
31 print("attributes",test['attributes'])
32 print("-------------------------------------------------------------------------------------")
33 for key in test.keys():
34 print(key,test[str(key)])
35
36 '''
37 try:
38 get_ipython().magic('matplotlib inline')
39 except:
40 pl.ion()
41 img, meta = wradlib.io.read_dx(filename)
42 print("Shape of polar array: %r\n" % (img.shape,))
43 print("Some meta data of the DX file:")
44 #print("\tdatetime: %r" % (meta["datetime"],))
45 #print("\tRadar ID: %s" % (meta["radarid"],))
46
47 img[200:250,:]= np.ones([50,img.shape[1]])*np.nan
48
49 img[300:360,:]= np.ones([60,img.shape[1]])*np.nan
50
51 cgax, pm= wradlib.vis.plot_ppi(img)
52 txt = pl.title('Simple PPI')
53 print("coordenada angular",img[:,0],len(img[:,0]))
54 print("COORDENADA 0",img[0],len(img[0]))
55 cbar = pl.gcf().colorbar(pm, pad=0.075)
56
57 #r = np.arange(40, 80)
58 #az = np.arange(200, 250)
59 #ax, pm = wradlib.vis.plot_ppi(img[200:250, 40:80], r, az, autoext=False)
60 #ax, pm = wradlib.vis.plot_ppi(img[200:250, 40:80], r, az)
61
62 #txt = pl.title('Sector PPI')
63 pl.show()
64 '''
Show file before | Show file after | Hide comments
@@ -1,707 +1,898
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
7 7 from schainpy.utils import log
8 8 # libreria wradlib
9 9 import wradlib as wrl
10 10
11 11 EARTH_RADIUS = 6.3710e3
12 12
13 13
14 14 def ll2xy(lat1, lon1, lat2, lon2):
15 15
16 16 p = 0.017453292519943295
17 17 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
18 18 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
19 19 r = 12742 * numpy.arcsin(numpy.sqrt(a))
20 20 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
21 21 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
22 22 theta = -theta + numpy.pi/2
23 23 return r*numpy.cos(theta), r*numpy.sin(theta)
24 24
25 25
26 26 def km2deg(km):
27 27 '''
28 28 Convert distance in km to degrees
29 29 '''
30 30
31 31 return numpy.rad2deg(km/EARTH_RADIUS)
32 32
33 33
34 34
35 35 class SpectralMomentsPlot(SpectraPlot):
36 36 '''
37 37 Plot for Spectral Moments
38 38 '''
39 39 CODE = 'spc_moments'
40 40 # colormap = 'jet'
41 41 # plot_type = 'pcolor'
42 42
43 43 class DobleGaussianPlot(SpectraPlot):
44 44 '''
45 45 Plot for Double Gaussian Plot
46 46 '''
47 47 CODE = 'gaussian_fit'
48 48 # colormap = 'jet'
49 49 # plot_type = 'pcolor'
50 50
51 51 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
52 52 '''
53 53 Plot SpectraCut with Double Gaussian Fit
54 54 '''
55 55 CODE = 'cut_gaussian_fit'
56 56
57 57 class SnrPlot(RTIPlot):
58 58 '''
59 59 Plot for SNR Data
60 60 '''
61 61
62 62 CODE = 'snr'
63 63 colormap = 'jet'
64 64
65 65 def update(self, dataOut):
66 66
67 67 data = {
68 68 'snr': 10*numpy.log10(dataOut.data_snr)
69 69 }
70 70
71 71 return data, {}
72 72
73 73 class DopplerPlot(RTIPlot):
74 74 '''
75 75 Plot for DOPPLER Data (1st moment)
76 76 '''
77 77
78 78 CODE = 'dop'
79 79 colormap = 'jet'
80 80
81 81 def update(self, dataOut):
82 82
83 83 data = {
84 84 'dop': 10*numpy.log10(dataOut.data_dop)
85 85 }
86 86
87 87 return data, {}
88 88
89 89 class PowerPlot(RTIPlot):
90 90 '''
91 91 Plot for Power Data (0 moment)
92 92 '''
93 93
94 94 CODE = 'pow'
95 95 colormap = 'jet'
96 96
97 97 def update(self, dataOut):
98 98 data = {
99 99 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
100 100 }
101 101 return data, {}
102 102
103 103 class SpectralWidthPlot(RTIPlot):
104 104 '''
105 105 Plot for Spectral Width Data (2nd moment)
106 106 '''
107 107
108 108 CODE = 'width'
109 109 colormap = 'jet'
110 110
111 111 def update(self, dataOut):
112 112
113 113 data = {
114 114 'width': dataOut.data_width
115 115 }
116 116
117 117 return data, {}
118 118
119 119 class SkyMapPlot(Plot):
120 120 '''
121 121 Plot for meteors detection data
122 122 '''
123 123
124 124 CODE = 'param'
125 125
126 126 def setup(self):
127 127
128 128 self.ncols = 1
129 129 self.nrows = 1
130 130 self.width = 7.2
131 131 self.height = 7.2
132 132 self.nplots = 1
133 133 self.xlabel = 'Zonal Zenith Angle (deg)'
134 134 self.ylabel = 'Meridional Zenith Angle (deg)'
135 135 self.polar = True
136 136 self.ymin = -180
137 137 self.ymax = 180
138 138 self.colorbar = False
139 139
140 140 def plot(self):
141 141
142 142 arrayParameters = numpy.concatenate(self.data['param'])
143 143 error = arrayParameters[:, -1]
144 144 indValid = numpy.where(error == 0)[0]
145 145 finalMeteor = arrayParameters[indValid, :]
146 146 finalAzimuth = finalMeteor[:, 3]
147 147 finalZenith = finalMeteor[:, 4]
148 148
149 149 x = finalAzimuth * numpy.pi / 180
150 150 y = finalZenith
151 151
152 152 ax = self.axes[0]
153 153
154 154 if ax.firsttime:
155 155 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
156 156 else:
157 157 ax.plot.set_data(x, y)
158 158
159 159 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
160 160 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
161 161 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
162 162 dt2,
163 163 len(x))
164 164 self.titles[0] = title
165 165
166 166
167 167 class GenericRTIPlot(Plot):
168 168 '''
169 169 Plot for data_xxxx object
170 170 '''
171 171
172 172 CODE = 'param'
173 173 colormap = 'viridis'
174 174 plot_type = 'pcolorbuffer'
175 175
176 176 def setup(self):
177 177 self.xaxis = 'time'
178 178 self.ncols = 1
179 179 self.nrows = self.data.shape('param')[0]
180 180 self.nplots = self.nrows
181 181 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
182 182
183 183 if not self.xlabel:
184 184 self.xlabel = 'Time'
185 185
186 186 self.ylabel = 'Range [km]'
187 187 if not self.titles:
188 188 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
189 189
190 190 def update(self, dataOut):
191 191
192 192 data = {
193 193 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
194 194 }
195 195
196 196 meta = {}
197 197
198 198 return data, meta
199 199
200 200 def plot(self):
201 201 # self.data.normalize_heights()
202 202 self.x = self.data.times
203 203 self.y = self.data.yrange
204 204 self.z = self.data['param']
205 205 self.z = 10*numpy.log10(self.z)
206 206 self.z = numpy.ma.masked_invalid(self.z)
207 207
208 208 if self.decimation is None:
209 209 x, y, z = self.fill_gaps(self.x, self.y, self.z)
210 210 else:
211 211 x, y, z = self.fill_gaps(*self.decimate())
212 212
213 213 for n, ax in enumerate(self.axes):
214 214
215 215 self.zmax = self.zmax if self.zmax is not None else numpy.max(
216 216 self.z[n])
217 217 self.zmin = self.zmin if self.zmin is not None else numpy.min(
218 218 self.z[n])
219 219
220 220 if ax.firsttime:
221 221 if self.zlimits is not None:
222 222 self.zmin, self.zmax = self.zlimits[n]
223 223
224 224 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
225 225 vmin=self.zmin,
226 226 vmax=self.zmax,
227 227 cmap=self.cmaps[n]
228 228 )
229 229 else:
230 230 if self.zlimits is not None:
231 231 self.zmin, self.zmax = self.zlimits[n]
232 232 ax.collections.remove(ax.collections[0])
233 233 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
234 234 vmin=self.zmin,
235 235 vmax=self.zmax,
236 236 cmap=self.cmaps[n]
237 237 )
238 238
239 239
240 240 class PolarMapPlot(Plot):
241 241 '''
242 242 Plot for weather radar
243 243 '''
244 244
245 245 CODE = 'param'
246 246 colormap = 'seismic'
247 247
248 248 def setup(self):
249 249 self.ncols = 1
250 250 self.nrows = 1
251 251 self.width = 9
252 252 self.height = 8
253 253 self.mode = self.data.meta['mode']
254 254 if self.channels is not None:
255 255 self.nplots = len(self.channels)
256 256 self.nrows = len(self.channels)
257 257 else:
258 258 self.nplots = self.data.shape(self.CODE)[0]
259 259 self.nrows = self.nplots
260 260 self.channels = list(range(self.nplots))
261 261 if self.mode == 'E':
262 262 self.xlabel = 'Longitude'
263 263 self.ylabel = 'Latitude'
264 264 else:
265 265 self.xlabel = 'Range (km)'
266 266 self.ylabel = 'Height (km)'
267 267 self.bgcolor = 'white'
268 268 self.cb_labels = self.data.meta['units']
269 269 self.lat = self.data.meta['latitude']
270 270 self.lon = self.data.meta['longitude']
271 271 self.xmin, self.xmax = float(
272 272 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
273 273 self.ymin, self.ymax = float(
274 274 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
275 275 # self.polar = True
276 276
277 277 def plot(self):
278 278
279 279 for n, ax in enumerate(self.axes):
280 280 data = self.data['param'][self.channels[n]]
281 281
282 282 zeniths = numpy.linspace(
283 283 0, self.data.meta['max_range'], data.shape[1])
284 284 if self.mode == 'E':
285 285 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
286 286 r, theta = numpy.meshgrid(zeniths, azimuths)
287 287 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
288 288 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
289 289 x = km2deg(x) + self.lon
290 290 y = km2deg(y) + self.lat
291 291 else:
292 292 azimuths = numpy.radians(self.data.yrange)
293 293 r, theta = numpy.meshgrid(zeniths, azimuths)
294 294 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
295 295 self.y = zeniths
296 296
297 297 if ax.firsttime:
298 298 if self.zlimits is not None:
299 299 self.zmin, self.zmax = self.zlimits[n]
300 300 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
301 301 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
302 302 vmin=self.zmin,
303 303 vmax=self.zmax,
304 304 cmap=self.cmaps[n])
305 305 else:
306 306 if self.zlimits is not None:
307 307 self.zmin, self.zmax = self.zlimits[n]
308 308 ax.collections.remove(ax.collections[0])
309 309 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
310 310 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
311 311 vmin=self.zmin,
312 312 vmax=self.zmax,
313 313 cmap=self.cmaps[n])
314 314
315 315 if self.mode == 'A':
316 316 continue
317 317
318 318 # plot district names
319 319 f = open('/data/workspace/schain_scripts/distrito.csv')
320 320 for line in f:
321 321 label, lon, lat = [s.strip() for s in line.split(',') if s]
322 322 lat = float(lat)
323 323 lon = float(lon)
324 324 # ax.plot(lon, lat, '.b', ms=2)
325 325 ax.text(lon, lat, label.decode('utf8'), ha='center',
326 326 va='bottom', size='8', color='black')
327 327
328 328 # plot limites
329 329 limites = []
330 330 tmp = []
331 331 for line in open('/data/workspace/schain_scripts/lima.csv'):
332 332 if '#' in line:
333 333 if tmp:
334 334 limites.append(tmp)
335 335 tmp = []
336 336 continue
337 337 values = line.strip().split(',')
338 338 tmp.append((float(values[0]), float(values[1])))
339 339 for points in limites:
340 340 ax.add_patch(
341 341 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
342 342
343 343 # plot Cuencas
344 344 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
345 345 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
346 346 values = [line.strip().split(',') for line in f]
347 347 points = [(float(s[0]), float(s[1])) for s in values]
348 348 ax.add_patch(Polygon(points, ec='b', fc='none'))
349 349
350 350 # plot grid
351 351 for r in (15, 30, 45, 60):
352 352 ax.add_artist(plt.Circle((self.lon, self.lat),
353 353 km2deg(r), color='0.6', fill=False, lw=0.2))
354 354 ax.text(
355 355 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
356 356 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
357 357 '{}km'.format(r),
358 358 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
359 359
360 360 if self.mode == 'E':
361 361 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
362 362 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
363 363 else:
364 364 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
365 365 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
366 366
367 367 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
368 368 self.titles = ['{} {}'.format(
369 369 self.data.parameters[x], title) for x in self.channels]
370 370
371 371 class WeatherPlot(Plot):
372 372 CODE = 'weather'
373 373 plot_name = 'weather'
374 374 plot_type = 'ppistyle'
375 375 buffering = False
376 376
377 377 def setup(self):
378 378 self.ncols = 1
379 379 self.nrows = 1
380 380 self.nplots= 1
381 381 self.ylabel= 'Range [Km]'
382 382 self.titles= ['Weather']
383 383 self.colorbar=False
384 384 self.width =8
385 385 self.height =8
386 386 self.ini =0
387 387 self.len_azi =0
388 388 self.buffer_ini = None
389 389 self.buffer_azi = None
390 390 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
391 391 self.flag =0
392 392 self.indicador= 0
393 393 self.last_data_azi = None
394 394 self.val_mean = None
395 395
396 396 def update(self, dataOut):
397 397
398 398 data = {}
399 399 meta = {}
400 400 if hasattr(dataOut, 'dataPP_POWER'):
401 401 factor = 1
402 402 if hasattr(dataOut, 'nFFTPoints'):
403 403 factor = dataOut.normFactor
404 404 data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor))
405 405 data['azi'] = dataOut.data_azi
406 406 data['ele'] = dataOut.data_ele
407 407 return data, meta
408 408
409 409 def get2List(self,angulos):
410 410 list1=[]
411 411 list2=[]
412 412 for i in reversed(range(len(angulos))):
413 413 diff_ = angulos[i]-angulos[i-1]
414 414 if diff_ >1.5:
415 415 list1.append(i-1)
416 416 list2.append(diff_)
417 417 return list(reversed(list1)),list(reversed(list2))
418 418
419 419 def fixData360(self,list_,ang_):
420 420 if list_[0]==-1:
421 421 vec = numpy.where(ang_<ang_[0])
422 422 ang_[vec] = ang_[vec]+360
423 423 return ang_
424 424 return ang_
425 425
426 426 def fixData360HL(self,angulos):
427 427 vec = numpy.where(angulos>=360)
428 428 angulos[vec]=angulos[vec]-360
429 429 return angulos
430 430
431 431 def search_pos(self,pos,list_):
432 432 for i in range(len(list_)):
433 433 if pos == list_[i]:
434 434 return True,i
435 435 i=None
436 436 return False,i
437 437
438 438 def fixDataComp(self,ang_,list1_,list2_):
439 439 size = len(ang_)
440 440 size2 = 0
441 441 for i in range(len(list2_)):
442 442 size2=size2+round(list2_[i])-1
443 443 new_size= size+size2
444 444 ang_new = numpy.zeros(new_size)
445 445 ang_new2 = numpy.zeros(new_size)
446 446
447 447 tmp = 0
448 448 c = 0
449 449 for i in range(len(ang_)):
450 450 ang_new[tmp +c] = ang_[i]
451 451 ang_new2[tmp+c] = ang_[i]
452 452 condition , value = self.search_pos(i,list1_)
453 453 if condition:
454 454 pos = tmp + c + 1
455 455 for k in range(round(list2_[value])-1):
456 456 ang_new[pos+k] = ang_new[pos+k-1]+1
457 457 ang_new2[pos+k] = numpy.nan
458 458 tmp = pos +k
459 459 c = 0
460 460 c=c+1
461 461 return ang_new,ang_new2
462 462
463 463 def globalCheckPED(self,angulos):
464 464 l1,l2 = self.get2List(angulos)
465 465 if len(l1)>0:
466 466 angulos2 = self.fixData360(list_=l1,ang_=angulos)
467 467 l1,l2 = self.get2List(angulos2)
468 468
469 469 ang1_,ang2_ = self.fixDataComp(ang_=angulos2,list1_=l1,list2_=l2)
470 470 ang1_ = self.fixData360HL(ang1_)
471 471 ang2_ = self.fixData360HL(ang2_)
472 472 else:
473 473 ang1_= angulos
474 474 ang2_= angulos
475 475 return ang1_,ang2_
476 476
477 477 def analizeDATA(self,data_azi):
478 478 list1 = []
479 479 list2 = []
480 480 dat = data_azi
481 481 for i in reversed(range(1,len(dat))):
482 482 if dat[i]>dat[i-1]:
483 483 diff = int(dat[i])-int(dat[i-1])
484 484 else:
485 485 diff = 360+int(dat[i])-int(dat[i-1])
486 486 if diff > 1:
487 487 list1.append(i-1)
488 488 list2.append(diff-1)
489 489 return list1,list2
490 490
491 491 def fixDATANEW(self,data_azi,data_weather):
492 492 list1,list2 = self.analizeDATA(data_azi)
493 493 if len(list1)== 0:
494 494 return data_azi,data_weather
495 495 else:
496 496 resize = 0
497 497 for i in range(len(list2)):
498 498 resize= resize + list2[i]
499 499 new_data_azi = numpy.resize(data_azi,resize)
500 500 new_data_weather= numpy.resize(date_weather,resize)
501 501
502 502 for i in range(len(list2)):
503 503 j=0
504 504 position=list1[i]+1
505 505 for j in range(list2[i]):
506 506 new_data_azi[position+j]=new_data_azi[position+j-1]+1
507 507 return new_data_azi
508 508
509 509 def fixDATA(self,data_azi):
510 510 data=data_azi
511 511 for i in range(len(data)):
512 512 if numpy.isnan(data[i]):
513 513 data[i]=data[i-1]+1
514 514 return data
515 515
516 516 def replaceNAN(self,data_weather,data_azi,val):
517 517 data= data_azi
518 518 data_T= data_weather
519 519 if data.shape[0]> data_T.shape[0]:
520 520 data_N = numpy.ones( [data.shape[0],data_T.shape[1]])
521 521 c = 0
522 522 for i in range(len(data)):
523 523 if numpy.isnan(data[i]):
524 524 data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
525 525 else:
526 526 data_N[i,:]=data_T[c,:]
527 527 sc=c+1
528 528 else:
529 529 for i in range(len(data)):
530 530 if numpy.isnan(data[i]):
531 531 data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
532 532 return data_T
533 533
534 534 def const_ploteo(self,data_weather,data_azi,step,res):
535 535 if self.ini==0:
536 536 #-------
537 537 n = (360/res)-len(data_azi)
538 538 #--------------------- new -------------------------
539 539 data_azi_new ,data_azi_old= self.globalCheckPED(data_azi)
540 540 #------------------------
541 541 start = data_azi_new[-1] + res
542 542 end = data_azi_new[0] - res
543 543 #------ new
544 544 self.last_data_azi = end
545 545 if start>end:
546 546 end = end + 360
547 547 azi_vacia = numpy.linspace(start,end,int(n))
548 548 azi_vacia = numpy.where(azi_vacia>360,azi_vacia-360,azi_vacia)
549 549 data_azi = numpy.hstack((data_azi_new,azi_vacia))
550 550 # RADAR
551 551 val_mean = numpy.mean(data_weather[:,-1])
552 552 self.val_mean = val_mean
553 553 data_weather_cmp = numpy.ones([(360-data_weather.shape[0]),data_weather.shape[1]])*val_mean
554 554 data_weather = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean)
555 555 data_weather = numpy.vstack((data_weather,data_weather_cmp))
556 556 else:
557 557 # azimuth
558 558 flag=0
559 559 start_azi = self.res_azi[0]
560 560 #-----------new------------
561 561 data_azi ,data_azi_old= self.globalCheckPED(data_azi)
562 562 data_weather = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean)
563 563 #--------------------------
564 564 start = data_azi[0]
565 565 end = data_azi[-1]
566 566 self.last_data_azi= end
567 567 if start< start_azi:
568 568 start = start +360
569 569 if end <start_azi:
570 570 end = end +360
571 571
572 572 pos_ini = int((start-start_azi)/res)
573 573 len_azi = len(data_azi)
574 574 if (360-pos_ini)<len_azi:
575 575 if pos_ini+1==360:
576 576 pos_ini=0
577 577 else:
578 578 flag=1
579 579 dif= 360-pos_ini
580 580 comp= len_azi-dif
581 581 #-----------------
582 582 if flag==0:
583 583 # AZIMUTH
584 584 self.res_azi[pos_ini:pos_ini+len_azi] = data_azi
585 585 # RADAR
586 586 self.res_weather[pos_ini:pos_ini+len_azi,:] = data_weather
587 587 else:
588 588 # AZIMUTH
589 589 self.res_azi[pos_ini:pos_ini+dif] = data_azi[0:dif]
590 590 self.res_azi[0:comp] = data_azi[dif:]
591 591 # RADAR
592 592 self.res_weather[pos_ini:pos_ini+dif,:] = data_weather[0:dif,:]
593 593 self.res_weather[0:comp,:] = data_weather[dif:,:]
594 594 flag=0
595 595 data_azi = self.res_azi
596 596 data_weather = self.res_weather
597 597
598 598 return data_weather,data_azi
599 599
600 600 def plot(self):
601 601 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
602 602 data = self.data[-1]
603 603 r = self.data.yrange
604 604 delta_height = r[1]-r[0]
605 605 r_mask = numpy.where(r>=0)[0]
606 606 r = numpy.arange(len(r_mask))*delta_height
607 607 self.y = 2*r
608 608 # RADAR
609 609 #data_weather = data['weather']
610 610 # PEDESTAL
611 611 #data_azi = data['azi']
612 612 res = 1
613 613 # STEP
614 614 step = (360/(res*data['weather'].shape[0]))
615 615
616 616 self.res_weather, self.res_azi = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_azi=data['azi'],step=step,res=res)
617 617 self.res_ele = numpy.mean(data['ele'])
618 618 ################# PLOTEO ###################
619 619
620 620 for i,ax in enumerate(self.axes):
621 621 if ax.firsttime:
622 622 plt.clf()
623 623 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=8, vmax=35)
624 624 else:
625 625 plt.clf()
626 626 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=8, vmax=35)
627 627 caax = cgax.parasites[0]
628 628 paax = cgax.parasites[1]
629 629 cbar = plt.gcf().colorbar(pm, pad=0.075)
630 630 caax.set_xlabel('x_range [km]')
631 631 caax.set_ylabel('y_range [km]')
632 632 plt.text(1.0, 1.05, 'Azimuth '+str(thisDatetime)+" Step "+str(self.ini)+ " Elev: "+str(round(self.res_ele,2)), transform=caax.transAxes, va='bottom',ha='right')
633 633
634 634 self.ini= self.ini+1
635 635
636 636
637 637 class WeatherRHIPlot(Plot):
638 638 CODE = 'weather'
639 639 plot_name = 'weather'
640 640 plot_type = 'rhistyle'
641 641 buffering = False
642 642
643 643 def setup(self):
644 644 self.ncols = 1
645 645 self.nrows = 1
646 646 self.nplots= 1
647 647 self.ylabel= 'Range [Km]'
648 648 self.titles= ['Weather']
649 649 self.colorbar=False
650 650 self.width =8
651 651 self.height =8
652 652 self.ini =0
653 653 self.len_azi =0
654 654 self.buffer_ini = None
655 self.buffer_azi = None
655 self.buffer_ele = None
656 656 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
657 657 self.flag =0
658 658 self.indicador= 0
659 self.last_data_azi = None
659 self.last_data_ele = None
660 660 self.val_mean = None
661 661
662 662 def update(self, dataOut):
663 663
664 664 data = {}
665 665 meta = {}
666 666 if hasattr(dataOut, 'dataPP_POWER'):
667 667 factor = 1
668 668 if hasattr(dataOut, 'nFFTPoints'):
669 669 factor = dataOut.normFactor
670 670 data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor))
671 671 data['azi'] = dataOut.data_azi
672 672 data['ele'] = dataOut.data_ele
673 673 return data, meta
674 674
675 def get2List(self,angulos):
676 list1=[]
677 list2=[]
678 for i in reversed(range(len(angulos))):
679 diff_ = angulos[i]-angulos[i-1]
680 if diff_ >1.5:
681 list1.append(i-1)
682 list2.append(diff_)
683 return list(reversed(list1)),list(reversed(list2))
684
685 def fixData180(self,list_,ang_):
686 if list_[0]==-1:
687 vec = numpy.where(ang_<ang_[0])
688 ang_[vec] = ang_[vec]+180
689 return ang_
690 return ang_
691
692 def fixData180HL(self,angulos):
693 vec = numpy.where(angulos>=180)
694 angulos[vec]=angulos[vec]-180
695 return angulos
696
697
698 def search_pos(self,pos,list_):
699 for i in range(len(list_)):
700 if pos == list_[i]:
701 return True,i
702 i=None
703 return False,i
704
705 def fixDataComp(self,ang_,list1_,list2_):
706 size = len(ang_)
707 size2 = 0
708 for i in range(len(list2_)):
709 size2=size2+round(list2_[i])-1
710 new_size= size+size2
711 ang_new = numpy.zeros(new_size)
712 ang_new2 = numpy.zeros(new_size)
713
714 tmp = 0
715 c = 0
716 for i in range(len(ang_)):
717 ang_new[tmp +c] = ang_[i]
718 ang_new2[tmp+c] = ang_[i]
719 condition , value = self.search_pos(i,list1_)
720 if condition:
721 pos = tmp + c + 1
722 for k in range(round(list2_[value])-1):
723 ang_new[pos+k] = ang_new[pos+k-1]+1
724 ang_new2[pos+k] = numpy.nan
725 tmp = pos +k
726 c = 0
727 c=c+1
728 return ang_new,ang_new2
729
730 def globalCheckPED(self,angulos):
731 l1,l2 = self.get2List(angulos)
732 if len(l1)>0:
733 angulos2 = self.fixData180(list_=l1,ang_=angulos)
734 l1,l2 = self.get2List(angulos2)
735
736 ang1_,ang2_ = self.fixDataComp(ang_=angulos2,list1_=l1,list2_=l2)
737 ang1_ = self.fixData180HL(ang1_)
738 ang2_ = self.fixData180HL(ang2_)
739 else:
740 ang1_= angulos
741 ang2_= angulos
742 return ang1_,ang2_
743
744
745 def replaceNAN(self,data_weather,data_ele,val):
746 data= data_ele
747 data_T= data_weather
748 if data.shape[0]> data_T.shape[0]:
749 data_N = numpy.ones( [data.shape[0],data_T.shape[1]])
750 c = 0
751 for i in range(len(data)):
752 if numpy.isnan(data[i]):
753 data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
754 else:
755 data_N[i,:]=data_T[c,:]
756 sc=c+1
757 else:
758 for i in range(len(data)):
759 if numpy.isnan(data[i]):
760 data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
761 return data_T
762
763 def const_ploteo(self,data_weather,data_ele,step,res):
764 if self.ini==0:
765 #-------
766 n = (180/res)-len(data_ele)
767 #--------------------- new -------------------------
768 data_ele_new ,data_ele_old= self.globalCheckPED(data_ele)
769 #------------------------
770 start = data_ele_new[-1] + res
771 end = data_ele_new[0] - res
772 #------ new
773 self.last_data_ele = end
774 if start>end:
775 end = end + 180
776 ele_vacia = numpy.linspace(start,end,int(n))
777 ele_vacia = numpy.where(ele_vacia>180,ele_vacia-180,ele_vacia)
778 data_ele = numpy.hstack((data_ele_new,ele_vacia))
779 # RADAR
780 val_mean = numpy.mean(data_weather[:,-1])
781 self.val_mean = val_mean
782 data_weather_cmp = numpy.ones([(180-data_weather.shape[0]),data_weather.shape[1]])*val_mean
783 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
784 data_weather = numpy.vstack((data_weather,data_weather_cmp))
785 else:
786 # azimuth
787 flag=0
788 start_ele = self.res_ele[0]
789 #-----------new------------
790 data_ele ,data_ele_old= self.globalCheckPED(data_ele)
791 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
792 #--------------------------
793 start = data_ele[0]
794 end = data_ele[-1]
795 self.last_data_ele= end
796 if start< start_ele:
797 start = start +180
798 if end <start_ele:
799 end = end +180
800
801 pos_ini = int((start-start_ele)/res)
802 len_ele = len(data_ele)
803 if (180-pos_ini)<len_ele:
804 if pos_ini+1==180:
805 pos_ini=0
806 else:
807 flag=1
808 dif= 180-pos_ini
809 comp= len_ele-dif
810 #-----------------
811 if flag==0:
812 # elevacion
813 self.res_ele[pos_ini:pos_ini+len_ele] = data_ele
814 # RADAR
815 self.res_weather[pos_ini:pos_ini+len_ele,:] = data_weather
816 else:
817 # elevacion
818 self.res_ele[pos_ini:pos_ini+dif] = data_ele[0:dif]
819 self.res_ele[0:comp] = data_ele[dif:]
820 # RADAR
821 self.res_weather[pos_ini:pos_ini+dif,:] = data_weather[0:dif,:]
822 self.res_weather[0:comp,:] = data_weather[dif:,:]
823 flag=0
824 data_ele = self.res_ele
825 data_weather = self.res_weather
826
827 return data_weather,data_ele
828
829
675 830 def plot(self):
676 831 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
677 832 data = self.data[-1]
678 833 r = self.data.yrange
679 834 delta_height = r[1]-r[0]
680 835 r_mask = numpy.where(r>=0)[0]
681 836 r = numpy.arange(len(r_mask))*delta_height
682 837 self.y = 2*r
838 '''
839 #-------------------------------------------------------------
840 # RADAR
841 #data_weather = data['weather']
842 # PEDESTAL
843 #data_azi = data['azi']
844 res = 1
845 # STEP
846 step = (180/(res*data['weather'].shape[0]))
847
848
849 self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_ele=data['ele'],step=step,res=res)
850 self.res_azi = numpy.mean(data['azi'])
851 print("self.res_ele------------------------------:",self.res_ele)
852 ################# PLOTEO ###################
853
854 for i,ax in enumerate(self.axes):
855 if ax.firsttime:
856 plt.clf()
857 cgax, pm = wrl.vis.plot_rhi(self.res_weather,r=r,th=self.res_ele,fig=self.figures[0], proj='cg', vmin=8, vmax=35)
858 else:
859 plt.clf()
860 cgax, pm = wrl.vis.plot_rhi(self.res_weather,r=r,th=self.res_ele,fig=self.figures[0], proj='cg', vmin=8, vmax=35)
861 caax = cgax.parasites[0]
862 paax = cgax.parasites[1]
863 cbar = plt.gcf().colorbar(pm, pad=0.075)
864 caax.set_xlabel('x_range [km]')
865 caax.set_ylabel('y_range [km]')
866 plt.text(1.0, 1.05, 'Elevacion '+str(thisDatetime)+" Step "+str(self.ini)+ " Azi: "+str(round(self.res_azi,2)), transform=caax.transAxes, va='bottom',ha='right')
867
868 self.ini= self.ini+1
869
870
871 '''
872 #--------------------------------------------------------
873
683 874 ###self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_azi=data['ele'],step=step,res=res)
684 875 ###self.res_azi = numpy.mean(data['azi'])
685 876 #-------------
686 877 # 90 angulos en el axis 0
687 878 # 1000 step en el axis 1
688 self.res_weather = numpy.ones([90,1000])
879 self.res_weather = numpy.ones([120,1000])
689 880 r = numpy.linspace(0,1999,1000)
690 self.res_ele = numpy.arange(0,90)
881 self.res_ele = numpy.arange(0,120)
691 882 self.res_azi = 240
692 883 #-------------
693 884 for i,ax in enumerate(self.axes):
694 885 if ax.firsttime:
695 886 plt.clf()
696 887 cgax, pm = wrl.vis.plot_rhi(self.res_weather,r=r,th=self.res_ele,fig=self.figures[0], proj='cg')
697 888 else:
698 889 plt.clf()
699 890 cgax, pm = wrl.vis.plot_rhi(self.res_weather,r=r,th=self.res_ele,fig=self.figures[0], proj='cg')
700 891 caax = cgax.parasites[0]
701 892 paax = cgax.parasites[1]
702 893 cbar = plt.gcf().colorbar(pm, pad=0.075)
703 894 caax.set_xlabel('x_range [km]')
704 895 caax.set_ylabel('y_range [km]')
705 896 plt.text(1.0, 1.05, 'Elevacion '+str(thisDatetime)+" Step "+str(self.ini)+ " Azi: "+str(round(self.res_azi,2)), transform=caax.transAxes, va='bottom',ha='right')
706 897
707 898 self.ini= self.ini+1
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
Show file before | Show file after | Hide comments
@@ -1,1633 +1,1633
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62 62 for channel in channelList:
63 63 if channel not in self.dataOut.channelList:
64 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
65 65
66 66 index = self.dataOut.channelList.index(channel)
67 67 channelIndexList.append(index)
68 68 self.selectChannelsByIndex(channelIndexList)
69 69 return self.dataOut
70 70
71 71 def selectChannelsByIndex(self, channelIndexList):
72 72 """
73 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
74 74
75 75 Input:
76 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
77 77
78 78 Affected:
79 79 self.dataOut.data
80 80 self.dataOut.channelIndexList
81 81 self.dataOut.nChannels
82 82 self.dataOut.m_ProcessingHeader.totalSpectra
83 83 self.dataOut.systemHeaderObj.numChannels
84 84 self.dataOut.m_ProcessingHeader.blockSize
85 85
86 86 Return:
87 87 None
88 88 """
89 89
90 90 for channelIndex in channelIndexList:
91 91 if channelIndex not in self.dataOut.channelIndexList:
92 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
93 93
94 94 if self.dataOut.type == 'Voltage':
95 95 if self.dataOut.flagDataAsBlock:
96 96 """
97 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
98 98 """
99 99 data = self.dataOut.data[channelIndexList,:,:]
100 100 else:
101 101 data = self.dataOut.data[channelIndexList,:]
102 102
103 103 self.dataOut.data = data
104 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
105 105 self.dataOut.channelList = range(len(channelIndexList))
106 106
107 107 elif self.dataOut.type == 'Spectra':
108 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
109 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
110 110
111 111 self.dataOut.data_spc = data_spc
112 112 self.dataOut.data_dc = data_dc
113 113
114 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 self.dataOut.channelList = range(len(channelIndexList))
116 116 self.__selectPairsByChannel(channelIndexList)
117 117
118 118 return 1
119 119
120 120 def __selectPairsByChannel(self, channelList=None):
121 121
122 122 if channelList == None:
123 123 return
124 124
125 125 pairsIndexListSelected = []
126 126 for pairIndex in self.dataOut.pairsIndexList:
127 127 # First pair
128 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
129 129 continue
130 130 # Second pair
131 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
132 132 continue
133 133
134 134 pairsIndexListSelected.append(pairIndex)
135 135
136 136 if not pairsIndexListSelected:
137 137 self.dataOut.data_cspc = None
138 138 self.dataOut.pairsList = []
139 139 return
140 140
141 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
142 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
143 143 for i in pairsIndexListSelected]
144 144
145 145 return
146 146
147 147 class selectHeights(Operation):
148 148
149 149 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
150 150 """
151 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
152 152 minHei <= height <= maxHei
153 153
154 154 Input:
155 155 minHei : valor minimo de altura a considerar
156 156 maxHei : valor maximo de altura a considerar
157 157
158 158 Affected:
159 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
160 160
161 161 Return:
162 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
163 163 """
164 164
165 165 self.dataOut = dataOut
166 166
167 167 if minHei and maxHei:
168 168
169 169 if (minHei < self.dataOut.heightList[0]):
170 170 minHei = self.dataOut.heightList[0]
171 171
172 172 if (maxHei > self.dataOut.heightList[-1]):
173 173 maxHei = self.dataOut.heightList[-1]
174 174
175 175 minIndex = 0
176 176 maxIndex = 0
177 177 heights = self.dataOut.heightList
178 178
179 179 inda = numpy.where(heights >= minHei)
180 180 indb = numpy.where(heights <= maxHei)
181 181
182 182 try:
183 183 minIndex = inda[0][0]
184 184 except:
185 185 minIndex = 0
186 186
187 187 try:
188 188 maxIndex = indb[0][-1]
189 189 except:
190 190 maxIndex = len(heights)
191 191
192 192 self.selectHeightsByIndex(minIndex, maxIndex)
193 193
194 194 return self.dataOut
195 195
196 196 def selectHeightsByIndex(self, minIndex, maxIndex):
197 197 """
198 198 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
199 199 minIndex <= index <= maxIndex
200 200
201 201 Input:
202 202 minIndex : valor de indice minimo de altura a considerar
203 203 maxIndex : valor de indice maximo de altura a considerar
204 204
205 205 Affected:
206 206 self.dataOut.data
207 207 self.dataOut.heightList
208 208
209 209 Return:
210 210 1 si el metodo se ejecuto con exito caso contrario devuelve 0
211 211 """
212 212
213 213 if self.dataOut.type == 'Voltage':
214 214 if (minIndex < 0) or (minIndex > maxIndex):
215 215 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
216 216
217 217 if (maxIndex >= self.dataOut.nHeights):
218 218 maxIndex = self.dataOut.nHeights
219 219
220 220 #voltage
221 221 if self.dataOut.flagDataAsBlock:
222 222 """
223 223 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
224 224 """
225 225 data = self.dataOut.data[:,:, minIndex:maxIndex]
226 226 else:
227 227 data = self.dataOut.data[:, minIndex:maxIndex]
228 228
229 229 # firstHeight = self.dataOut.heightList[minIndex]
230 230
231 231 self.dataOut.data = data
232 232 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
233 233
234 234 if self.dataOut.nHeights <= 1:
235 235 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
236 236 elif self.dataOut.type == 'Spectra':
237 237 if (minIndex < 0) or (minIndex > maxIndex):
238 238 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
239 239 minIndex, maxIndex))
240 240
241 241 if (maxIndex >= self.dataOut.nHeights):
242 242 maxIndex = self.dataOut.nHeights - 1
243 243
244 244 # Spectra
245 245 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
246 246
247 247 data_cspc = None
248 248 if self.dataOut.data_cspc is not None:
249 249 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
250 250
251 251 data_dc = None
252 252 if self.dataOut.data_dc is not None:
253 253 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
254 254
255 255 self.dataOut.data_spc = data_spc
256 256 self.dataOut.data_cspc = data_cspc
257 257 self.dataOut.data_dc = data_dc
258 258
259 259 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
260 260
261 261 return 1
262 262
263 263
264 264 class filterByHeights(Operation):
265 265
266 266 def run(self, dataOut, window):
267 267
268 268 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
269 269
270 270 if window == None:
271 271 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
272 272
273 273 newdelta = deltaHeight * window
274 274 r = dataOut.nHeights % window
275 275 newheights = (dataOut.nHeights-r)/window
276 276
277 277 if newheights <= 1:
278 278 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
279 279
280 280 if dataOut.flagDataAsBlock:
281 281 """
282 282 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
283 283 """
284 284 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
285 285 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
286 286 buffer = numpy.sum(buffer,3)
287 287
288 288 else:
289 289 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
290 290 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
291 291 buffer = numpy.sum(buffer,2)
292 292
293 293 dataOut.data = buffer
294 294 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
295 295 dataOut.windowOfFilter = window
296 296
297 297 return dataOut
298 298
299 299
300 300 class setH0(Operation):
301 301
302 302 def run(self, dataOut, h0, deltaHeight = None):
303 303
304 304 if not deltaHeight:
305 305 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
306 306
307 307 nHeights = dataOut.nHeights
308 308
309 309 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
310 310
311 311 dataOut.heightList = newHeiRange
312 312
313 313 return dataOut
314 314
315 315
316 316 class deFlip(Operation):
317 317
318 318 def run(self, dataOut, channelList = []):
319 319
320 320 data = dataOut.data.copy()
321 321
322 322 if dataOut.flagDataAsBlock:
323 323 flip = self.flip
324 324 profileList = list(range(dataOut.nProfiles))
325 325
326 326 if not channelList:
327 327 for thisProfile in profileList:
328 328 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
329 329 flip *= -1.0
330 330 else:
331 331 for thisChannel in channelList:
332 332 if thisChannel not in dataOut.channelList:
333 333 continue
334 334
335 335 for thisProfile in profileList:
336 336 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
337 337 flip *= -1.0
338 338
339 339 self.flip = flip
340 340
341 341 else:
342 342 if not channelList:
343 343 data[:,:] = data[:,:]*self.flip
344 344 else:
345 345 for thisChannel in channelList:
346 346 if thisChannel not in dataOut.channelList:
347 347 continue
348 348
349 349 data[thisChannel,:] = data[thisChannel,:]*self.flip
350 350
351 351 self.flip *= -1.
352 352
353 353 dataOut.data = data
354 354
355 355 return dataOut
356 356
357 357
358 358 class setAttribute(Operation):
359 359 '''
360 360 Set an arbitrary attribute(s) to dataOut
361 361 '''
362 362
363 363 def __init__(self):
364 364
365 365 Operation.__init__(self)
366 366 self._ready = False
367 367
368 368 def run(self, dataOut, **kwargs):
369 369
370 370 for key, value in kwargs.items():
371 371 setattr(dataOut, key, value)
372 372
373 373 return dataOut
374 374
375 375
376 376 @MPDecorator
377 377 class printAttribute(Operation):
378 378 '''
379 379 Print an arbitrary attribute of dataOut
380 380 '''
381 381
382 382 def __init__(self):
383 383
384 384 Operation.__init__(self)
385 385
386 386 def run(self, dataOut, attributes):
387 387
388 388 if isinstance(attributes, str):
389 389 attributes = [attributes]
390 390 for attr in attributes:
391 391 if hasattr(dataOut, attr):
392 392 log.log(getattr(dataOut, attr), attr)
393 393
394 394
395 395 class interpolateHeights(Operation):
396 396
397 397 def run(self, dataOut, topLim, botLim):
398 398 #69 al 72 para julia
399 399 #82-84 para meteoros
400 400 if len(numpy.shape(dataOut.data))==2:
401 401 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
402 402 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
403 403 #dataOut.data[:,botLim:limSup+1] = sampInterp
404 404 dataOut.data[:,botLim:topLim+1] = sampInterp
405 405 else:
406 406 nHeights = dataOut.data.shape[2]
407 407 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
408 408 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
409 409 f = interpolate.interp1d(x, y, axis = 2)
410 410 xnew = numpy.arange(botLim,topLim+1)
411 411 ynew = f(xnew)
412 412 dataOut.data[:,:,botLim:topLim+1] = ynew
413 413
414 414 return dataOut
415 415
416 416
417 417 class CohInt(Operation):
418 418
419 419 isConfig = False
420 420 __profIndex = 0
421 421 __byTime = False
422 422 __initime = None
423 423 __lastdatatime = None
424 424 __integrationtime = None
425 425 __buffer = None
426 426 __bufferStride = []
427 427 __dataReady = False
428 428 __profIndexStride = 0
429 429 __dataToPutStride = False
430 430 n = None
431 431
432 432 def __init__(self, **kwargs):
433 433
434 434 Operation.__init__(self, **kwargs)
435 435
436 436 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
437 437 """
438 438 Set the parameters of the integration class.
439 439
440 440 Inputs:
441 441
442 442 n : Number of coherent integrations
443 443 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
444 444 overlapping :
445 445 """
446 446
447 447 self.__initime = None
448 448 self.__lastdatatime = 0
449 449 self.__buffer = None
450 450 self.__dataReady = False
451 451 self.byblock = byblock
452 452 self.stride = stride
453 453
454 454 if n == None and timeInterval == None:
455 455 raise ValueError("n or timeInterval should be specified ...")
456 456
457 457 if n != None:
458 458 self.n = n
459 459 self.__byTime = False
460 460 else:
461 461 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
462 462 self.n = 9999
463 463 self.__byTime = True
464 464
465 465 if overlapping:
466 466 self.__withOverlapping = True
467 467 self.__buffer = None
468 468 else:
469 469 self.__withOverlapping = False
470 470 self.__buffer = 0
471 471
472 472 self.__profIndex = 0
473 473
474 474 def putData(self, data):
475 475
476 476 """
477 477 Add a profile to the __buffer and increase in one the __profileIndex
478 478
479 479 """
480 480
481 481 if not self.__withOverlapping:
482 482 self.__buffer += data.copy()
483 483 self.__profIndex += 1
484 484 return
485 485
486 486 #Overlapping data
487 487 nChannels, nHeis = data.shape
488 488 data = numpy.reshape(data, (1, nChannels, nHeis))
489 489
490 490 #If the buffer is empty then it takes the data value
491 491 if self.__buffer is None:
492 492 self.__buffer = data
493 493 self.__profIndex += 1
494 494 return
495 495
496 496 #If the buffer length is lower than n then stakcing the data value
497 497 if self.__profIndex < self.n:
498 498 self.__buffer = numpy.vstack((self.__buffer, data))
499 499 self.__profIndex += 1
500 500 return
501 501
502 502 #If the buffer length is equal to n then replacing the last buffer value with the data value
503 503 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
504 504 self.__buffer[self.n-1] = data
505 505 self.__profIndex = self.n
506 506 return
507 507
508 508
509 509 def pushData(self):
510 510 """
511 511 Return the sum of the last profiles and the profiles used in the sum.
512 512
513 513 Affected:
514 514
515 515 self.__profileIndex
516 516
517 517 """
518 518
519 519 if not self.__withOverlapping:
520 520 data = self.__buffer
521 521 n = self.__profIndex
522 522
523 523 self.__buffer = 0
524 524 self.__profIndex = 0
525 525
526 526 return data, n
527 527
528 528 #Integration with Overlapping
529 529 data = numpy.sum(self.__buffer, axis=0)
530 530 # print data
531 531 # raise
532 532 n = self.__profIndex
533 533
534 534 return data, n
535 535
536 536 def byProfiles(self, data):
537 537
538 538 self.__dataReady = False
539 539 avgdata = None
540 540 # n = None
541 541 # print data
542 542 # raise
543 543 self.putData(data)
544 544
545 545 if self.__profIndex == self.n:
546 546 avgdata, n = self.pushData()
547 547 self.__dataReady = True
548 548
549 549 return avgdata
550 550
551 551 def byTime(self, data, datatime):
552 552
553 553 self.__dataReady = False
554 554 avgdata = None
555 555 n = None
556 556
557 557 self.putData(data)
558 558
559 559 if (datatime - self.__initime) >= self.__integrationtime:
560 560 avgdata, n = self.pushData()
561 561 self.n = n
562 562 self.__dataReady = True
563 563
564 564 return avgdata
565 565
566 566 def integrateByStride(self, data, datatime):
567 567 # print data
568 568 if self.__profIndex == 0:
569 569 self.__buffer = [[data.copy(), datatime]]
570 570 else:
571 571 self.__buffer.append([data.copy(),datatime])
572 572 self.__profIndex += 1
573 573 self.__dataReady = False
574 574
575 575 if self.__profIndex == self.n * self.stride :
576 576 self.__dataToPutStride = True
577 577 self.__profIndexStride = 0
578 578 self.__profIndex = 0
579 579 self.__bufferStride = []
580 580 for i in range(self.stride):
581 581 current = self.__buffer[i::self.stride]
582 582 data = numpy.sum([t[0] for t in current], axis=0)
583 583 avgdatatime = numpy.average([t[1] for t in current])
584 584 # print data
585 585 self.__bufferStride.append((data, avgdatatime))
586 586
587 587 if self.__dataToPutStride:
588 588 self.__dataReady = True
589 589 self.__profIndexStride += 1
590 590 if self.__profIndexStride == self.stride:
591 591 self.__dataToPutStride = False
592 592 # print self.__bufferStride[self.__profIndexStride - 1]
593 593 # raise
594 594 return self.__bufferStride[self.__profIndexStride - 1]
595 595
596 596
597 597 return None, None
598 598
599 599 def integrate(self, data, datatime=None):
600 600
601 601 if self.__initime == None:
602 602 self.__initime = datatime
603 603
604 604 if self.__byTime:
605 605 avgdata = self.byTime(data, datatime)
606 606 else:
607 607 avgdata = self.byProfiles(data)
608 608
609 609
610 610 self.__lastdatatime = datatime
611 611
612 612 if avgdata is None:
613 613 return None, None
614 614
615 615 avgdatatime = self.__initime
616 616
617 617 deltatime = datatime - self.__lastdatatime
618 618
619 619 if not self.__withOverlapping:
620 620 self.__initime = datatime
621 621 else:
622 622 self.__initime += deltatime
623 623
624 624 return avgdata, avgdatatime
625 625
626 626 def integrateByBlock(self, dataOut):
627 627
628 628 times = int(dataOut.data.shape[1]/self.n)
629 629 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
630 630
631 631 id_min = 0
632 632 id_max = self.n
633 633
634 634 for i in range(times):
635 635 junk = dataOut.data[:,id_min:id_max,:]
636 636 avgdata[:,i,:] = junk.sum(axis=1)
637 637 id_min += self.n
638 638 id_max += self.n
639 639
640 640 timeInterval = dataOut.ippSeconds*self.n
641 641 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
642 642 self.__dataReady = True
643 643 return avgdata, avgdatatime
644 644
645 645 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
646 646
647 647 if not self.isConfig:
648 648 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
649 649 self.isConfig = True
650 650
651 651 if dataOut.flagDataAsBlock:
652 652 """
653 653 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
654 654 """
655 655 avgdata, avgdatatime = self.integrateByBlock(dataOut)
656 656 dataOut.nProfiles /= self.n
657 657 else:
658 658 if stride is None:
659 659 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
660 660 else:
661 661 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
662 662
663 663
664 664 # dataOut.timeInterval *= n
665 665 dataOut.flagNoData = True
666 666
667 667 if self.__dataReady:
668 668 dataOut.data = avgdata
669 669 if not dataOut.flagCohInt:
670 670 dataOut.nCohInt *= self.n
671 671 dataOut.flagCohInt = True
672 672 dataOut.utctime = avgdatatime
673 673 # print avgdata, avgdatatime
674 674 # raise
675 675 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
676 676 dataOut.flagNoData = False
677 677 return dataOut
678 678
679 679 class Decoder(Operation):
680 680
681 681 isConfig = False
682 682 __profIndex = 0
683 683
684 684 code = None
685 685
686 686 nCode = None
687 687 nBaud = None
688 688
689 689 def __init__(self, **kwargs):
690 690
691 691 Operation.__init__(self, **kwargs)
692 692
693 693 self.times = None
694 694 self.osamp = None
695 695 # self.__setValues = False
696 696 self.isConfig = False
697 697 self.setupReq = False
698 698 def setup(self, code, osamp, dataOut):
699 699
700 700 self.__profIndex = 0
701 701
702 702 self.code = code
703 703
704 704 self.nCode = len(code)
705 705 self.nBaud = len(code[0])
706 706
707 707 if (osamp != None) and (osamp >1):
708 708 self.osamp = osamp
709 709 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
710 710 self.nBaud = self.nBaud*self.osamp
711 711
712 712 self.__nChannels = dataOut.nChannels
713 713 self.__nProfiles = dataOut.nProfiles
714 714 self.__nHeis = dataOut.nHeights
715 715
716 716 if self.__nHeis < self.nBaud:
717 717 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
718 718
719 719 #Frequency
720 720 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
721 721
722 722 __codeBuffer[:,0:self.nBaud] = self.code
723 723
724 724 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
725 725
726 726 if dataOut.flagDataAsBlock:
727 727
728 728 self.ndatadec = self.__nHeis #- self.nBaud + 1
729 729
730 730 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
731 731
732 732 else:
733 733
734 734 #Time
735 735 self.ndatadec = self.__nHeis #- self.nBaud + 1
736 736
737 737 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
738 738
739 739 def __convolutionInFreq(self, data):
740 740
741 741 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
742 742
743 743 fft_data = numpy.fft.fft(data, axis=1)
744 744
745 745 conv = fft_data*fft_code
746 746
747 747 data = numpy.fft.ifft(conv,axis=1)
748 748
749 749 return data
750 750
751 751 def __convolutionInFreqOpt(self, data):
752 752
753 753 raise NotImplementedError
754 754
755 755 def __convolutionInTime(self, data):
756 756
757 757 code = self.code[self.__profIndex]
758 758 for i in range(self.__nChannels):
759 759 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
760 760
761 761 return self.datadecTime
762 762
763 763 def __convolutionByBlockInTime(self, data):
764 764
765 765 repetitions = int(self.__nProfiles / self.nCode)
766 766 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
767 767 junk = junk.flatten()
768 768 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
769 769 profilesList = range(self.__nProfiles)
770 770
771 771 for i in range(self.__nChannels):
772 772 for j in profilesList:
773 773 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
774 774 return self.datadecTime
775 775
776 776 def __convolutionByBlockInFreq(self, data):
777 777
778 778 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
779 779
780 780
781 781 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
782 782
783 783 fft_data = numpy.fft.fft(data, axis=2)
784 784
785 785 conv = fft_data*fft_code
786 786
787 787 data = numpy.fft.ifft(conv,axis=2)
788 788
789 789 return data
790 790
791 791
792 792 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
793 793
794 794 if dataOut.flagDecodeData:
795 795 print("This data is already decoded, recoding again ...")
796 796
797 797 if not self.isConfig:
798 798
799 799 if code is None:
800 800 if dataOut.code is None:
801 801 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
802 802
803 803 code = dataOut.code
804 804 else:
805 805 code = numpy.array(code).reshape(nCode,nBaud)
806 806 self.setup(code, osamp, dataOut)
807 807
808 808 self.isConfig = True
809 809
810 810 if mode == 3:
811 811 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
812 812
813 813 if times != None:
814 814 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
815 815
816 816 if self.code is None:
817 817 print("Fail decoding: Code is not defined.")
818 818 return
819 819
820 820 self.__nProfiles = dataOut.nProfiles
821 821 datadec = None
822 822
823 823 if mode == 3:
824 824 mode = 0
825 825
826 826 if dataOut.flagDataAsBlock:
827 827 """
828 828 Decoding when data have been read as block,
829 829 """
830 830
831 831 if mode == 0:
832 832 datadec = self.__convolutionByBlockInTime(dataOut.data)
833 833 if mode == 1:
834 834 datadec = self.__convolutionByBlockInFreq(dataOut.data)
835 835 else:
836 836 """
837 837 Decoding when data have been read profile by profile
838 838 """
839 839 if mode == 0:
840 840 datadec = self.__convolutionInTime(dataOut.data)
841 841
842 842 if mode == 1:
843 843 datadec = self.__convolutionInFreq(dataOut.data)
844 844
845 845 if mode == 2:
846 846 datadec = self.__convolutionInFreqOpt(dataOut.data)
847 847
848 848 if datadec is None:
849 849 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
850 850
851 851 dataOut.code = self.code
852 852 dataOut.nCode = self.nCode
853 853 dataOut.nBaud = self.nBaud
854 854
855 855 dataOut.data = datadec
856 856
857 857 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
858 858
859 859 dataOut.flagDecodeData = True #asumo q la data esta decodificada
860 860
861 861 if self.__profIndex == self.nCode-1:
862 862 self.__profIndex = 0
863 863 return dataOut
864 864
865 865 self.__profIndex += 1
866 866
867 867 return dataOut
868 868 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
869 869
870 870
871 871 class ProfileConcat(Operation):
872 872
873 873 isConfig = False
874 874 buffer = None
875 875
876 876 def __init__(self, **kwargs):
877 877
878 878 Operation.__init__(self, **kwargs)
879 879 self.profileIndex = 0
880 880
881 881 def reset(self):
882 882 self.buffer = numpy.zeros_like(self.buffer)
883 883 self.start_index = 0
884 884 self.times = 1
885 885
886 886 def setup(self, data, m, n=1):
887 887 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
888 888 self.nHeights = data.shape[1]#.nHeights
889 889 self.start_index = 0
890 890 self.times = 1
891 891
892 892 def concat(self, data):
893 893
894 894 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
895 895 self.start_index = self.start_index + self.nHeights
896 896
897 897 def run(self, dataOut, m):
898 898 dataOut.flagNoData = True
899 899
900 900 if not self.isConfig:
901 901 self.setup(dataOut.data, m, 1)
902 902 self.isConfig = True
903 903
904 904 if dataOut.flagDataAsBlock:
905 905 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
906 906
907 907 else:
908 908 self.concat(dataOut.data)
909 909 self.times += 1
910 910 if self.times > m:
911 911 dataOut.data = self.buffer
912 912 self.reset()
913 913 dataOut.flagNoData = False
914 914 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
915 915 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
916 916 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
917 917 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
918 918 dataOut.ippSeconds *= m
919 919 return dataOut
920 920
921 921 class ProfileSelector(Operation):
922 922
923 923 profileIndex = None
924 924 # Tamanho total de los perfiles
925 925 nProfiles = None
926 926
927 927 def __init__(self, **kwargs):
928 928
929 929 Operation.__init__(self, **kwargs)
930 930 self.profileIndex = 0
931 931
932 932 def incProfileIndex(self):
933 933
934 934 self.profileIndex += 1
935 935
936 936 if self.profileIndex >= self.nProfiles:
937 937 self.profileIndex = 0
938 938
939 939 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
940 940
941 941 if profileIndex < minIndex:
942 942 return False
943 943
944 944 if profileIndex > maxIndex:
945 945 return False
946 946
947 947 return True
948 948
949 949 def isThisProfileInList(self, profileIndex, profileList):
950 950
951 951 if profileIndex not in profileList:
952 952 return False
953 953
954 954 return True
955 955
956 956 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
957 957
958 958 """
959 959 ProfileSelector:
960 960
961 961 Inputs:
962 962 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
963 963
964 964 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
965 965
966 966 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
967 967
968 968 """
969 969
970 970 if rangeList is not None:
971 971 if type(rangeList[0]) not in (tuple, list):
972 972 rangeList = [rangeList]
973 973
974 974 dataOut.flagNoData = True
975 975
976 976 if dataOut.flagDataAsBlock:
977 977 """
978 978 data dimension = [nChannels, nProfiles, nHeis]
979 979 """
980 980 if profileList != None:
981 981 dataOut.data = dataOut.data[:,profileList,:]
982 982
983 983 if profileRangeList != None:
984 984 minIndex = profileRangeList[0]
985 985 maxIndex = profileRangeList[1]
986 986 profileList = list(range(minIndex, maxIndex+1))
987 987
988 988 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
989 989
990 990 if rangeList != None:
991 991
992 992 profileList = []
993 993
994 994 for thisRange in rangeList:
995 995 minIndex = thisRange[0]
996 996 maxIndex = thisRange[1]
997 997
998 998 profileList.extend(list(range(minIndex, maxIndex+1)))
999 999
1000 1000 dataOut.data = dataOut.data[:,profileList,:]
1001 1001
1002 1002 dataOut.nProfiles = len(profileList)
1003 1003 dataOut.profileIndex = dataOut.nProfiles - 1
1004 1004 dataOut.flagNoData = False
1005 1005
1006 1006 return dataOut
1007 1007
1008 1008 """
1009 1009 data dimension = [nChannels, nHeis]
1010 1010 """
1011 1011
1012 1012 if profileList != None:
1013 1013
1014 1014 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1015 1015
1016 1016 self.nProfiles = len(profileList)
1017 1017 dataOut.nProfiles = self.nProfiles
1018 1018 dataOut.profileIndex = self.profileIndex
1019 1019 dataOut.flagNoData = False
1020 1020
1021 1021 self.incProfileIndex()
1022 1022 return dataOut
1023 1023
1024 1024 if profileRangeList != None:
1025 1025
1026 1026 minIndex = profileRangeList[0]
1027 1027 maxIndex = profileRangeList[1]
1028 1028
1029 1029 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1030 1030
1031 1031 self.nProfiles = maxIndex - minIndex + 1
1032 1032 dataOut.nProfiles = self.nProfiles
1033 1033 dataOut.profileIndex = self.profileIndex
1034 1034 dataOut.flagNoData = False
1035 1035
1036 1036 self.incProfileIndex()
1037 1037 return dataOut
1038 1038
1039 1039 if rangeList != None:
1040 1040
1041 1041 nProfiles = 0
1042 1042
1043 1043 for thisRange in rangeList:
1044 1044 minIndex = thisRange[0]
1045 1045 maxIndex = thisRange[1]
1046 1046
1047 1047 nProfiles += maxIndex - minIndex + 1
1048 1048
1049 1049 for thisRange in rangeList:
1050 1050
1051 1051 minIndex = thisRange[0]
1052 1052 maxIndex = thisRange[1]
1053 1053
1054 1054 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1055 1055
1056 1056 self.nProfiles = nProfiles
1057 1057 dataOut.nProfiles = self.nProfiles
1058 1058 dataOut.profileIndex = self.profileIndex
1059 1059 dataOut.flagNoData = False
1060 1060
1061 1061 self.incProfileIndex()
1062 1062
1063 1063 break
1064 1064
1065 1065 return dataOut
1066 1066
1067 1067
1068 1068 if beam != None: #beam is only for AMISR data
1069 1069 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1070 1070 dataOut.flagNoData = False
1071 1071 dataOut.profileIndex = self.profileIndex
1072 1072
1073 1073 self.incProfileIndex()
1074 1074
1075 1075 return dataOut
1076 1076
1077 1077 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1078 1078
1079 1079
1080 1080 class Reshaper(Operation):
1081 1081
1082 1082 def __init__(self, **kwargs):
1083 1083
1084 1084 Operation.__init__(self, **kwargs)
1085 1085
1086 1086 self.__buffer = None
1087 1087 self.__nitems = 0
1088 1088
1089 1089 def __appendProfile(self, dataOut, nTxs):
1090 1090
1091 1091 if self.__buffer is None:
1092 1092 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1093 1093 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1094 1094
1095 1095 ini = dataOut.nHeights * self.__nitems
1096 1096 end = ini + dataOut.nHeights
1097 1097
1098 1098 self.__buffer[:, ini:end] = dataOut.data
1099 1099
1100 1100 self.__nitems += 1
1101 1101
1102 1102 return int(self.__nitems*nTxs)
1103 1103
1104 1104 def __getBuffer(self):
1105 1105
1106 1106 if self.__nitems == int(1./self.__nTxs):
1107 1107
1108 1108 self.__nitems = 0
1109 1109
1110 1110 return self.__buffer.copy()
1111 1111
1112 1112 return None
1113 1113
1114 1114 def __checkInputs(self, dataOut, shape, nTxs):
1115 1115
1116 1116 if shape is None and nTxs is None:
1117 1117 raise ValueError("Reshaper: shape of factor should be defined")
1118 1118
1119 1119 if nTxs:
1120 1120 if nTxs < 0:
1121 1121 raise ValueError("nTxs should be greater than 0")
1122 1122
1123 1123 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1124 1124 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1125 1125
1126 1126 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1127 1127
1128 1128 return shape, nTxs
1129 1129
1130 1130 if len(shape) != 2 and len(shape) != 3:
1131 1131 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1132 1132
1133 1133 if len(shape) == 2:
1134 1134 shape_tuple = [dataOut.nChannels]
1135 1135 shape_tuple.extend(shape)
1136 1136 else:
1137 1137 shape_tuple = list(shape)
1138 1138
1139 1139 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1140 1140
1141 1141 return shape_tuple, nTxs
1142 1142
1143 1143 def run(self, dataOut, shape=None, nTxs=None):
1144 1144
1145 1145 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1146 1146
1147 1147 dataOut.flagNoData = True
1148 1148 profileIndex = None
1149 1149
1150 1150 if dataOut.flagDataAsBlock:
1151 1151
1152 1152 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1153 1153 dataOut.flagNoData = False
1154 1154
1155 1155 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1156 1156
1157 1157 else:
1158 1158
1159 1159 if self.__nTxs < 1:
1160 1160
1161 1161 self.__appendProfile(dataOut, self.__nTxs)
1162 1162 new_data = self.__getBuffer()
1163 1163
1164 1164 if new_data is not None:
1165 1165 dataOut.data = new_data
1166 1166 dataOut.flagNoData = False
1167 1167
1168 1168 profileIndex = dataOut.profileIndex*nTxs
1169 1169
1170 1170 else:
1171 1171 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1172 1172
1173 1173 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1174 1174
1175 1175 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1176 1176
1177 1177 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1178 1178
1179 1179 dataOut.profileIndex = profileIndex
1180 1180
1181 1181 dataOut.ippSeconds /= self.__nTxs
1182 1182
1183 1183 return dataOut
1184 1184
1185 1185 class SplitProfiles(Operation):
1186 1186
1187 1187 def __init__(self, **kwargs):
1188 1188
1189 1189 Operation.__init__(self, **kwargs)
1190 1190
1191 1191 def run(self, dataOut, n):
1192 1192
1193 1193 dataOut.flagNoData = True
1194 1194 profileIndex = None
1195 1195
1196 1196 if dataOut.flagDataAsBlock:
1197 1197
1198 1198 #nchannels, nprofiles, nsamples
1199 1199 shape = dataOut.data.shape
1200 1200
1201 1201 if shape[2] % n != 0:
1202 1202 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1203 1203
1204 1204 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1205 1205
1206 1206 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1207 1207 dataOut.flagNoData = False
1208 1208
1209 1209 profileIndex = int(dataOut.nProfiles/n) - 1
1210 1210
1211 1211 else:
1212 1212
1213 1213 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1214 1214
1215 1215 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1216 1216
1217 1217 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1218 1218
1219 1219 dataOut.nProfiles = int(dataOut.nProfiles*n)
1220 1220
1221 1221 dataOut.profileIndex = profileIndex
1222 1222
1223 1223 dataOut.ippSeconds /= n
1224 1224
1225 1225 return dataOut
1226 1226
1227 1227 class CombineProfiles(Operation):
1228 1228 def __init__(self, **kwargs):
1229 1229
1230 1230 Operation.__init__(self, **kwargs)
1231 1231
1232 1232 self.__remData = None
1233 1233 self.__profileIndex = 0
1234 1234
1235 1235 def run(self, dataOut, n):
1236 1236
1237 1237 dataOut.flagNoData = True
1238 1238 profileIndex = None
1239 1239
1240 1240 if dataOut.flagDataAsBlock:
1241 1241
1242 1242 #nchannels, nprofiles, nsamples
1243 1243 shape = dataOut.data.shape
1244 1244 new_shape = shape[0], shape[1]/n, shape[2]*n
1245 1245
1246 1246 if shape[1] % n != 0:
1247 1247 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1248 1248
1249 1249 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1250 1250 dataOut.flagNoData = False
1251 1251
1252 1252 profileIndex = int(dataOut.nProfiles*n) - 1
1253 1253
1254 1254 else:
1255 1255
1256 1256 #nchannels, nsamples
1257 1257 if self.__remData is None:
1258 1258 newData = dataOut.data
1259 1259 else:
1260 1260 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1261 1261
1262 1262 self.__profileIndex += 1
1263 1263
1264 1264 if self.__profileIndex < n:
1265 1265 self.__remData = newData
1266 1266 #continue
1267 1267 return
1268 1268
1269 1269 self.__profileIndex = 0
1270 1270 self.__remData = None
1271 1271
1272 1272 dataOut.data = newData
1273 1273 dataOut.flagNoData = False
1274 1274
1275 1275 profileIndex = dataOut.profileIndex/n
1276 1276
1277 1277
1278 1278 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1279 1279
1280 1280 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1281 1281
1282 1282 dataOut.nProfiles = int(dataOut.nProfiles/n)
1283 1283
1284 1284 dataOut.profileIndex = profileIndex
1285 1285
1286 1286 dataOut.ippSeconds *= n
1287 1287
1288 1288 return dataOut
1289 1289
1290 1290 class PulsePair(Operation):
1291 1291 '''
1292 1292 Function PulsePair(Signal Power, Velocity)
1293 1293 The real component of Lag[0] provides Intensity Information
1294 1294 The imag component of Lag[1] Phase provides Velocity Information
1295 1295
1296 1296 Configuration Parameters:
1297 1297 nPRF = Number of Several PRF
1298 1298 theta = Degree Azimuth angel Boundaries
1299 1299
1300 1300 Input:
1301 1301 self.dataOut
1302 1302 lag[N]
1303 1303 Affected:
1304 1304 self.dataOut.spc
1305 1305 '''
1306 1306 isConfig = False
1307 1307 __profIndex = 0
1308 1308 __initime = None
1309 1309 __lastdatatime = None
1310 1310 __buffer = None
1311 1311 noise = None
1312 1312 __dataReady = False
1313 1313 n = None
1314 1314 __nch = 0
1315 1315 __nHeis = 0
1316 1316 removeDC = False
1317 1317 ipp = None
1318 1318 lambda_ = 0
1319 1319
1320 1320 def __init__(self,**kwargs):
1321 1321 Operation.__init__(self,**kwargs)
1322 1322
1323 1323 def setup(self, dataOut, n = None, removeDC=False):
1324 1324 '''
1325 1325 n= Numero de PRF's de entrada
1326 1326 '''
1327 1327 self.__initime = None
1328 1328 ####print("[INICIO]-setup del METODO PULSE PAIR")
1329 1329 self.__lastdatatime = 0
1330 1330 self.__dataReady = False
1331 1331 self.__buffer = 0
1332 1332 self.__profIndex = 0
1333 1333 self.noise = None
1334 1334 self.__nch = dataOut.nChannels
1335 1335 self.__nHeis = dataOut.nHeights
1336 1336 self.removeDC = removeDC
1337 1337 self.lambda_ = 3.0e8/(9345.0e6)
1338 1338 self.ippSec = dataOut.ippSeconds
1339 1339 self.nCohInt = dataOut.nCohInt
1340 1340 ####print("IPPseconds",dataOut.ippSeconds)
1341 1341 ####print("ELVALOR DE n es:", n)
1342 1342 if n == None:
1343 1343 raise ValueError("n should be specified.")
1344 1344
1345 1345 if n != None:
1346 1346 if n<2:
1347 1347 raise ValueError("n should be greater than 2")
1348 1348
1349 1349 self.n = n
1350 1350 self.__nProf = n
1351 1351
1352 1352 self.__buffer = numpy.zeros((dataOut.nChannels,
1353 1353 n,
1354 1354 dataOut.nHeights),
1355 1355 dtype='complex')
1356 1356
1357 1357 def putData(self,data):
1358 1358 '''
1359 1359 Add a profile to he __buffer and increase in one the __profiel Index
1360 1360 '''
1361 1361 self.__buffer[:,self.__profIndex,:]= data
1362 1362 self.__profIndex += 1
1363 1363 return
1364 1364
1365 1365 def pushData(self,dataOut):
1366 1366 '''
1367 1367 Return the PULSEPAIR and the profiles used in the operation
1368 1368 Affected : self.__profileIndex
1369 1369 '''
1370 1370 #----------------- Remove DC-----------------------------------
1371 1371 if self.removeDC==True:
1372 1372 mean = numpy.mean(self.__buffer,1)
1373 1373 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1374 1374 dc= numpy.tile(tmp,[1,self.__nProf,1])
1375 1375 self.__buffer = self.__buffer - dc
1376 1376 #------------------Calculo de Potencia ------------------------
1377 1377 pair0 = self.__buffer*numpy.conj(self.__buffer)
1378 1378 pair0 = pair0.real
1379 1379 lag_0 = numpy.sum(pair0,1)
1380 1380 #-----------------Calculo de Cscp------------------------------ New
1381 cspc_pair01 = self.__buffer[0]*__self.buffer[1]
1381 cspc_pair01 = self.__buffer[0]*self.__buffer[1]
1382 1382 #------------------Calculo de Ruido x canal--------------------
1383 1383 self.noise = numpy.zeros(self.__nch)
1384 1384 for i in range(self.__nch):
1385 1385 daux = numpy.sort(pair0[i,:,:],axis= None)
1386 1386 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1387 1387
1388 1388 self.noise = self.noise.reshape(self.__nch,1)
1389 1389 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1390 1390 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1391 1391 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1392 1392 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1393 1393 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1394 1394 #-------------------- Power --------------------------------------------------
1395 1395 data_power = lag_0/(self.n*self.nCohInt)
1396 1396 #--------------------CCF------------------------------------------------------
1397 1397 data_ccf =numpy.sum(cspc_pair01,axis=0)/(self.n*self.nCohInt)
1398 1398 #------------------ Senal --------------------------------------------------
1399 1399 data_intensity = pair0 - noise_buffer
1400 1400 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1401 1401 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1402 1402 for i in range(self.__nch):
1403 1403 for j in range(self.__nHeis):
1404 1404 if data_intensity[i][j] < 0:
1405 1405 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1406 1406
1407 1407 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1408 1408 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1409 1409 lag_1 = numpy.sum(pair1,1)
1410 1410 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1411 1411 data_velocity = (self.lambda_/2.0)*data_freq
1412 1412
1413 1413 #---------------- Potencia promedio estimada de la Senal-----------
1414 1414 lag_0 = lag_0/self.n
1415 1415 S = lag_0-self.noise
1416 1416
1417 1417 #---------------- Frecuencia Doppler promedio ---------------------
1418 1418 lag_1 = lag_1/(self.n-1)
1419 1419 R1 = numpy.abs(lag_1)
1420 1420
1421 1421 #---------------- Calculo del SNR----------------------------------
1422 1422 data_snrPP = S/self.noise
1423 1423 for i in range(self.__nch):
1424 1424 for j in range(self.__nHeis):
1425 1425 if data_snrPP[i][j] < 1.e-20:
1426 1426 data_snrPP[i][j] = 1.e-20
1427 1427
1428 1428 #----------------- Calculo del ancho espectral ----------------------
1429 1429 L = S/R1
1430 1430 L = numpy.where(L<0,1,L)
1431 1431 L = numpy.log(L)
1432 1432 tmp = numpy.sqrt(numpy.absolute(L))
1433 1433 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1434 1434 n = self.__profIndex
1435 1435
1436 1436 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1437 1437 self.__profIndex = 0
1438 1438 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,data_ccf,n
1439 1439
1440 1440
1441 1441 def pulsePairbyProfiles(self,dataOut):
1442 1442
1443 1443 self.__dataReady = False
1444 1444 data_power = None
1445 1445 data_intensity = None
1446 1446 data_velocity = None
1447 1447 data_specwidth = None
1448 1448 data_snrPP = None
1449 1449 data_ccf = None
1450 1450 self.putData(data=dataOut.data)
1451 1451 if self.__profIndex == self.n:
1452 1452 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, n = self.pushData(dataOut=dataOut)
1453 1453 self.__dataReady = True
1454 1454
1455 1455 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf
1456 1456
1457 1457
1458 1458 def pulsePairOp(self, dataOut, datatime= None):
1459 1459
1460 1460 if self.__initime == None:
1461 1461 self.__initime = datatime
1462 1462 data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf = self.pulsePairbyProfiles(dataOut)
1463 1463 self.__lastdatatime = datatime
1464 1464
1465 1465 if data_power is None:
1466 return None, None, None,None,None,None
1466 return None, None, None,None,None,None,None
1467 1467
1468 1468 avgdatatime = self.__initime
1469 1469 deltatime = datatime - self.__lastdatatime
1470 1470 self.__initime = datatime
1471 1471
1472 1472 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf, avgdatatime
1473 1473
1474 1474 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1475 1475
1476 1476 if not self.isConfig:
1477 1477 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1478 1478 self.isConfig = True
1479 1479 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1480 1480 dataOut.flagNoData = True
1481 1481
1482 1482 if self.__dataReady:
1483 1483 ###print("READY ----------------------------------")
1484 1484 dataOut.nCohInt *= self.n
1485 1485 dataOut.dataPP_POW = data_intensity # S
1486 1486 dataOut.dataPP_POWER = data_power # P valor que corresponde a POTENCIA MOMENTO
1487 1487 dataOut.dataPP_DOP = data_velocity
1488 1488 dataOut.dataPP_SNR = data_snrPP
1489 1489 dataOut.dataPP_WIDTH = data_specwidth
1490 1490 dataOut.dataPP_CCF = data_ccf
1491 1491 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1492 1492 dataOut.nProfiles = int(dataOut.nProfiles/n)
1493 1493 dataOut.utctime = avgdatatime
1494 1494 dataOut.flagNoData = False
1495 1495 return dataOut
1496 1496
1497 1497
1498 1498
1499 1499 # import collections
1500 1500 # from scipy.stats import mode
1501 1501 #
1502 1502 # class Synchronize(Operation):
1503 1503 #
1504 1504 # isConfig = False
1505 1505 # __profIndex = 0
1506 1506 #
1507 1507 # def __init__(self, **kwargs):
1508 1508 #
1509 1509 # Operation.__init__(self, **kwargs)
1510 1510 # # self.isConfig = False
1511 1511 # self.__powBuffer = None
1512 1512 # self.__startIndex = 0
1513 1513 # self.__pulseFound = False
1514 1514 #
1515 1515 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1516 1516 #
1517 1517 # #Read data
1518 1518 #
1519 1519 # powerdB = dataOut.getPower(channel = channel)
1520 1520 # noisedB = dataOut.getNoise(channel = channel)[0]
1521 1521 #
1522 1522 # self.__powBuffer.extend(powerdB.flatten())
1523 1523 #
1524 1524 # dataArray = numpy.array(self.__powBuffer)
1525 1525 #
1526 1526 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1527 1527 #
1528 1528 # maxValue = numpy.nanmax(filteredPower)
1529 1529 #
1530 1530 # if maxValue < noisedB + 10:
1531 1531 # #No se encuentra ningun pulso de transmision
1532 1532 # return None
1533 1533 #
1534 1534 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1535 1535 #
1536 1536 # if len(maxValuesIndex) < 2:
1537 1537 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1538 1538 # return None
1539 1539 #
1540 1540 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1541 1541 #
1542 1542 # #Seleccionar solo valores con un espaciamiento de nSamples
1543 1543 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1544 1544 #
1545 1545 # if len(pulseIndex) < 2:
1546 1546 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1547 1547 # return None
1548 1548 #
1549 1549 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1550 1550 #
1551 1551 # #remover senales que se distancien menos de 10 unidades o muestras
1552 1552 # #(No deberian existir IPP menor a 10 unidades)
1553 1553 #
1554 1554 # realIndex = numpy.where(spacing > 10 )[0]
1555 1555 #
1556 1556 # if len(realIndex) < 2:
1557 1557 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1558 1558 # return None
1559 1559 #
1560 1560 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1561 1561 # realPulseIndex = pulseIndex[realIndex]
1562 1562 #
1563 1563 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1564 1564 #
1565 1565 # print "IPP = %d samples" %period
1566 1566 #
1567 1567 # self.__newNSamples = dataOut.nHeights #int(period)
1568 1568 # self.__startIndex = int(realPulseIndex[0])
1569 1569 #
1570 1570 # return 1
1571 1571 #
1572 1572 #
1573 1573 # def setup(self, nSamples, nChannels, buffer_size = 4):
1574 1574 #
1575 1575 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1576 1576 # maxlen = buffer_size*nSamples)
1577 1577 #
1578 1578 # bufferList = []
1579 1579 #
1580 1580 # for i in range(nChannels):
1581 1581 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1582 1582 # maxlen = buffer_size*nSamples)
1583 1583 #
1584 1584 # bufferList.append(bufferByChannel)
1585 1585 #
1586 1586 # self.__nSamples = nSamples
1587 1587 # self.__nChannels = nChannels
1588 1588 # self.__bufferList = bufferList
1589 1589 #
1590 1590 # def run(self, dataOut, channel = 0):
1591 1591 #
1592 1592 # if not self.isConfig:
1593 1593 # nSamples = dataOut.nHeights
1594 1594 # nChannels = dataOut.nChannels
1595 1595 # self.setup(nSamples, nChannels)
1596 1596 # self.isConfig = True
1597 1597 #
1598 1598 # #Append new data to internal buffer
1599 1599 # for thisChannel in range(self.__nChannels):
1600 1600 # bufferByChannel = self.__bufferList[thisChannel]
1601 1601 # bufferByChannel.extend(dataOut.data[thisChannel])
1602 1602 #
1603 1603 # if self.__pulseFound:
1604 1604 # self.__startIndex -= self.__nSamples
1605 1605 #
1606 1606 # #Finding Tx Pulse
1607 1607 # if not self.__pulseFound:
1608 1608 # indexFound = self.__findTxPulse(dataOut, channel)
1609 1609 #
1610 1610 # if indexFound == None:
1611 1611 # dataOut.flagNoData = True
1612 1612 # return
1613 1613 #
1614 1614 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1615 1615 # self.__pulseFound = True
1616 1616 # self.__startIndex = indexFound
1617 1617 #
1618 1618 # #If pulse was found ...
1619 1619 # for thisChannel in range(self.__nChannels):
1620 1620 # bufferByChannel = self.__bufferList[thisChannel]
1621 1621 # #print self.__startIndex
1622 1622 # x = numpy.array(bufferByChannel)
1623 1623 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1624 1624 #
1625 1625 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1626 1626 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1627 1627 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1628 1628 #
1629 1629 # dataOut.data = self.__arrayBuffer
1630 1630 #
1631 1631 # self.__startIndex += self.__newNSamples
1632 1632 #
1633 1633 # return
Show file before | Show file after | Hide comments
@@ -1,275 +1,303
1 1 #!python
2 2 '''
3 3 '''
4 4
5 5 import os, sys
6 6 import datetime
7 7 import time
8 8
9 9 #path = os.path.dirname(os.getcwd())
10 10 #path = os.path.dirname(path)
11 11 #sys.path.insert(0, path)
12 12
13 13 from schainpy.controller import Project
14 14
15 15 desc = "USRP_test"
16 16 filename = "USRP_processing.xml"
17 17 controllerObj = Project()
18 18 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
19 19
20 20 ############## USED TO PLOT IQ VOLTAGE, POWER AND SPECTRA #############
21 21
22 22 #######################################################################
23 23 ######PATH DE LECTURA, ESCRITURA, GRAFICOS Y ENVIO WEB#################
24 24 #######################################################################
25 25 #path = '/media/data/data/vientos/57.2063km/echoes/NCO_Woodman'
26 26 #path = '/DATA_RM/TEST_INTEGRACION'
27 27 #path = '/DATA_RM/PRUEBA_USRP_RP'
28 path = '/DATA_RM/PRUEBA_USRP_RP'
28 #path = '/DATA_RM/PRUEBA_USRP_RP'
29 29
30 figpath = '/home/soporte/Pictures/TEST_RP_0001'
31 figpath = '/home/soporte/Pictures/TEST_RP_6000'
32 figpath = '/home/soporte/Pictures/USRP'
30 path = '/DATA_RM/TEST_2M'
31 path = '/DATA_RM/TEST_2M_UD'
32 path = '/DATA_RM/2MHZ17022022'
33 path = '/DATA_RM/10MHZTEST/'
34 path = '/DATA_RM/10MHZDRONE/'
35
36 #figpath = '/home/soporte/Pictures/TEST_RP_0001'
37 #figpath = '/home/soporte/Pictures/TEST_RP_6000'
38 figpath = '/home/soporte/Pictures/USRP_TEST_2M'
39 figpath = '/home/soporte/Pictures/USRP_TEST_2M_UD'
40 figpaht = '/home/soporte/Pictures/10MHZDRONE'
33 41 #remotefolder = "/home/wmaster/graficos"
34 42 #######################################################################
35 43 ################# RANGO DE PLOTEO######################################
36 44 #######################################################################
37 dBmin = '-5'
38 dBmax = '20'
45 dBmin = '20'
46 dBmax = '60'
39 47 xmin = '0'
40 48 xmax ='24'
41 49 ymin = '0'
42 50 ymax = '600'
43 51 #######################################################################
44 52 ########################FECHA##########################################
45 53 #######################################################################
46 54 str = datetime.date.today()
47 55 today = str.strftime("%Y/%m/%d")
48 56 str2 = str - datetime.timedelta(days=1)
49 57 yesterday = str2.strftime("%Y/%m/%d")
50 58 #######################################################################
51 59 ######################## UNIDAD DE LECTURA#############################
52 60 #######################################################################
53 61 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
54 62 path=path,
55 startDate="2021/07/02",#today,
56 endDate="2021/07/02",#today,
57 startTime='14:50:00',# inicio libre
63 startDate="2022/02/18",#today,
64 endDate="2022/02/18",#today,
65 startTime='00:00:00',# inicio libre
58 66 #startTime='00:00:00',
59 endTime='14:55:59',
67 endTime='23:59:59',
60 68 delay=0,
61 69 #set=0,
62 online=0,
70 online=1,
63 71 walk=1,
64 ippKm = 6000)
72 ippKm = 60)
65 73
66 74 opObj11 = readUnitConfObj.addOperation(name='printInfo')
67 75 #opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
68 76 #######################################################################
69 77 ################ OPERACIONES DOMINIO DEL TIEMPO########################
70 78 #######################################################################
71 79
72 80 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
73 81
82 '''
83 # OJO SCOPE
84 opObj10 = procUnitConfObjA.addOperation(name='ScopePlot', optype='external')
85 opObj10.addParameter(name='id', value='10', format='int')
86 ##opObj10.addParameter(name='xmin', value='0', format='int')
87 ##opObj10.addParameter(name='xmax', value='50', format='int')
88 opObj10.addParameter(name='type', value='iq')
89 opObj10.addParameter(name='ymin', value='-1200', format='int')
90 opObj10.addParameter(name='ymax', value='1200', format='int')
91 opObj10.addParameter(name='save', value=figpath, format='str')
92 opObj10.addParameter(name='save_period', value=10, format='int')
93 '''
94 '''
74 95 opObj11 = procUnitConfObjA.addOperation(name='selectHeights')
75 96 opObj11.addParameter(name='minIndex', value='1', format='int')
76 97 # opObj11.addParameter(name='maxIndex', value='10000', format='int')
77 98 opObj11.addParameter(name='maxIndex', value='39980', format='int')
78
99 '''
79 100 #
80 101 # codigo64='1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1,1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0,'+\
81 102 # '1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1,0,0,0,1,0,0,1,0,0,0,0,1,1,1,0,1,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1'
82 103
83 104 #opObj11 = procUnitConfObjA.addOperation(name='setRadarFrequency')
84 105 #opObj11.addParameter(name='frequency', value='49920000')
85 106
86 107 '''
87 108 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
88 109 opObj11.addParameter(name='n', value='625', format='int')#10
89 110 opObj11.addParameter(name='removeDC', value=1, format='int')
90 111 '''
91 112
92 113 # Ploteo TEST
93 114 '''
94 115 opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')
95 116 opObj11 = procUnitConfObjA.addOperation(name='PulsepairSignalPlot', optype='other')
96 117 opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')
97 118 #opObj11.addParameter(name='xmax', value=8)
98 119 opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')
99 120 '''
100 121 # OJO SCOPE
101 122 #opObj10 = procUnitConfObjA.addOperation(name='ScopePlot', optype='external')
102 123 #opObj10.addParameter(name='id', value='10', format='int')
103 124 ##opObj10.addParameter(name='xmin', value='0', format='int')
104 125 ##opObj10.addParameter(name='xmax', value='50', format='int')
105 126 #opObj10.addParameter(name='type', value='iq')
106 127 ##opObj10.addParameter(name='ymin', value='-5000', format='int')
107 128 ##opObj10.addParameter(name='ymax', value='8500', format='int')
108 129 #opObj11.addParameter(name='save', value=figpath, format='str')
109 130 #opObj11.addParameter(name='save_period', value=10, format='int')
110 131
111 132 #opObj10 = procUnitConfObjA.addOperation(name='setH0')
112 133 #opObj10.addParameter(name='h0', value='-5000', format='float')
113 134
114 135 #opObj11 = procUnitConfObjA.addOperation(name='filterByHeights')
115 136 #opObj11.addParameter(name='window', value='1', format='int')
116 137
117 138 #codigo='1,1,-1,1,1,-1,1,-1,-1,1,-1,-1,-1,1,-1,-1,-1,1,-1,-1,-1,1,1,1,1,-1,-1,-1'
118 139 #opObj11 = procUnitConfObjSousy.addOperation(name='Decoder', optype='other')
119 140 #opObj11.addParameter(name='code', value=codigo, format='floatlist')
120 141 #opObj11.addParameter(name='nCode', value='1', format='int')
121 142 #opObj11.addParameter(name='nBaud', value='28', format='int')
122 143
123 144 #opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other')
124 145 #opObj11.addParameter(name='n', value='100', format='int')
125 146
126 147 #######################################################################
127 148 ########## OPERACIONES ParametersProc########################
128 149 #######################################################################
129 150 ###procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
130 151 '''
131 152
132 153 opObj11 = procUnitConfObjA.addOperation(name='PedestalInformation')
133 154 opObj11.addParameter(name='path_ped', value=path_ped)
134 155 opObj11.addParameter(name='path_adq', value=path_adq)
135 156 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
136 157 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
137 158 opObj11.addParameter(name='blocksPerfile', value='100', format='int')
138 159 opObj11.addParameter(name='f_a_p', value='25', format='int')
139 160 opObj11.addParameter(name='online', value='0', format='int')
140 161
141 162 opObj11 = procUnitConfObjA.addOperation(name='Block360')
142 163 opObj11.addParameter(name='n', value='40', format='int')
143 164
144 165 opObj11= procUnitConfObjA.addOperation(name='WeatherPlot',optype='other')
145 166 opObj11.addParameter(name='save', value=figpath)
146 167 opObj11.addParameter(name='save_period', value=1)
147 168
148 169 8
149 170 '''
150 171
172
173
174 opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other')
175 opObj11.addParameter(name='n', value='250', format='int')
176
151 177 #######################################################################
152 178 ########## OPERACIONES DOMINIO DE LA FRECUENCIA########################
153 179 #######################################################################
154 180
155 #procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
156 #procUnitConfObjB.addParameter(name='nFFTPoints', value='32', format='int')
157 #procUnitConfObjB.addParameter(name='nProfiles', value='32', format='int')
158
181 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
182 procUnitConfObjB.addParameter(name='nFFTPoints', value='32', format='int')
183 procUnitConfObjB.addParameter(name='nProfiles', value='32', format='int')
184 '''
159 185 procUnitConfObjC = controllerObj.addProcUnit(datatype='SpectraHeisProc', inputId=procUnitConfObjA.getId())
160 186 #procUnitConfObjB.addParameter(name='nFFTPoints', value='64', format='int')
161 187 #procUnitConfObjB.addParameter(name='nProfiles', value='64', format='int')
162 188 opObj11 = procUnitConfObjC.addOperation(name='IncohInt4SpectraHeis', optype='other')
163 189 #opObj11.addParameter(name='timeInterval', value='4', format='int')
164 190 opObj11.addParameter(name='n', value='100', format='int')
165 191
166 192 #procUnitConfObjB.addParameter(name='pairsList', value='(0,0),(1,1),(0,1)', format='pairsList')
167 193
168 194 #opObj13 = procUnitConfObjB.addOperation(name='removeDC')
169 195 #opObj13.addParameter(name='mode', value='2', format='int')
170 196
171 197 #opObj11 = procUnitConfObjB.addOperation(name='IncohInt', optype='other')
172 198 #opObj11.addParameter(name='n', value='8', format='float')
173 199 #######################################################################
174 200 ########## PLOTEO DOMINIO DE LA FRECUENCIA#############################
175 201 #######################################################################
176 202 #----
177
203 '''
204 '''
178 205 opObj11 = procUnitConfObjC.addOperation(name='SpectraHeisPlot')
179 206 opObj11.addParameter(name='id', value='10', format='int')
180 207 opObj11.addParameter(name='wintitle', value='Spectra_Alturas', format='str')
181 208 #opObj11.addParameter(name='xmin', value=-100000, format='float')
182 209 #opObj11.addParameter(name='xmax', value=100000, format='float')
183 210 opObj11.addParameter(name='oneFigure', value=False,format='bool')
184 211 #opObj11.addParameter(name='zmin', value=-10, format='int')
185 212 #opObj11.addParameter(name='zmax', value=40, format='int')
186 213 opObj11.addParameter(name='ymin', value=10, format='int')
187 214 opObj11.addParameter(name='ymax', value=55, format='int')
188 215 opObj11.addParameter(name='grid', value=True, format='bool')
189 216 #opObj11.addParameter(name='showprofile', value='1', format='int')
190 217 opObj11.addParameter(name='save', value=figpath, format='str')
191 218 #opObj11.addParameter(name='save_period', value=10, format='int')
192
219 '''
193 220 '''
194 221 opObj11 = procUnitConfObjC.addOperation(name='RTIHeisPlot')
195 222 opObj11.addParameter(name='id', value='10', format='int')
196 223 opObj11.addParameter(name='wintitle', value='RTI_Alturas', format='str')
197 224 opObj11.addParameter(name='xmin', value=11.0, format='float')
198 225 opObj11.addParameter(name='xmax', value=18.0, format='float')
199 226 opObj11.addParameter(name='zmin', value=10, format='int')
200 227 opObj11.addParameter(name='zmax', value=30, format='int')
201 228 opObj11.addParameter(name='ymin', value=5, format='int')
202 229 opObj11.addParameter(name='ymax', value=28, format='int')
203 230 opObj11.addParameter(name='showprofile', value='1', format='int')
204 231 opObj11.addParameter(name='save', value=figpath, format='str')
205 232 opObj11.addParameter(name='save_period', value=10, format='int')
206 233 '''
207 '''
234
208 235 #SpectraPlot
209 236
210 237 opObj11 = procUnitConfObjB.addOperation(name='SpectraPlot', optype='external')
211 238 opObj11.addParameter(name='id', value='1', format='int')
212 239 opObj11.addParameter(name='wintitle', value='Spectra', format='str')
213 240 #opObj11.addParameter(name='xmin', value=-0.01, format='float')
214 241 #opObj11.addParameter(name='xmax', value=0.01, format='float')
215 242 opObj11.addParameter(name='zmin', value=dBmin, format='int')
216 243 opObj11.addParameter(name='zmax', value=dBmax, format='int')
217 244 #opObj11.addParameter(name='ymin', value=ymin, format='int')
218 245 #opObj11.addParameter(name='ymax', value=ymax, format='int')
219 246 opObj11.addParameter(name='showprofile', value='1', format='int')
220 247 opObj11.addParameter(name='save', value=figpath, format='str')
221 248 opObj11.addParameter(name='save_period', value=10, format='int')
222 249
250
223 251 #RTIPLOT
224 252
225 253 opObj11 = procUnitConfObjB.addOperation(name='RTIPlot', optype='external')
226 254 opObj11.addParameter(name='id', value='2', format='int')
227 255 opObj11.addParameter(name='wintitle', value='RTIPlot', format='str')
228 256 opObj11.addParameter(name='zmin', value=dBmin, format='int')
229 257 opObj11.addParameter(name='zmax', value=dBmax, format='int')
230 258 #opObj11.addParameter(name='ymin', value=ymin, format='int')
231 259 #opObj11.addParameter(name='ymax', value=ymax, format='int')
232 260 #opObj11.addParameter(name='xmin', value=15, format='int')
233 261 #opObj11.addParameter(name='xmax', value=16, format='int')
234 262
235 263 opObj11.addParameter(name='showprofile', value='1', format='int')
236 264 opObj11.addParameter(name='save', value=figpath, format='str')
237 265 opObj11.addParameter(name='save_period', value=10, format='int')
238 266
239 267 '''
240 268 # opObj11 = procUnitConfObjB.addOperation(name='CrossSpectraPlot', optype='other')
241 269 # opObj11.addParameter(name='id', value='3', format='int')
242 270 # opObj11.addParameter(name='wintitle', value='CrossSpectraPlot', format='str')
243 271 # opObj11.addParameter(name='ymin', value=ymin, format='int')
244 272 # opObj11.addParameter(name='ymax', value=ymax, format='int')
245 273 # opObj11.addParameter(name='phase_cmap', value='jet', format='str')
246 274 # opObj11.addParameter(name='zmin', value=dBmin, format='int')
247 275 # opObj11.addParameter(name='zmax', value=dBmax, format='int')
248 276 # opObj11.addParameter(name='figpath', value=figures_path, format='str')
249 277 # opObj11.addParameter(name='save', value=0, format='bool')
250 278 # opObj11.addParameter(name='pairsList', value='(0,1)', format='pairsList')
251 279 # #
252 280 # opObj11 = procUnitConfObjB.addOperation(name='CoherenceMap', optype='other')
253 281 # opObj11.addParameter(name='id', value='4', format='int')
254 282 # opObj11.addParameter(name='wintitle', value='Coherence', format='str')
255 283 # opObj11.addParameter(name='phase_cmap', value='jet', format='str')
256 284 # opObj11.addParameter(name='xmin', value=xmin, format='float')
257 285 # opObj11.addParameter(name='xmax', value=xmax, format='float')
258 286 # opObj11.addParameter(name='figpath', value=figures_path, format='str')
259 287 # opObj11.addParameter(name='save', value=0, format='bool')
260 288 # opObj11.addParameter(name='pairsList', value='(0,1)', format='pairsList')
261 289 #
262
290 '''
263 291 '''
264 292 #######################################################################
265 293 ############### UNIDAD DE ESCRITURA ###################################
266 294 #######################################################################
267 295 #opObj11 = procUnitConfObjB.addOperation(name='SpectraWriter', optype='other')
268 296 #opObj11.addParameter(name='path', value=wr_path)
269 297 #opObj11.addParameter(name='blocksPerFile', value='50', format='int')
270 298 print ("Escribiendo el archivo XML")
271 299 print ("Leyendo el archivo XML")
272 300 '''
273 301
274 302
275 303 controllerObj.start()
Show file before | Show file after | Hide comments
@@ -1,184 +1,185
1 1 # Ing. AVP
2 2 # 04/01/2022
3 3 # ARCHIVO DE LECTURA
4 4 import os, sys
5 5 import datetime
6 6 import time
7 7 import numpy
8 8 import json
9 9 from ext_met import getfirstFilefromPath,getDatavaluefromDirFilename
10 10 from schainpy.controller import Project
11 11 #-----------------------------------------------------------------------------------------
12 12 # path_ped = "/DATA_RM/TEST_PEDESTAL/P20211110-171003"
13 13 ## print("PATH PEDESTAL :",path_ped)
14 14
15 15 print("[SETUP]-RADAR METEOROLOGICO-")
16 16 path_ped = "/DATA_RM/TEST_PEDESTAL/P20211111-173856"
17 17 print("PATH PEDESTAL :",path_ped)
18 18 path_adq = "/DATA_RM/11"
19 path_adq = "/DATA_RM/10MHZDRONE/"
19 20 print("PATH DATA :",path_adq)
20 21
21 22
22 23 figpath_pp_rti = "/home/soporte/Pictures/TEST_PP_RTI"
23 24 print("PATH PP RTI :",figpath_pp_rti)
24 25 figpath_pp_ppi = "/home/soporte/Pictures/TEST_PP_PPI"
25 26 print("PATH PP PPI :",figpath_pp_ppi)
26 27 path_pp_save_int = "/DATA_RM/TEST_NEW_FORMAT"
27 28 print("PATH SAVE PP INT :",path_pp_save_int)
28 29 print(" ")
29 30 #-------------------------------------------------------------------------------------------
30 31 print("SELECCIONAR MODO: PPI (0) O RHI (1)")
31 32 mode_wr = 0
32 33 if mode_wr==0:
33 34 print("[ ON ] MODE PPI")
34 35 list_ped = getfirstFilefromPath(path=path_ped,meta="PE",ext=".hdf5")
35 36 ff_pedestal = list_ped[2]
36 37 azi_vel = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="azi_vel")
37 38 V = round(azi_vel[0])
38 39 print("VELOCIDAD AZI :", int(numpy.mean(azi_vel)),"Β°/seg")
39 40 else:
40 41 print("[ ON ] MODE RHI")
41 42 list_ped = getfirstFilefromPath(path=path_ped,meta="PE",ext=".hdf5")
42 43 ff_pedestal = list_ped[2]
43 44 V = round(ele_vel[0])
44 45 ele_vel = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="ele_vel")
45 46 print("VELOCIDAD ELE :", int(numpy.mean(ele_vel)),"Β°/seg")
46 47 print(" ")
47 48 #---------------------------------------------------------------------------------------
48 49 print("SELECCIONAR MODO: PULSE PAIR (0) O FREQUENCY (1)")
49 50 mode_proc = 0
50 51 if mode_proc==0:
51 52 print("[ ON ] MODE PULSEPAIR")
52 53 else:
53 54 print("[ ON ] MODE FREQUENCY")
54 55 ipp = 60.0
55 56 print("IPP(Km.) : %1.2f"%ipp)
56 57 ipp_sec = (ipp*1.0e3/150.0)*1.0e-6
57 58 print("IPP(useg.) : %1.2f"%(ipp_sec*(1.0e6)))
58 59 VEL=V
59 60 n= int(1/(VEL*ipp_sec))
60 61 print("NΒ° Profiles : ", n)
61 62 #---------------------------------------------------------------------------------------
62 63 plot_rti = 0
63 plot_ppi = 0
64 plot_ppi = 1
64 65 integration = 1
65 save = 1
66 save = 0
66 67 #---------------------------RANGO DE PLOTEO----------------------------------
67 68 dBmin = '1'
68 69 dBmax = '85'
69 70 xmin = '14'
70 71 xmax = '16'
71 72 ymin = '0'
72 73 ymax = '600'
73 74 #----------------------------------------------------------------------------
74 75 time.sleep(3)
75 76 #---------------------SIGNAL CHAIN ------------------------------------
76 77 desc_wr= {
77 78 'Data': {
78 79 'dataPP_POW': 'Power',
79 80 'utctime': 'Time',
80 81 'azimuth': 'az',
81 82 'elevation':'el'
82 83 },
83 84 'Metadata': {
84 85 'heightList': 'range',
85 86 'channelList': 'Channels'
86 87 }
87 88 }
88 89
89 90
90 91 desc = "USRP_WEATHER_RADAR"
91 92 filename = "USRP_processing.xml"
92 93 controllerObj = Project()
93 94 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
94 95 #---------------------UNIDAD DE LECTURA--------------------------------
95 96 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
96 97 path=path_adq,
97 startDate="2021/11/11",#today,
98 endDate="2021/12/30",#today,
99 startTime='17:39:17',
98 startDate="2022/02/19",#today,
99 endDate="2022/02/18",#today,
100 startTime='00:00:00',
100 101 endTime='23:59:59',
101 102 delay=0,
102 103 #set=0,
103 104 online=0,
104 105 walk=1,
105 106 ippKm=ipp)
106 107
107 108 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc',inputId=readUnitConfObj.getId())
108 109
109 110 opObj11 = procUnitConfObjA.addOperation(name='selectHeights')
110 111 opObj11.addParameter(name='minIndex', value='1', format='int')
111 112 # opObj11.addParameter(name='maxIndex', value='10000', format='int')
112 113 opObj11.addParameter(name='maxIndex', value='400', format='int')
113 114
114 115 if mode_proc==0:
115 116 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
116 117 opObj11.addParameter(name='n', value=int(n), format='int')
117 118 procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
118 119 # REVISAR EL test_sim00013.py
119 120 if plot_rti==1:
120 121 opObj11 = procUnitConfObjB.addOperation(name='GenericRTIPlot',optype='external')
121 122 opObj11.addParameter(name='attr_data', value='dataPP_POW')
122 123 opObj11.addParameter(name='colormap', value='jet')
123 124 opObj11.addParameter(name='xmin', value=xmin)
124 125 opObj11.addParameter(name='xmax', value=xmax)
125 126 opObj11.addParameter(name='zmin', value=dBmin)
126 127 opObj11.addParameter(name='zmax', value=dBmax)
127 128 opObj11.addParameter(name='save', value=figpath_pp_rti)
128 129 opObj11.addParameter(name='showprofile', value=0)
129 130 opObj11.addParameter(name='save_period', value=50)
130 131 if integration==1:
131 132 opObj11 = procUnitConfObjB.addOperation(name='PedestalInformation')
132 133 opObj11.addParameter(name='path_ped', value=path_ped)
133 134 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
134 135 opObj11.addParameter(name='wr_exp', value='PPI')
135 136 #------------------------------------------------------------------------------
136 137 '''
137 138 opObj11.addParameter(name='Datatype', value='RadialSet')
138 139 opObj11.addParameter(name='Scantype', value='PPI')
139 140 opObj11.addParameter(name='Latitude', value='-11.96')
140 141 opObj11.addParameter(name='Longitud', value='-76.54')
141 142 opObj11.addParameter(name='Heading', value='293')
142 143 opObj11.addParameter(name='Height', value='293')
143 144 opObj11.addParameter(name='Waveform', value='OFM')
144 145 opObj11.addParameter(name='PRF', value='2000')
145 146 opObj11.addParameter(name='CreatedBy', value='WeatherRadarJROTeam')
146 147 opObj11.addParameter(name='ContactInformation', value='avaldez@igp.gob.pe')
147 148 '''
148 149 if plot_ppi==1:
149 150 opObj11 = procUnitConfObjB.addOperation(name='Block360')
150 151 opObj11.addParameter(name='n', value='10', format='int')
151 152 opObj11.addParameter(name='mode', value=mode_proc, format='int')
152 153 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
153 154 opObj11= procUnitConfObjB.addOperation(name='WeatherPlot',optype='other')
154 155 opObj11.addParameter(name='save', value=figpath_pp_ppi)
155 156 opObj11.addParameter(name='save_period', value=1)
156 157
157 158 if save==1:
158 159 opObj10 = procUnitConfObjB.addOperation(name='HDFWriter')
159 160 opObj10.addParameter(name='path',value=path_pp_save_int)
160 161 opObj10.addParameter(name='mode',value="weather")
161 162 opObj10.addParameter(name='type_data',value='F')
162 163 opObj10.addParameter(name='blocksPerFile',value='360',format='int')
163 164 #opObj10.addParameter(name='metadataList',value='utctimeInit,paramInterval,channelList,heightList,flagDataAsBlock',format='list')
164 165 opObj10.addParameter(name='metadataList',value='heightList,channelList,Typename,Datatype,Scantype,Latitude,Longitud,Heading,Height,Waveform,PRF,CreatedBy,ContactInformation',format='list')
165 166 #--------------------
166 167 opObj10.addParameter(name='Typename', value='Differential_Reflectivity')
167 168 opObj10.addParameter(name='Datatype', value='RadialSet')
168 169 opObj10.addParameter(name='Scantype', value='PPI')
169 170 opObj10.addParameter(name='Latitude', value='-11.96')
170 171 opObj10.addParameter(name='Longitud', value='-76.54')
171 172 opObj10.addParameter(name='Heading', value='293')
172 173 opObj10.addParameter(name='Height', value='293')
173 174 opObj10.addParameter(name='Waveform', value='OFM')
174 175 opObj10.addParameter(name='PRF', value='2000')
175 176 opObj10.addParameter(name='CreatedBy', value='WeatherRadarJROTeam')
176 177 opObj10.addParameter(name='ContactInformation', value='avaldez@igp.gob.pe')
177 178 #---------------------------------------------------
178 179 #opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,azimuth,elevation,utctime',format='list')#,format='list'
179 180 #opObj10.addParameter(name='metadataList',value='utctimeInit,timeZone,paramInterval,profileIndex,channelList,heightList,flagDataAsBlock',format='list')
180 181
181 182 opObj10.addParameter(name='dataList',value='dataPP_POW,azimuth,elevation,utctime',format='list')#,format='list'
182 183 opObj10.addParameter(name='description',value=json.dumps(desc_wr))
183 184
184 185 controllerObj.start()
Show file before | Show file after | Hide comments
@@ -1,122 +1,125
1 1 # Ing. AVP
2 2 # 04/01/2022
3 3 # ARCHIVO DE LECTURA
4 4 #---- DATA RHI --- 23 DE NOVIEMBRE DEL 2021 --- 23/11/2021---
5 5 #---- PEDESTAL ----------------------------------------------
6 6 #------- HORA 143826 /DATA_RM/TEST_PEDESTAL/P20211123-143826 14:38-15:10
7 7 #---- RADAR ----------------------------------------------
8 8 #------- 14:26-15:00
9 9 #------- /DATA_RM/DRONE/2MHZ_5V_ELEVACION/
10 10 #------- /DATA_RM/DRONE/2MHZ_5V_ELEVACION/ch0/2021-11-23T19-00-00
11 11
12 12 import os, sys
13 13 import datetime
14 14 import time
15 15 import numpy
16 16 from ext_met import getfirstFilefromPath,getDatavaluefromDirFilename
17 17 from schainpy.controller import Project
18 18 #-----------------------------------------------------------------------------------------
19 19 print("[SETUP]-RADAR METEOROLOGICO-")
20 20 path_ped = "/DATA_RM/TEST_PEDESTAL/P20211123-143826"
21 path_ped = "/DATA_RM/TEST_PEDESTAL/P20220217-172216"
21 22 print("PATH PEDESTAL :",path_ped)
22 23 path_adq = "/DATA_RM/DRONE/2MHZ_5V_ELEVACION/"
24 path_adq = "/DATA_RM/2MHZTEST/"
23 25 print("PATH DATA :",path_adq)
24 26 figpath_pp_rti = "/home/soporte/Pictures/TEST_PP_RHI"
25 27 print("PATH PP RTI :",figpath_pp_rti)
26 28 figpath_pp_rhi = "/home/soporte/Pictures/TEST_PP_RHI"
27 29 print("PATH PP RHI :",figpath_pp_rhi)
28 30 path_pp_save_int = "/DATA_RM/TEST_SAVE_PP_INT_RHI"
29 31 print("PATH SAVE PP INT :",path_pp_save_int)
30 32 print(" ")
31 33 #-------------------------------------------------------------------------------------------
32 34 print("SELECCIONAR MODO: PPI (0) O RHI (1)")
33 35 mode_wr = 1
34 36 if mode_wr==0:
35 37 print("[ ON ] MODE PPI")
36 38 list_ped = getfirstFilefromPath(path=path_ped,meta="PE",ext=".hdf5")
37 39 ff_pedestal = list_ped[2]
38 40 azi_vel = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="azi_vel")
39 41 V = round(azi_vel[0])
40 42 print("VELOCIDAD AZI :", int(numpy.mean(azi_vel)),"Β°/seg")
41 43 else:
42 44 print("[ ON ] MODE RHI")
43 45 list_ped = getfirstFilefromPath(path=path_ped,meta="PE",ext=".hdf5")
44 46 ff_pedestal = list_ped[2]
45 47 ele_vel = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="ele_vel")
46 48 V = round(ele_vel[0])
47 V = 10.0
49 V = 8.0#10.0
48 50 print("VELOCIDAD ELE :", int(numpy.mean(ele_vel)),"Β°/seg")
49 51 print(" ")
50 52 #---------------------------------------------------------------------------------------
51 53 print("SELECCIONAR MODO: PULSE PAIR (0) O FREQUENCY (1)")
52 54 mode_proc = 0
53 55 if mode_proc==0:
54 56 print("[ ON ] MODE PULSEPAIR")
55 57 else:
56 58 print("[ ON ] MODE FREQUENCY")
57 59 ipp = 60.0
58 60 print("IPP(Km.) : %1.2f"%ipp)
59 61 ipp_sec = (ipp*1.0e3/150.0)*1.0e-6
60 62 print("IPP(useg.) : %1.2f"%(ipp_sec*(1.0e6)))
61 63 VEL=V
62 64 n= int(1/(VEL*ipp_sec))
63 65 print("NΒ° Profiles : ", n)
64 66 #--------------------------------------------
65 67 plot_rti = 0
66 68 plot_rhi = 1
67 69 integration = 1
68 70 save = 0
69 71 #---------------------------RANGO DE PLOTEO----------------------------------
70 72 dBmin = '1'
71 73 dBmax = '85'
72 74 xmin = '17'
73 75 xmax = '17.25'
74 76 ymin = '0'
75 77 ymax = '600'
76 78 #----------------------------------------------------------------------------
77 79 time.sleep(3)
78 80 #---------------------SIGNAL CHAIN ------------------------------------
79 81 desc = "USRP_WEATHER_RADAR"
80 82 filename = "USRP_processing.xml"
81 83 controllerObj = Project()
82 84 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
83 85 #---------------------UNIDAD DE LECTURA--------------------------------
84 86 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
85 87 path=path_adq,
86 startDate="2021/11/23",#today,
87 endDate="2021/12/30",#today,
88 startTime='14:38:23',
88 startDate="2022/02/17",#today,
89 endDate="2022/02/17",#today,
90 startTime='00:00:00',
89 91 endTime='23:59:59',
90 92 delay=0,
91 93 #set=0,
92 94 online=0,
93 95 walk=1,
94 96 ippKm=ipp)
95 97
96 98 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc',inputId=readUnitConfObj.getId())
97 99
98 100 opObj11 = procUnitConfObjA.addOperation(name='selectHeights')
99 101 opObj11.addParameter(name='minIndex', value='1', format='int')
100 102 # opObj11.addParameter(name='maxIndex', value='10000', format='int')
101 103 opObj11.addParameter(name='maxIndex', value='400', format='int')
102 104
103 105 if mode_proc==0:
104 106 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
105 107 opObj11.addParameter(name='n', value=int(n), format='int')
106 108 procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
107 109
108 110 if integration==1:
109 111 opObj11 = procUnitConfObjB.addOperation(name='PedestalInformation')
110 112 opObj11.addParameter(name='path_ped', value=path_ped)
111 113 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
114 opObj11.addParameter(name='wr_exp', value='RHI')
112 115
113 116 if plot_rhi==1:
114 117 opObj11 = procUnitConfObjB.addOperation(name='Block360')
115 118 opObj11.addParameter(name='n', value='10', format='int')
116 119 opObj11.addParameter(name='mode', value=mode_proc, format='int')
117 120 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
118 121 opObj11= procUnitConfObjB.addOperation(name='WeatherRHIPlot',optype='other')
119 122 opObj11.addParameter(name='save', value=figpath_pp_rhi)
120 123 opObj11.addParameter(name='save_period', value=1)
121 124
122 125 controllerObj.start()
Show file before | Show file after | Hide comments
@@ -1,123 +1,126
1 1 # Ing-AlexanderValdez
2 2 # Monitoreo de Pedestal
3 3
4 4 ############## IMPORTA LIBRERIAS ###################
5 5 import os,numpy,h5py
6 6 import sys,time
7 7 import matplotlib.pyplot as plt
8 8 ####################################################
9 9 #################################################################
10 10 # LA FECHA 21-10-20 CORRESPONDE A LAS PRUEBAS DEL DIA MIERCOLES
11 11 # 1:15:51 pm hasta 3:49:32 pm
12 12 #################################################################
13 13
14 14 #path_ped = '/DATA_RM/TEST_PEDESTAL/P20211012-082745'
15 15 path_ped = '/DATA_RM/TEST_PEDESTAL/P20211020-131248'
16 16 path_ped = '/DATA_RM/TEST_PEDESTAL/P20211110-171003'
17 17 path_ped = '/DATA_RM/TEST_PEDESTAL/P20211111-173856'
18 18 path_ped = '/DATA_RM/TEST_PEDESTAL/P20211123-143826'
19 path_ped = "/DATA_RM/TEST_PEDESTAL/P20220217-172216"
19 20 #path_ped = '/DATA_RM/TEST_PEDESTAL/P20211111-173409'
20 21 # Metodo para verificar numero
21 22 def isNumber(str):
22 23 try:
23 24 float(str)
24 25 return True
25 26 except:
26 27 return False
27 28 # Metodo para extraer el arreglo
28 29 def getDatavaluefromDirFilename(path,file,value):
29 30 dir_file= path+"/"+file
30 31 fp = h5py.File(dir_file,'r')
31 32 array = fp['Data'].get(value)[()]
32 33 fp.close()
33 34 return array
34 35
35 36 # LISTA COMPLETA DE ARCHIVOS HDF5 Pedestal
36 37 LIST= sorted(os.listdir(path_ped))
37 38 m=len(LIST)
38 39 print("TOTAL DE ARCHIVOS DE PEDESTAL:",m)
39 40 # Contadores temporales
40 41 k= 0
41 42 l= 0
42 43 t= 0
43 44 # Marca de tiempo temporal
44 45 time_ = numpy.zeros([m])
45 46 # creacion de
46 47 for i in range(m):
47 48 print("order:",i)
48 49 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_pos")
49 50 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="ele_pos")
50 51 tmp_azi_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_vel")
51 tmp_ele_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_vel")# nuevo :D
52 tmp_ele_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="ele_vel")# nuevo :D
52 53
53 54 time_[i] = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="utc")
54 55
55 56 k=k +tmp_azi_pos.shape[0]
56 57 l=l +tmp_ele_pos.shape[0]
57 58 t=t +tmp_azi_vel.shape[0]
58 59
59 60 print("TOTAL DE MUESTRAS, ARCHIVOS X100:",k)
60 61 time.sleep(5)
61 62 ######CREACION DE ARREGLOS CANTIDAD DE VALORES POR MUESTRA#################
62 63 azi_pos = numpy.zeros([k])
63 64 ele_pos = numpy.zeros([l])
64 65 time_azi_pos= numpy.zeros([k])
65 66 # Contadores temporales
66 67 p=0
67 68 r=0
68 69 z=0
69 70 # VARIABLES TMP para almacenar azimuth, elevacion y tiempo
70 71
71 72 #for filename in sorted(os.listdir(path_ped)):
72 73 # CONDICION POR LEER EN TIEMPO REAL NO OFFLINE
73 74
74 75 for filename in LIST:
75 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="azi_pos")
76 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="ele_pos")
76 #tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="azi_pos")
77 #tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="ele_pos")
78 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="ele_vel")
79 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="azi_vel")
77 80 # CONDICION POR LEER EN TIEMPO REAL NO OFFLINE
78 81
79 82 if z==(m-1):
80 83 tmp_azi_time=numpy.arange(time_[z],time_[z]+1,1/(tmp_azi_pos.shape[0]))
81 84 else:
82 85 tmp_azi_time=numpy.arange(time_[z],time_[z+1],(time_[z+1]-time_[z])/(tmp_azi_pos.shape[0]))
83 86
84 87 print(filename,time_[z])
85 88 print(z,tmp_azi_pos.shape[0])
86 89
87 90 i=0
88 91 for i in range(tmp_azi_pos.shape[0]):
89 92 index=p+i
90 93 azi_pos[index]=tmp_azi_pos[i]
91 94 time_azi_pos[index]=tmp_azi_time[i]
92 95 p=p+tmp_azi_pos.shape[0]
93 96 i=0
94 97 for i in range(tmp_ele_pos.shape[0]):
95 98 index=r+i
96 99 ele_pos[index]=tmp_ele_pos[i]
97 100 r=r+tmp_ele_pos.shape[0]
98 101
99 102
100 103 z+=1
101 104
102 105
103 106 ######## GRAFIQUEMOS Y VEAMOS LOS DATOS DEL Pedestal
104 107 fig, ax = plt.subplots(figsize=(16,8))
105 108 print(time_azi_pos.shape)
106 109 print(azi_pos.shape)
107 110 t=numpy.arange(time_azi_pos.shape[0])*0.01/(60.0)
108 111 plt.plot(t,azi_pos,label='AZIMUTH_POS',color='blue')
109 112
110 113 # AQUI ESTOY ADICIONANDO LA POSICION EN elevaciont=numpy.arange(len(ele_pos))*0.01/60.0
111 114 t=numpy.arange(len(ele_pos))*0.01/60.0
112 115 plt.plot(t,ele_pos,label='ELEVATION_POS',color='red')#*10
113 116
114 117 ax.set_xlim(0, 4)
115 118 ax.set_ylim(-5, 50)
116 119 plt.ylabel("Azimuth Position")
117 120 plt.xlabel("Muestra")
118 121 plt.title('Azimuth Position vs Muestra ', fontsize=20)
119 122 axes = plt.gca()
120 123 axes.yaxis.grid()
121 124 plt.xticks(fontsize=16)
122 125 plt.yticks(fontsize=16)
123 126 plt.show()
General Comments 0
You need to be logged in to leave comments. Login now