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