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