##// END OF EJS Templates
schain
RSS
ax.plt.remove(), se añaden 3 lineas de codigo
avaldez -
r1785:6f1a660b8ab2
parent child
Show More
Side by Side Unified
Context file:
r1785:6f1a660b8ab2 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,494 +1,499
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
7 7 from schainpy.utils import log
8 8
9 9 EARTH_RADIUS = 6.3710e3
10 10
11 11
12 12 def ll2xy(lat1, lon1, lat2, lon2):
13 13
14 14 p = 0.017453292519943295
15 15 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
16 16 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
17 17 r = 12742 * numpy.arcsin(numpy.sqrt(a))
18 18 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
19 19 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
20 20 theta = -theta + numpy.pi/2
21 21 return r*numpy.cos(theta), r*numpy.sin(theta)
22 22
23 23
24 24 def km2deg(km):
25 25 '''
26 26 Convert distance in km to degrees
27 27 '''
28 28
29 29 return numpy.rad2deg(km/EARTH_RADIUS)
30 30
31 31
32 32
33 33 class SpectralMomentsPlot(SpectraPlot):
34 34 '''
35 35 Plot for Spectral Moments
36 36 '''
37 37 CODE = 'spc_moments'
38 38 # colormap = 'jet'
39 39 # plot_type = 'pcolor'
40 40
41 41 class DobleGaussianPlot(SpectraPlot):
42 42 '''
43 43 Plot for Double Gaussian Plot
44 44 '''
45 45 CODE = 'gaussian_fit'
46 46 # colormap = 'jet'
47 47 # plot_type = 'pcolor'
48 48
49 49
50 50 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
51 51 '''
52 52 Plot SpectraCut with Double Gaussian Fit
53 53 '''
54 54 CODE = 'cut_gaussian_fit'
55 55
56 56
57 57 class SpectralFitObliquePlot(SpectraPlot):
58 58 '''
59 59 Plot for Spectral Oblique
60 60 '''
61 61 CODE = 'spc_moments'
62 62 colormap = 'jet'
63 63 plot_type = 'pcolor'
64 64
65 65
66 66
67 67 class SnrPlot(RTIPlot):
68 68 '''
69 69 Plot for SNR Data
70 70 '''
71 71
72 72 CODE = 'snr'
73 73 colormap = 'jet'
74 74
75 75 def update(self, dataOut):
76 76
77 77 data = {
78 78 'snr': 10*numpy.log10(dataOut.data_snr)
79 79 }
80 80
81 81 return data, {}
82 82
83 83 class DopplerPlot(RTIPlot):
84 84 '''
85 85 Plot for DOPPLER Data (1st moment)
86 86 '''
87 87
88 88 CODE = 'dop'
89 89 colormap = 'jet'
90 90
91 91 def update(self, dataOut):
92 92
93 93 data = {
94 94 'dop': 10*numpy.log10(dataOut.data_dop)
95 95 }
96 96
97 97 return data, {}
98 98
99 99 class DopplerEEJPlot_V0(RTIPlot):
100 100 '''
101 101 Written by R. Flores
102 102 '''
103 103 '''
104 104 Plot for EEJ
105 105 '''
106 106
107 107 CODE = 'dop'
108 108 colormap = 'RdBu_r'
109 109 colormap = 'jet'
110 110
111 111 def setup(self):
112 112
113 113 self.xaxis = 'time'
114 114 self.ncols = 1
115 115 self.nrows = len(self.data.channels)
116 116 self.nplots = len(self.data.channels)
117 117 self.ylabel = 'Range [km]'
118 118 self.xlabel = 'Time'
119 119 self.cb_label = '(m/s)'
120 120 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
121 121 self.titles = ['{} Channel {}'.format(
122 122 self.CODE.upper(), x) for x in range(self.nrows)]
123 123
124 124 def update(self, dataOut):
125 125 #print(self.EEJtype)
126 126
127 127 if self.EEJtype == 1:
128 128 data = {
129 129 'dop': dataOut.Oblique_params[:,-2,:]
130 130 }
131 131 elif self.EEJtype == 2:
132 132 data = {
133 133 'dop': dataOut.Oblique_params[:,-1,:]
134 134 }
135 135
136 136 return data, {}
137 137
138 138 class DopplerEEJPlot(RTIPlot):
139 139 '''
140 140 Written by R. Flores
141 141 '''
142 142 '''
143 143 Plot for Doppler Shift EEJ
144 144 '''
145 145
146 146 CODE = 'dop'
147 147 colormap = 'RdBu_r'
148 148 #colormap = 'jet'
149 149
150 150 def setup(self):
151 151
152 152 self.xaxis = 'time'
153 153 self.ncols = 1
154 154 self.nrows = 2
155 155 self.nplots = 2
156 156 self.ylabel = 'Range [km]'
157 157 self.xlabel = 'Time'
158 158 self.cb_label = '(m/s)'
159 159 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
160 160 self.titles = ['{} EJJ Type {} /'.format(
161 161 self.CODE.upper(), x) for x in range(1,1+self.nrows)]
162 162
163 163 def update(self, dataOut):
164 164
165 165 if dataOut.mode == 11: #Double Gaussian
166 166 doppler = numpy.append(dataOut.Oblique_params[:,1,:],dataOut.Oblique_params[:,4,:],axis=0)
167 167 elif dataOut.mode == 9: #Double Skew Gaussian
168 168 doppler = numpy.append(dataOut.Oblique_params[:,-2,:],dataOut.Oblique_params[:,-1,:],axis=0)
169 169 data = {
170 170 'dop': doppler
171 171 }
172 172
173 173 return data, {}
174 174
175 175 class SpcWidthEEJPlot(RTIPlot):
176 176 '''
177 177 Written by R. Flores
178 178 '''
179 179 '''
180 180 Plot for EEJ Spectral Width
181 181 '''
182 182
183 183 CODE = 'width'
184 184 colormap = 'RdBu_r'
185 185 colormap = 'jet'
186 186
187 187 def setup(self):
188 188
189 189 self.xaxis = 'time'
190 190 self.ncols = 1
191 191 self.nrows = 2
192 192 self.nplots = 2
193 193 self.ylabel = 'Range [km]'
194 194 self.xlabel = 'Time'
195 195 self.cb_label = '(m/s)'
196 196 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
197 197 self.titles = ['{} EJJ Type {} /'.format(
198 198 self.CODE.upper(), x) for x in range(1,1+self.nrows)]
199 199
200 200 def update(self, dataOut):
201 201
202 202 if dataOut.mode == 11: #Double Gaussian
203 203 width = numpy.append(dataOut.Oblique_params[:,2,:],dataOut.Oblique_params[:,5,:],axis=0)
204 204 elif dataOut.mode == 9: #Double Skew Gaussian
205 205 width = numpy.append(dataOut.Oblique_params[:,2,:],dataOut.Oblique_params[:,6,:],axis=0)
206 206 data = {
207 207 'width': width
208 208 }
209 209
210 210 return data, {}
211 211
212 212 class PowerPlot(RTIPlot):
213 213 '''
214 214 Plot for Power Data (0 moment)
215 215 '''
216 216
217 217 CODE = 'pow'
218 218 colormap = 'jet'
219 219
220 220 def update(self, dataOut):
221 221
222 222 data = {
223 223 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
224 224 }
225 225
226 226 return data, {}
227 227
228 228 class SpectralWidthPlot(RTIPlot):
229 229 '''
230 230 Plot for Spectral Width Data (2nd moment)
231 231 '''
232 232
233 233 CODE = 'width'
234 234 colormap = 'jet'
235 235
236 236 def update(self, dataOut):
237 237
238 238 data = {
239 239 'width': dataOut.data_width
240 240 }
241 241
242 242 return data, {}
243 243
244 244 class SkyMapPlot(Plot):
245 245 '''
246 246 Plot for meteors detection data
247 247 '''
248 248
249 249 CODE = 'param'
250 250
251 251 def setup(self):
252 252
253 253 self.ncols = 1
254 254 self.nrows = 1
255 255 self.width = 7.2
256 256 self.height = 7.2
257 257 self.nplots = 1
258 258 self.xlabel = 'Zonal Zenith Angle (deg)'
259 259 self.ylabel = 'Meridional Zenith Angle (deg)'
260 260 self.polar = True
261 261 self.ymin = -180
262 262 self.ymax = 180
263 263 self.colorbar = False
264 264
265 265 def plot(self):
266 266
267 267 arrayParameters = numpy.concatenate(self.data['param'])
268 268 error = arrayParameters[:, -1]
269 269 indValid = numpy.where(error == 0)[0]
270 270 finalMeteor = arrayParameters[indValid, :]
271 271 finalAzimuth = finalMeteor[:, 3]
272 272 finalZenith = finalMeteor[:, 4]
273 273
274 274 x = finalAzimuth * numpy.pi / 180
275 275 y = finalZenith
276 276
277 277 ax = self.axes[0]
278 278
279 279 if ax.firsttime:
280 280 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
281 281 else:
282 282 ax.plot.set_data(x, y)
283 283
284 284 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
285 285 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
286 286 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
287 287 dt2,
288 288 len(x))
289 289 self.titles[0] = title
290 290
291 291
292 292 class GenericRTIPlot(Plot):
293 293 '''
294 294 Plot for data_xxxx object
295 295 '''
296 296
297 297 CODE = 'param'
298 298 colormap = 'viridis'
299 299 plot_type = 'pcolorbuffer'
300 300
301 301 def setup(self):
302 302 self.xaxis = 'time'
303 303 self.ncols = 1
304 304 self.nrows = self.data.shape('param')[0]
305 305 self.nplots = self.nrows
306 306 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
307 307
308 308 if not self.xlabel:
309 309 self.xlabel = 'Time'
310 310
311 311 self.ylabel = 'Range [km]'
312 312 if not self.titles:
313 313 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
314 314
315 315 def update(self, dataOut):
316 316
317 317 data = {
318 318 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
319 319 }
320 320
321 321 meta = {}
322 322
323 323 return data, meta
324 324
325 325 def plot(self):
326 326 # self.data.normalize_heights()
327 327 self.x = self.data.times
328 328 self.y = self.data.yrange
329 329 self.z = self.data['param']
330 330
331 331 self.z = numpy.ma.masked_invalid(self.z)
332 332
333 333 if self.decimation is None:
334 334 x, y, z = self.fill_gaps(self.x, self.y, self.z)
335 335 else:
336 336 x, y, z = self.fill_gaps(*self.decimate())
337 337
338 338 for n, ax in enumerate(self.axes):
339 339
340 340 self.zmax = self.zmax if self.zmax is not None else numpy.max(
341 341 self.z[n])
342 342 self.zmin = self.zmin if self.zmin is not None else numpy.min(
343 343 self.z[n])
344 344
345 345 if ax.firsttime:
346 346 if self.zlimits is not None:
347 347 self.zmin, self.zmax = self.zlimits[n]
348 348
349 349 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
350 350 vmin=self.zmin,
351 351 vmax=self.zmax,
352 352 cmap=self.cmaps[n]
353 353 )
354 354 else:
355 355 if self.zlimits is not None:
356 356 self.zmin, self.zmax = self.zlimits[n]
357 ax.plt.remove()
357
358 try:
359 ax.collections.remove(ax.collections[0])
360 except:
361 pass
362 # ax.plt.remove()
358 363 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
359 364 vmin=self.zmin,
360 365 vmax=self.zmax,
361 366 cmap=self.cmaps[n]
362 367 )
363 368
364 369
365 370 class PolarMapPlot(Plot):
366 371 '''
367 372 Plot for weather radar
368 373 '''
369 374
370 375 CODE = 'param'
371 376 colormap = 'seismic'
372 377
373 378 def setup(self):
374 379 self.ncols = 1
375 380 self.nrows = 1
376 381 self.width = 9
377 382 self.height = 8
378 383 self.mode = self.data.meta['mode']
379 384 if self.channels is not None:
380 385 self.nplots = len(self.channels)
381 386 self.nrows = len(self.channels)
382 387 else:
383 388 self.nplots = self.data.shape(self.CODE)[0]
384 389 self.nrows = self.nplots
385 390 self.channels = list(range(self.nplots))
386 391 if self.mode == 'E':
387 392 self.xlabel = 'Longitude'
388 393 self.ylabel = 'Latitude'
389 394 else:
390 395 self.xlabel = 'Range (km)'
391 396 self.ylabel = 'Height (km)'
392 397 self.bgcolor = 'white'
393 398 self.cb_labels = self.data.meta['units']
394 399 self.lat = self.data.meta['latitude']
395 400 self.lon = self.data.meta['longitude']
396 401 self.xmin, self.xmax = float(
397 402 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
398 403 self.ymin, self.ymax = float(
399 404 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
400 405 # self.polar = True
401 406
402 407 def plot(self):
403 408
404 409 for n, ax in enumerate(self.axes):
405 410 data = self.data['param'][self.channels[n]]
406 411
407 412 zeniths = numpy.linspace(
408 413 0, self.data.meta['max_range'], data.shape[1])
409 414 if self.mode == 'E':
410 415 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
411 416 r, theta = numpy.meshgrid(zeniths, azimuths)
412 417 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
413 418 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
414 419 x = km2deg(x) + self.lon
415 420 y = km2deg(y) + self.lat
416 421 else:
417 422 azimuths = numpy.radians(self.data.yrange)
418 423 r, theta = numpy.meshgrid(zeniths, azimuths)
419 424 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
420 425 self.y = zeniths
421 426
422 427 if ax.firsttime:
423 428 if self.zlimits is not None:
424 429 self.zmin, self.zmax = self.zlimits[n]
425 430 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
426 431 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
427 432 vmin=self.zmin,
428 433 vmax=self.zmax,
429 434 cmap=self.cmaps[n])
430 435 else:
431 436 if self.zlimits is not None:
432 437 self.zmin, self.zmax = self.zlimits[n]
433 438 ax.plt.remove()
434 439 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
435 440 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
436 441 vmin=self.zmin,
437 442 vmax=self.zmax,
438 443 cmap=self.cmaps[n])
439 444
440 445 if self.mode == 'A':
441 446 continue
442 447
443 448 # plot district names
444 449 f = open('/data/workspace/schain_scripts/distrito.csv')
445 450 for line in f:
446 451 label, lon, lat = [s.strip() for s in line.split(',') if s]
447 452 lat = float(lat)
448 453 lon = float(lon)
449 454 # ax.plot(lon, lat, '.b', ms=2)
450 455 ax.text(lon, lat, label.decode('utf8'), ha='center',
451 456 va='bottom', size='8', color='black')
452 457
453 458 # plot limites
454 459 limites = []
455 460 tmp = []
456 461 for line in open('/data/workspace/schain_scripts/lima.csv'):
457 462 if '#' in line:
458 463 if tmp:
459 464 limites.append(tmp)
460 465 tmp = []
461 466 continue
462 467 values = line.strip().split(',')
463 468 tmp.append((float(values[0]), float(values[1])))
464 469 for points in limites:
465 470 ax.add_patch(
466 471 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
467 472
468 473 # plot Cuencas
469 474 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
470 475 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
471 476 values = [line.strip().split(',') for line in f]
472 477 points = [(float(s[0]), float(s[1])) for s in values]
473 478 ax.add_patch(Polygon(points, ec='b', fc='none'))
474 479
475 480 # plot grid
476 481 for r in (15, 30, 45, 60):
477 482 ax.add_artist(plt.Circle((self.lon, self.lat),
478 483 km2deg(r), color='0.6', fill=False, lw=0.2))
479 484 ax.text(
480 485 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
481 486 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
482 487 '{}km'.format(r),
483 488 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
484 489
485 490 if self.mode == 'E':
486 491 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
487 492 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
488 493 else:
489 494 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
490 495 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
491 496
492 497 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
493 498 self.titles = ['{} {}'.format(
494 499 self.data.parameters[x], title) for x in self.channels]
General Comments 0
You need to be logged in to leave comments. Login now