##// END OF EJS Templates
schain
RSS
Use function elaz to latlon instead of wradlib
Juan C. Espinoza -
r1530:69c85e5f5aa8
parent child
Show More
Side by Side Unified
Context file:
r1530:69c85e5f5aa8 Loading diff...:

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

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