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