##// END OF EJS Templates
schain
RSS
Fix weatherplot for two channels
jespinoza -
r1544:24f1bea5f47d
parent child
Show More
Side by Side Unified
Context file:
r1544:24f1bea5f47d Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,691 +1,697
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
7 7 from schainpy.model.graphics.jroplot_base import Plot, plt
8 8 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
9 9 from schainpy.utils import log
10 10
11 11
12 12 EARTH_RADIUS = 6.3710e3
13 13
14 14
15 15 def antenna_to_cartesian(ranges, azimuths, elevations):
16 16 """
17 17 Return Cartesian coordinates from antenna coordinates.
18 18
19 19 Parameters
20 20 ----------
21 21 ranges : array
22 22 Distances to the center of the radar gates (bins) in kilometers.
23 23 azimuths : array
24 24 Azimuth angle of the radar in degrees.
25 25 elevations : array
26 26 Elevation angle of the radar in degrees.
27 27
28 28 Returns
29 29 -------
30 30 x, y, z : array
31 31 Cartesian coordinates in meters from the radar.
32 32
33 33 Notes
34 34 -----
35 35 The calculation for Cartesian coordinate is adapted from equations
36 36 2.28(b) and 2.28(c) of Doviak and Zrnic [1]_ assuming a
37 37 standard atmosphere (4/3 Earth's radius model).
38 38
39 39 .. math::
40 40
41 41 z = \\sqrt{r^2+R^2+2*r*R*sin(\\theta_e)} - R
42 42
43 43 s = R * arcsin(\\frac{r*cos(\\theta_e)}{R+z})
44 44
45 45 x = s * sin(\\theta_a)
46 46
47 47 y = s * cos(\\theta_a)
48 48
49 49 Where r is the distance from the radar to the center of the gate,
50 50 :math:`\\theta_a` is the azimuth angle, :math:`\\theta_e` is the
51 51 elevation angle, s is the arc length, and R is the effective radius
52 52 of the earth, taken to be 4/3 the mean radius of earth (6371 km).
53 53
54 54 References
55 55 ----------
56 56 .. [1] Doviak and Zrnic, Doppler Radar and Weather Observations, Second
57 57 Edition, 1993, p. 21.
58 58
59 59 """
60 60 theta_e = numpy.deg2rad(elevations) # elevation angle in radians.
61 61 theta_a = numpy.deg2rad(azimuths) # azimuth angle in radians.
62 62 R = 6371.0 * 1000.0 * 4.0 / 3.0 # effective radius of earth in meters.
63 63 r = ranges * 1000.0 # distances to gates in meters.
64 64
65 65 z = (r ** 2 + R ** 2 + 2.0 * r * R * numpy.sin(theta_e)) ** 0.5 - R
66 66 s = R * numpy.arcsin(r * numpy.cos(theta_e) / (R + z)) # arc length in m.
67 67 x = s * numpy.sin(theta_a)
68 68 y = s * numpy.cos(theta_a)
69 69 return x, y, z
70 70
71 71 def cartesian_to_geographic_aeqd(x, y, lon_0, lat_0, R=EARTH_RADIUS):
72 72 """
73 73 Azimuthal equidistant Cartesian to geographic coordinate transform.
74 74
75 75 Transform a set of Cartesian/Cartographic coordinates (x, y) to
76 76 geographic coordinate system (lat, lon) using a azimuthal equidistant
77 77 map projection [1]_.
78 78
79 79 .. math::
80 80
81 81 lat = \\arcsin(\\cos(c) * \\sin(lat_0) +
82 82 (y * \\sin(c) * \\cos(lat_0) / \\rho))
83 83
84 84 lon = lon_0 + \\arctan2(
85 85 x * \\sin(c),
86 86 \\rho * \\cos(lat_0) * \\cos(c) - y * \\sin(lat_0) * \\sin(c))
87 87
88 88 \\rho = \\sqrt(x^2 + y^2)
89 89
90 90 c = \\rho / R
91 91
92 92 Where x, y are the Cartesian position from the center of projection;
93 93 lat, lon the corresponding latitude and longitude; lat_0, lon_0 are the
94 94 latitude and longitude of the center of the projection; R is the radius of
95 95 the earth (defaults to ~6371 km). lon is adjusted to be between -180 and
96 96 180.
97 97
98 98 Parameters
99 99 ----------
100 100 x, y : array-like
101 101 Cartesian coordinates in the same units as R, typically meters.
102 102 lon_0, lat_0 : float
103 103 Longitude and latitude, in degrees, of the center of the projection.
104 104 R : float, optional
105 105 Earth radius in the same units as x and y. The default value is in
106 106 units of meters.
107 107
108 108 Returns
109 109 -------
110 110 lon, lat : array
111 111 Longitude and latitude of Cartesian coordinates in degrees.
112 112
113 113 References
114 114 ----------
115 115 .. [1] Snyder, J. P. Map Projections--A Working Manual. U. S. Geological
116 116 Survey Professional Paper 1395, 1987, pp. 191-202.
117 117
118 118 """
119 119 x = numpy.atleast_1d(numpy.asarray(x))
120 120 y = numpy.atleast_1d(numpy.asarray(y))
121 121
122 122 lat_0_rad = numpy.deg2rad(lat_0)
123 123 lon_0_rad = numpy.deg2rad(lon_0)
124 124
125 125 rho = numpy.sqrt(x*x + y*y)
126 126 c = rho / R
127 127
128 128 with warnings.catch_warnings():
129 129 # division by zero may occur here but is properly addressed below so
130 130 # the warnings can be ignored
131 131 warnings.simplefilter("ignore", RuntimeWarning)
132 132 lat_rad = numpy.arcsin(numpy.cos(c) * numpy.sin(lat_0_rad) +
133 133 y * numpy.sin(c) * numpy.cos(lat_0_rad) / rho)
134 134 lat_deg = numpy.rad2deg(lat_rad)
135 135 # fix cases where the distance from the center of the projection is zero
136 136 lat_deg[rho == 0] = lat_0
137 137
138 138 x1 = x * numpy.sin(c)
139 139 x2 = rho*numpy.cos(lat_0_rad)*numpy.cos(c) - y*numpy.sin(lat_0_rad)*numpy.sin(c)
140 140 lon_rad = lon_0_rad + numpy.arctan2(x1, x2)
141 141 lon_deg = numpy.rad2deg(lon_rad)
142 142 # Longitudes should be from -180 to 180 degrees
143 143 lon_deg[lon_deg > 180] -= 360.
144 144 lon_deg[lon_deg < -180] += 360.
145 145
146 146 return lon_deg, lat_deg
147 147
148 148 def antenna_to_geographic(ranges, azimuths, elevations, site):
149 149
150 150 x, y, z = antenna_to_cartesian(numpy.array(ranges), numpy.array(azimuths), numpy.array(elevations))
151 151 lon, lat = cartesian_to_geographic_aeqd(x, y, site[0], site[1], R=6370997.)
152 152
153 153 return lon, lat
154 154
155 155 def ll2xy(lat1, lon1, lat2, lon2):
156 156
157 157 p = 0.017453292519943295
158 158 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
159 159 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
160 160 r = 12742 * numpy.arcsin(numpy.sqrt(a))
161 161 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
162 162 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
163 163 theta = -theta + numpy.pi/2
164 164 return r*numpy.cos(theta), r*numpy.sin(theta)
165 165
166 166
167 167 def km2deg(km):
168 168 '''
169 169 Convert distance in km to degrees
170 170 '''
171 171
172 172 return numpy.rad2deg(km/EARTH_RADIUS)
173 173
174 174
175 175
176 176 class SpectralMomentsPlot(SpectraPlot):
177 177 '''
178 178 Plot for Spectral Moments
179 179 '''
180 180 CODE = 'spc_moments'
181 181 # colormap = 'jet'
182 182 # plot_type = 'pcolor'
183 183
184 184 class DobleGaussianPlot(SpectraPlot):
185 185 '''
186 186 Plot for Double Gaussian Plot
187 187 '''
188 188 CODE = 'gaussian_fit'
189 189 # colormap = 'jet'
190 190 # plot_type = 'pcolor'
191 191
192 192 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
193 193 '''
194 194 Plot SpectraCut with Double Gaussian Fit
195 195 '''
196 196 CODE = 'cut_gaussian_fit'
197 197
198 198 class SnrPlot(RTIPlot):
199 199 '''
200 200 Plot for SNR Data
201 201 '''
202 202
203 203 CODE = 'snr'
204 204 colormap = 'jet'
205 205
206 206 def update(self, dataOut):
207 207
208 208 data = {
209 209 'snr': 10*numpy.log10(dataOut.data_snr)
210 210 }
211 211
212 212 return data, {}
213 213
214 214 class DopplerPlot(RTIPlot):
215 215 '''
216 216 Plot for DOPPLER Data (1st moment)
217 217 '''
218 218
219 219 CODE = 'dop'
220 220 colormap = 'jet'
221 221
222 222 def update(self, dataOut):
223 223
224 224 data = {
225 225 'dop': 10*numpy.log10(dataOut.data_dop)
226 226 }
227 227
228 228 return data, {}
229 229
230 230 class PowerPlot(RTIPlot):
231 231 '''
232 232 Plot for Power Data (0 moment)
233 233 '''
234 234
235 235 CODE = 'pow'
236 236 colormap = 'jet'
237 237
238 238 def update(self, dataOut):
239 239 data = {
240 240 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
241 241 }
242 242 return data, {}
243 243
244 244 class SpectralWidthPlot(RTIPlot):
245 245 '''
246 246 Plot for Spectral Width Data (2nd moment)
247 247 '''
248 248
249 249 CODE = 'width'
250 250 colormap = 'jet'
251 251
252 252 def update(self, dataOut):
253 253
254 254 data = {
255 255 'width': dataOut.data_width
256 256 }
257 257
258 258 return data, {}
259 259
260 260 class SkyMapPlot(Plot):
261 261 '''
262 262 Plot for meteors detection data
263 263 '''
264 264
265 265 CODE = 'param'
266 266
267 267 def setup(self):
268 268
269 269 self.ncols = 1
270 270 self.nrows = 1
271 271 self.width = 7.2
272 272 self.height = 7.2
273 273 self.nplots = 1
274 274 self.xlabel = 'Zonal Zenith Angle (deg)'
275 275 self.ylabel = 'Meridional Zenith Angle (deg)'
276 276 self.polar = True
277 277 self.ymin = -180
278 278 self.ymax = 180
279 279 self.colorbar = False
280 280
281 281 def plot(self):
282 282
283 283 arrayParameters = numpy.concatenate(self.data['param'])
284 284 error = arrayParameters[:, -1]
285 285 indValid = numpy.where(error == 0)[0]
286 286 finalMeteor = arrayParameters[indValid, :]
287 287 finalAzimuth = finalMeteor[:, 3]
288 288 finalZenith = finalMeteor[:, 4]
289 289
290 290 x = finalAzimuth * numpy.pi / 180
291 291 y = finalZenith
292 292
293 293 ax = self.axes[0]
294 294
295 295 if ax.firsttime:
296 296 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
297 297 else:
298 298 ax.plot.set_data(x, y)
299 299
300 300 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
301 301 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
302 302 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
303 303 dt2,
304 304 len(x))
305 305 self.titles[0] = title
306 306
307 307
308 308 class GenericRTIPlot(Plot):
309 309 '''
310 310 Plot for data_xxxx object
311 311 '''
312 312
313 313 CODE = 'param'
314 314 colormap = 'viridis'
315 315 plot_type = 'pcolorbuffer'
316 316
317 317 def setup(self):
318 318 self.xaxis = 'time'
319 319 self.ncols = 1
320 320 self.nrows = self.data.shape('param')[0]
321 321 self.nplots = self.nrows
322 322 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
323 323
324 324 if not self.xlabel:
325 325 self.xlabel = 'Time'
326 326
327 327 self.ylabel = 'Range [km]'
328 328 if not self.titles:
329 329 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
330 330
331 331 def update(self, dataOut):
332 332
333 333 data = {
334 334 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
335 335 }
336 336
337 337 meta = {}
338 338
339 339 return data, meta
340 340
341 341 def plot(self):
342 342 # self.data.normalize_heights()
343 343 self.x = self.data.times
344 344 self.y = self.data.yrange
345 345 self.z = self.data['param']
346 346 self.z = 10*numpy.log10(self.z)
347 347 self.z = numpy.ma.masked_invalid(self.z)
348 348
349 349 if self.decimation is None:
350 350 x, y, z = self.fill_gaps(self.x, self.y, self.z)
351 351 else:
352 352 x, y, z = self.fill_gaps(*self.decimate())
353 353
354 354 for n, ax in enumerate(self.axes):
355 355
356 356 self.zmax = self.zmax if self.zmax is not None else numpy.max(
357 357 self.z[n])
358 358 self.zmin = self.zmin if self.zmin is not None else numpy.min(
359 359 self.z[n])
360 360
361 361 if ax.firsttime:
362 362 if self.zlimits is not None:
363 363 self.zmin, self.zmax = self.zlimits[n]
364 364
365 365 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
366 366 vmin=self.zmin,
367 367 vmax=self.zmax,
368 368 cmap=self.cmaps[n]
369 369 )
370 370 else:
371 371 if self.zlimits is not None:
372 372 self.zmin, self.zmax = self.zlimits[n]
373 373 ax.collections.remove(ax.collections[0])
374 374 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
375 375 vmin=self.zmin,
376 376 vmax=self.zmax,
377 377 cmap=self.cmaps[n]
378 378 )
379 379
380 380
381 381 class PolarMapPlot(Plot):
382 382 '''
383 383 Plot for weather radar
384 384 '''
385 385
386 386 CODE = 'param'
387 387 colormap = 'seismic'
388 388
389 389 def setup(self):
390 390 self.ncols = 1
391 391 self.nrows = 1
392 392 self.width = 9
393 393 self.height = 8
394 394 self.mode = self.data.meta['mode']
395 395 if self.channels is not None:
396 396 self.nplots = len(self.channels)
397 397 self.nrows = len(self.channels)
398 398 else:
399 399 self.nplots = self.data.shape(self.CODE)[0]
400 400 self.nrows = self.nplots
401 401 self.channels = list(range(self.nplots))
402 402 if self.mode == 'E':
403 403 self.xlabel = 'Longitude'
404 404 self.ylabel = 'Latitude'
405 405 else:
406 406 self.xlabel = 'Range (km)'
407 407 self.ylabel = 'Height (km)'
408 408 self.bgcolor = 'white'
409 409 self.cb_labels = self.data.meta['units']
410 410 self.lat = self.data.meta['latitude']
411 411 self.lon = self.data.meta['longitude']
412 412 self.xmin, self.xmax = float(
413 413 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
414 414 self.ymin, self.ymax = float(
415 415 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
416 416 # self.polar = True
417 417
418 418 def plot(self):
419 419
420 420 for n, ax in enumerate(self.axes):
421 421 data = self.data['param'][self.channels[n]]
422 422
423 423 zeniths = numpy.linspace(
424 424 0, self.data.meta['max_range'], data.shape[1])
425 425 if self.mode == 'E':
426 426 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
427 427 r, theta = numpy.meshgrid(zeniths, azimuths)
428 428 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
429 429 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
430 430 x = km2deg(x) + self.lon
431 431 y = km2deg(y) + self.lat
432 432 else:
433 433 azimuths = numpy.radians(self.data.yrange)
434 434 r, theta = numpy.meshgrid(zeniths, azimuths)
435 435 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
436 436 self.y = zeniths
437 437
438 438 if ax.firsttime:
439 439 if self.zlimits is not None:
440 440 self.zmin, self.zmax = self.zlimits[n]
441 441 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
442 442 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
443 443 vmin=self.zmin,
444 444 vmax=self.zmax,
445 445 cmap=self.cmaps[n])
446 446 else:
447 447 if self.zlimits is not None:
448 448 self.zmin, self.zmax = self.zlimits[n]
449 449 ax.collections.remove(ax.collections[0])
450 450 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
451 451 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
452 452 vmin=self.zmin,
453 453 vmax=self.zmax,
454 454 cmap=self.cmaps[n])
455 455
456 456 if self.mode == 'A':
457 457 continue
458 458
459 459 # plot district names
460 460 f = open('/data/workspace/schain_scripts/distrito.csv')
461 461 for line in f:
462 462 label, lon, lat = [s.strip() for s in line.split(',') if s]
463 463 lat = float(lat)
464 464 lon = float(lon)
465 465 # ax.plot(lon, lat, '.b', ms=2)
466 466 ax.text(lon, lat, label.decode('utf8'), ha='center',
467 467 va='bottom', size='8', color='black')
468 468
469 469 # plot limites
470 470 limites = []
471 471 tmp = []
472 472 for line in open('/data/workspace/schain_scripts/lima.csv'):
473 473 if '#' in line:
474 474 if tmp:
475 475 limites.append(tmp)
476 476 tmp = []
477 477 continue
478 478 values = line.strip().split(',')
479 479 tmp.append((float(values[0]), float(values[1])))
480 480 for points in limites:
481 481 ax.add_patch(
482 482 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
483 483
484 484 # plot Cuencas
485 485 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
486 486 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
487 487 values = [line.strip().split(',') for line in f]
488 488 points = [(float(s[0]), float(s[1])) for s in values]
489 489 ax.add_patch(Polygon(points, ec='b', fc='none'))
490 490
491 491 # plot grid
492 492 for r in (15, 30, 45, 60):
493 493 ax.add_artist(plt.Circle((self.lon, self.lat),
494 494 km2deg(r), color='0.6', fill=False, lw=0.2))
495 495 ax.text(
496 496 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
497 497 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
498 498 '{}km'.format(r),
499 499 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
500 500
501 501 if self.mode == 'E':
502 502 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
503 503 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
504 504 else:
505 505 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
506 506 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
507 507
508 508 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
509 509 self.titles = ['{} {}'.format(
510 510 self.data.parameters[x], title) for x in self.channels]
511 511
512 512 class WeatherParamsPlot(Plot):
513 513 #CODE = 'RHI'
514 514 #plot_name = 'RHI'
515 515 plot_type = 'scattermap'
516 516 buffering = False
517 517
518 518 def setup(self):
519 519
520 520 self.ncols = 1
521 521 self.nrows = 1
522 522 self.nplots= 1
523 523 self.ylabel= 'Range [km]'
524 524 self.xlabel= 'Range [km]'
525 525 self.polar = True
526 526 self.grid = True
527 527 if self.channels is not None:
528 528 self.nplots = len(self.channels)
529 self.nrows = len(self.channels)
529 self.ncols = len(self.channels)
530 530 else:
531 531 self.nplots = self.data.shape(self.CODE)[0]
532 self.nrows = self.nplots
532 self.ncols = self.nplots
533 533 self.channels = list(range(self.nplots))
534 534
535 535 self.colorbar=True
536 if len(self.channels)>1:
537 self.width = 12
538 else:
536 539 self.width =8
537 540 self.height =8
538 541 self.ini =0
539 542 self.len_azi =0
540 543 self.buffer_ini = None
541 544 self.buffer_ele = None
542 545 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
543 546 self.flag =0
544 547 self.indicador= 0
545 548 self.last_data_ele = None
546 549 self.val_mean = None
547 550
548 551 def update(self, dataOut):
549 552
550 553 vars = {
551 554 'S' : 0,
552 555 'V' : 1,
553 556 'W' : 2,
554 557 'SNR' : 3,
555 558 'Z' : 4,
556 559 'D' : 5,
557 560 'P' : 6,
558 561 'R' : 7,
559 562 }
560 563
561 564 data = {}
562 565 meta = {}
563 566
564 567 if hasattr(dataOut, 'nFFTPoints'):
565 568 factor = dataOut.normFactor
566 569 else:
567 570 factor = 1
568 571
569 572 if hasattr(dataOut, 'dparam'):
570 573 tmp = getattr(dataOut, 'data_param')
571 574 else:
572 575
573 576 if 'S' in self.attr_data[0]:
574 577 tmp = 10*numpy.log10(10.0*getattr(dataOut, 'data_param')[:,0,:]/(factor))
575 578 else:
576 579 tmp = getattr(dataOut, 'data_param')[:,vars[self.attr_data[0]],:]
577 580
578 581 if self.mask:
579 582 mask = dataOut.data_param[:,3,:] < self.mask
580 583 tmp = numpy.ma.masked_array(tmp, mask=mask)
581 584
582 585 r = dataOut.heightList
583 586 delta_height = r[1]-r[0]
584 587 valid = numpy.where(r>=0)[0]
585 588 data['r'] = numpy.arange(len(valid))*delta_height
586 589
587 try:
588 data['data'] = tmp[self.channels[0]][:,valid]
589 except:
590 data['data'] = tmp[0][:,valid]
590 data['data'] = [0, 0]
591
592 #try:
593 data['data'][0] = tmp[0][:,valid]
594 data['data'][1] = tmp[1][:,valid]
595 #except:
596 # data['data'] = tmp[0][:,valid]
591 597
592 598 if dataOut.mode_op == 'PPI':
593 599 self.CODE = 'PPI'
594 600 self.title = self.CODE
595 601 elif dataOut.mode_op == 'RHI':
596 602 self.CODE = 'RHI'
597 603 self.title = self.CODE
598 604
599 605 data['azi'] = dataOut.data_azi
600 606 data['ele'] = dataOut.data_ele
601 607 data['mode_op'] = dataOut.mode_op
602 608 self.mode = dataOut.mode_op
603 var = data['data'].flatten()
604 r = numpy.tile(data['r'], data['data'].shape[0])
605 az = numpy.repeat(data['azi'], data['data'].shape[1])
606 el = numpy.repeat(data['ele'], data['data'].shape[1])
609 var = data['data'][0].flatten()
610 r = numpy.tile(data['r'], data['data'][0].shape[0])
611 az = numpy.repeat(data['azi'], data['data'][0].shape[1])
612 el = numpy.repeat(data['ele'], data['data'][0].shape[1])
607 613
608 614 # lla = georef.spherical_to_proj(r, data['azi'], data['ele'], (-75.295893, -12.040436, 3379.2147))
609 615
610 616 latlon = antenna_to_geographic(r, az, el, (-75.295893, -12.040436))
611 617
612 618 if self.mask:
613 619 meta['lat'] = latlon[1][var.mask==False]
614 620 meta['lon'] = latlon[0][var.mask==False]
615 621 data['var'] = numpy.array([var[var.mask==False]])
616 622 else:
617 623 meta['lat'] = latlon[1]
618 624 meta['lon'] = latlon[0]
619 625 data['var'] = numpy.array([var])
620 626
621 627 return data, meta
622 628
623 629 def plot(self):
624 630 data = self.data[-1]
625 631 z = data['data']
626 632 r = data['r']
627 633 self.titles = []
628 634
629 635 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
630 636 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
631 637 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
632 638 self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z)
633 639
634 640 if isinstance(data['mode_op'], bytes):
635 641 data['mode_op'] = data['mode_op'].decode()
636 642
637 643 if data['mode_op'] == 'RHI':
638 644 try:
639 645 if self.data['mode_op'][-2] == 'PPI':
640 646 self.ang_min = None
641 647 self.ang_max = None
642 648 except:
643 649 pass
644 650 self.ang_min = self.ang_min if self.ang_min else 0
645 651 self.ang_max = self.ang_max if self.ang_max else 90
646 652 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )
647 653 elif data['mode_op'] == 'PPI':
648 654 try:
649 655 if self.data['mode_op'][-2] == 'RHI':
650 656 self.ang_min = None
651 657 self.ang_max = None
652 658 except:
653 659 pass
654 660 self.ang_min = self.ang_min if self.ang_min else 0
655 661 self.ang_max = self.ang_max if self.ang_max else 360
656 662 r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )
657 663
658 664 self.clear_figures()
659 665
660 666 for i,ax in enumerate(self.axes):
661 667
662 668 if ax.firsttime:
663 669 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
664 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
670 ax.plt = ax.pcolormesh(theta, r, z[i], cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
665 671 if data['mode_op'] == 'PPI':
666 672 ax.set_theta_direction(-1)
667 673 ax.set_theta_offset(numpy.pi/2)
668 674
669 675 else:
670 676 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
671 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
677 ax.plt = ax.pcolormesh(theta, r, z[i], cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
672 678 if data['mode_op'] == 'PPI':
673 679 ax.set_theta_direction(-1)
674 680 ax.set_theta_offset(numpy.pi/2)
675 681
676 682 ax.grid(True)
677 683 if data['mode_op'] == 'RHI':
678 684 len_aux = int(data['azi'].shape[0]/4)
679 685 mean = numpy.mean(data['azi'][len_aux:-len_aux])
680 686 if len(self.channels) !=1:
681 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
687 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels]
682 688 else:
683 689 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
684 690 elif data['mode_op'] == 'PPI':
685 691 len_aux = int(data['ele'].shape[0]/4)
686 692 mean = numpy.mean(data['ele'][len_aux:-len_aux])
687 693 if len(self.channels) !=1:
688 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
694 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels]
689 695 else:
690 696 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
691 697 self.mode_value = round(mean,1) No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now