@@ -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 |
|
|
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_ |
|
726 | reader_d = shpreader.BasicReader(shape_d, encoding='latin1') | |
725 |
reader_p = shpreader.BasicReader(shape_ |
|
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