##// END OF EJS Templates
schain
RSS
jroplot_parameters.py jroplot_spectra.py jroIO_digitalRF.py jroproc_parameters.py sophy_proc.py
eynilupu -
r1653:f661d02a9309
parent child
Show More
Side by Side Unified
Context file:
r1653:f661d02a9309 Loading diff...:

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

Show file before | Show file after | Hide comments
@@ -1,725 +1,736
1 1 import os
2 2 import datetime
3 3 import warnings
4 4 import numpy
5 5 from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
6 6 from matplotlib.patches import Circle
7 7 import cartopy.crs as ccrs
8 8 from cartopy.feature import ShapelyFeature
9 9 import cartopy.io.shapereader as shpreader
10 10
11 11 from schainpy.model.graphics.jroplot_base import Plot, plt
12 12 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
13 13 from schainpy.utils import log
14 14 from schainpy.model.graphics.plotting_codes import cb_tables
15 15
16 16
17 17 EARTH_RADIUS = 6.3710e3
18 18
19 19
20 20 def antenna_to_cartesian(ranges, azimuths, elevations):
21 21 """
22 22 Return Cartesian coordinates from antenna coordinates.
23 23
24 24 Parameters
25 25 ----------
26 26 ranges : array
27 27 Distances to the center of the radar gates (bins) in kilometers.
28 28 azimuths : array
29 29 Azimuth angle of the radar in degrees.
30 30 elevations : array
31 31 Elevation angle of the radar in degrees.
32 32
33 33 Returns
34 34 -------
35 35 x, y, z : array
36 36 Cartesian coordinates in meters from the radar.
37 37
38 38 Notes
39 39 -----
40 40 The calculation for Cartesian coordinate is adapted from equations
41 41 2.28(b) and 2.28(c) of Doviak and Zrnic [1]_ assuming a
42 42 standard atmosphere (4/3 Earth's radius model).
43 43
44 44 .. math::
45 45
46 46 z = \\sqrt{r^2+R^2+2*r*R*sin(\\theta_e)} - R
47 47
48 48 s = R * arcsin(\\frac{r*cos(\\theta_e)}{R+z})
49 49
50 50 x = s * sin(\\theta_a)
51 51
52 52 y = s * cos(\\theta_a)
53 53
54 54 Where r is the distance from the radar to the center of the gate,
55 55 :math:`\\theta_a` is the azimuth angle, :math:`\\theta_e` is the
56 56 elevation angle, s is the arc length, and R is the effective radius
57 57 of the earth, taken to be 4/3 the mean radius of earth (6371 km).
58 58
59 59 References
60 60 ----------
61 61 .. [1] Doviak and Zrnic, Doppler Radar and Weather Observations, Second
62 62 Edition, 1993, p. 21.
63 63
64 64 """
65 65 theta_e = numpy.deg2rad(elevations) # elevation angle in radians.
66 66 theta_a = numpy.deg2rad(azimuths) # azimuth angle in radians.
67 67 R = 6371.0 * 1000.0 * 4.0 / 3.0 # effective radius of earth in meters.
68 68 r = ranges * 1000.0 # distances to gates in meters.
69 69
70 70 z = (r ** 2 + R ** 2 + 2.0 * r * R * numpy.sin(theta_e)) ** 0.5 - R
71 71 s = R * numpy.arcsin(r * numpy.cos(theta_e) / (R + z)) # arc length in m.
72 72 x = s * numpy.sin(theta_a)
73 73 y = s * numpy.cos(theta_a)
74 74 return x, y, z
75 75
76 76 def cartesian_to_geographic_aeqd(x, y, lon_0, lat_0, R=EARTH_RADIUS):
77 77 """
78 78 Azimuthal equidistant Cartesian to geographic coordinate transform.
79 79
80 80 Transform a set of Cartesian/Cartographic coordinates (x, y) to
81 81 geographic coordinate system (lat, lon) using a azimuthal equidistant
82 82 map projection [1]_.
83 83
84 84 .. math::
85 85
86 86 lat = \\arcsin(\\cos(c) * \\sin(lat_0) +
87 87 (y * \\sin(c) * \\cos(lat_0) / \\rho))
88 88
89 89 lon = lon_0 + \\arctan2(
90 90 x * \\sin(c),
91 91 \\rho * \\cos(lat_0) * \\cos(c) - y * \\sin(lat_0) * \\sin(c))
92 92
93 93 \\rho = \\sqrt(x^2 + y^2)
94 94
95 95 c = \\rho / R
96 96
97 97 Where x, y are the Cartesian position from the center of projection;
98 98 lat, lon the corresponding latitude and longitude; lat_0, lon_0 are the
99 99 latitude and longitude of the center of the projection; R is the radius of
100 100 the earth (defaults to ~6371 km). lon is adjusted to be between -180 and
101 101 180.
102 102
103 103 Parameters
104 104 ----------
105 105 x, y : array-like
106 106 Cartesian coordinates in the same units as R, typically meters.
107 107 lon_0, lat_0 : float
108 108 Longitude and latitude, in degrees, of the center of the projection.
109 109 R : float, optional
110 110 Earth radius in the same units as x and y. The default value is in
111 111 units of meters.
112 112
113 113 Returns
114 114 -------
115 115 lon, lat : array
116 116 Longitude and latitude of Cartesian coordinates in degrees.
117 117
118 118 References
119 119 ----------
120 120 .. [1] Snyder, J. P. Map Projections--A Working Manual. U. S. Geological
121 121 Survey Professional Paper 1395, 1987, pp. 191-202.
122 122
123 123 """
124 124 x = numpy.atleast_1d(numpy.asarray(x))
125 125 y = numpy.atleast_1d(numpy.asarray(y))
126 126
127 127 lat_0_rad = numpy.deg2rad(lat_0)
128 128 lon_0_rad = numpy.deg2rad(lon_0)
129 129
130 130 rho = numpy.sqrt(x*x + y*y)
131 131 c = rho / R
132 132
133 133 with warnings.catch_warnings():
134 134 # division by zero may occur here but is properly addressed below so
135 135 # the warnings can be ignored
136 136 warnings.simplefilter("ignore", RuntimeWarning)
137 137 lat_rad = numpy.arcsin(numpy.cos(c) * numpy.sin(lat_0_rad) +
138 138 y * numpy.sin(c) * numpy.cos(lat_0_rad) / rho)
139 139 lat_deg = numpy.rad2deg(lat_rad)
140 140 # fix cases where the distance from the center of the projection is zero
141 141 lat_deg[rho == 0] = lat_0
142 142
143 143 x1 = x * numpy.sin(c)
144 144 x2 = rho*numpy.cos(lat_0_rad)*numpy.cos(c) - y*numpy.sin(lat_0_rad)*numpy.sin(c)
145 145 lon_rad = lon_0_rad + numpy.arctan2(x1, x2)
146 146 lon_deg = numpy.rad2deg(lon_rad)
147 147 # Longitudes should be from -180 to 180 degrees
148 148 lon_deg[lon_deg > 180] -= 360.
149 149 lon_deg[lon_deg < -180] += 360.
150 150
151 151 return lon_deg, lat_deg
152 152
153 153 def antenna_to_geographic(ranges, azimuths, elevations, site):
154 154
155 155 x, y, z = antenna_to_cartesian(numpy.array(ranges), numpy.array(azimuths), numpy.array(elevations))
156 156 lon, lat = cartesian_to_geographic_aeqd(x, y, site[0], site[1], R=6370997.)
157 157
158 158 return lon, lat
159 159
160 160 def ll2xy(lat1, lon1, lat2, lon2):
161 161
162 162 p = 0.017453292519943295
163 163 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
164 164 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
165 165 r = 12742 * numpy.arcsin(numpy.sqrt(a))
166 166 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
167 167 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
168 168 theta = -theta + numpy.pi/2
169 169 return r*numpy.cos(theta), r*numpy.sin(theta)
170 170
171 171
172 172 def km2deg(km):
173 173 '''
174 174 Convert distance in km to degrees
175 175 '''
176 176
177 177 return numpy.rad2deg(km/EARTH_RADIUS)
178 178
179 179
180 180
181 181 class SpectralMomentsPlot(SpectraPlot):
182 182 '''
183 183 Plot for Spectral Moments
184 184 '''
185 185 CODE = 'spc_moments'
186 186 # colormap = 'jet'
187 187 # plot_type = 'pcolor'
188 188
189 189 class DobleGaussianPlot(SpectraPlot):
190 190 '''
191 191 Plot for Double Gaussian Plot
192 192 '''
193 193 CODE = 'gaussian_fit'
194 194 # colormap = 'jet'
195 195 # plot_type = 'pcolor'
196 196
197 197 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
198 198 '''
199 199 Plot SpectraCut with Double Gaussian Fit
200 200 '''
201 201 CODE = 'cut_gaussian_fit'
202 202
203 203 class SnrPlot(RTIPlot):
204 204 '''
205 205 Plot for SNR Data
206 206 '''
207 207
208 208 CODE = 'snr'
209 209 colormap = 'jet'
210 210
211 211 def update(self, dataOut):
212 212
213 213 data = {
214 214 'snr': 10*numpy.log10(dataOut.data_snr)
215 215 }
216 216
217 217 return data, {}
218 218
219 219 class DopplerPlot(RTIPlot):
220 220 '''
221 221 Plot for DOPPLER Data (1st moment)
222 222 '''
223 223
224 224 CODE = 'dop'
225 225 colormap = 'jet'
226 226
227 227 def update(self, dataOut):
228 228
229 229 data = {
230 230 'dop': 10*numpy.log10(dataOut.data_dop)
231 231 }
232 232
233 233 return data, {}
234 234
235 235 class PowerPlot(RTIPlot):
236 236 '''
237 237 Plot for Power Data (0 moment)
238 238 '''
239 239
240 240 CODE = 'pow'
241 241 colormap = 'jet'
242 242
243 243 def update(self, dataOut):
244 244 data = {
245 245 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
246 246 }
247 247 return data, {}
248 248
249 249 class SpectralWidthPlot(RTIPlot):
250 250 '''
251 251 Plot for Spectral Width Data (2nd moment)
252 252 '''
253 253
254 254 CODE = 'width'
255 255 colormap = 'jet'
256 256
257 257 def update(self, dataOut):
258 258
259 259 data = {
260 260 'width': dataOut.data_width
261 261 }
262 262
263 263 return data, {}
264 264
265 265 class SkyMapPlot(Plot):
266 266 '''
267 267 Plot for meteors detection data
268 268 '''
269 269
270 270 CODE = 'param'
271 271
272 272 def setup(self):
273 273
274 274 self.ncols = 1
275 275 self.nrows = 1
276 276 self.width = 7.2
277 277 self.height = 7.2
278 278 self.nplots = 1
279 279 self.xlabel = 'Zonal Zenith Angle (deg)'
280 280 self.ylabel = 'Meridional Zenith Angle (deg)'
281 281 self.polar = True
282 282 self.ymin = -180
283 283 self.ymax = 180
284 284 self.colorbar = False
285 285
286 286 def plot(self):
287 287
288 288 arrayParameters = numpy.concatenate(self.data['param'])
289 289 error = arrayParameters[:, -1]
290 290 indValid = numpy.where(error == 0)[0]
291 291 finalMeteor = arrayParameters[indValid, :]
292 292 finalAzimuth = finalMeteor[:, 3]
293 293 finalZenith = finalMeteor[:, 4]
294 294
295 295 x = finalAzimuth * numpy.pi / 180
296 296 y = finalZenith
297 297
298 298 ax = self.axes[0]
299 299
300 300 if ax.firsttime:
301 301 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
302 302 else:
303 303 ax.plot.set_data(x, y)
304 304
305 305 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
306 306 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
307 307 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
308 308 dt2,
309 309 len(x))
310 310 self.titles[0] = title
311 311
312 312
313 313 class GenericRTIPlot(Plot):
314 314 '''
315 315 Plot for data_xxxx object
316 316 '''
317 317
318 318 CODE = 'param'
319 319 colormap = 'viridis'
320 320 plot_type = 'pcolorbuffer'
321 321
322 322 def setup(self):
323 323 self.xaxis = 'time'
324 324 self.ncols = 1
325 325 self.nrows = self.data.shape('param')[0]
326 326 self.nplots = self.nrows
327 327 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
328 328
329 329 if not self.xlabel:
330 330 self.xlabel = 'Time'
331 331
332 332 self.ylabel = 'Range [km]'
333 333 if not self.titles:
334 334 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
335 335
336 336 def update(self, dataOut):
337 337
338 338 data = {
339 339 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
340 340 }
341 341
342 342 meta = {}
343 343
344 344 return data, meta
345 345
346 346 def plot(self):
347 347 # self.data.normalize_heights()
348 348 self.x = self.data.times
349 349 self.y = self.data.yrange
350 350 self.z = self.data['param']
351 351 self.z = 10*numpy.log10(self.z)
352 352 self.z = numpy.ma.masked_invalid(self.z)
353 353
354 354 if self.decimation is None:
355 355 x, y, z = self.fill_gaps(self.x, self.y, self.z)
356 356 else:
357 357 x, y, z = self.fill_gaps(*self.decimate())
358 358
359 359 for n, ax in enumerate(self.axes):
360 360
361 361 self.zmax = self.zmax if self.zmax is not None else numpy.max(
362 362 self.z[n])
363 363 self.zmin = self.zmin if self.zmin is not None else numpy.min(
364 364 self.z[n])
365 365
366 366 if ax.firsttime:
367 367 if self.zlimits is not None:
368 368 self.zmin, self.zmax = self.zlimits[n]
369 369
370 370 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
371 371 vmin=self.zmin,
372 372 vmax=self.zmax,
373 373 cmap=self.cmaps[n]
374 374 )
375 375 else:
376 376 if self.zlimits is not None:
377 377 self.zmin, self.zmax = self.zlimits[n]
378 378 ax.collections.remove(ax.collections[0])
379 379 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
380 380 vmin=self.zmin,
381 381 vmax=self.zmax,
382 382 cmap=self.cmaps[n]
383 383 )
384 384
385 385
386 386 class PolarMapPlot(Plot):
387 387 '''
388 388 Plot for weather radar
389 389 '''
390 390
391 391 CODE = 'param'
392 392 colormap = 'seismic'
393 393
394 394 def setup(self):
395 395 self.ncols = 1
396 396 self.nrows = 1
397 397 self.width = 9
398 398 self.height = 8
399 399 self.mode = self.data.meta['mode']
400 400 if self.channels is not None:
401 401 self.nplots = len(self.channels)
402 402 self.nrows = len(self.channels)
403 403 else:
404 404 self.nplots = self.data.shape(self.CODE)[0]
405 405 self.nrows = self.nplots
406 406 self.channels = list(range(self.nplots))
407 407 if self.mode == 'E':
408 408 self.xlabel = 'Longitude'
409 409 self.ylabel = 'Latitude'
410 410 else:
411 411 self.xlabel = 'Range (km)'
412 412 self.ylabel = 'Height (km)'
413 413 self.bgcolor = 'white'
414 414 self.cb_labels = self.data.meta['units']
415 415 self.lat = self.data.meta['latitude']
416 416 self.lon = self.data.meta['longitude']
417 417 self.xmin, self.xmax = float(
418 418 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
419 419 self.ymin, self.ymax = float(
420 420 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
421 421 # self.polar = True
422 422
423 423 def plot(self):
424 424
425 425 for n, ax in enumerate(self.axes):
426 426 data = self.data['param'][self.channels[n]]
427 427
428 428 zeniths = numpy.linspace(
429 429 0, self.data.meta['max_range'], data.shape[1])
430 430 if self.mode == 'E':
431 431 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
432 432 r, theta = numpy.meshgrid(zeniths, azimuths)
433 433 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
434 434 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
435 435 x = km2deg(x) + self.lon
436 436 y = km2deg(y) + self.lat
437 437 else:
438 438 azimuths = numpy.radians(self.data.yrange)
439 439 r, theta = numpy.meshgrid(zeniths, azimuths)
440 440 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
441 441 self.y = zeniths
442 442
443 443 if ax.firsttime:
444 444 if self.zlimits is not None:
445 445 self.zmin, self.zmax = self.zlimits[n]
446 446 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
447 447 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
448 448 vmin=self.zmin,
449 449 vmax=self.zmax,
450 450 cmap=self.cmaps[n])
451 451 else:
452 452 if self.zlimits is not None:
453 453 self.zmin, self.zmax = self.zlimits[n]
454 454 ax.collections.remove(ax.collections[0])
455 455 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
456 456 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
457 457 vmin=self.zmin,
458 458 vmax=self.zmax,
459 459 cmap=self.cmaps[n])
460 460
461 461 if self.mode == 'A':
462 462 continue
463 463
464 464 # plot district names
465 465 f = open('/data/workspace/schain_scripts/distrito.csv')
466 466 for line in f:
467 467 label, lon, lat = [s.strip() for s in line.split(',') if s]
468 468 lat = float(lat)
469 469 lon = float(lon)
470 470 # ax.plot(lon, lat, '.b', ms=2)
471 471 ax.text(lon, lat, label.decode('utf8'), ha='center',
472 472 va='bottom', size='8', color='black')
473 473
474 474 # plot limites
475 475 limites = []
476 476 tmp = []
477 477 for line in open('/data/workspace/schain_scripts/lima.csv'):
478 478 if '#' in line:
479 479 if tmp:
480 480 limites.append(tmp)
481 481 tmp = []
482 482 continue
483 483 values = line.strip().split(',')
484 484 tmp.append((float(values[0]), float(values[1])))
485 485 for points in limites:
486 486 ax.add_patch(
487 487 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
488 488
489 489 # plot Cuencas
490 490 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
491 491 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
492 492 values = [line.strip().split(',') for line in f]
493 493 points = [(float(s[0]), float(s[1])) for s in values]
494 494 ax.add_patch(Polygon(points, ec='b', fc='none'))
495 495
496 496 # plot grid
497 497 for r in (15, 30, 45, 60):
498 498 ax.add_artist(plt.Circle((self.lon, self.lat),
499 499 km2deg(r), color='0.6', fill=False, lw=0.2))
500 500 ax.text(
501 501 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
502 502 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
503 503 '{}km'.format(r),
504 504 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
505 505
506 506 if self.mode == 'E':
507 507 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
508 508 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
509 509 else:
510 510 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
511 511 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
512 512
513 513 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
514 514 self.titles = ['{} {}'.format(
515 515 self.data.parameters[x], title) for x in self.channels]
516 516
517 517 class WeatherParamsPlot(Plot):
518 518 #CODE = 'RHI'
519 519 #plot_name = 'RHI'
520 520 plot_type = 'scattermap'
521 521 buffering = False
522 522 projection = ccrs.PlateCarree()
523 523
524 524 def setup(self):
525 525
526 526 self.ncols = 1
527 527 self.nrows = 1
528 528 self.nplots= 1
529 529 self.ylabel= 'Height [km]'
530 530 self.xlabel= 'Distance from radar [km]'
531
531
532 532 if self.channels is not None:
533 533 self.nplots = len(self.channels)
534 534 self.ncols = len(self.channels)
535 535 else:
536 536 self.nplots = self.data.shape(self.CODE)[0]
537 537 self.ncols = self.nplots
538 538 self.channels = list(range(self.nplots))
539 539
540 540 self.colorbar=True
541 541 if len(self.channels)>1:
542 542 self.width = 12
543 543 else:
544 544 self.width =8
545 545 self.height =7
546 546 self.ini =0
547 547 self.len_azi =0
548 548 self.buffer_ini = None
549 549 self.buffer_ele = None
550 550 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.1})
551 551 self.flag =0
552 552 self.indicador= 0
553 553 self.last_data_ele = None
554 554 self.val_mean = None
555 555
556 556 def update(self, dataOut):
557 557
558 558 vars = {
559 559 'S' : 0,
560 560 'V' : 1,
561 561 'W' : 2,
562 562 'SNR' : 3,
563 563 'Z' : 4,
564 564 'D' : 5,
565 565 'P' : 6,
566 566 'R' : 7,
567 567 }
568 568
569 569 data = {}
570 570 meta = {}
571 571
572 572 if hasattr(dataOut, 'nFFTPoints'):
573 573 factor = dataOut.normFactor
574 574 else:
575 575 factor = 1
576 576
577 577 if hasattr(dataOut, 'dparam'):
578 578 tmp = getattr(dataOut, 'data_param')
579 579 else:
580
580 #print("-------------------self.attr_data[0]",self.attr_data[0])
581 581 if 'S' in self.attr_data[0]:
582 tmp = 10*numpy.log10(10.0*getattr(dataOut, 'data_param')[:,0,:]/(factor))
582 if self.attr_data[0]=='S':
583 tmp = 10*numpy.log10(10.0*getattr(dataOut, 'data_param')[:,0,:]/(factor))
584 if self.attr_data[0]=='SNR':
585 tmp = 10*numpy.log10(getattr(dataOut, 'data_param')[:,3,:])
583 586 else:
584 587 tmp = getattr(dataOut, 'data_param')[:,vars[self.attr_data[0]],:]
585 588
586 589 if self.mask:
587 590 mask = dataOut.data_param[:,3,:] < self.mask
588 591 tmp = numpy.ma.masked_array(tmp, mask=mask)
589 592
590 593 r = dataOut.heightList
591 594 delta_height = r[1]-r[0]
592 595 valid = numpy.where(r>=0)[0]
593 596 data['r'] = numpy.arange(len(valid))*delta_height
594 597
595 598 data['data'] = [0, 0]
596 599
597 600 try:
598 601 data['data'][0] = tmp[0][:,valid]
599 602 data['data'][1] = tmp[1][:,valid]
600 603 except:
601 604 data['data'][0] = tmp[0][:,valid]
602 605 data['data'][1] = tmp[0][:,valid]
603 606
604 607 if dataOut.mode_op == 'PPI':
605 608 self.CODE = 'PPI'
606 609 self.title = self.CODE
607 610 elif dataOut.mode_op == 'RHI':
608 611 self.CODE = 'RHI'
609 612 self.title = self.CODE
610 613
611 614 data['azi'] = dataOut.data_azi
612 615 data['ele'] = dataOut.data_ele
616
617 if isinstance(dataOut.mode_op, bytes):
618 try:
619 dataOut.mode_op = dataOut.mode_op.decode()
620 except:
621 dataOut.mode_op = str(dataOut.mode_op, 'utf-8')
613 622 data['mode_op'] = dataOut.mode_op
614 623 self.mode = dataOut.mode_op
615 624
616 625 return data, meta
617 626
618 627 def plot(self):
619 628 data = self.data[-1]
620 629 z = data['data']
621 630 r = data['r']
622 631 self.titles = []
623 632
624 633 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
625 634 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
626 635 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
627 636 self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z)
628 637
629 638 if isinstance(data['mode_op'], bytes):
630 639 data['mode_op'] = data['mode_op'].decode()
631 640
632 641 if data['mode_op'] == 'RHI':
633 642 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']))
634 643 len_aux = int(data['azi'].shape[0]/4)
635 644 mean = numpy.mean(data['azi'][len_aux:-len_aux])
636 645 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
637 elif data['mode_op'] == 'PPI':
646 elif data['mode_op'] == 'PPI':
638 647 r, theta = numpy.meshgrid(r, -numpy.radians(data['azi'])+numpy.pi/2)
639 648 len_aux = int(data['ele'].shape[0]/4)
640 649 mean = numpy.mean(data['ele'][len_aux:-len_aux])
641 650 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(mean)), r*numpy.sin(
642 651 theta)*numpy.cos(numpy.radians(mean))
643 652 x = km2deg(x) + -75.295893
644 653 y = km2deg(y) + -12.040436
645 654
646 655 self.clear_figures()
647 656
648 657 if data['mode_op'] == 'PPI':
649 658 axes = self.axes['PPI']
650 659 else:
651 660 axes = self.axes['RHI']
652 661
653 662 if self.colormap in cb_tables:
654 663 norm = cb_tables[self.colormap]['norm']
655 664 else:
656 665 norm = None
657
666
658 667 for i, ax in enumerate(axes):
659 668 if data['mode_op'] == 'PPI':
660 669 ax.set_extent([-75.745893, -74.845893, -12.490436, -11.590436])
661 670
662 671 if norm is None:
663 672 ax.plt = ax.pcolormesh(x, y, z[i], cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
664 673 else:
665 674 ax.plt = ax.pcolormesh(x, y, z[i], cmap=self.colormap, norm=norm)
666 675
667 676 if data['mode_op'] == 'RHI':
668 677 len_aux = int(data['azi'].shape[0]/4)
669 678 mean = numpy.mean(data['azi'][len_aux:-len_aux])
670 679 if len(self.channels) !=1:
671 680 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels]
672 681 else:
673 682 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
674 683 elif data['mode_op'] == 'PPI':
675 684 len_aux = int(data['ele'].shape[0]/4)
676 685 mean = numpy.mean(data['ele'][len_aux:-len_aux])
677 686 if len(self.channels) !=1:
678 687 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels]
679 688 else:
680 689 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
681 690 self.mode_value = round(mean,1)
682 691
683 692 if data['mode_op'] == 'PPI':
684 693 gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
685 694 linewidth=1, color='gray', alpha=0.5, linestyle='--')
686 695 gl.xlabel_style = {'size': 8}
687 696 gl.ylabel_style = {'size': 8}
688 697 gl.xlabels_top = False
689 698 gl.ylabels_right = False
699 #self.shapes="/home/soporte/workspace/sirm/volumes/schain/shapes/"
700 #print("self.shapes",self.shapes)
690 701 shape_p = os.path.join(self.shapes,'PER_ADM2/PER_ADM2.shp')
691 702 shape_d = os.path.join(self.shapes,'PER_ADM1/PER_ADM1.shp')
692 703 capitales = os.path.join(self.shapes,'CAPITALES/cap_provincia.shp')
693 704 vias = os.path.join(self.shapes,'Carreteras/VIAS_NACIONAL_250000.shp')
694 705 reader_d = shpreader.BasicReader(shape_p, encoding='latin1')
695 706 reader_p = shpreader.BasicReader(shape_d, encoding='latin1')
696 707 reader_c = shpreader.BasicReader(capitales, encoding='latin1')
697 708 reader_v = shpreader.BasicReader(vias, encoding='latin1')
698 caps = [x for x in reader_c.records() if x.attributes["Departa"] in ("JUNIN", "LIMA", "AYACUCHO", "HUANCAVELICA")]
709 caps = [x for x in reader_c.records() if x.attributes["Departa"] in ("JUNIN", "LIMA", "AYACUCHO", "HUANCAVELICA")]
699 710 districts = [x for x in reader_d.records() if x.attributes["Name"] in ("JUNÍN", "CHANCHAMAYO", "CHUPACA", "CONCEPCIÓN", "HUANCAYO", "JAUJA", "SATIPO", "TARMA", "YAUYOS", "HUAROCHIRÍ", "CANTA", "HUANTA", "TAYACAJA")]
700 711 provs = [x for x in reader_p.records() if x.attributes["NAME"] in ("Junín", "Lima")]
701 712 vias = [x for x in reader_v.records() if x.attributes["DEP"] in ("JUNIN", "LIMA")]
702 713
703 714 # Display limits and streets
704 715 shape_feature = ShapelyFeature([x.geometry for x in districts], ccrs.PlateCarree(), facecolor="none", edgecolor='grey', lw=0.5)
705 716 ax.add_feature(shape_feature)
706 717 shape_feature = ShapelyFeature([x.geometry for x in provs], ccrs.PlateCarree(), facecolor="none", edgecolor='white', lw=1)
707 718 ax.add_feature(shape_feature)
708 719 shape_feature = ShapelyFeature([x.geometry for x in vias], ccrs.PlateCarree(), facecolor="none", edgecolor='yellow', lw=1)
709 720 ax.add_feature(shape_feature)
710 721
711 722 for cap in caps:
712 723 if cap.attributes['Nombre'] in ("LA OROYA", "CONCEPCIÓN", "HUANCAYO", "JAUJA", "CHUPACA", "YAUYOS", "HUANTA", "PAMPAS"):
713 724 ax.text(cap.attributes['X'], cap.attributes['Y'], cap.attributes['Nombre'].title(), size=7, color='white')
714 725 ax.text(-75.052003, -11.915552, 'Huaytapallana', size=7, color='cyan')
715 726 ax.plot(-75.052003, -11.915552, '*')
716
727
717 728 for R in (10, 20, 30 , 40, 50):
718 729 circle = Circle((-75.295893, -12.040436), km2deg(R), facecolor='none',
719 730 edgecolor='skyblue', linewidth=1, alpha=0.5)
720 731 ax.add_patch(circle)
721 ax.text(km2deg(R)*numpy.cos(numpy.radians(45))-75.295893,
722 km2deg(R)*numpy.sin(numpy.radians(45))-12.040436,
732 ax.text(km2deg(R)*numpy.cos(numpy.radians(45))-75.295893,
733 km2deg(R)*numpy.sin(numpy.radians(45))-12.040436,
723 734 '{}km'.format(R), color='skyblue', size=7)
724 735 elif data['mode_op'] == 'RHI':
725 736 ax.grid(color='grey', alpha=0.5, linestyle='--', linewidth=1)
Show file before | Show file after | Hide comments
@@ -1,743 +1,743
1 1 # Copyright (c) 2012-2021 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13
14 14
15 15 class SpectraPlot(Plot):
16 16 '''
17 17 Plot for Spectra data
18 18 '''
19 19
20 20 CODE = 'spc'
21 21 colormap = 'jet'
22 22 plot_type = 'pcolor'
23 23 buffering = False
24 24
25 25 def setup(self):
26 26 self.nplots = len(self.data.channels)
27 27 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
28 28 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
29 29 self.height = 2.6 * self.nrows
30 30 self.cb_label = 'dB'
31 31 if self.showprofile:
32 32 self.width = 4 * self.ncols
33 33 else:
34 34 self.width = 3.5 * self.ncols
35 35 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
36 36 self.ylabel = 'Range [km]'
37 37
38 38 def update(self, dataOut):
39 39
40 40 data = {}
41 41 meta = {}
42 42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 43 data['spc'] = spc
44 44 data['rti'] = dataOut.getPower()
45 45 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
46 meta['xrange'] = (dataOut.getFreqRange(0)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(0))
47 47
48 48 if self.CODE == 'spc_moments':
49 49 data['moments'] = dataOut.moments
50 50 # data['spc'] = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
51 51 if self.CODE == 'gaussian_fit':
52 52 # data['moments'] = dataOut.moments
53 53 data['gaussfit'] = dataOut.DGauFitParams
54 54 # data['spc'] = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
55 55
56 56 return data, meta
57 57
58 58 def plot(self):
59 59 if self.xaxis == "frequency":
60 60 x = self.data.xrange[0]
61 61 self.xlabel = "Frequency (kHz)"
62 62 elif self.xaxis == "time":
63 63 x = self.data.xrange[1]
64 64 self.xlabel = "Time (ms)"
65 65 else:
66 66 x = self.data.xrange[2]
67 67 self.xlabel = "Velocity (m/s)"
68 68
69 69 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
70 70 x = self.data.xrange[2]
71 71 self.xlabel = "Velocity (m/s)"
72 72
73 73 self.titles = []
74 74
75 75 y = self.data.yrange
76 76 self.y = y
77 77
78 78 data = self.data[-1]
79 79 z = data['spc']
80 80
81 81 for n, ax in enumerate(self.axes):
82 82 noise = data['noise'][n]
83 83 if self.CODE == 'spc_moments':
84 84 mean = data['moments'][n, 1]
85 85 if self.CODE == 'gaussian_fit':
86 86 # mean = data['moments'][n, 1]
87 87 gau0 = data['gaussfit'][n][2,:,0]
88 88 gau1 = data['gaussfit'][n][2,:,1]
89 89 if ax.firsttime:
90 90 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
91 91 self.xmin = self.xmin if self.xmin else -self.xmax
92 92 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
93 93 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
94 94 ax.plt = ax.pcolormesh(x, y, z[n].T,
95 95 vmin=self.zmin,
96 96 vmax=self.zmax,
97 97 cmap=plt.get_cmap(self.colormap)
98 98 )
99 99
100 100 if self.showprofile:
101 101 ax.plt_profile = self.pf_axes[n].plot(
102 102 data['rti'][n], y)[0]
103 103 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
104 104 color="k", linestyle="dashed", lw=1)[0]
105 105 if self.CODE == 'spc_moments':
106 106 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
107 107 if self.CODE == 'gaussian_fit':
108 108 # ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
109 109 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
110 110 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
111 111 else:
112 112 ax.plt.set_array(z[n].T.ravel())
113 113 if self.showprofile:
114 114 ax.plt_profile.set_data(data['rti'][n], y)
115 115 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
116 116 if self.CODE == 'spc_moments':
117 117 ax.plt_mean.set_data(mean, y)
118 118 if self.CODE == 'gaussian_fit':
119 119 # ax.plt_mean.set_data(mean, y)
120 120 ax.plt_gau0.set_data(gau0, y)
121 121 ax.plt_gau1.set_data(gau1, y)
122 122 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
123 123
124 124
125 125 class CrossSpectraPlot(Plot):
126 126
127 127 CODE = 'cspc'
128 128 colormap = 'jet'
129 129 plot_type = 'pcolor'
130 130 zmin_coh = None
131 131 zmax_coh = None
132 132 zmin_phase = None
133 133 zmax_phase = None
134 134
135 135 def setup(self):
136 136
137 137 self.ncols = 4
138 138 self.nplots = len(self.data.pairs) * 2
139 139 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
140 140 self.width = 3.1 * self.ncols
141 141 self.height = 2.6 * self.nrows
142 142 self.ylabel = 'Range [km]'
143 143 self.showprofile = False
144 144 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
145 145
146 146 def update(self, dataOut):
147 147
148 148 data = {}
149 149 meta = {}
150 150
151 151 spc = dataOut.data_spc
152 152 cspc = dataOut.data_cspc
153 153 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
154 154 meta['pairs'] = dataOut.pairsList
155 155
156 156 tmp = []
157 157
158 158 for n, pair in enumerate(meta['pairs']):
159 159 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
160 160 coh = numpy.abs(out)
161 161 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
162 162 tmp.append(coh)
163 163 tmp.append(phase)
164 164
165 165 data['cspc'] = numpy.array(tmp)
166 166
167 167 return data, meta
168 168
169 169 def plot(self):
170 170
171 171 if self.xaxis == "frequency":
172 172 x = self.data.xrange[0]
173 173 self.xlabel = "Frequency (kHz)"
174 174 elif self.xaxis == "time":
175 175 x = self.data.xrange[1]
176 176 self.xlabel = "Time (ms)"
177 177 else:
178 178 x = self.data.xrange[2]
179 179 self.xlabel = "Velocity (m/s)"
180 180
181 181 self.titles = []
182 182
183 183 y = self.data.yrange
184 184 self.y = y
185 185
186 186 data = self.data[-1]
187 187 cspc = data['cspc']
188 188
189 189 for n in range(len(self.data.pairs)):
190 190 pair = self.data.pairs[n]
191 191 coh = cspc[n*2]
192 192 phase = cspc[n*2+1]
193 193 ax = self.axes[2 * n]
194 194 if ax.firsttime:
195 195 ax.plt = ax.pcolormesh(x, y, coh.T,
196 196 vmin=0,
197 197 vmax=1,
198 198 cmap=plt.get_cmap(self.colormap_coh)
199 199 )
200 200 else:
201 201 ax.plt.set_array(coh.T.ravel())
202 202 self.titles.append(
203 203 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
204 204
205 205 ax = self.axes[2 * n + 1]
206 206 if ax.firsttime:
207 207 ax.plt = ax.pcolormesh(x, y, phase.T,
208 208 vmin=-180,
209 209 vmax=180,
210 210 cmap=plt.get_cmap(self.colormap_phase)
211 211 )
212 212 else:
213 213 ax.plt.set_array(phase.T.ravel())
214 214 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
215 215
216 216
217 217 class RTIPlot(Plot):
218 218 '''
219 219 Plot for RTI data
220 220 '''
221 221
222 222 CODE = 'rti'
223 223 colormap = 'jet'
224 224 plot_type = 'pcolorbuffer'
225 225
226 226 def setup(self):
227 227 self.xaxis = 'time'
228 228 self.ncols = 1
229 229 self.nrows = len(self.data.channels)
230 230 self.nplots = len(self.data.channels)
231 231 self.ylabel = 'Range [km]'
232 232 self.xlabel = 'Time'
233 233 self.cb_label = 'dB'
234 234 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
235 235 self.titles = ['{} Channel {}'.format(
236 236 self.CODE.upper(), x) for x in range(self.nrows)]
237 237
238 238 def update(self, dataOut):
239 239
240 240 data = {}
241 241 meta = {}
242 242 data['rti'] = dataOut.getPower()
243 243 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
244 244
245 245 return data, meta
246 246
247 247 def plot(self):
248 248 self.x = self.data.times
249 249 self.y = self.data.yrange
250 250 self.z = self.data[self.CODE]
251 251 self.z = numpy.ma.masked_invalid(self.z)
252 252
253 253 if self.decimation is None:
254 254 x, y, z = self.fill_gaps(self.x, self.y, self.z)
255 255 else:
256 256 x, y, z = self.fill_gaps(*self.decimate())
257 257
258 258 for n, ax in enumerate(self.axes):
259 259 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
260 260 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
261 261 data = self.data[-1]
262 262 if ax.firsttime:
263 263 ax.plt = ax.pcolormesh(x, y, z[n].T,
264 264 vmin=self.zmin,
265 265 vmax=self.zmax,
266 266 cmap=plt.get_cmap(self.colormap)
267 267 )
268 268 if self.showprofile:
269 269 ax.plot_profile = self.pf_axes[n].plot(
270 270 data['rti'][n], self.y)[0]
271 271 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
272 272 color="k", linestyle="dashed", lw=1)[0]
273 273 else:
274 274 ax.collections.remove(ax.collections[0])
275 275 ax.plt = ax.pcolormesh(x, y, z[n].T,
276 276 vmin=self.zmin,
277 277 vmax=self.zmax,
278 278 cmap=plt.get_cmap(self.colormap)
279 279 )
280 280 if self.showprofile:
281 281 ax.plot_profile.set_data(data['rti'][n], self.y)
282 282 ax.plot_noise.set_data(numpy.repeat(
283 283 data['noise'][n], len(self.y)), self.y)
284 284
285 285
286 286 class CoherencePlot(RTIPlot):
287 287 '''
288 288 Plot for Coherence data
289 289 '''
290 290
291 291 CODE = 'coh'
292 292
293 293 def setup(self):
294 294 self.xaxis = 'time'
295 295 self.ncols = 1
296 296 self.nrows = len(self.data.pairs)
297 297 self.nplots = len(self.data.pairs)
298 298 self.ylabel = 'Range [km]'
299 299 self.xlabel = 'Time'
300 300 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
301 301 if self.CODE == 'coh':
302 302 self.cb_label = ''
303 303 self.titles = [
304 304 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
305 305 else:
306 306 self.cb_label = 'Degrees'
307 307 self.titles = [
308 308 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
309 309
310 310 def update(self, dataOut):
311 311
312 312 data = {}
313 313 meta = {}
314 314 data['coh'] = dataOut.getCoherence()
315 315 meta['pairs'] = dataOut.pairsList
316 316
317 317 return data, meta
318 318
319 319 class PhasePlot(CoherencePlot):
320 320 '''
321 321 Plot for Phase map data
322 322 '''
323 323
324 324 CODE = 'phase'
325 325 colormap = 'seismic'
326 326
327 327 def update(self, dataOut):
328 328
329 329 data = {}
330 330 meta = {}
331 331 data['phase'] = dataOut.getCoherence(phase=True)
332 332 meta['pairs'] = dataOut.pairsList
333 333
334 334 return data, meta
335 335
336 336 class NoisePlot(Plot):
337 337 '''
338 338 Plot for noise
339 339 '''
340 340
341 341 CODE = 'noise'
342 342 plot_type = 'scatterbuffer'
343 343
344 344 def setup(self):
345 345 self.xaxis = 'time'
346 346 self.ncols = 1
347 347 self.nrows = 1
348 348 self.nplots = 1
349 349 self.ylabel = 'Intensity [dB]'
350 350 self.xlabel = 'Time'
351 351 self.titles = ['Noise']
352 352 self.colorbar = False
353 353 self.plots_adjust.update({'right': 0.85 })
354 354
355 355 def update(self, dataOut):
356 356
357 357 data = {}
358 358 meta = {}
359 359 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
360 360 meta['yrange'] = numpy.array([])
361 361
362 362 return data, meta
363 363
364 364 def plot(self):
365 365
366 366 x = self.data.times
367 367 xmin = self.data.min_time
368 368 xmax = xmin + self.xrange * 60 * 60
369 369 Y = self.data['noise']
370 370
371 371 if self.axes[0].firsttime:
372 372 self.ymin = numpy.nanmin(Y) - 5
373 373 self.ymax = numpy.nanmax(Y) + 5
374 374 for ch in self.data.channels:
375 375 y = Y[ch]
376 376 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
377 377 plt.legend(bbox_to_anchor=(1.18, 1.0))
378 378 else:
379 379 for ch in self.data.channels:
380 380 y = Y[ch]
381 381 self.axes[0].lines[ch].set_data(x, y)
382 382
383 383
384 384 class PowerProfilePlot(Plot):
385 385
386 386 CODE = 'pow_profile'
387 387 plot_type = 'scatter'
388 388
389 389 def setup(self):
390 390
391 391 self.ncols = 1
392 392 self.nrows = 1
393 393 self.nplots = 1
394 394 self.height = 4
395 395 self.width = 3
396 396 self.ylabel = 'Range [km]'
397 397 self.xlabel = 'Intensity [dB]'
398 398 self.titles = ['Power Profile']
399 399 self.colorbar = False
400 400
401 401 def update(self, dataOut):
402 402
403 403 data = {}
404 404 meta = {}
405 405 data[self.CODE] = dataOut.getPower()
406 406
407 407 return data, meta
408 408
409 409 def plot(self):
410 410
411 411 y = self.data.yrange
412 412 self.y = y
413 413
414 414 x = self.data[-1][self.CODE]
415 415
416 416 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
417 417 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
418 418
419 419 if self.axes[0].firsttime:
420 420 for ch in self.data.channels:
421 421 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
422 422 plt.legend()
423 423 else:
424 424 for ch in self.data.channels:
425 425 self.axes[0].lines[ch].set_data(x[ch], y)
426 426
427 427
428 428 class SpectraCutPlot(Plot):
429 429
430 430 CODE = 'spc_cut'
431 431 plot_type = 'scatter'
432 432 buffering = False
433 433
434 434 def setup(self):
435 435
436 436 self.nplots = len(self.data.channels)
437 437 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
438 438 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
439 439 self.width = 3.4 * self.ncols + 1.5
440 440 self.height = 3 * self.nrows
441 441 self.ylabel = 'Power [dB]'
442 442 self.colorbar = False
443 443 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
444 444
445 445 def update(self, dataOut):
446 446
447 447 data = {}
448 448 meta = {}
449 449 spc = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
450 450 data['spc'] = spc
451 451 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
452 452 if self.CODE == 'cut_gaussian_fit':
453 453 data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
454 454 data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
455 455 return data, meta
456 456
457 457 def plot(self):
458 458 if self.xaxis == "frequency":
459 459 x = self.data.xrange[0][1:]
460 460 self.xlabel = "Frequency (kHz)"
461 461 elif self.xaxis == "time":
462 462 x = self.data.xrange[1]
463 463 self.xlabel = "Time (ms)"
464 464 else:
465 465 x = self.data.xrange[2][:-1]
466 466 self.xlabel = "Velocity (m/s)"
467 467
468 468 if self.CODE == 'cut_gaussian_fit':
469 469 x = self.data.xrange[2][:-1]
470 470 self.xlabel = "Velocity (m/s)"
471 471
472 472 self.titles = []
473 473
474 474 y = self.data.yrange
475 475 data = self.data[-1]
476 476 z = data['spc']
477 477
478 478 if self.height_index:
479 479 index = numpy.array(self.height_index)
480 480 else:
481 481 index = numpy.arange(0, len(y), int((len(y))/9))
482 482
483 483 for n, ax in enumerate(self.axes):
484 484 if self.CODE == 'cut_gaussian_fit':
485 485 gau0 = data['gauss_fit0']
486 486 gau1 = data['gauss_fit1']
487 487 if ax.firsttime:
488 488 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
489 489 self.xmin = self.xmin if self.xmin else -self.xmax
490 490 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
491 491 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
492 492 ax.plt = ax.plot(x, z[n, :, index].T, lw=0.25)
493 493 if self.CODE == 'cut_gaussian_fit':
494 494 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
495 495 for i, line in enumerate(ax.plt_gau0):
496 496 line.set_color(ax.plt[i].get_color())
497 497 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
498 498 for i, line in enumerate(ax.plt_gau1):
499 499 line.set_color(ax.plt[i].get_color())
500 500 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
501 501 self.figures[0].legend(ax.plt, labels, loc='center right')
502 502 else:
503 503 for i, line in enumerate(ax.plt):
504 504 line.set_data(x, z[n, :, index[i]].T)
505 505 for i, line in enumerate(ax.plt_gau0):
506 506 line.set_data(x, gau0[n, :, index[i]].T)
507 507 line.set_color(ax.plt[i].get_color())
508 508 for i, line in enumerate(ax.plt_gau1):
509 509 line.set_data(x, gau1[n, :, index[i]].T)
510 510 line.set_color(ax.plt[i].get_color())
511 511 self.titles.append('CH {}'.format(n))
512 512
513 513
514 514 class BeaconPhase(Plot):
515 515
516 516 __isConfig = None
517 517 __nsubplots = None
518 518
519 519 PREFIX = 'beacon_phase'
520 520
521 521 def __init__(self):
522 522 Plot.__init__(self)
523 523 self.timerange = 24*60*60
524 524 self.isConfig = False
525 525 self.__nsubplots = 1
526 526 self.counter_imagwr = 0
527 527 self.WIDTH = 800
528 528 self.HEIGHT = 400
529 529 self.WIDTHPROF = 120
530 530 self.HEIGHTPROF = 0
531 531 self.xdata = None
532 532 self.ydata = None
533 533
534 534 self.PLOT_CODE = BEACON_CODE
535 535
536 536 self.FTP_WEI = None
537 537 self.EXP_CODE = None
538 538 self.SUB_EXP_CODE = None
539 539 self.PLOT_POS = None
540 540
541 541 self.filename_phase = None
542 542
543 543 self.figfile = None
544 544
545 545 self.xmin = None
546 546 self.xmax = None
547 547
548 548 def getSubplots(self):
549 549
550 550 ncol = 1
551 551 nrow = 1
552 552
553 553 return nrow, ncol
554 554
555 555 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
556 556
557 557 self.__showprofile = showprofile
558 558 self.nplots = nplots
559 559
560 560 ncolspan = 7
561 561 colspan = 6
562 562 self.__nsubplots = 2
563 563
564 564 self.createFigure(id = id,
565 565 wintitle = wintitle,
566 566 widthplot = self.WIDTH+self.WIDTHPROF,
567 567 heightplot = self.HEIGHT+self.HEIGHTPROF,
568 568 show=show)
569 569
570 570 nrow, ncol = self.getSubplots()
571 571
572 572 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
573 573
574 574 def save_phase(self, filename_phase):
575 575 f = open(filename_phase,'w+')
576 576 f.write('\n\n')
577 577 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
578 578 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
579 579 f.close()
580 580
581 581 def save_data(self, filename_phase, data, data_datetime):
582 582 f=open(filename_phase,'a')
583 583 timetuple_data = data_datetime.timetuple()
584 584 day = str(timetuple_data.tm_mday)
585 585 month = str(timetuple_data.tm_mon)
586 586 year = str(timetuple_data.tm_year)
587 587 hour = str(timetuple_data.tm_hour)
588 588 minute = str(timetuple_data.tm_min)
589 589 second = str(timetuple_data.tm_sec)
590 590 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
591 591 f.close()
592 592
593 593 def plot(self):
594 594 log.warning('TODO: Not yet implemented...')
595 595
596 596 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
597 597 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
598 598 timerange=None,
599 599 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
600 600 server=None, folder=None, username=None, password=None,
601 601 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
602 602
603 603 if dataOut.flagNoData:
604 604 return dataOut
605 605
606 606 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
607 607 return
608 608
609 609 if pairsList == None:
610 610 pairsIndexList = dataOut.pairsIndexList[:10]
611 611 else:
612 612 pairsIndexList = []
613 613 for pair in pairsList:
614 614 if pair not in dataOut.pairsList:
615 615 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
616 616 pairsIndexList.append(dataOut.pairsList.index(pair))
617 617
618 618 if pairsIndexList == []:
619 619 return
620 620
621 621 # if len(pairsIndexList) > 4:
622 622 # pairsIndexList = pairsIndexList[0:4]
623 623
624 624 hmin_index = None
625 625 hmax_index = None
626 626
627 627 if hmin != None and hmax != None:
628 628 indexes = numpy.arange(dataOut.nHeights)
629 629 hmin_list = indexes[dataOut.heightList >= hmin]
630 630 hmax_list = indexes[dataOut.heightList <= hmax]
631 631
632 632 if hmin_list.any():
633 633 hmin_index = hmin_list[0]
634 634
635 635 if hmax_list.any():
636 636 hmax_index = hmax_list[-1]+1
637 637
638 638 x = dataOut.getTimeRange()
639 639
640 640 thisDatetime = dataOut.datatime
641 641
642 642 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
643 643 xlabel = "Local Time"
644 644 ylabel = "Phase (degrees)"
645 645
646 646 update_figfile = False
647 647
648 648 nplots = len(pairsIndexList)
649 649 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
650 650 phase_beacon = numpy.zeros(len(pairsIndexList))
651 651 for i in range(nplots):
652 652 pair = dataOut.pairsList[pairsIndexList[i]]
653 653 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
654 654 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
655 655 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
656 656 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
657 657 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
658 658
659 659 if dataOut.beacon_heiIndexList:
660 660 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
661 661 else:
662 662 phase_beacon[i] = numpy.average(phase)
663 663
664 664 if not self.isConfig:
665 665
666 666 nplots = len(pairsIndexList)
667 667
668 668 self.setup(id=id,
669 669 nplots=nplots,
670 670 wintitle=wintitle,
671 671 showprofile=showprofile,
672 672 show=show)
673 673
674 674 if timerange != None:
675 675 self.timerange = timerange
676 676
677 677 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
678 678
679 679 if ymin == None: ymin = 0
680 680 if ymax == None: ymax = 360
681 681
682 682 self.FTP_WEI = ftp_wei
683 683 self.EXP_CODE = exp_code
684 684 self.SUB_EXP_CODE = sub_exp_code
685 685 self.PLOT_POS = plot_pos
686 686
687 687 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
688 688 self.isConfig = True
689 689 self.figfile = figfile
690 690 self.xdata = numpy.array([])
691 691 self.ydata = numpy.array([])
692 692
693 693 update_figfile = True
694 694
695 695 #open file beacon phase
696 696 path = '%s%03d' %(self.PREFIX, self.id)
697 697 beacon_file = os.path.join(path,'%s.txt'%self.name)
698 698 self.filename_phase = os.path.join(figpath,beacon_file)
699 699 #self.save_phase(self.filename_phase)
700 700
701 701
702 702 #store data beacon phase
703 703 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
704 704
705 705 self.setWinTitle(title)
706 706
707 707
708 708 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
709 709
710 710 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
711 711
712 712 axes = self.axesList[0]
713 713
714 714 self.xdata = numpy.hstack((self.xdata, x[0:1]))
715 715
716 716 if len(self.ydata)==0:
717 717 self.ydata = phase_beacon.reshape(-1,1)
718 718 else:
719 719 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
720 720
721 721
722 722 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
723 723 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
724 724 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
725 725 XAxisAsTime=True, grid='both'
726 726 )
727 727
728 728 self.draw()
729 729
730 730 if dataOut.ltctime >= self.xmax:
731 731 self.counter_imagwr = wr_period
732 732 self.isConfig = False
733 733 update_figfile = True
734 734
735 735 self.save(figpath=figpath,
736 736 figfile=figfile,
737 737 save=save,
738 738 ftp=ftp,
739 739 wr_period=wr_period,
740 740 thisDatetime=thisDatetime,
741 741 update_figfile=update_figfile)
742 742
743 743 return dataOut
Show file before | Show file after | Hide comments
@@ -1,844 +1,863
1 1 '''
2 2 Created on Jul 3, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 # SUBCHANNELS EN VEZ DE CHANNELS
7 7 # BENCHMARKS -> PROBLEMAS CON ARCHIVOS GRANDES -> INCONSTANTE EN EL TIEMPO
8 8 # ACTUALIZACION DE VERSION
9 9 # HEADERS
10 10 # MODULO DE ESCRITURA
11 11 # METADATA
12 12
13 13 import os
14 14 import time
15 15 import datetime
16 16 import numpy
17 17 import timeit
18 18 from fractions import Fraction
19 19 from time import time
20 20 from time import sleep
21 21
22 22 import schainpy.admin
23 23 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
24 24 from schainpy.model.data.jrodata import Voltage
25 25 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
26 26
27 27 import pickle
28 28 try:
29 os.environ["HDF5_USE_FILE_LOCKING"] = "FALSE"
29 30 import digital_rf
30 31 except:
31 32 pass
32 33
33 34
34 35 class DigitalRFReader(ProcessingUnit):
35 36 '''
36 37 classdocs
37 38 '''
38 39
39 40 def __init__(self):
40 41 '''
41 42 Constructor
42 43 '''
43 44
44 45 ProcessingUnit.__init__(self)
45 46
46 47 self.dataOut = Voltage()
47 48 self.__printInfo = True
48 49 self.__flagDiscontinuousBlock = False
49 50 self.__bufferIndex = 9999999
50 51 self.__codeType = 0
51 52 self.__ippKm = None
52 53 self.__nCode = None
53 54 self.__nBaud = None
54 55 self.__code = None
55 56 self.dtype = None
56 57 self.oldAverage = None
57 58 self.path = None
58 59
59 60 def close(self):
60 61 print('Average of writing to digital rf format is ', self.oldAverage * 1000)
61 62 return
62 63
63 64 def __getCurrentSecond(self):
64 65
65 66 return self.__thisUnixSample / self.__sample_rate
66 67
67 68 thisSecond = property(__getCurrentSecond, "I'm the 'thisSecond' property.")
68 69
69 70 def __setFileHeader(self):
70 71 '''
71 72 In this method will be initialized every parameter of dataOut object (header, no data)
72 73 '''
73 74 ippSeconds = 1.0 * self.__nSamples / self.__sample_rate
74 75 if not self.getByBlock:
75 76 nProfiles = 1.0 / ippSeconds # Number of profiles in one second
76 77 else:
77 78 nProfiles = self.nProfileBlocks # Number of profiles in one block
78 79
79 80 try:
80 81 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
81 82 self.__radarControllerHeader)
82 83 except:
83 84 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
84 85 txA=0,
85 86 txB=0,
86 87 nWindows=1,
87 88 nHeights=self.__nSamples,
88 89 firstHeight=self.__firstHeigth,
89 90 deltaHeight=self.__deltaHeigth,
90 91 codeType=self.__codeType,
91 92 nCode=self.__nCode, nBaud=self.__nBaud,
92 93 code=self.__code)
93 94
94 95 try:
95 96 self.dataOut.systemHeaderObj = SystemHeader(self.__systemHeader)
96 97 except:
97 98 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
98 99 nProfiles=nProfiles,
99 100 nChannels=len(
100 101 self.__channelList),
101 102 adcResolution=14)
102 103 self.dataOut.type = "Voltage"
103 104
104 105 self.dataOut.data = None
105 106
106 107 self.dataOut.dtype = self.dtype
107 108
108 109 # self.dataOut.nChannels = 0
109 110
110 111 # self.dataOut.nHeights = 0
111 112
112 113 self.dataOut.nProfiles = int(nProfiles)
113 114
114 115 self.dataOut.heightList = self.__firstHeigth + \
115 116 numpy.arange(self.__nSamples, dtype=numpy.float) * \
116 117 self.__deltaHeigth
117 118
118 119 #self.dataOut.channelList = list(range(self.__num_subchannels))
119 120 self.dataOut.channelList = list(range(len(self.__channelList)))
120 121 if not self.getByBlock:
121 122
122 123 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
123 124 else:
124 125 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights*self.nProfileBlocks
125 126
126 127 # self.dataOut.channelIndexList = None
127 128
128 129 self.dataOut.flagNoData = True
129 130 if not self.getByBlock:
130 131 self.dataOut.flagDataAsBlock = False
131 132 else:
132 133 self.dataOut.flagDataAsBlock = True
133 134 # Set to TRUE if the data is discontinuous
134 135 self.dataOut.flagDiscontinuousBlock = False
135 136
136 137 self.dataOut.utctime = None
137 138
138 139 # timezone like jroheader, difference in minutes between UTC and localtime
139 140 self.dataOut.timeZone = self.__timezone / 60
140 141
141 142 self.dataOut.dstFlag = 0
142 143
143 144 self.dataOut.errorCount = 0
144 145
145 146 try:
146 147 self.dataOut.nCohInt = self.fixed_metadata_dict.get(
147 148 'nCohInt', self.nCohInt)
148 149
149 150 # asumo que la data esta decodificada
150 151 self.dataOut.flagDecodeData = self.fixed_metadata_dict.get(
151 152 'flagDecodeData', self.flagDecodeData)
152 153
153 154 # asumo que la data esta sin flip
154 155 self.dataOut.flagDeflipData = self.fixed_metadata_dict['flagDeflipData']
155 156
156 157 self.dataOut.flagShiftFFT = self.fixed_metadata_dict['flagShiftFFT']
157 158
158 159 self.dataOut.useLocalTime = self.fixed_metadata_dict['useLocalTime']
159 160 except:
160 161 pass
161 162
162 163 self.dataOut.ippSeconds = ippSeconds
163 164
164 165 # Time interval between profiles
165 166 # self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
166 167
167 168 self.dataOut.frequency = self.__frequency
168 169
169 170 self.dataOut.realtime = self.__online
170 171
171 172 def findDatafiles(self, path, startDate=None, endDate=None):
172 173
173 174 if not os.path.isdir(path):
174 175 return []
175 176
176 177 try:
177 178 digitalReadObj = digital_rf.DigitalRFReader(
178 179 path, load_all_metadata=True)
179 180 except:
180 181 digitalReadObj = digital_rf.DigitalRFReader(path)
181 182
182 183 channelNameList = digitalReadObj.get_channels()
183 184
184 185 if not channelNameList:
185 186 return []
186 187
187 188 metadata_dict = digitalReadObj.get_rf_file_metadata(channelNameList[0])
188 189
189 190 sample_rate = metadata_dict['sample_rate'][0]
190 191
191 192 this_metadata_file = digitalReadObj.get_metadata(channelNameList[0])
192 193
193 194 try:
194 195 timezone = this_metadata_file['timezone'].value
195 196 except:
196 197 timezone = 0
197 198
198 199 startUTCSecond, endUTCSecond = digitalReadObj.get_bounds(
199 200 channelNameList[0]) / sample_rate - timezone
200 201
201 202 startDatetime = datetime.datetime.utcfromtimestamp(startUTCSecond)
202 203 endDatatime = datetime.datetime.utcfromtimestamp(endUTCSecond)
203 204
204 205 if not startDate:
205 206 startDate = startDatetime.date()
206 207
207 208 if not endDate:
208 209 endDate = endDatatime.date()
209 210
210 211 dateList = []
211 212
212 213 thisDatetime = startDatetime
213 214
214 215 while(thisDatetime <= endDatatime):
215 216
216 217 thisDate = thisDatetime.date()
217 218
218 219 if thisDate < startDate:
219 220 continue
220 221
221 222 if thisDate > endDate:
222 223 break
223 224
224 225 dateList.append(thisDate)
225 226 thisDatetime += datetime.timedelta(1)
226 227
227 228 return dateList
228 229
229 230 def setup(self, path=None,
230 231 startDate=None,
231 232 endDate=None,
232 233 startTime=datetime.time(0, 0, 0),
233 234 endTime=datetime.time(23, 59, 59),
234 235 channelList=None,
235 236 nSamples=None,
236 237 online=False,
237 238 delay=60,
238 239 buffer_size=1024,
239 240 ippKm=None,
240 241 nCohInt=1,
241 242 nCode=1,
242 243 nBaud=1,
243 244 flagDecodeData=False,
244 245 code=numpy.ones((1, 1), dtype=numpy.int),
245 246 getByBlock=0,
246 247 nProfileBlocks=1,
247 248 **kwargs):
248 249 '''
249 250 In this method we should set all initial parameters.
250 251
251 252 Inputs:
252 253 path
253 254 startDate
254 255 endDate
255 256 startTime
256 257 endTime
257 258 set
258 259 expLabel
259 260 ext
260 261 online
261 262 delay
262 263 '''
263 264 self.path = path
264 265 self.nCohInt = nCohInt
265 266 self.flagDecodeData = flagDecodeData
266 267 self.i = 0
267 268
268 269 self.getByBlock = getByBlock
269 270 self.nProfileBlocks = nProfileBlocks
271 if online:
272 print('Waiting for RF data..')
273 sleep(40)
274
270 275 if not os.path.isdir(path):
271 276 raise ValueError("[Reading] Directory %s does not exist" % path)
272 277
278 #print("path",path)
273 279 try:
274 280 self.digitalReadObj = digital_rf.DigitalRFReader(
275 281 path, load_all_metadata=True)
276 282 except:
277 283 self.digitalReadObj = digital_rf.DigitalRFReader(path)
278 284
279 285 channelNameList = self.digitalReadObj.get_channels()
280 286
281 287 if not channelNameList:
282 288 raise ValueError("[Reading] Directory %s does not have any files" % path)
283 289
284 290 if not channelList:
285 291 channelList = list(range(len(channelNameList)))
286 292
287 293 ########## Reading metadata ######################
288 294
289 295 top_properties = self.digitalReadObj.get_properties(
290 296 channelNameList[channelList[0]])
291 297
292 298 self.__num_subchannels = top_properties['num_subchannels']
293 299 self.__sample_rate = 1.0 * \
294 300 top_properties['sample_rate_numerator'] / \
295 301 top_properties['sample_rate_denominator']
296 302 # self.__samples_per_file = top_properties['samples_per_file'][0]
297 303 self.__deltaHeigth = 1e6 * 0.15 / self.__sample_rate # why 0.15?
298 304
299 305 this_metadata_file = self.digitalReadObj.get_digital_metadata(
300 306 channelNameList[channelList[0]])
301 307 metadata_bounds = this_metadata_file.get_bounds()
302 308 self.fixed_metadata_dict = this_metadata_file.read(
303 309 metadata_bounds[0])[metadata_bounds[0]] # GET FIRST HEADER
304 310
305 311 try:
306 312 self.__processingHeader = self.fixed_metadata_dict['processingHeader']
307 313 self.__radarControllerHeader = self.fixed_metadata_dict['radarControllerHeader']
308 314 self.__systemHeader = self.fixed_metadata_dict['systemHeader']
309 315 self.dtype = pickle.loads(self.fixed_metadata_dict['dtype'])
310 316 except:
311 317 pass
312 318
313 319 self.__frequency = None
314 320
315 321 self.__frequency = self.fixed_metadata_dict.get('frequency', 1)
316 322
317 323 self.__timezone = self.fixed_metadata_dict.get('timezone', 18000)
318 324
319 325 try:
320 326 nSamples = self.fixed_metadata_dict['nSamples']
321 327 except:
322 328 nSamples = None
323 329
324 330 self.__firstHeigth = 0
325 331
326 332 try:
327 333 codeType = self.__radarControllerHeader['codeType']
328 334 except:
329 335 codeType = 0
330 336
331 337 try:
332 338 if codeType:
333 339 nCode = self.__radarControllerHeader['nCode']
334 340 nBaud = self.__radarControllerHeader['nBaud']
335 341 code = self.__radarControllerHeader['code']
336 342 except:
337 343 pass
338 344
339 345 if not ippKm:
340 346 try:
341 347 # seconds to km
342 348 ippKm = self.__radarControllerHeader['ipp']
343 349 except:
344 350 ippKm = None
345 351 ####################################################
346 352 self.__ippKm = ippKm
347 353 startUTCSecond = None
348 354 endUTCSecond = None
349 355
350 356 if startDate:
351 357 startDatetime = datetime.datetime.combine(startDate, startTime)
352 358 startUTCSecond = (
353 startDatetime - datetime.datetime(1970, 1, 1)).total_seconds() + self.__timezone
359 startDatetime - datetime.datetime(1970, 1, 1)).total_seconds()# + self.__timezone
354 360
355 361 if endDate:
356 362 endDatetime = datetime.datetime.combine(endDate, endTime)
357 363 endUTCSecond = (endDatetime - datetime.datetime(1970,
358 1, 1)).total_seconds() + self.__timezone
359
360
361 #print(startUTCSecond,endUTCSecond)
362 start_index, end_index = self.digitalReadObj.get_bounds(
363 channelNameList[channelList[0]])
364
365 #print("*****",start_index,end_index)
364 1, 1)).total_seconds()# + self.__timezone
365 start_index, end_index = self.digitalReadObj.get_bounds(channelNameList[channelList[0]])
366 if start_index==None or end_index==None:
367 print("Check error No data, start_index: ",start_index,",end_index: ",end_index)
368 #return 0
366 369 if not startUTCSecond:
367 370 startUTCSecond = start_index / self.__sample_rate
368
369 371 if start_index > startUTCSecond * self.__sample_rate:
370 372 startUTCSecond = start_index / self.__sample_rate
371 373
372 374 if not endUTCSecond:
373 375 endUTCSecond = end_index / self.__sample_rate
376
374 377 if end_index < endUTCSecond * self.__sample_rate:
375 378 endUTCSecond = end_index / self.__sample_rate #Check UTC and LT time
379
376 380 if not nSamples:
377 381 if not ippKm:
378 382 raise ValueError("[Reading] nSamples or ippKm should be defined")
379 383 nSamples = int(ippKm / (1e6 * 0.15 / self.__sample_rate))
380 384
381 385 channelBoundList = []
382 386 channelNameListFiltered = []
383 387
384 388 for thisIndexChannel in channelList:
385 389 thisChannelName = channelNameList[thisIndexChannel]
386 390 start_index, end_index = self.digitalReadObj.get_bounds(
387 391 thisChannelName)
388 392 channelBoundList.append((start_index, end_index))
389 393 channelNameListFiltered.append(thisChannelName)
390 394
391 395 self.profileIndex = 0
392 396 self.i = 0
393 397 self.__delay = delay
394 398
395 399 self.__codeType = codeType
396 400 self.__nCode = nCode
397 401 self.__nBaud = nBaud
398 402 self.__code = code
399 403
400 404 self.__datapath = path
401 405 self.__online = online
402 406 self.__channelList = channelList
403 407 self.__channelNameList = channelNameListFiltered
404 408 self.__channelBoundList = channelBoundList
405 409 self.__nSamples = nSamples
406 410 if self.getByBlock:
407 411 nSamples = nSamples*nProfileBlocks
408 412
409 413
410 414 self.__samples_to_read = int(nSamples) # FIJO: AHORA 40
411 415 self.__nChannels = len(self.__channelList)
412 416 #print("------------------------------------------")
413 417 #print("self.__samples_to_read",self.__samples_to_read)
414 418 #print("self.__nSamples",self.__nSamples)
415 419 # son iguales y el buffer_index da 0
416 420 self.__startUTCSecond = startUTCSecond
417 421 self.__endUTCSecond = endUTCSecond
418 422
419 423 self.__timeInterval = 1.0 * self.__samples_to_read / \
420 424 self.__sample_rate # Time interval
421 425
422 426 if online:
423 427 # self.__thisUnixSample = int(endUTCSecond*self.__sample_rate - 4*self.__samples_to_read)
424 428 startUTCSecond = numpy.floor(endUTCSecond)
425 429
426 430 # por que en el otro metodo lo primero q se hace es sumar samplestoread
427 431 self.__thisUnixSample = int(startUTCSecond * self.__sample_rate) - self.__samples_to_read
428 432
429 433 #self.__data_buffer = numpy.zeros(
430 434 # (self.__num_subchannels, self.__samples_to_read), dtype=numpy.complex)
435 print("samplestoread",self.__samples_to_read)
431 436 self.__data_buffer = numpy.zeros((int(len(channelList)), self.__samples_to_read), dtype=numpy.complex)
432 437
433 438
434 439 self.__setFileHeader()
435 440 self.isConfig = True
436 441
437 442 print("[Reading] Digital RF Data was found from %s to %s " % (
438 443 datetime.datetime.utcfromtimestamp(
439 444 self.__startUTCSecond - self.__timezone),
440 445 datetime.datetime.utcfromtimestamp(
441 446 self.__endUTCSecond - self.__timezone)
442 447 ))
443 448
444 449 print("[Reading] Starting process from %s to %s" % (datetime.datetime.utcfromtimestamp(startUTCSecond - self.__timezone),
445 datetime.datetime.utcfromtimestamp(
446 endUTCSecond - self.__timezone)
447 ))
450 datetime.datetime.utcfromtimestamp(endUTCSecond - self.__timezone)))
448 451 self.oldAverage = None
449 452 self.count = 0
450 453 self.executionTime = 0
451 454
452 455 def __reload(self):
453 456 # print
454 457 # print "%s not in range [%s, %s]" %(
455 458 # datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
456 459 # datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone),
457 460 # datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone)
458 461 # )
459 462 print("[Reading] reloading metadata ...")
460 463
461 464 try:
462 465 self.digitalReadObj.reload(complete_update=True)
463 466 except:
464 467 self.digitalReadObj = digital_rf.DigitalRFReader(self.path)
465 468
466 469 start_index, end_index = self.digitalReadObj.get_bounds(
467 470 self.__channelNameList[self.__channelList[0]])
468 471
469 472 if start_index > self.__startUTCSecond * self.__sample_rate:
470 473 self.__startUTCSecond = 1.0 * start_index / self.__sample_rate
471 474
472 475 if end_index > self.__endUTCSecond * self.__sample_rate:
473 476 self.__endUTCSecond = 1.0 * end_index / self.__sample_rate
474 477 print()
475 478 print("[Reading] New timerange found [%s, %s] " % (
476 479 datetime.datetime.utcfromtimestamp(
477 480 self.__startUTCSecond - self.__timezone),
478 481 datetime.datetime.utcfromtimestamp(
479 482 self.__endUTCSecond - self.__timezone)
480 483 ))
481 484
482 485 return True
483 486
484 487 return False
485 488
486 489 def timeit(self, toExecute):
487 490 t0 = time.time()
488 491 toExecute()
489 492 self.executionTime = time.time() - t0
490 493 if self.oldAverage is None:
491 494 self.oldAverage = self.executionTime
492 495 self.oldAverage = (self.executionTime + self.count *
493 496 self.oldAverage) / (self.count + 1.0)
494 497 self.count = self.count + 1.0
495 498 return
496 499
497 500 def __readNextBlock(self, seconds=30, volt_scale=1/20000.0):
498 501 '''
499 502 NOTA: APLICACION RADAR METEOROLOGICO
500 503 VALORES OBTENIDOS CON LA USRP, volt_scale = 1,conexion directa al Ch Rx.
501 504
502 505 MAXIMO
503 506 9886 -> 0.980 Voltiospp
504 507 4939 -> 0.480 Voltiospp
505 508 14825 -> 1.440 Voltiospp
506 509 18129 -> 1.940 Voltiospp
507 510 Para llevar al valor correspondiente de Voltaje, debemos dividir por 20000
508 511 y obtenemos la Amplitud correspondiente de entrada IQ.
509 512 volt_scale = (1/20000.0)
510 513 '''
511 514 # Set the next data
512 515 self.__flagDiscontinuousBlock = False
513 516 self.__thisUnixSample += self.__samples_to_read
514 517
515 518 if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
516 519 print ("[Reading] There are no more data into selected time-range")
517 520 if self.__online:
518 521 sleep(3)
519 522 self.__reload()
520 523 else:
521 524 return False
522 525
523 526 if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
524 527 return False
525 528 self.__thisUnixSample -= self.__samples_to_read
526 529
527 530 indexChannel = 0
528 531
529 532 dataOk = False
530 533
531 534 for thisChannelName in self.__channelNameList: # TODO VARIOS CHANNELS?
532 535 for indexSubchannel in range(self.__num_subchannels):
533 536 try:
534 537 t0 = time()
538 #print("thisUNixSample",self.__thisUnixSample)
535 539 result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample,
536 540 self.__samples_to_read,
537 541 thisChannelName, sub_channel=indexSubchannel)
542 #print("result--------------",result)
538 543 self.executionTime = time() - t0
539 544 if self.oldAverage is None:
540 545 self.oldAverage = self.executionTime
541 546 self.oldAverage = (
542 547 self.executionTime + self.count * self.oldAverage) / (self.count + 1.0)
543 548 self.count = self.count + 1.0
544 549
545 550 except IOError as e:
546 551 # read next profile
547 552 self.__flagDiscontinuousBlock = True
548 553 print("[Reading] %s" % datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), e)
549 break
554 bot = 0
555 while(self.__flagDiscontinuousBlock):
556 bot +=1
557 self.__thisUnixSample += self.__sample_rate
558 try:
559 result = result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample,self.__samples_to_read,thisChannelName, sub_channel=indexSubchannel)
560 self.__flagDiscontinuousBlock=False
561 print("Searching.. N°: ",bot,"Success",self.__thisUnixSample)
562 except:
563 print("Searching...N°: ",bot,"Fail", self.__thisUnixSample)
564 if self.__flagDiscontinuousBlock==True:
565 break
566 else:
567 print("New data index found...",self.__thisUnixSample)
568 #break
550 569
551 570 if result.shape[0] != self.__samples_to_read:
552 571 self.__flagDiscontinuousBlock = True
553 572 print("[Reading] %s: Too few samples were found, just %d/%d samples" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
554 573 result.shape[0],
555 574 self.__samples_to_read))
556 575 break
557 576
558 577 self.__data_buffer[indexChannel, :] = result * volt_scale
559 578 indexChannel+=1
560 579
561 580 dataOk = True
562 581
563 582 self.__utctime = self.__thisUnixSample / self.__sample_rate
564 583
565 584 if not dataOk:
566 585 return False
567 586
568 587 print("[Reading] %s: %d samples <> %f sec" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
569 588 self.__samples_to_read,
570 589 self.__timeInterval))
571 590
572 591 self.__bufferIndex = 0
573 592
574 593 return True
575 594
576 595 def __isBufferEmpty(self):
577 596
578 597 return self.__bufferIndex > self.__samples_to_read - self.__nSamples # 40960 - 40
579 598
580 599 def getData(self, seconds=30, nTries=5):
581 600 '''
582 601 This method gets the data from files and put the data into the dataOut object
583 602
584 603 In addition, increase el the buffer counter in one.
585 604
586 605 Return:
587 606 data : retorna un perfil de voltages (alturas * canales) copiados desde el
588 607 buffer. Si no hay mas archivos a leer retorna None.
589 608
590 609 Affected:
591 610 self.dataOut
592 611 self.profileIndex
593 612 self.flagDiscontinuousBlock
594 613 self.flagIsNewBlock
595 614 '''
596 615 #print("getdata")
597 616 err_counter = 0
598 617 self.dataOut.flagNoData = True
599 618
600 619
601 620 if self.__isBufferEmpty():
602 621 #print("hi")
603 622 self.__flagDiscontinuousBlock = False
604 623
605 624 while True:
606 625 if self.__readNextBlock():
607 626 break
608 627 if self.__thisUnixSample > self.__endUTCSecond * self.__sample_rate:
609 628 raise schainpy.admin.SchainError('Error')
610 629 return
611 630
612 631 if self.__flagDiscontinuousBlock:
613 632 raise schainpy.admin.SchainError('discontinuous block found')
614 633 return
615 634
616 635 if not self.__online:
617 636 raise schainpy.admin.SchainError('Online?')
618 637 return
619 638
620 639 err_counter += 1
621 640 if err_counter > nTries:
622 641 raise schainpy.admin.SchainError('Max retrys reach')
623 642 return
624 643
625 644 print('[Reading] waiting %d seconds to read a new block' % seconds)
626 645 sleep(seconds)
627 646
628 647
629 648 if not self.getByBlock:
630 649
631 650 #print("self.__bufferIndex",self.__bufferIndex)# este valor siempre es cero aparentemente
632 651 self.dataOut.data = self.__data_buffer[:, self.__bufferIndex:self.__bufferIndex + self.__nSamples]
633 652 self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
634 653 self.dataOut.flagNoData = False
635 654 self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
636 655 self.dataOut.profileIndex = self.profileIndex
637 656
638 657 self.__bufferIndex += self.__nSamples
639 658 self.profileIndex += 1
640 659
641 660 if self.profileIndex == self.dataOut.nProfiles:
642 661 self.profileIndex = 0
643 662 else:
644 663 # ojo debo anadir el readNextBLock y el __isBufferEmpty(
645 664 self.dataOut.flagNoData = False
646 665 buffer = self.__data_buffer[:,self.__bufferIndex:self.__bufferIndex + self.__samples_to_read]
647 666 buffer = buffer.reshape((self.__nChannels, self.nProfileBlocks, int(self.__samples_to_read/self.nProfileBlocks)))
648 667 self.dataOut.nProfileBlocks = self.nProfileBlocks
649 668 self.dataOut.data = buffer
650 669 self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
651 670 self.profileIndex += self.__samples_to_read
652 671 self.__bufferIndex += self.__samples_to_read
653 672 self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
654 673 return True
655 674
656 675
657 676 def printInfo(self):
658 677 '''
659 678 '''
660 679 if self.__printInfo == False:
661 680 return
662 681
663 682 # self.systemHeaderObj.printInfo()
664 683 # self.radarControllerHeaderObj.printInfo()
665 684
666 685 self.__printInfo = False
667 686
668 687 def printNumberOfBlock(self):
669 688 '''
670 689 '''
671 690 return
672 691 # print self.profileIndex
673 692
674 693 def run(self, **kwargs):
675 694 '''
676 695 This method will be called many times so here you should put all your code
677 696 '''
678 697
679 698 if not self.isConfig:
680 699 self.setup(**kwargs)
681 700
682 701 self.getData(seconds=self.__delay)
683 702
684 703 return
685 704
686 705 @MPDecorator
687 706 class DigitalRFWriter(Operation):
688 707 '''
689 708 classdocs
690 709 '''
691 710
692 711 def __init__(self, **kwargs):
693 712 '''
694 713 Constructor
695 714 '''
696 715 Operation.__init__(self, **kwargs)
697 716 self.metadata_dict = {}
698 717 self.dataOut = None
699 718 self.dtype = None
700 719 self.oldAverage = 0
701 720
702 721 def setHeader(self):
703 722
704 723 self.metadata_dict['frequency'] = self.dataOut.frequency
705 724 self.metadata_dict['timezone'] = self.dataOut.timeZone
706 725 self.metadata_dict['dtype'] = pickle.dumps(self.dataOut.dtype)
707 726 self.metadata_dict['nProfiles'] = self.dataOut.nProfiles
708 727 self.metadata_dict['heightList'] = self.dataOut.heightList
709 728 self.metadata_dict['channelList'] = self.dataOut.channelList
710 729 self.metadata_dict['flagDecodeData'] = self.dataOut.flagDecodeData
711 730 self.metadata_dict['flagDeflipData'] = self.dataOut.flagDeflipData
712 731 self.metadata_dict['flagShiftFFT'] = self.dataOut.flagShiftFFT
713 732 self.metadata_dict['useLocalTime'] = self.dataOut.useLocalTime
714 733 self.metadata_dict['nCohInt'] = self.dataOut.nCohInt
715 734 self.metadata_dict['type'] = self.dataOut.type
716 735 self.metadata_dict['flagDataAsBlock']= getattr(
717 736 self.dataOut, 'flagDataAsBlock', None) # chequear
718 737
719 738 def setup(self, dataOut, path, frequency, fileCadence, dirCadence, metadataCadence, set=0, metadataFile='metadata', ext='.h5'):
720 739 '''
721 740 In this method we should set all initial parameters.
722 741 Input:
723 742 dataOut: Input data will also be outputa data
724 743 '''
725 744 self.setHeader()
726 745 self.__ippSeconds = dataOut.ippSeconds
727 746 self.__deltaH = dataOut.getDeltaH()
728 747 self.__sample_rate = 1e6 * 0.15 / self.__deltaH
729 748 self.__dtype = dataOut.dtype
730 749 if len(dataOut.dtype) == 2:
731 750 self.__dtype = dataOut.dtype[0]
732 751 self.__nSamples = dataOut.systemHeaderObj.nSamples
733 752 self.__nProfiles = dataOut.nProfiles
734 753
735 754 if self.dataOut.type != 'Voltage':
736 755 raise 'Digital RF cannot be used with this data type'
737 756 self.arr_data = numpy.ones((1, dataOut.nFFTPoints * len(
738 757 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
739 758 else:
740 759 self.arr_data = numpy.ones((self.__nSamples, len(
741 760 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
742 761
743 762 file_cadence_millisecs = 1000
744 763
745 764 sample_rate_fraction = Fraction(self.__sample_rate).limit_denominator()
746 765 sample_rate_numerator = int(sample_rate_fraction.numerator)
747 766 sample_rate_denominator = int(sample_rate_fraction.denominator)
748 767 start_global_index = dataOut.utctime * self.__sample_rate
749 768
750 769 uuid = 'prueba'
751 770 compression_level = 0
752 771 checksum = False
753 772 is_complex = True
754 773 num_subchannels = len(dataOut.channelList)
755 774 is_continuous = True
756 775 marching_periods = False
757 776
758 777 self.digitalWriteObj = digital_rf.DigitalRFWriter(path, self.__dtype, dirCadence,
759 778 fileCadence, start_global_index,
760 779 sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum,
761 780 is_complex, num_subchannels, is_continuous, marching_periods)
762 781 metadata_dir = os.path.join(path, 'metadata')
763 782 os.system('mkdir %s' % (metadata_dir))
764 783 self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, # 236, file_cadence_millisecs / 1000
765 784 sample_rate_numerator, sample_rate_denominator,
766 785 metadataFile)
767 786 self.isConfig = True
768 787 self.currentSample = 0
769 788 self.oldAverage = 0
770 789 self.count = 0
771 790 return
772 791
773 792 def writeMetadata(self):
774 793 start_idx = self.__sample_rate * self.dataOut.utctime
775 794
776 795 self.metadata_dict['processingHeader'] = self.dataOut.processingHeaderObj.getAsDict(
777 796 )
778 797 self.metadata_dict['radarControllerHeader'] = self.dataOut.radarControllerHeaderObj.getAsDict(
779 798 )
780 799 self.metadata_dict['systemHeader'] = self.dataOut.systemHeaderObj.getAsDict(
781 800 )
782 801 self.digitalMetadataWriteObj.write(start_idx, self.metadata_dict)
783 802 return
784 803
785 804 def timeit(self, toExecute):
786 805 t0 = time()
787 806 toExecute()
788 807 self.executionTime = time() - t0
789 808 if self.oldAverage is None:
790 809 self.oldAverage = self.executionTime
791 810 self.oldAverage = (self.executionTime + self.count *
792 811 self.oldAverage) / (self.count + 1.0)
793 812 self.count = self.count + 1.0
794 813 return
795 814
796 815 def writeData(self):
797 816 if self.dataOut.type != 'Voltage':
798 817 raise 'Digital RF cannot be used with this data type'
799 818 for channel in self.dataOut.channelList:
800 819 for i in range(self.dataOut.nFFTPoints):
801 820 self.arr_data[1][channel * self.dataOut.nFFTPoints +
802 821 i]['r'] = self.dataOut.data[channel][i].real
803 822 self.arr_data[1][channel * self.dataOut.nFFTPoints +
804 823 i]['i'] = self.dataOut.data[channel][i].imag
805 824 else:
806 825 for i in range(self.dataOut.systemHeaderObj.nSamples):
807 826 for channel in self.dataOut.channelList:
808 827 self.arr_data[i][channel]['r'] = self.dataOut.data[channel][i].real
809 828 self.arr_data[i][channel]['i'] = self.dataOut.data[channel][i].imag
810 829
811 830 def f(): return self.digitalWriteObj.rf_write(self.arr_data)
812 831 self.timeit(f)
813 832
814 833 return
815 834
816 835 def run(self, dataOut, frequency=49.92e6, path=None, fileCadence=1000, dirCadence=36000, metadataCadence=1, **kwargs):
817 836 '''
818 837 This method will be called many times so here you should put all your code
819 838 Inputs:
820 839 dataOut: object with the data
821 840 '''
822 841 # print dataOut.__dict__
823 842 self.dataOut = dataOut
824 843 if not self.isConfig:
825 844 self.setup(dataOut, path, frequency, fileCadence,
826 845 dirCadence, metadataCadence, **kwargs)
827 846 self.writeMetadata()
828 847
829 848 self.writeData()
830 849
831 850 ## self.currentSample += 1
832 851 # if self.dataOut.flagDataAsBlock or self.currentSample == 1:
833 852 # self.writeMetadata()
834 853 ## if self.currentSample == self.__nProfiles: self.currentSample = 0
835 854
836 855 return dataOut# en la version 2.7 no aparece este return
837 856
838 857 def close(self):
839 858 print('[Writing] - Closing files ')
840 859 print('Average of writing to digital rf format is ', self.oldAverage * 1000)
841 860 try:
842 861 self.digitalWriteObj.close()
843 862 except:
844 863 pass
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,375 +1,375
1 1 # SOPHY PROC script
2 2 import os, sys, json, argparse
3 3 import datetime
4 4 import time
5 5
6 6 PATH = '/DATA_RM/DATA'
7 7 PATH = '/media/jespinoza/Elements'
8 8 PATH = '/media/jespinoza/data/SOPHY'
9 9 PATH = '/home/soporte/Documents/EVENTO'
10 10
11 11 PARAM = {
12 12 'S': {'zmin': -45, 'zmax': -25, 'colormap': 'jet', 'label': 'Power', 'wrname': 'power','cb_label': 'dBm', 'ch':0},
13 13 'SNR': {'zmin': -40, 'zmax': -20, 'colormap': 'jet', 'label': 'SNR', 'wrname': 'snr','cb_label': 'dB', 'ch':0},
14 14 'V': {'zmin': -12, 'zmax': 12, 'colormap': 'sophy_v', 'label': 'Velocity', 'wrname': 'velocity', 'cb_label': 'm/s', 'ch':0},
15 15 'R': {'zmin': 0, 'zmax': 1, 'colormap': 'jet', 'label': 'RhoHV', 'wrname':'rhoHV', 'cb_label': '*', 'ch':0},
16 16 'P': {'zmin': -180,'zmax': 180,'colormap': 'RdBu_r', 'label': 'PhiDP', 'wrname':'phiDP' , 'cb_label': 'º', 'ch':0},
17 17 'D': {'zmin': -30, 'zmax': 80, 'colormap': 'sophy_r','label': 'ZDR','wrname':'differential_reflectivity' , 'cb_label': 'dBz','ch':0},
18 18 'Z': {'zmin': -30, 'zmax': 80, 'colormap': 'sophy_r','label': 'Reflectivity ', 'wrname':'reflectivity', 'cb_label': 'dBz','ch':0},
19 19 'W': {'zmin': 0, 'zmax': 15, 'colormap': 'sophy_w','label': 'Spectral Width', 'wrname':'spectral_width', 'cb_label': 'm/s', 'ch':0}
20 20 }
21 21
22 22 def max_index(r, sample_rate, ipp):
23 23
24 24 return int(sample_rate*ipp*1e6 * r / 60) + int(sample_rate*ipp*1e6 * 1.2 / 60)
25 25
26 26 def main(args):
27 27
28 28 experiment = args.experiment
29 29 fp = open(os.path.join(PATH, experiment, 'experiment.conf'))
30 30 conf = json.loads(fp.read())
31 31
32 32 ipp_km = conf['usrp_tx']['ipp']
33 33 ipp = ipp_km * 2 /300000
34 34 sample_rate = conf['usrp_rx']['sample_rate']
35 35 axis = ['0' if x=='elevation' else '1' for x in conf['pedestal']['axis']] # AZIMUTH 1 ELEVACION 0
36 36 speed_axis = conf['pedestal']['speed']
37 37 steps = conf['pedestal']['table']
38 38 time_offset = args.time_offset
39 39 parameters = args.parameters
40 40 start_date = experiment.split('@')[1].split('T')[0].replace('-', '/')
41 41 end_date = start_date
42 42 if args.start_time:
43 43 start_time = args.start_time
44 44 else:
45 45 start_time = experiment.split('@')[1].split('T')[1].replace('-', ':')
46 46 end_time = '23:59:59'
47 47 N = int(1/(speed_axis[0]*ipp)) # 1 GRADO DE RESOLUCION
48 48 path = os.path.join(PATH, experiment, 'rawdata')
49 49 path_ped = os.path.join(PATH, experiment, 'position')
50 50 path_plots = os.path.join(PATH, experiment, 'plotsC0N'+str(args.range))
51 51 path_save = os.path.join(PATH, experiment, 'paramC0N'+str(args.range))
52 52 RMIX = 1.62
53 53 H0 = -1.68
54 54 MASK = 0.3
55 55
56 56 from schainpy.controller import Project
57 57
58 58 project = Project()
59 59 project.setup(id='1', name='Sophy', description='sophy proc')
60 60
61 61 reader = project.addReadUnit(datatype='DigitalRFReader',
62 62 path=path,
63 63 startDate=start_date,
64 64 endDate=end_date,
65 65 startTime=start_time,
66 66 endTime=end_time,
67 67 delay=30,
68 68 online=args.online,
69 69 walk=1,
70 70 ippKm = ipp_km,
71 71 getByBlock = 1,
72 72 nProfileBlocks = N,
73 73 )
74 74
75 75 if not conf['usrp_tx']['enable_2']: # One Pulse
76 76 voltage = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
77 77
78 78 if conf['usrp_tx']['code_type_1'] != 'None':
79 79 codes = [ c.strip() for c in conf['usrp_tx']['code_1'].split(',')]
80 80 code = []
81 81 for c in codes:
82 82 code.append([int(x) for x in c])
83 83 op = voltage.addOperation(name='Decoder', optype='other')
84 84 op.addParameter(name='code', value=code)
85 85 op.addParameter(name='nCode', value=len(code), format='int')
86 86 op.addParameter(name='nBaud', value=len(code[0]), format='int')
87 87
88 88 op = voltage.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
89 89 op.addParameter(name='n', value=len(code), format='int')
90 90 ncode = len(code)
91 91 else:
92 92 ncode = 1
93 93 code = ['0']
94 94
95 95 op = voltage.addOperation(name='setH0')
96 96 op.addParameter(name='h0', value=H0)
97 97
98 98 if args.range > 0:
99 99 op = voltage.addOperation(name='selectHeights')
100 100 op.addParameter(name='minIndex', value='0', format='int')
101 101 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
102 102
103 103 op = voltage.addOperation(name='PulsePair_vRF', optype='other')
104 104 op.addParameter(name='n', value=int(N)/ncode, format='int')
105 105 #op.addParameter(name='removeDC', value=1, format='int')
106 106
107 107
108 108 proc = project.addProcUnit(datatype='ParametersProc', inputId=voltage.getId())
109 109
110 110 opObj10 = proc.addOperation(name="WeatherRadar")
111 111 opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0]))
112 112 opObj10.addParameter(name='Pt',value=((1e-6/sample_rate)*len(code[0])/ipp)*200)
113 113
114 114 op = proc.addOperation(name='PedestalInformation')
115 115 op.addParameter(name='path', value=path_ped, format='str')
116 116 op.addParameter(name='interval', value='0.04')
117 117 op.addParameter(name='time_offset', value=time_offset)
118 118 op.addParameter(name='mode', value='PPI')
119 119
120 120 for param in parameters:
121 121 op = proc.addOperation(name='Block360')
122 122 op.addParameter(name='runNextOp', value=True)
123 123
124 124 op= proc.addOperation(name='WeatherParamsPlot')
125 125 if args.save: op.addParameter(name='save', value=path_plots, format='str')
126 126 op.addParameter(name='save_period', value=-1)
127 127 op.addParameter(name='show', value=args.show)
128 128 op.addParameter(name='channels', value='1,')
129 129 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
130 130 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
131 131 op.addParameter(name='attr_data', value=param, format='str')
132 132 op.addParameter(name='labels', value=[PARAM[param]['label']])
133 133 op.addParameter(name='save_code', value=param)
134 134 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
135 135 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
136 136 op.addParameter(name='bgcolor', value='black')
137 137 if MASK: op.addParameter(name='mask', value=MASK, format='float')
138 138 if args.server:
139 139 op.addParameter(name='server', value='0.0.0.0:4444')
140 140 op.addParameter(name='exp_code', value='400')
141 141
142 142 desc = {
143 143 'Data': {
144 144 param: PARAM[param]['wrname'],
145 145 'utctime': 'time'
146 146 },
147 147 'Metadata': {
148 148 'heightList': 'range',
149 149 'data_azi': 'azimuth',
150 150 'data_ele': 'elevation',
151 151 'mode_op': 'scan_type',
152 152 'h0': 'range_correction',
153 153 }
154 154 }
155 155
156 156 if args.save:
157 157 opObj10 = proc.addOperation(name='HDFWriter')
158 158 writer.addParameter(name='path', value=path_save, format='str')
159 159 writer.addParameter(name='Reset', value=True)
160 160 writer.addParameter(name='setType', value='weather')
161 161 writer.addParameter(name='description', value=json.dumps(desc))
162 162 writer.addParameter(name='blocksPerFile', value='1',format='int')
163 163 writer.addParameter(name='metadataList', value='heightList,data_azi,data_ele,mode_op,latitude,longitude,altitude,heading,radar_name,institution,contact,h0,range_unit')
164 164 writer.addParameter(name='dataList', value='{},utctime'.format(param))
165 165 writer.addParameter(name='mask', value=MASK, format='float')
166 166 # meta
167 167 writer.addParameter(name='latitude', value='-12.040436')
168 168 writer.addParameter(name='longitude', value='-75.295893')
169 169 writer.addParameter(name='altitude', value='3379.2147')
170 170 writer.addParameter(name='heading', value='0')
171 171 writer.addParameter(name='radar_name', value='SOPHy')
172 172 writer.addParameter(name='institution', value='IGP')
173 173 writer.addParameter(name='contact', value='dscipion@igp.gob.pe')
174 174 writer.addParameter(name='created_by', value='Signal Chain (https://pypi.org/project/schainpy/)')
175 175 writer.addParameter(name='range_unit', value='km')
176 176
177 177 else: #Two pulses
178 178
179 179 voltage1 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
180 180
181 181 op = voltage1.addOperation(name='ProfileSelector')
182 182 op.addParameter(name='profileRangeList', value='0,{}'.format(conf['usrp_tx']['repetitions_1']-1))
183 183
184 184 if conf['usrp_tx']['code_type_1'] != 'None':
185 185 codes = [ c.strip() for c in conf['usrp_tx']['code_1'].split(',')]
186 186 code = []
187 187 for c in codes:
188 188 code.append([int(x) for x in c])
189 189 op = voltage1.addOperation(name='Decoder', optype='other')
190 190 op.addParameter(name='code', value=code)
191 191 op.addParameter(name='nCode', value=len(code), format='int')
192 192 op.addParameter(name='nBaud', value=len(code[0]), format='int')
193 193 else:
194 194 code = ['0']
195 195
196 196 op = voltage1.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
197 197 op.addParameter(name='n', value=2, format='int')
198 198
199 199 if args.range > 0:
200 200 op = voltage1.addOperation(name='selectHeights')
201 201 op.addParameter(name='minIndex', value='0', format='int')
202 202 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
203 203
204 204 op = voltage1.addOperation(name='setH0')
205 205 op.addParameter(name='h0', value=H0, format='float')
206 206
207 207 op = voltage1.addOperation(name='PulsePair_vRF', optype='other')
208 208 op.addParameter(name='n', value=int(conf['usrp_tx']['repetitions_1'])/2, format='int')
209 209 #op.addParameter(name='removeDC', value=1, format='int')
210 210
211 211
212 212 proc1 = project.addProcUnit(datatype='ParametersProc', inputId=voltage1.getId())
213 213 proc1.addParameter(name='runNextUnit', value=True)
214 214
215 215 opObj10 = proc1.addOperation(name="WeatherRadar")
216 216 opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0]))
217 opObj10.addParameter(name='Pt',value=((1e-6/sample_rate)*len(code[0])/ipp)*200)
217 opObj10.addParameter(name='Pt',value=200)
218 218
219 219 op = proc1.addOperation(name='PedestalInformation')
220 220 op.addParameter(name='path', value=path_ped, format='str')
221 221 op.addParameter(name='interval', value='0.04')
222 222 op.addParameter(name='time_offset', value=time_offset)
223 223 op.addParameter(name='mode', value='PPI')
224 224
225 225 op = proc1.addOperation(name='Block360')
226 226 op.addParameter(name='attr_data', value='data_param')
227 227 op.addParameter(name='runNextOp', value=True)
228 228
229 229
230 230 voltage2 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
231 231
232 232 op = voltage2.addOperation(name='ProfileSelector')
233 233 op.addParameter(name='profileRangeList', value='{},{}'.format(conf['usrp_tx']['repetitions_1'], conf['usrp_tx']['repetitions_1']+conf['usrp_tx']['repetitions_2']-1))
234 234
235 235 if conf['usrp_tx']['code_type_2']:
236 236 codes = [ c.strip() for c in conf['usrp_tx']['code_2'].split(',')]
237 237 code = []
238 238 for c in codes:
239 239 code.append([int(x) for x in c])
240 240 op = voltage2.addOperation(name='Decoder', optype='other')
241 241 op.addParameter(name='code', value=code)
242 242 op.addParameter(name='nCode', value=len(code), format='int')
243 243 op.addParameter(name='nBaud', value=len(code[0]), format='int')
244 244
245 245 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
246 246 op.addParameter(name='n', value=len(code), format='int')
247 247 ncode = len(code)
248 248 else:
249 249 ncode = 1
250 250
251 251 if args.range > 0:
252 252 op = voltage2.addOperation(name='selectHeights')
253 253 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
254 254 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
255 255
256 256 op = voltage2.addOperation(name='setH0')
257 257 op.addParameter(name='h0', value=H0, format='float')
258 258
259 259 op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
260 260 op.addParameter(name='n', value=int(conf['usrp_tx']['repetitions_2'])/ncode, format='int')
261 261 #op.addParameter(name='removeDC', value=1, format='int')
262 262
263 263
264 264 proc2 = project.addProcUnit(datatype='ParametersProc', inputId=voltage2.getId())
265 265 proc2.addParameter(name='runNextUnit', value=True)
266 266
267 267 opObj10 = proc2.addOperation(name="WeatherRadar")
268 268 opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0]))
269 opObj10.addParameter(name='Pt',value=((1e-6/sample_rate)*len(code[0])/ipp)*200)
269 opObj10.addParameter(name='Pt',value=200)
270 270
271 271 op = proc2.addOperation(name='PedestalInformation')
272 272 op.addParameter(name='path', value=path_ped, format='str')
273 273 op.addParameter(name='interval', value='0.04')
274 274 op.addParameter(name='time_offset', value=time_offset)
275 275 op.addParameter(name='mode', value='PPI')
276 276
277 277 op = proc2.addOperation(name='Block360')
278 278 op.addParameter(name='attr_data', value='data_param')
279 279 op.addParameter(name='runNextOp', value=True)
280 280
281 281 merge = project.addProcUnit(datatype='MergeProc', inputId=[proc1.getId(), proc2.getId()])
282 282 merge.addParameter(name='attr_data', value='data_param')
283 283 merge.addParameter(name='mode', value='7') #RM
284 284
285 285
286 286 for param in parameters:
287 287
288 288 if args.plot:
289 289 op= merge.addOperation(name='WeatherParamsPlot')
290 290 if args.save:
291 291 op.addParameter(name='save', value=path_plots, format='str')
292 292 op.addParameter(name='save_period', value=-1)
293 293 op.addParameter(name='show', value=args.show)
294 294 op.addParameter(name='channels', value='0,')
295 295 op.addParameter(name='zmin', value=PARAM[param]['zmin'], format='int')
296 296 op.addParameter(name='zmax', value=PARAM[param]['zmax'], format='int')
297 297 op.addParameter(name='attr_data', value=param, format='str')
298 298 op.addParameter(name='labels', value=[PARAM[param]['label']])
299 299 op.addParameter(name='save_code', value=param)
300 300 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
301 301 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
302 302 op.addParameter(name='bgcolor', value='black')
303 303 if MASK: op.addParameter(name='mask', value=MASK, format='float')
304 304 if args.server:
305 305 op.addParameter(name='server', value='0.0.0.0:4444')
306 306 op.addParameter(name='exp_code', value='400')
307 307
308 308 desc = {
309 309 'Data': {
310 310 'data_param': {PARAM[param]['wrname']: ['H', 'V']},
311 311 'utctime': 'time'
312 312 },
313 313 'Metadata': {
314 314 'heightList': 'range',
315 315 'data_azi': 'azimuth',
316 316 'data_ele': 'elevation',
317 317 'mode_op': 'scan_type',
318 318 'h0': 'range_correction',
319 319 }
320 320 }
321 321
322 322 if args.save:
323 323 writer = merge.addOperation(name='HDFWriter')
324 324 writer.addParameter(name='path', value=path_save, format='str')
325 325 writer.addParameter(name='Reset', value=True)
326 326 writer.addParameter(name='setType', value='weather')
327 327 writer.addParameter(name='description', value=json.dumps(desc))
328 328 writer.addParameter(name='blocksPerFile', value='1',format='int')
329 329 writer.addParameter(name='metadataList', value='heightList,data_azi,data_ele,mode_op,latitude,longitude,altitude,heading,radar_name,institution,contact,h0,range_unit')
330 330 writer.addParameter(name='dataList', value='data_param,utctime')
331 331 writer.addParameter(name='weather_var', value=param)
332 332 writer.addParameter(name='mask', value=MASK, format='float')
333 333 # meta
334 334 writer.addParameter(name='latitude', value='-12.040436')
335 335 writer.addParameter(name='longitude', value='-75.295893')
336 336 writer.addParameter(name='altitude', value='3379.2147')
337 337 writer.addParameter(name='heading', value='0')
338 338 writer.addParameter(name='radar_name', value='SOPHy')
339 339 writer.addParameter(name='institution', value='IGP')
340 340 writer.addParameter(name='contact', value='dscipion@igp.gob.pe')
341 341 writer.addParameter(name='created_by', value='Signal Chain (https://pypi.org/project/schainpy/)')
342 342 writer.addParameter(name='range_unit', value='km')
343 343
344 344 project.start()
345 345
346 346 if __name__ == '__main__':
347 347
348 348 parser = argparse.ArgumentParser(description='Script to process SOPHy data.')
349 349 parser.add_argument('experiment',
350 350 help='Experiment name')
351 351 parser.add_argument('--parameters', nargs='*', default=['S'],
352 352 help='Variables to process: P, Z, V')
353 353 parser.add_argument('--time_offset', default=0,
354 354 help='Fix time offset')
355 355 parser.add_argument('--range', default=0, type=float,
356 356 help='Max range to plot')
357 357 parser.add_argument('--save', action='store_true',
358 358 help='Create output files')
359 359 parser.add_argument('--plot', action='store_true',
360 360 help='Create plot files')
361 361 parser.add_argument('--show', action='store_true',
362 362 help='Show matplotlib plot.')
363 363 parser.add_argument('--online', action='store_true',
364 364 help='Set online mode.')
365 365 parser.add_argument('--server', action='store_true',
366 366 help='Send to realtime')
367 367 parser.add_argument('--start_time', default='',
368 368 help='Set start time.')
369 369
370 370
371 371 args = parser.parse_args()
372 372
373 373 main(args)
374 374
375 375 # python sophy_proc.py HYO_PM@2022-06-09T15-05-12 --parameters V --plot --save --show --range 36S
General Comments 0
You need to be logged in to leave comments. Login now