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