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