##// END OF EJS Templates
schain
RSS
act online
joabAM -
r1280:3e5818298b5e
parent child
Show More
Side by Side Unified
Context file:
r1280:3e5818298b5e Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,779 +1,779
1 1 '''
2 2 New Plots Operations
3 3
4 4 @author: juan.espinoza@jro.igp.gob.pe
5 5 '''
6 6
7 7
8 8 import time
9 9 import datetime
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt
13 13 from schainpy.utils import log
14 14
15 15 EARTH_RADIUS = 6.3710e3
16 16
17 17
18 18 def ll2xy(lat1, lon1, lat2, lon2):
19 19
20 20 p = 0.017453292519943295
21 21 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
22 22 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
23 23 r = 12742 * numpy.arcsin(numpy.sqrt(a))
24 24 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
25 25 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
26 26 theta = -theta + numpy.pi/2
27 27 return r*numpy.cos(theta), r*numpy.sin(theta)
28 28
29 29
30 30 def km2deg(km):
31 31 '''
32 32 Convert distance in km to degrees
33 33 '''
34 34
35 35 return numpy.rad2deg(km/EARTH_RADIUS)
36 36
37 37
38 38 class SpectraPlot(Plot):
39 39 '''
40 40 Plot for Spectra data
41 41 '''
42 42
43 43 CODE = 'spc'
44 44 colormap = 'jro'
45 45 plot_name = 'Spectra'
46 46 plot_type = 'pcolor'
47 47
48 48 def setup(self):
49 49 self.nplots = len(self.data.channels)
50 50 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
51 51 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
52 52 self.width = 3.4 * self.ncols
53 53 self.height = 3 * self.nrows
54 54 self.cb_label = 'dB'
55 55 if self.showprofile:
56 56 self.width += 0.8 * self.ncols
57 57
58 58 self.ylabel = 'Range [km]'
59 59
60 60 def plot(self):
61 61 if self.xaxis == "frequency":
62 62 x = self.data.xrange[0]
63 63 self.xlabel = "Frequency (kHz)"
64 64 elif self.xaxis == "time":
65 65 x = self.data.xrange[1]
66 66 self.xlabel = "Time (ms)"
67 67 else:
68 68 x = self.data.xrange[2]
69 69 self.xlabel = "Velocity (m/s)"
70 70
71 71 if self.CODE == 'spc_moments':
72 72 x = self.data.xrange[2]
73 73 self.xlabel = "Velocity (m/s)"
74 74
75 75 self.titles = []
76 76
77 77 y = self.data.heights
78 78 self.y = y
79 79 z = self.data['spc']
80 80
81 81 for n, ax in enumerate(self.axes):
82 82 noise = self.data['noise'][n][-1]
83 83 if self.CODE == 'spc_moments':
84 84 mean = self.data['moments'][n, :, 1, :][-1]
85 85 if ax.firsttime:
86 86 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
87 87 self.xmin = self.xmin if self.xmin else -self.xmax
88 88 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
89 89 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
90 90 ax.plt = ax.pcolormesh(x, y, z[n].T,
91 91 vmin=self.zmin,
92 92 vmax=self.zmax,
93 93 cmap=plt.get_cmap(self.colormap)
94 94 )
95 95
96 96 if self.showprofile:
97 97 ax.plt_profile = self.pf_axes[n].plot(
98 98 self.data['rti'][n][-1], y)[0]
99 99 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
100 100 color="k", linestyle="dashed", lw=1)[0]
101 101 if self.CODE == 'spc_moments':
102 102 ax.plt_mean = ax.plot(mean, y, color='k')[0]
103 103 else:
104 104 ax.plt.set_array(z[n].T.ravel())
105 105 if self.showprofile:
106 106 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
107 107 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
108 108 if self.CODE == 'spc_moments':
109 109 ax.plt_mean.set_data(mean, y)
110 110 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
111 111
112 112
113 113 class CrossSpectraPlot(Plot):
114 114
115 115 CODE = 'cspc'
116 116 colormap = 'jet'
117 117 plot_name = 'CrossSpectra'
118 118 plot_type = 'pcolor'
119 119 zmin_coh = None
120 120 zmax_coh = None
121 121 zmin_phase = None
122 122 zmax_phase = None
123 123
124 124 def setup(self):
125 125
126 126 self.ncols = 4
127 127 self.nrows = len(self.data.pairs)
128 128 self.nplots = self.nrows * 4
129 129 self.width = 3.4 * self.ncols
130 130 self.height = 3 * self.nrows
131 131 self.ylabel = 'Range [km]'
132 132 self.showprofile = False
133 133
134 134 def plot(self):
135 135
136 136 if self.xaxis == "frequency":
137 137 x = self.data.xrange[0]
138 138 self.xlabel = "Frequency (kHz)"
139 139 elif self.xaxis == "time":
140 140 x = self.data.xrange[1]
141 141 self.xlabel = "Time (ms)"
142 142 else:
143 143 x = self.data.xrange[2]
144 144 self.xlabel = "Velocity (m/s)"
145
145
146 146 self.titles = []
147 147
148 148 y = self.data.heights
149 149 self.y = y
150 150 spc = self.data['spc']
151 151 cspc = self.data['cspc']
152 152
153 153 for n in range(self.nrows):
154 154 noise = self.data['noise'][n][-1]
155 155 pair = self.data.pairs[n]
156 156 ax = self.axes[4 * n]
157 157 spc0 = 10.*numpy.log10(spc[pair[0]]/self.data.factor)
158 if ax.firsttime:
159 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
160 self.xmin = self.xmin if self.xmin else -self.xmax
161 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
162 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
158 if ax.firsttime:
159 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
160 self.xmin = self.xmin if self.xmin else -self.xmax
161 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
162 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
163 163 ax.plt = ax.pcolormesh(x , y , spc0.T,
164 164 vmin=self.zmin,
165 165 vmax=self.zmax,
166 166 cmap=plt.get_cmap(self.colormap)
167 )
168 else:
167 )
168 else:
169 169 ax.plt.set_array(spc0.T.ravel())
170 170 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise))
171 171
172 172 ax = self.axes[4 * n + 1]
173 173 spc1 = 10.*numpy.log10(spc[pair[1]]/self.data.factor)
174 174 if ax.firsttime:
175 175 ax.plt = ax.pcolormesh(x , y, spc1.T,
176 176 vmin=self.zmin,
177 177 vmax=self.zmax,
178 178 cmap=plt.get_cmap(self.colormap)
179 179 )
180 else:
180 else:
181 181 ax.plt.set_array(spc1.T.ravel())
182 182 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise))
183
183
184 184 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
185 185 coh = numpy.abs(out)
186 186 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
187 187
188 188 ax = self.axes[4 * n + 2]
189 189 if ax.firsttime:
190 190 ax.plt = ax.pcolormesh(x, y, coh.T,
191 191 vmin=0,
192 192 vmax=1,
193 193 cmap=plt.get_cmap(self.colormap_coh)
194 194 )
195 195 else:
196 196 ax.plt.set_array(coh.T.ravel())
197 197 self.titles.append(
198 198 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
199
199
200 200 ax = self.axes[4 * n + 3]
201 201 if ax.firsttime:
202 202 ax.plt = ax.pcolormesh(x, y, phase.T,
203 203 vmin=-180,
204 204 vmax=180,
205 cmap=plt.get_cmap(self.colormap_phase)
205 cmap=plt.get_cmap(self.colormap_phase)
206 206 )
207 207 else:
208 208 ax.plt.set_array(phase.T.ravel())
209 209 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
210 210
211 211
212 212 class SpectralMomentsPlot(SpectraPlot):
213 213 '''
214 214 Plot for Spectral Moments
215 215 '''
216 216 CODE = 'spc_moments'
217 217 colormap = 'jro'
218 218 plot_name = 'SpectralMoments'
219 219 plot_type = 'pcolor'
220 220
221 221
222 222 class RTIPlot(Plot):
223 223 '''
224 224 Plot for RTI data
225 225 '''
226 226
227 227 CODE = 'rti'
228 228 colormap = 'jro'
229 229 plot_name = 'RTI'
230 230 plot_type = 'pcolorbuffer'
231 231
232 232 def setup(self):
233 233 self.xaxis = 'time'
234 234 self.ncols = 1
235 235 self.nrows = len(self.data.channels)
236 236 self.nplots = len(self.data.channels)
237 237 self.ylabel = 'Range [km]'
238 238 self.cb_label = 'dB'
239 239 self.titles = ['{} Channel {}'.format(
240 240 self.CODE.upper(), x) for x in range(self.nrows)]
241 241
242 242 def plot(self):
243 243 self.x = self.data.times
244 244 self.y = self.data.heights
245 245 self.z = self.data[self.CODE]
246 246 self.z = numpy.ma.masked_invalid(self.z)
247 247
248 248 if self.decimation is None:
249 249 x, y, z = self.fill_gaps(self.x, self.y, self.z)
250 250 else:
251 251 x, y, z = self.fill_gaps(*self.decimate())
252 252
253 253 for n, ax in enumerate(self.axes):
254 254 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
255 255 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
256 256 if ax.firsttime:
257 257 ax.plt = ax.pcolormesh(x, y, z[n].T,
258 258 vmin=self.zmin,
259 259 vmax=self.zmax,
260 260 cmap=plt.get_cmap(self.colormap)
261 261 )
262 262 if self.showprofile:
263 263 ax.plot_profile = self.pf_axes[n].plot(
264 264 self.data['rti'][n][-1], self.y)[0]
265 265 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
266 266 color="k", linestyle="dashed", lw=1)[0]
267 267 else:
268 268 ax.collections.remove(ax.collections[0])
269 269 ax.plt = ax.pcolormesh(x, y, z[n].T,
270 270 vmin=self.zmin,
271 271 vmax=self.zmax,
272 272 cmap=plt.get_cmap(self.colormap)
273 273 )
274 274 if self.showprofile:
275 275 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
276 276 ax.plot_noise.set_data(numpy.repeat(
277 277 self.data['noise'][n][-1], len(self.y)), self.y)
278 278
279 279
280 280 class CoherencePlot(RTIPlot):
281 281 '''
282 282 Plot for Coherence data
283 283 '''
284 284
285 285 CODE = 'coh'
286 286 plot_name = 'Coherence'
287 287
288 288 def setup(self):
289 289 self.xaxis = 'time'
290 290 self.ncols = 1
291 291 self.nrows = len(self.data.pairs)
292 292 self.nplots = len(self.data.pairs)
293 293 self.ylabel = 'Range [km]'
294 294 if self.CODE == 'coh':
295 295 self.cb_label = ''
296 296 self.titles = [
297 297 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
298 298 else:
299 299 self.cb_label = 'Degrees'
300 300 self.titles = [
301 301 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
302 302
303 303
304 304 class PhasePlot(CoherencePlot):
305 305 '''
306 306 Plot for Phase map data
307 307 '''
308 308
309 309 CODE = 'phase'
310 310 colormap = 'seismic'
311 311 plot_name = 'Phase'
312 312
313 313
314 314 class NoisePlot(Plot):
315 315 '''
316 Plot for noise
316 Plot for noise
317 317 '''
318 318
319 319 CODE = 'noise'
320 320 plot_name = 'Noise'
321 321 plot_type = 'scatterbuffer'
322 322
323 323
324 324 def setup(self):
325 325 self.xaxis = 'time'
326 326 self.ncols = 1
327 327 self.nrows = 1
328 328 self.nplots = 1
329 329 self.ylabel = 'Intensity [dB]'
330 330 self.titles = ['Noise']
331 331 self.colorbar = False
332 332
333 333 def plot(self):
334 334
335 335 x = self.data.times
336 336 xmin = self.data.min_time
337 337 xmax = xmin + self.xrange * 60 * 60
338 338 Y = self.data[self.CODE]
339 339
340 340 if self.axes[0].firsttime:
341 341 for ch in self.data.channels:
342 342 y = Y[ch]
343 343 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
344 344 plt.legend()
345 345 else:
346 346 for ch in self.data.channels:
347 347 y = Y[ch]
348 348 self.axes[0].lines[ch].set_data(x, y)
349 349
350 350 self.ymin = numpy.nanmin(Y) - 5
351 351 self.ymax = numpy.nanmax(Y) + 5
352 352
353 353
354 354 class SnrPlot(RTIPlot):
355 355 '''
356 356 Plot for SNR Data
357 357 '''
358 358
359 359 CODE = 'snr'
360 360 colormap = 'jet'
361 361 plot_name = 'SNR'
362 362
363 363
364 364 class DopplerPlot(RTIPlot):
365 365 '''
366 366 Plot for DOPPLER Data (1st moment)
367 367 '''
368 368
369 369 CODE = 'dop'
370 370 colormap = 'jet'
371 371 plot_name = 'DopplerShift'
372 372
373 373
374 374 class PowerPlot(RTIPlot):
375 375 '''
376 376 Plot for Power Data (0 moment)
377 377 '''
378 378
379 379 CODE = 'pow'
380 380 colormap = 'jet'
381 381 plot_name = 'TotalPower'
382 382
383 383
384 384 class SpectralWidthPlot(RTIPlot):
385 385 '''
386 386 Plot for Spectral Width Data (2nd moment)
387 387 '''
388 388
389 389 CODE = 'width'
390 390 colormap = 'jet'
391 391 plot_name = 'SpectralWidth'
392 392
393 393
394 394 class SkyMapPlot(Plot):
395 395 '''
396 396 Plot for meteors detection data
397 397 '''
398 398
399 399 CODE = 'param'
400 400
401 401 def setup(self):
402 402
403 403 self.ncols = 1
404 404 self.nrows = 1
405 405 self.width = 7.2
406 406 self.height = 7.2
407 407 self.nplots = 1
408 408 self.xlabel = 'Zonal Zenith Angle (deg)'
409 409 self.ylabel = 'Meridional Zenith Angle (deg)'
410 410 self.polar = True
411 411 self.ymin = -180
412 412 self.ymax = 180
413 413 self.colorbar = False
414 414
415 415 def plot(self):
416 416
417 417 arrayParameters = numpy.concatenate(self.data['param'])
418 418 error = arrayParameters[:, -1]
419 419 indValid = numpy.where(error == 0)[0]
420 420 finalMeteor = arrayParameters[indValid, :]
421 421 finalAzimuth = finalMeteor[:, 3]
422 422 finalZenith = finalMeteor[:, 4]
423 423
424 424 x = finalAzimuth * numpy.pi / 180
425 425 y = finalZenith
426 426
427 427 ax = self.axes[0]
428 428
429 429 if ax.firsttime:
430 430 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
431 431 else:
432 432 ax.plot.set_data(x, y)
433 433
434 434 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
435 435 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
436 436 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
437 437 dt2,
438 438 len(x))
439 439 self.titles[0] = title
440 440
441 441
442 442 class ParametersPlot(RTIPlot):
443 443 '''
444 444 Plot for data_param object
445 445 '''
446 446
447 447 CODE = 'param'
448 448 colormap = 'seismic'
449 449 plot_name = 'Parameters'
450 450
451 451 def setup(self):
452 452 self.xaxis = 'time'
453 453 self.ncols = 1
454 454 self.nrows = self.data.shape(self.CODE)[0]
455 455 self.nplots = self.nrows
456 456 if self.showSNR:
457 457 self.nrows += 1
458 458 self.nplots += 1
459 459
460 460 self.ylabel = 'Height [km]'
461 461 if not self.titles:
462 462 self.titles = self.data.parameters \
463 463 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
464 464 if self.showSNR:
465 465 self.titles.append('SNR')
466 466
467 467 def plot(self):
468 468 self.data.normalize_heights()
469 469 self.x = self.data.times
470 470 self.y = self.data.heights
471 471 if self.showSNR:
472 472 self.z = numpy.concatenate(
473 473 (self.data[self.CODE], self.data['snr'])
474 474 )
475 475 else:
476 476 self.z = self.data[self.CODE]
477 477
478 478 self.z = numpy.ma.masked_invalid(self.z)
479 479
480 480 if self.decimation is None:
481 481 x, y, z = self.fill_gaps(self.x, self.y, self.z)
482 482 else:
483 483 x, y, z = self.fill_gaps(*self.decimate())
484 484
485 485 for n, ax in enumerate(self.axes):
486 486
487 487 self.zmax = self.zmax if self.zmax is not None else numpy.max(
488 488 self.z[n])
489 489 self.zmin = self.zmin if self.zmin is not None else numpy.min(
490 490 self.z[n])
491 491
492 492 if ax.firsttime:
493 493 if self.zlimits is not None:
494 494 self.zmin, self.zmax = self.zlimits[n]
495 495
496 496 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
497 497 vmin=self.zmin,
498 498 vmax=self.zmax,
499 499 cmap=self.cmaps[n]
500 500 )
501 501 else:
502 502 if self.zlimits is not None:
503 503 self.zmin, self.zmax = self.zlimits[n]
504 504 ax.collections.remove(ax.collections[0])
505 505 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
506 506 vmin=self.zmin,
507 507 vmax=self.zmax,
508 508 cmap=self.cmaps[n]
509 509 )
510 510
511 511
512 512 class OutputPlot(ParametersPlot):
513 513 '''
514 514 Plot data_output object
515 515 '''
516 516
517 517 CODE = 'output'
518 518 colormap = 'seismic'
519 519 plot_name = 'Output'
520 520
521 521
522 522 class PolarMapPlot(Plot):
523 523 '''
524 524 Plot for weather radar
525 525 '''
526 526
527 527 CODE = 'param'
528 528 colormap = 'seismic'
529 529
530 530 def setup(self):
531 531 self.ncols = 1
532 532 self.nrows = 1
533 533 self.width = 9
534 534 self.height = 8
535 535 self.mode = self.data.meta['mode']
536 536 if self.channels is not None:
537 537 self.nplots = len(self.channels)
538 538 self.nrows = len(self.channels)
539 539 else:
540 540 self.nplots = self.data.shape(self.CODE)[0]
541 541 self.nrows = self.nplots
542 542 self.channels = list(range(self.nplots))
543 543 if self.mode == 'E':
544 544 self.xlabel = 'Longitude'
545 545 self.ylabel = 'Latitude'
546 546 else:
547 547 self.xlabel = 'Range (km)'
548 548 self.ylabel = 'Height (km)'
549 549 self.bgcolor = 'white'
550 550 self.cb_labels = self.data.meta['units']
551 551 self.lat = self.data.meta['latitude']
552 552 self.lon = self.data.meta['longitude']
553 553 self.xmin, self.xmax = float(
554 554 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
555 555 self.ymin, self.ymax = float(
556 556 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
557 557 # self.polar = True
558 558
559 559 def plot(self):
560 560
561 561 for n, ax in enumerate(self.axes):
562 562 data = self.data['param'][self.channels[n]]
563 563
564 564 zeniths = numpy.linspace(
565 565 0, self.data.meta['max_range'], data.shape[1])
566 566 if self.mode == 'E':
567 567 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
568 568 r, theta = numpy.meshgrid(zeniths, azimuths)
569 569 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
570 570 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
571 571 x = km2deg(x) + self.lon
572 572 y = km2deg(y) + self.lat
573 573 else:
574 574 azimuths = numpy.radians(self.data.heights)
575 575 r, theta = numpy.meshgrid(zeniths, azimuths)
576 576 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
577 577 self.y = zeniths
578 578
579 579 if ax.firsttime:
580 580 if self.zlimits is not None:
581 581 self.zmin, self.zmax = self.zlimits[n]
582 582 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
583 583 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
584 584 vmin=self.zmin,
585 585 vmax=self.zmax,
586 586 cmap=self.cmaps[n])
587 587 else:
588 588 if self.zlimits is not None:
589 589 self.zmin, self.zmax = self.zlimits[n]
590 590 ax.collections.remove(ax.collections[0])
591 591 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
592 592 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
593 593 vmin=self.zmin,
594 594 vmax=self.zmax,
595 595 cmap=self.cmaps[n])
596 596
597 597 if self.mode == 'A':
598 598 continue
599 599
600 600 # plot district names
601 601 f = open('/data/workspace/schain_scripts/distrito.csv')
602 602 for line in f:
603 603 label, lon, lat = [s.strip() for s in line.split(',') if s]
604 604 lat = float(lat)
605 605 lon = float(lon)
606 606 # ax.plot(lon, lat, '.b', ms=2)
607 607 ax.text(lon, lat, label.decode('utf8'), ha='center',
608 608 va='bottom', size='8', color='black')
609 609
610 610 # plot limites
611 611 limites = []
612 612 tmp = []
613 613 for line in open('/data/workspace/schain_scripts/lima.csv'):
614 614 if '#' in line:
615 615 if tmp:
616 616 limites.append(tmp)
617 617 tmp = []
618 618 continue
619 619 values = line.strip().split(',')
620 620 tmp.append((float(values[0]), float(values[1])))
621 621 for points in limites:
622 622 ax.add_patch(
623 623 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
624 624
625 625 # plot Cuencas
626 626 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
627 627 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
628 628 values = [line.strip().split(',') for line in f]
629 629 points = [(float(s[0]), float(s[1])) for s in values]
630 630 ax.add_patch(Polygon(points, ec='b', fc='none'))
631 631
632 632 # plot grid
633 633 for r in (15, 30, 45, 60):
634 634 ax.add_artist(plt.Circle((self.lon, self.lat),
635 635 km2deg(r), color='0.6', fill=False, lw=0.2))
636 636 ax.text(
637 637 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
638 638 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
639 639 '{}km'.format(r),
640 640 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
641 641
642 642 if self.mode == 'E':
643 643 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
644 644 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
645 645 else:
646 646 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
647 647 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
648 648
649 649 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
650 650 self.titles = ['{} {}'.format(
651 651 self.data.parameters[x], title) for x in self.channels]
652 652
653 653
654 654 class ScopePlot(Plot):
655 655
656 656 '''
657 657 Plot for Scope
658 '''
658 '''
659 659
660 660 CODE = 'scope'
661 661 plot_name = 'Scope'
662 662 plot_type = 'scatter'
663
663
664 664 def setup(self):
665 665
666 666 self.xaxis = 'Range (Km)'
667 667 self.ncols = 1
668 668 self.nrows = 1
669 669 self.nplots = 1
670 670 self.ylabel = 'Intensity [dB]'
671 671 self.titles = ['Scope']
672 672 self.colorbar = False
673 673 colspan = 3
674 674 rowspan = 1
675 675
676 676 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
677
677
678 678 yreal = y[channelIndexList,:].real
679 679 yimag = y[channelIndexList,:].imag
680 680 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
681 681 self.xlabel = "Range (Km)"
682 682 self.ylabel = "Intensity - IQ"
683
683
684 684 self.y = yreal
685 685 self.x = x
686 686 self.xmin = min(x)
687 687 self.xmax = max(x)
688
689 688
690 self.titles[0] = title
689
690 self.titles[0] = title
691 691
692 692 for i,ax in enumerate(self.axes):
693 693 title = "Channel %d" %(i)
694 694 if ax.firsttime:
695 695 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
696 696 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
697 697 else:
698 698 #pass
699 699 ax.plt_r.set_data(x, yreal[i,:])
700 700 ax.plt_i.set_data(x, yimag[i,:])
701
701
702 702 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
703 703 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
704 704 yreal = y.real
705 705 self.y = yreal
706 706 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
707 707 self.xlabel = "Range (Km)"
708 708 self.ylabel = "Intensity"
709 709 self.xmin = min(x)
710 710 self.xmax = max(x)
711
712
711
712
713 713 self.titles[0] = title
714 714
715 715 for i,ax in enumerate(self.axes):
716 716 title = "Channel %d" %(i)
717
717
718 718 ychannel = yreal[i,:]
719
720 if ax.firsttime:
719
720 if ax.firsttime:
721 721 ax.plt_r = ax.plot(x, ychannel)[0]
722 722 else:
723 723 #pass
724 724 ax.plt_r.set_data(x, ychannel)
725
725
726 726
727 727 def plot(self):
728
728
729 729 if self.channels:
730 730 channels = self.channels
731 731 else:
732 732 channels = self.data.channels
733 733
734
735
734
735
736 736 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
737
737
738 738 scope = self.data['scope']
739
740
739
740
741 741 if self.data.flagDataAsBlock:
742
742
743 743 for i in range(self.data.nProfiles):
744 744
745 745 wintitle1 = " [Profile = %d] " %i
746 746
747 747 if self.type == "power":
748 self.plot_power(self.data.heights,
749 scope[:,i,:],
748 self.plot_power(self.data.heights,
749 scope[:,i,:],
750 750 channels,
751 thisDatetime,
751 thisDatetime,
752 752 wintitle1
753 753 )
754 754
755 755 if self.type == "iq":
756 self.plot_iq(self.data.heights,
756 self.plot_iq(self.data.heights,
757 757 scope[:,i,:],
758 758 channels,
759 thisDatetime,
760 wintitle1
759 thisDatetime,
760 wintitle1
761 761 )
762 762 else:
763 763 wintitle = " [Profile = %d] " %self.data.profileIndex
764
764
765 765 if self.type == "power":
766 self.plot_power(self.data.heights,
766 self.plot_power(self.data.heights,
767 767 scope,
768 768 channels,
769 769 thisDatetime,
770 770 wintitle
771 771 )
772
772
773 773 if self.type == "iq":
774 self.plot_iq(self.data.heights,
774 self.plot_iq(self.data.heights,
775 775 scope,
776 776 channels,
777 777 thisDatetime,
778 778 wintitle
779 779 )
Show file before | Show file after | Hide comments
@@ -1,649 +1,650
1 1 '''
2 2 Created on Set 9, 2015
3 3
4 4 @author: roj-idl71 Karim Kuyeng
5 5 '''
6 6
7 7 import os
8 8 import sys
9 9 import glob
10 10 import fnmatch
11 11 import datetime
12 12 import time
13 13 import re
14 14 import h5py
15 15 import numpy
16 16
17 17 try:
18 18 from gevent import sleep
19 19 except:
20 20 from time import sleep
21 21
22 22 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
23 23 from schainpy.model.data.jrodata import Voltage
24 24 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
25 25 from numpy import imag
26 26
27 27 @MPDecorator
28 28 class AMISRReader(ProcessingUnit):
29 29 '''
30 30 classdocs
31 31 '''
32 32
33 33 def __init__(self):
34 34 '''
35 35 Constructor
36 36 '''
37 37
38 38 ProcessingUnit.__init__(self)
39 39
40 40 self.set = None
41 41 self.subset = None
42 42 self.extension_file = '.h5'
43 43 self.dtc_str = 'dtc'
44 44 self.dtc_id = 0
45 45 self.status = True
46 46 self.isConfig = False
47 47 self.dirnameList = []
48 48 self.filenameList = []
49 49 self.fileIndex = None
50 50 self.flagNoMoreFiles = False
51 51 self.flagIsNewFile = 0
52 52 self.filename = ''
53 53 self.amisrFilePointer = None
54 54
55 55
56 56 #self.dataset = None
57 57
58 58
59 59
60 60
61 61 self.profileIndex = 0
62 62
63 63
64 64 self.beamCodeByFrame = None
65 65 self.radacTimeByFrame = None
66 66
67 67 self.dataset = None
68 68
69 69
70 70
71 71
72 72 self.__firstFile = True
73 73
74 74 self.buffer = None
75 75
76 76
77 77 self.timezone = 'ut'
78 78
79 79 self.__waitForNewFile = 20
80 80 self.__filename_online = None
81 81 #Is really necessary create the output object in the initializer
82 82 self.dataOut = Voltage()
83 83 self.dataOut.error=False
84 84
85 85 def setup(self,path=None,
86 86 startDate=None,
87 87 endDate=None,
88 88 startTime=None,
89 89 endTime=None,
90 90 walk=True,
91 91 timezone='ut',
92 92 all=0,
93 93 code = None,
94 94 nCode = 0,
95 95 nBaud = 0,
96 96 online=False):
97 97
98 98 #print ("T",path)
99 99
100 100 self.timezone = timezone
101 101 self.all = all
102 102 self.online = online
103 103
104 104 self.code = code
105 105 self.nCode = int(nCode)
106 106 self.nBaud = int(nBaud)
107 107
108 108
109 109
110 110 #self.findFiles()
111 111 if not(online):
112 112 #Busqueda de archivos offline
113 113 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk)
114 114 else:
115 115 self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
116 116
117 117 if not(self.filenameList):
118 118 print("There is no files into the folder: %s"%(path))
119 119 sys.exit(-1)
120 120
121 121 self.fileIndex = -1
122 122
123 123 self.readNextFile(online)
124 124
125 125 '''
126 126 Add code
127 127 '''
128 128 self.isConfig = True
129 129
130 130 pass
131 131
132 132
133 133 def readAMISRHeader(self,fp):
134 134 header = 'Raw11/Data/RadacHeader'
135 135 self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
136 136 self.beamCode = fp.get('Raw11/Data/Beamcodes') # NUMBER OF CHANNELS AND IDENTIFY POSITION TO CREATE A FILE WITH THAT INFO
137 137 #self.code = fp.get(header+'/Code') # NOT USE FOR THIS
138 138 self.frameCount = fp.get(header+'/FrameCount')# NOT USE FOR THIS
139 139 self.modeGroup = fp.get(header+'/ModeGroup')# NOT USE FOR THIS
140 140 self.nsamplesPulse = fp.get(header+'/NSamplesPulse')# TO GET NSA OR USING DATA FOR THAT
141 141 self.pulseCount = fp.get(header+'/PulseCount')# NOT USE FOR THIS
142 142 self.radacTime = fp.get(header+'/RadacTime')# 1st TIME ON FILE ANDE CALCULATE THE REST WITH IPP*nindexprofile
143 143 self.timeCount = fp.get(header+'/TimeCount')# NOT USE FOR THIS
144 144 self.timeStatus = fp.get(header+'/TimeStatus')# NOT USE FOR THIS
145 145 self.rangeFromFile = fp.get('Raw11/Data/Samples/Range')
146 146 self.frequency = fp.get('Rx/Frequency')
147 147 txAus = fp.get('Raw11/Data/Pulsewidth')
148 148
149 149
150 150 self.nblocks = self.pulseCount.shape[0] #nblocks
151 151
152 152 self.nprofiles = self.pulseCount.shape[1] #nprofile
153 153 self.nsa = self.nsamplesPulse[0,0] #ngates
154 154 self.nchannels = self.beamCode.shape[1]
155 155 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
156 156 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
157 157 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
158 158
159 159 #filling radar controller header parameters
160 160 self.__ippKm = self.ippSeconds *.15*1e6 # in km
161 161 self.__txA = (txAus.value)*.15 #(ipp[us]*.15km/1us) in km
162 162 self.__txB = 0
163 163 nWindows=1
164 164 self.__nSamples = self.nsa
165 165 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
166 166 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
167 167
168 168 #for now until understand why the code saved is different (code included even though code not in tuf file)
169 169 #self.__codeType = 0
170 170 # self.__nCode = None
171 171 # self.__nBaud = None
172 172 self.__code = self.code
173 173 self.__codeType = 0
174 174 if self.code != None:
175 175 self.__codeType = 1
176 176 self.__nCode = self.nCode
177 177 self.__nBaud = self.nBaud
178 178 #self.__code = 0
179 179
180 180 #filling system header parameters
181 181 self.__nSamples = self.nsa
182 182 self.newProfiles = self.nprofiles/self.nchannels
183 183 self.__channelList = list(range(self.nchannels))
184 184
185 185 self.__frequency = self.frequency[0][0]
186 186
187 187
188 188
189 189 def createBuffers(self):
190 190
191 191 pass
192 192
193 193 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
194 194 self.path = path
195 195 self.startDate = startDate
196 196 self.endDate = endDate
197 197 self.startTime = startTime
198 198 self.endTime = endTime
199 199 self.walk = walk
200 200
201 201 def __checkPath(self):
202 202 if os.path.exists(self.path):
203 203 self.status = 1
204 204 else:
205 205 self.status = 0
206 206 print('Path:%s does not exists'%self.path)
207 207
208 208 return
209 209
210 210
211 211 def __selDates(self, amisr_dirname_format):
212 212 try:
213 213 year = int(amisr_dirname_format[0:4])
214 214 month = int(amisr_dirname_format[4:6])
215 215 dom = int(amisr_dirname_format[6:8])
216 216 thisDate = datetime.date(year,month,dom)
217 217
218 218 if (thisDate>=self.startDate and thisDate <= self.endDate):
219 219 return amisr_dirname_format
220 220 except:
221 221 return None
222 222
223 223
224 224 def __findDataForDates(self,online=False):
225 225
226 226 if not(self.status):
227 227 return None
228 228
229 229 pat = '\d+.\d+'
230 230 dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
231 231 dirnameList = [x for x in dirnameList if x!=None]
232 232 dirnameList = [x.string for x in dirnameList]
233 233 if not(online):
234 234 dirnameList = [self.__selDates(x) for x in dirnameList]
235 235 dirnameList = [x for x in dirnameList if x!=None]
236 236 if len(dirnameList)>0:
237 237 self.status = 1
238 238 self.dirnameList = dirnameList
239 239 self.dirnameList.sort()
240 240 else:
241 241 self.status = 0
242 242 return None
243 243
244 244 def __getTimeFromData(self):
245 245 startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
246 246 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
247 247
248 248 print('Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader))
249 249 print('........................................')
250 250 filter_filenameList = []
251 251 self.filenameList.sort()
252 252 #for i in range(len(self.filenameList)-1):
253 253 for i in range(len(self.filenameList)):
254 254 filename = self.filenameList[i]
255 255 fp = h5py.File(filename,'r')
256 256 time_str = fp.get('Time/RadacTimeString')
257 257
258 258 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
259 259 #startDateTimeStr_File = "2019-12-16 09:21:11"
260 260 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
261 261 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
262 262
263 263 #endDateTimeStr_File = "2019-12-16 11:10:11"
264 264 endDateTimeStr_File = time_str[-1][-1].decode('UTF-8').split('.')[0]
265 265 junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
266 266 endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
267 267
268 268 fp.close()
269 269
270 270 #print("check time", startDateTime_File)
271 271 if self.timezone == 'lt':
272 272 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
273 273 endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
274 274 if (endDateTime_File>=startDateTime_Reader and endDateTime_File<endDateTime_Reader):
275 275 #self.filenameList.remove(filename)
276 276 filter_filenameList.append(filename)
277 277
278 278 if (endDateTime_File>=endDateTime_Reader):
279 279 break
280 280
281 281
282 282 filter_filenameList.sort()
283 283 self.filenameList = filter_filenameList
284 284 return 1
285 285
286 286 def __filterByGlob1(self, dirName):
287 287 filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
288 288 filter_files.sort()
289 289 filterDict = {}
290 290 filterDict.setdefault(dirName)
291 291 filterDict[dirName] = filter_files
292 292 return filterDict
293 293
294 294 def __getFilenameList(self, fileListInKeys, dirList):
295 295 for value in fileListInKeys:
296 296 dirName = list(value.keys())[0]
297 297 for file in value[dirName]:
298 298 filename = os.path.join(dirName, file)
299 299 self.filenameList.append(filename)
300 300
301 301
302 302 def __selectDataForTimes(self, online=False):
303 303 #aun no esta implementado el filtro for tiempo
304 304 if not(self.status):
305 305 return None
306 306
307 307 dirList = [os.path.join(self.path,x) for x in self.dirnameList]
308 308
309 309 fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
310 310
311 311 self.__getFilenameList(fileListInKeys, dirList)
312 312 if not(online):
313 313 #filtro por tiempo
314 314 if not(self.all):
315 315 self.__getTimeFromData()
316 316
317 317 if len(self.filenameList)>0:
318 318 self.status = 1
319 319 self.filenameList.sort()
320 320 else:
321 321 self.status = 0
322 322 return None
323 323
324 324 else:
325 325 #get the last file - 1
326 326 self.filenameList = [self.filenameList[-2]]
327
328 327 new_dirnameList = []
329 328 for dirname in self.dirnameList:
330 329 junk = numpy.array([dirname in x for x in self.filenameList])
331 330 junk_sum = junk.sum()
332 331 if junk_sum > 0:
333 332 new_dirnameList.append(dirname)
334 333 self.dirnameList = new_dirnameList
335 334 return 1
336 335
337 336 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
338 337 endTime=datetime.time(23,59,59),walk=True):
339 338
340 339 if endDate ==None:
341 340 startDate = datetime.datetime.utcnow().date()
342 341 endDate = datetime.datetime.utcnow().date()
343 342
344 343 self.__setParameters(path=path, startDate=startDate, endDate=endDate,startTime = startTime,endTime=endTime, walk=walk)
345 344
346 345 self.__checkPath()
347 346
348 347 self.__findDataForDates(online=True)
349 348
350 349 self.dirnameList = [self.dirnameList[-1]]
351 350
352 351 self.__selectDataForTimes(online=True)
353 352
354 353 return
355 354
356 355
357 356 def searchFilesOffLine(self,
358 357 path,
359 358 startDate,
360 359 endDate,
361 360 startTime=datetime.time(0,0,0),
362 361 endTime=datetime.time(23,59,59),
363 362 walk=True):
364 363
365 364 self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
366 365
367 366 self.__checkPath()
368 367
369 368 self.__findDataForDates()
370 369
371 370 self.__selectDataForTimes()
372 371
373 372 for i in range(len(self.filenameList)):
374 373 print("%s" %(self.filenameList[i]))
375 374
376 375 return
377 376
378 377 def __setNextFileOffline(self):
379 378 idFile = self.fileIndex
380 379
381 380 while (True):
382 381 idFile += 1
383 382 if not(idFile < len(self.filenameList)):
384 383 self.flagNoMoreFiles = 1
385 384 print("No more Files")
386 self.dataOut.error = True
387 385 return 0
388 386
389 387 filename = self.filenameList[idFile]
390 388
391 389 amisrFilePointer = h5py.File(filename,'r')
392 390
393 391 break
394 392
395 393 self.flagIsNewFile = 1
396 394 self.fileIndex = idFile
397 395 self.filename = filename
398 396
399 397 self.amisrFilePointer = amisrFilePointer
400 398
401 399 print("Setting the file: %s"%self.filename)
402 400
403 401 return 1
404 402
405 403
406 404 def __setNextFileOnline(self):
407 405 filename = self.filenameList[0]
408 406 if self.__filename_online != None:
409 407 self.__selectDataForTimes(online=True)
410 408 filename = self.filenameList[0]
411 409 wait = 0
410 #self.__waitForNewFile=5 ## DEBUG:
412 411 while self.__filename_online == filename:
413 412 print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
414 413 if wait == 5:
414 self.flagNoMoreFiles = 1
415 415 return 0
416 416 sleep(self.__waitForNewFile)
417 417 self.__selectDataForTimes(online=True)
418 418 filename = self.filenameList[0]
419 419 wait += 1
420 420
421 421 self.__filename_online = filename
422 422
423 423 self.amisrFilePointer = h5py.File(filename,'r')
424 424 self.flagIsNewFile = 1
425 425 self.filename = filename
426 426 print("Setting the file: %s"%self.filename)
427 427 return 1
428 428
429 429
430 430 def readData(self):
431 431 buffer = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
432 432 re = buffer[:,:,:,0]
433 433 im = buffer[:,:,:,1]
434 434 dataset = re + im*1j
435 435
436 436 self.radacTime = self.amisrFilePointer.get('Raw11/Data/RadacHeader/RadacTime')
437 437 timeset = self.radacTime[:,0]
438 438
439 439 return dataset,timeset
440 440
441 441 def reshapeData(self):
442 442 #self.beamCodeByPulse, self.beamCode, self.nblocks, self.nprofiles, self.nsa,
443 443 channels = self.beamCodeByPulse[0,:]
444 444 nchan = self.nchannels
445 445 #self.newProfiles = self.nprofiles/nchan #must be defined on filljroheader
446 446 nblocks = self.nblocks
447 447 nsamples = self.nsa
448 448
449 449 #Dimensions : nChannels, nProfiles, nSamples
450 450 new_block = numpy.empty((nblocks, nchan, numpy.int_(self.newProfiles), nsamples), dtype="complex64")
451 451 ############################################
452 452
453 453 for thisChannel in range(nchan):
454 454 new_block[:,thisChannel,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[0][thisChannel])[0],:]
455 455
456 456
457 457 new_block = numpy.transpose(new_block, (1,0,2,3))
458 458 new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
459 459
460 460 return new_block
461 461
462 462 def updateIndexes(self):
463 463
464 464 pass
465 465
466 466 def fillJROHeader(self):
467 467
468 468 #fill radar controller header
469 469 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(ipp=self.__ippKm,
470 470 txA=self.__txA,
471 471 txB=0,
472 472 nWindows=1,
473 473 nHeights=self.__nSamples,
474 474 firstHeight=self.__firstHeight,
475 475 deltaHeight=self.__deltaHeight,
476 476 codeType=self.__codeType,
477 477 nCode=self.__nCode, nBaud=self.__nBaud,
478 478 code = self.__code,
479 479 fClock=1)
480 480
481 481 #fill system header
482 482 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
483 483 nProfiles=self.newProfiles,
484 484 nChannels=len(self.__channelList),
485 485 adcResolution=14,
486 486 pciDioBusWidth=32)
487 487
488 488 self.dataOut.type = "Voltage"
489 489
490 490 self.dataOut.data = None
491 491
492 492 self.dataOut.dtype = numpy.dtype([('real','<i8'),('imag','<i8')])
493 493
494 494 # self.dataOut.nChannels = 0
495 495
496 496 # self.dataOut.nHeights = 0
497 497
498 498 self.dataOut.nProfiles = self.newProfiles*self.nblocks
499 499
500 500 #self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
501 501 ranges = numpy.reshape(self.rangeFromFile.value,(-1))
502 502 self.dataOut.heightList = ranges/1000.0 #km
503 503
504 504
505 505 self.dataOut.channelList = self.__channelList
506 506
507 507 self.dataOut.blocksize = self.dataOut.getNChannels() * self.dataOut.getNHeights()
508 508
509 509 # self.dataOut.channelIndexList = None
510 510
511 511 self.dataOut.flagNoData = True
512 512
513 513 #Set to TRUE if the data is discontinuous
514 514 self.dataOut.flagDiscontinuousBlock = False
515 515
516 516 self.dataOut.utctime = None
517 517
518 518 #self.dataOut.timeZone = -5 #self.__timezone/60 #timezone like jroheader, difference in minutes between UTC and localtime
519 519 if self.timezone == 'lt':
520 520 self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
521 521 else:
522 522 self.dataOut.timeZone = 0 #by default time is UTC
523 523
524 524 self.dataOut.dstFlag = 0
525 525
526 526 self.dataOut.errorCount = 0
527 527
528 528 self.dataOut.nCohInt = 1
529 529
530 530 self.dataOut.flagDecodeData = False #asumo que la data esta decodificada
531 531
532 532 self.dataOut.flagDeflipData = False #asumo que la data esta sin flip
533 533
534 534 self.dataOut.flagShiftFFT = False
535 535
536 536 self.dataOut.ippSeconds = self.ippSeconds
537 537
538 538 #Time interval between profiles
539 539 #self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
540 540
541 541 self.dataOut.frequency = self.__frequency
542 542 self.dataOut.realtime = self.online
543 543 pass
544 544
545 545 def readNextFile(self,online=False):
546 546
547 547 if not(online):
548 548 newFile = self.__setNextFileOffline()
549 549 else:
550 550 newFile = self.__setNextFileOnline()
551 551
552 552 if not(newFile):
553 self.dataOut.error = True
553 554 return 0
554 555 #if self.__firstFile:
555 556 self.readAMISRHeader(self.amisrFilePointer)
556 557
557 558 self.createBuffers()
558 559
559 560 self.fillJROHeader()
560 561
561 562 #self.__firstFile = False
562 563
563 564
564 565
565 566 self.dataset,self.timeset = self.readData()
566 567
567 568 if self.endDate!=None:
568 569 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
569 570 time_str = self.amisrFilePointer.get('Time/RadacTimeString')
570 571 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
571 572 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
572 573 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
573 574 if self.timezone == 'lt':
574 575 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
575 576 if (startDateTime_File>endDateTime_Reader):
576 577 return 0
577 578
578 579 self.jrodataset = self.reshapeData()
579 580 #----self.updateIndexes()
580 581 self.profileIndex = 0
581 582
582 583 return 1
583 584
584 585
585 586 def __hasNotDataInBuffer(self):
586 587 if self.profileIndex >= (self.newProfiles*self.nblocks):
587 588 return 1
588 589 return 0
589 590
590 591
591 592 def getData(self):
592 593
593 594 if self.flagNoMoreFiles:
594 595 self.dataOut.flagNoData = True
595 596 return 0
596 597
597 598 if self.__hasNotDataInBuffer():
598 599 if not (self.readNextFile(self.online)):
599 600 return 0
600 601
601 602
602 603 if self.dataset is None: # setear esta condicion cuando no hayan datos por leer
603 604 self.dataOut.flagNoData = True
604 605 return 0
605 606
606 607 #self.dataOut.data = numpy.reshape(self.jrodataset[self.profileIndex,:],(1,-1))
607 608
608 609 self.dataOut.data = self.jrodataset[:,self.profileIndex,:]
609 610
610 611 #print("R_t",self.timeset)
611 612
612 613 #self.dataOut.utctime = self.jrotimeset[self.profileIndex]
613 614 #verificar basic header de jro data y ver si es compatible con este valor
614 615 #self.dataOut.utctime = self.timeset + (self.profileIndex * self.ippSeconds * self.nchannels)
615 616 indexprof = numpy.mod(self.profileIndex, self.newProfiles)
616 617 indexblock = self.profileIndex/self.newProfiles
617 618 #print (indexblock, indexprof)
618 619 diffUTC = 1.8e4 #UTC diference from peru in seconds --Joab
619 620 diffUTC = 0
620 621 t_comp = (indexprof * self.ippSeconds * self.nchannels) + diffUTC #
621 622 #cambio posible 18/02/2020
622 623
623 624
624 625
625 626 #print("utc :",indexblock," __ ",t_comp)
626 627 #print(numpy.shape(self.timeset))
627 628 self.dataOut.utctime = self.timeset[numpy.int_(indexblock)] + t_comp
628 629 #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
629 630 #print(self.dataOut.utctime)
630 631 self.dataOut.profileIndex = self.profileIndex
631 632 self.dataOut.flagNoData = False
632 633 # if indexprof == 0:
633 634 # print self.dataOut.utctime
634 635
635 636 self.profileIndex += 1
636 637
637 638 return self.dataOut.data
638 639
639 640
640 641 def run(self, **kwargs):
641 642 '''
642 643 This method will be called many times so here you should put all your code
643 644 '''
644 645
645 646 if not self.isConfig:
646 647 self.setup(**kwargs)
647 648 self.isConfig = True
648 649
649 650 self.getData()
General Comments 0
You need to be logged in to leave comments. Login now