##// END OF EJS Templates
schain
RSS
Se inclute SpectralMoments y DoubleGaussianPlot de tipo SpectraPlot
Danny Scipión -
r1358:fd32f6159d97
parent child
Show More
Side by Side Unified
Context file:
r1358:fd32f6159d97 Loading diff...:
Show file before | Show file after | Hide comments
1 NO CONTENT: new file 100644, binary diff hidden
Show file before | Show file after | Hide comments
@@ -1,358 +1,371
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot
6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
7 7 from schainpy.utils import log
8 8
9 9 EARTH_RADIUS = 6.3710e3
10 10
11 11
12 12 def ll2xy(lat1, lon1, lat2, lon2):
13 13
14 14 p = 0.017453292519943295
15 15 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
16 16 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
17 17 r = 12742 * numpy.arcsin(numpy.sqrt(a))
18 18 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
19 19 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
20 20 theta = -theta + numpy.pi/2
21 21 return r*numpy.cos(theta), r*numpy.sin(theta)
22 22
23 23
24 24 def km2deg(km):
25 25 '''
26 26 Convert distance in km to degrees
27 27 '''
28 28
29 29 return numpy.rad2deg(km/EARTH_RADIUS)
30 30
31 31
32 32
33 33 class SpectralMomentsPlot(SpectraPlot):
34 34 '''
35 35 Plot for Spectral Moments
36 36 '''
37 37 CODE = 'spc_moments'
38 colormap = 'jet'
39 plot_type = 'pcolor'
38 # colormap = 'jet'
39 # plot_type = 'pcolor'
40
41 class DobleGaussianPlot(SpectraPlot):
42 '''
43 Plot for Double Gaussian Plot
44 '''
45 CODE = 'gaussian_fit'
46 # colormap = 'jet'
47 # plot_type = 'pcolor'
40 48
49 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
50 '''
51 Plot SpectraCut with Double Gaussian Fit
52 '''
53 CODE = 'cut_gaussian_fit'
41 54
42 55 class SnrPlot(RTIPlot):
43 56 '''
44 57 Plot for SNR Data
45 58 '''
46 59
47 60 CODE = 'snr'
48 61 colormap = 'jet'
49 62
50 63 def update(self, dataOut):
51 64
52 65 data = {
53 66 'snr': 10*numpy.log10(dataOut.data_snr)
54 67 }
55 68
56 69 return data, {}
57 70
58 71 class DopplerPlot(RTIPlot):
59 72 '''
60 73 Plot for DOPPLER Data (1st moment)
61 74 '''
62 75
63 76 CODE = 'dop'
64 77 colormap = 'jet'
65 78
66 79 def update(self, dataOut):
67 80
68 81 data = {
69 82 'dop': 10*numpy.log10(dataOut.data_dop)
70 83 }
71 84
72 85 return data, {}
73 86
74 87 class PowerPlot(RTIPlot):
75 88 '''
76 89 Plot for Power Data (0 moment)
77 90 '''
78 91
79 92 CODE = 'pow'
80 93 colormap = 'jet'
81 94
82 95 def update(self, dataOut):
83 96
84 97 data = {
85 'pow': 10*numpy.log10(dataOut.data_pow)
98 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
86 99 }
87 100
88 101 return data, {}
89 102
90 103 class SpectralWidthPlot(RTIPlot):
91 104 '''
92 105 Plot for Spectral Width Data (2nd moment)
93 106 '''
94 107
95 108 CODE = 'width'
96 109 colormap = 'jet'
97 110
98 111 def update(self, dataOut):
99 112
100 113 data = {
101 114 'width': dataOut.data_width
102 115 }
103 116
104 117 return data, {}
105 118
106 119 class SkyMapPlot(Plot):
107 120 '''
108 121 Plot for meteors detection data
109 122 '''
110 123
111 124 CODE = 'param'
112 125
113 126 def setup(self):
114 127
115 128 self.ncols = 1
116 129 self.nrows = 1
117 130 self.width = 7.2
118 131 self.height = 7.2
119 132 self.nplots = 1
120 133 self.xlabel = 'Zonal Zenith Angle (deg)'
121 134 self.ylabel = 'Meridional Zenith Angle (deg)'
122 135 self.polar = True
123 136 self.ymin = -180
124 137 self.ymax = 180
125 138 self.colorbar = False
126 139
127 140 def plot(self):
128 141
129 142 arrayParameters = numpy.concatenate(self.data['param'])
130 143 error = arrayParameters[:, -1]
131 144 indValid = numpy.where(error == 0)[0]
132 145 finalMeteor = arrayParameters[indValid, :]
133 146 finalAzimuth = finalMeteor[:, 3]
134 147 finalZenith = finalMeteor[:, 4]
135 148
136 149 x = finalAzimuth * numpy.pi / 180
137 150 y = finalZenith
138 151
139 152 ax = self.axes[0]
140 153
141 154 if ax.firsttime:
142 155 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
143 156 else:
144 157 ax.plot.set_data(x, y)
145 158
146 159 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
147 160 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
148 161 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
149 162 dt2,
150 163 len(x))
151 164 self.titles[0] = title
152 165
153 166
154 167 class GenericRTIPlot(Plot):
155 168 '''
156 169 Plot for data_xxxx object
157 170 '''
158 171
159 172 CODE = 'param'
160 173 colormap = 'viridis'
161 174 plot_type = 'pcolorbuffer'
162 175
163 176 def setup(self):
164 177 self.xaxis = 'time'
165 178 self.ncols = 1
166 179 self.nrows = self.data.shape(self.attr_data)[0]
167 180 self.nplots = self.nrows
168 181 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
169 182
170 183 if not self.xlabel:
171 184 self.xlabel = 'Time'
172 185
173 self.ylabel = 'Height [km]'
186 self.ylabel = 'Range [km]'
174 187 if not self.titles:
175 188 self.titles = self.data.parameters \
176 189 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
177 190
178 191 def update(self, dataOut):
179 192
180 193 data = {
181 194 self.attr_data : getattr(dataOut, self.attr_data)
182 195 }
183 196
184 197 meta = {}
185 198
186 199 return data, meta
187 200
188 201 def plot(self):
189 202 # self.data.normalize_heights()
190 203 self.x = self.data.times
191 204 self.y = self.data.yrange
192 205 self.z = self.data[self.attr_data]
193 206
194 207 self.z = numpy.ma.masked_invalid(self.z)
195 208
196 209 if self.decimation is None:
197 210 x, y, z = self.fill_gaps(self.x, self.y, self.z)
198 211 else:
199 212 x, y, z = self.fill_gaps(*self.decimate())
200 213
201 214 for n, ax in enumerate(self.axes):
202 215
203 216 self.zmax = self.zmax if self.zmax is not None else numpy.max(
204 217 self.z[n])
205 218 self.zmin = self.zmin if self.zmin is not None else numpy.min(
206 219 self.z[n])
207 220
208 221 if ax.firsttime:
209 222 if self.zlimits is not None:
210 223 self.zmin, self.zmax = self.zlimits[n]
211 224
212 225 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
213 226 vmin=self.zmin,
214 227 vmax=self.zmax,
215 228 cmap=self.cmaps[n]
216 229 )
217 230 else:
218 231 if self.zlimits is not None:
219 232 self.zmin, self.zmax = self.zlimits[n]
220 233 ax.collections.remove(ax.collections[0])
221 234 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
222 235 vmin=self.zmin,
223 236 vmax=self.zmax,
224 237 cmap=self.cmaps[n]
225 238 )
226 239
227 240
228 241 class PolarMapPlot(Plot):
229 242 '''
230 243 Plot for weather radar
231 244 '''
232 245
233 246 CODE = 'param'
234 247 colormap = 'seismic'
235 248
236 249 def setup(self):
237 250 self.ncols = 1
238 251 self.nrows = 1
239 252 self.width = 9
240 253 self.height = 8
241 254 self.mode = self.data.meta['mode']
242 255 if self.channels is not None:
243 256 self.nplots = len(self.channels)
244 257 self.nrows = len(self.channels)
245 258 else:
246 259 self.nplots = self.data.shape(self.CODE)[0]
247 260 self.nrows = self.nplots
248 261 self.channels = list(range(self.nplots))
249 262 if self.mode == 'E':
250 263 self.xlabel = 'Longitude'
251 264 self.ylabel = 'Latitude'
252 265 else:
253 266 self.xlabel = 'Range (km)'
254 267 self.ylabel = 'Height (km)'
255 268 self.bgcolor = 'white'
256 269 self.cb_labels = self.data.meta['units']
257 270 self.lat = self.data.meta['latitude']
258 271 self.lon = self.data.meta['longitude']
259 272 self.xmin, self.xmax = float(
260 273 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
261 274 self.ymin, self.ymax = float(
262 275 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
263 276 # self.polar = True
264 277
265 278 def plot(self):
266 279
267 280 for n, ax in enumerate(self.axes):
268 281 data = self.data['param'][self.channels[n]]
269 282
270 283 zeniths = numpy.linspace(
271 284 0, self.data.meta['max_range'], data.shape[1])
272 285 if self.mode == 'E':
273 286 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
274 287 r, theta = numpy.meshgrid(zeniths, azimuths)
275 288 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
276 289 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
277 290 x = km2deg(x) + self.lon
278 291 y = km2deg(y) + self.lat
279 292 else:
280 293 azimuths = numpy.radians(self.data.yrange)
281 294 r, theta = numpy.meshgrid(zeniths, azimuths)
282 295 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
283 296 self.y = zeniths
284 297
285 298 if ax.firsttime:
286 299 if self.zlimits is not None:
287 300 self.zmin, self.zmax = self.zlimits[n]
288 301 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
289 302 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
290 303 vmin=self.zmin,
291 304 vmax=self.zmax,
292 305 cmap=self.cmaps[n])
293 306 else:
294 307 if self.zlimits is not None:
295 308 self.zmin, self.zmax = self.zlimits[n]
296 309 ax.collections.remove(ax.collections[0])
297 310 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
298 311 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
299 312 vmin=self.zmin,
300 313 vmax=self.zmax,
301 314 cmap=self.cmaps[n])
302 315
303 316 if self.mode == 'A':
304 317 continue
305 318
306 319 # plot district names
307 320 f = open('/data/workspace/schain_scripts/distrito.csv')
308 321 for line in f:
309 322 label, lon, lat = [s.strip() for s in line.split(',') if s]
310 323 lat = float(lat)
311 324 lon = float(lon)
312 325 # ax.plot(lon, lat, '.b', ms=2)
313 326 ax.text(lon, lat, label.decode('utf8'), ha='center',
314 327 va='bottom', size='8', color='black')
315 328
316 329 # plot limites
317 330 limites = []
318 331 tmp = []
319 332 for line in open('/data/workspace/schain_scripts/lima.csv'):
320 333 if '#' in line:
321 334 if tmp:
322 335 limites.append(tmp)
323 336 tmp = []
324 337 continue
325 338 values = line.strip().split(',')
326 339 tmp.append((float(values[0]), float(values[1])))
327 340 for points in limites:
328 341 ax.add_patch(
329 342 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
330 343
331 344 # plot Cuencas
332 345 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
333 346 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
334 347 values = [line.strip().split(',') for line in f]
335 348 points = [(float(s[0]), float(s[1])) for s in values]
336 349 ax.add_patch(Polygon(points, ec='b', fc='none'))
337 350
338 351 # plot grid
339 352 for r in (15, 30, 45, 60):
340 353 ax.add_artist(plt.Circle((self.lon, self.lat),
341 354 km2deg(r), color='0.6', fill=False, lw=0.2))
342 355 ax.text(
343 356 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
344 357 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
345 358 '{}km'.format(r),
346 359 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
347 360
348 361 if self.mode == 'E':
349 362 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
350 363 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
351 364 else:
352 365 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
353 366 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
354 367
355 368 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
356 369 self.titles = ['{} {}'.format(
357 370 self.data.parameters[x], title) for x in self.channels]
358 371
Show file before | Show file after | Hide comments
@@ -1,702 +1,743
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2021 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13
14 14
15 15 class SpectraPlot(Plot):
16 16 '''
17 17 Plot for Spectra data
18 18 '''
19 19
20 20 CODE = 'spc'
21 21 colormap = 'jet'
22 22 plot_type = 'pcolor'
23 23 buffering = False
24 24
25 25 def setup(self):
26 26 self.nplots = len(self.data.channels)
27 27 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
28 28 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
29 29 self.height = 2.6 * self.nrows
30 30 self.cb_label = 'dB'
31 31 if self.showprofile:
32 32 self.width = 4 * self.ncols
33 33 else:
34 34 self.width = 3.5 * self.ncols
35 35 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
36 36 self.ylabel = 'Range [km]'
37 37
38 38 def update(self, dataOut):
39 39
40 40 data = {}
41 41 meta = {}
42 42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 43 data['spc'] = spc
44 44 data['rti'] = dataOut.getPower()
45 45 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 46 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
47
47 48 if self.CODE == 'spc_moments':
48 49 data['moments'] = dataOut.moments
49
50 # data['spc'] = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
51 if self.CODE == 'gaussian_fit':
52 # data['moments'] = dataOut.moments
53 data['gaussfit'] = dataOut.DGauFitParams
54 # data['spc'] = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
55
50 56 return data, meta
51 57
52 58 def plot(self):
53 59 if self.xaxis == "frequency":
54 60 x = self.data.xrange[0]
55 61 self.xlabel = "Frequency (kHz)"
56 62 elif self.xaxis == "time":
57 63 x = self.data.xrange[1]
58 64 self.xlabel = "Time (ms)"
59 65 else:
60 66 x = self.data.xrange[2]
61 67 self.xlabel = "Velocity (m/s)"
62 68
63 if self.CODE == 'spc_moments':
69 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
64 70 x = self.data.xrange[2]
65 71 self.xlabel = "Velocity (m/s)"
66 72
67 73 self.titles = []
68 74
69 75 y = self.data.yrange
70 76 self.y = y
71 77
72 78 data = self.data[-1]
73 79 z = data['spc']
74 80
75 81 for n, ax in enumerate(self.axes):
76 82 noise = data['noise'][n]
77 83 if self.CODE == 'spc_moments':
78 mean = data['moments'][n, 2]
84 mean = data['moments'][n, 1]
85 if self.CODE == 'gaussian_fit':
86 # mean = data['moments'][n, 1]
87 gau0 = data['gaussfit'][n][2,:,0]
88 gau1 = data['gaussfit'][n][2,:,1]
79 89 if ax.firsttime:
80 90 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
81 91 self.xmin = self.xmin if self.xmin else -self.xmax
82 92 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
83 93 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
84 94 ax.plt = ax.pcolormesh(x, y, z[n].T,
85 95 vmin=self.zmin,
86 96 vmax=self.zmax,
87 97 cmap=plt.get_cmap(self.colormap)
88 98 )
89 99
90 100 if self.showprofile:
91 101 ax.plt_profile = self.pf_axes[n].plot(
92 102 data['rti'][n], y)[0]
93 103 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
94 104 color="k", linestyle="dashed", lw=1)[0]
95 105 if self.CODE == 'spc_moments':
96 ax.plt_mean = ax.plot(mean, y, color='k')[0]
106 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
107 if self.CODE == 'gaussian_fit':
108 # ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
109 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
110 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
97 111 else:
98 112 ax.plt.set_array(z[n].T.ravel())
99 113 if self.showprofile:
100 114 ax.plt_profile.set_data(data['rti'][n], y)
101 115 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
102 116 if self.CODE == 'spc_moments':
103 117 ax.plt_mean.set_data(mean, y)
118 if self.CODE == 'gaussian_fit':
119 # ax.plt_mean.set_data(mean, y)
120 ax.plt_gau0.set_data(gau0, y)
121 ax.plt_gau1.set_data(gau1, y)
104 122 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
105 123
106 124
107 125 class CrossSpectraPlot(Plot):
108 126
109 127 CODE = 'cspc'
110 128 colormap = 'jet'
111 129 plot_type = 'pcolor'
112 130 zmin_coh = None
113 131 zmax_coh = None
114 132 zmin_phase = None
115 133 zmax_phase = None
116 134
117 135 def setup(self):
118 136
119 137 self.ncols = 4
120 138 self.nplots = len(self.data.pairs) * 2
121 139 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
122 140 self.width = 3.1 * self.ncols
123 141 self.height = 2.6 * self.nrows
124 142 self.ylabel = 'Range [km]'
125 143 self.showprofile = False
126 144 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
127 145
128 146 def update(self, dataOut):
129 147
130 148 data = {}
131 149 meta = {}
132 150
133 151 spc = dataOut.data_spc
134 152 cspc = dataOut.data_cspc
135 153 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
136 154 meta['pairs'] = dataOut.pairsList
137 155
138 156 tmp = []
139 157
140 158 for n, pair in enumerate(meta['pairs']):
141 159 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
142 160 coh = numpy.abs(out)
143 161 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
144 162 tmp.append(coh)
145 163 tmp.append(phase)
146 164
147 165 data['cspc'] = numpy.array(tmp)
148 166
149 167 return data, meta
150 168
151 169 def plot(self):
152 170
153 171 if self.xaxis == "frequency":
154 172 x = self.data.xrange[0]
155 173 self.xlabel = "Frequency (kHz)"
156 174 elif self.xaxis == "time":
157 175 x = self.data.xrange[1]
158 176 self.xlabel = "Time (ms)"
159 177 else:
160 178 x = self.data.xrange[2]
161 179 self.xlabel = "Velocity (m/s)"
162 180
163 181 self.titles = []
164 182
165 183 y = self.data.yrange
166 184 self.y = y
167 185
168 186 data = self.data[-1]
169 187 cspc = data['cspc']
170 188
171 189 for n in range(len(self.data.pairs)):
172 190 pair = self.data.pairs[n]
173 191 coh = cspc[n*2]
174 192 phase = cspc[n*2+1]
175 193 ax = self.axes[2 * n]
176 194 if ax.firsttime:
177 195 ax.plt = ax.pcolormesh(x, y, coh.T,
178 196 vmin=0,
179 197 vmax=1,
180 198 cmap=plt.get_cmap(self.colormap_coh)
181 199 )
182 200 else:
183 201 ax.plt.set_array(coh.T.ravel())
184 202 self.titles.append(
185 203 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
186 204
187 205 ax = self.axes[2 * n + 1]
188 206 if ax.firsttime:
189 207 ax.plt = ax.pcolormesh(x, y, phase.T,
190 208 vmin=-180,
191 209 vmax=180,
192 210 cmap=plt.get_cmap(self.colormap_phase)
193 211 )
194 212 else:
195 213 ax.plt.set_array(phase.T.ravel())
196 214 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
197 215
198 216
199 217 class RTIPlot(Plot):
200 218 '''
201 219 Plot for RTI data
202 220 '''
203 221
204 222 CODE = 'rti'
205 223 colormap = 'jet'
206 224 plot_type = 'pcolorbuffer'
207 225
208 226 def setup(self):
209 227 self.xaxis = 'time'
210 228 self.ncols = 1
211 229 self.nrows = len(self.data.channels)
212 230 self.nplots = len(self.data.channels)
213 231 self.ylabel = 'Range [km]'
214 232 self.xlabel = 'Time'
215 233 self.cb_label = 'dB'
216 234 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
217 235 self.titles = ['{} Channel {}'.format(
218 236 self.CODE.upper(), x) for x in range(self.nrows)]
219 237
220 238 def update(self, dataOut):
221 239
222 240 data = {}
223 241 meta = {}
224 242 data['rti'] = dataOut.getPower()
225 243 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
226 244
227 245 return data, meta
228 246
229 247 def plot(self):
230 248 self.x = self.data.times
231 249 self.y = self.data.yrange
232 250 self.z = self.data[self.CODE]
233 251 self.z = numpy.ma.masked_invalid(self.z)
234 252
235 253 if self.decimation is None:
236 254 x, y, z = self.fill_gaps(self.x, self.y, self.z)
237 255 else:
238 256 x, y, z = self.fill_gaps(*self.decimate())
239 257
240 258 for n, ax in enumerate(self.axes):
241 259 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
242 260 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
243 261 data = self.data[-1]
244 262 if ax.firsttime:
245 263 ax.plt = ax.pcolormesh(x, y, z[n].T,
246 264 vmin=self.zmin,
247 265 vmax=self.zmax,
248 266 cmap=plt.get_cmap(self.colormap)
249 267 )
250 268 if self.showprofile:
251 269 ax.plot_profile = self.pf_axes[n].plot(
252 270 data['rti'][n], self.y)[0]
253 271 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
254 272 color="k", linestyle="dashed", lw=1)[0]
255 273 else:
256 274 ax.collections.remove(ax.collections[0])
257 275 ax.plt = ax.pcolormesh(x, y, z[n].T,
258 276 vmin=self.zmin,
259 277 vmax=self.zmax,
260 278 cmap=plt.get_cmap(self.colormap)
261 279 )
262 280 if self.showprofile:
263 281 ax.plot_profile.set_data(data['rti'][n], self.y)
264 282 ax.plot_noise.set_data(numpy.repeat(
265 283 data['noise'][n], len(self.y)), self.y)
266 284
267 285
268 286 class CoherencePlot(RTIPlot):
269 287 '''
270 288 Plot for Coherence data
271 289 '''
272 290
273 291 CODE = 'coh'
274 292
275 293 def setup(self):
276 294 self.xaxis = 'time'
277 295 self.ncols = 1
278 296 self.nrows = len(self.data.pairs)
279 297 self.nplots = len(self.data.pairs)
280 298 self.ylabel = 'Range [km]'
281 299 self.xlabel = 'Time'
282 300 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
283 301 if self.CODE == 'coh':
284 302 self.cb_label = ''
285 303 self.titles = [
286 304 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
287 305 else:
288 306 self.cb_label = 'Degrees'
289 307 self.titles = [
290 308 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
291 309
292 310 def update(self, dataOut):
293 311
294 312 data = {}
295 313 meta = {}
296 314 data['coh'] = dataOut.getCoherence()
297 315 meta['pairs'] = dataOut.pairsList
298 316
299 317 return data, meta
300 318
301 319 class PhasePlot(CoherencePlot):
302 320 '''
303 321 Plot for Phase map data
304 322 '''
305 323
306 324 CODE = 'phase'
307 325 colormap = 'seismic'
308 326
309 327 def update(self, dataOut):
310 328
311 329 data = {}
312 330 meta = {}
313 331 data['phase'] = dataOut.getCoherence(phase=True)
314 332 meta['pairs'] = dataOut.pairsList
315 333
316 334 return data, meta
317 335
318 336 class NoisePlot(Plot):
319 337 '''
320 338 Plot for noise
321 339 '''
322 340
323 341 CODE = 'noise'
324 342 plot_type = 'scatterbuffer'
325 343
326 344 def setup(self):
327 345 self.xaxis = 'time'
328 346 self.ncols = 1
329 347 self.nrows = 1
330 348 self.nplots = 1
331 349 self.ylabel = 'Intensity [dB]'
332 350 self.xlabel = 'Time'
333 351 self.titles = ['Noise']
334 352 self.colorbar = False
335 353 self.plots_adjust.update({'right': 0.85 })
336 354
337 355 def update(self, dataOut):
338 356
339 357 data = {}
340 358 meta = {}
341 359 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
342 360 meta['yrange'] = numpy.array([])
343 361
344 362 return data, meta
345 363
346 364 def plot(self):
347 365
348 366 x = self.data.times
349 367 xmin = self.data.min_time
350 368 xmax = xmin + self.xrange * 60 * 60
351 369 Y = self.data['noise']
352 370
353 371 if self.axes[0].firsttime:
354 372 self.ymin = numpy.nanmin(Y) - 5
355 373 self.ymax = numpy.nanmax(Y) + 5
356 374 for ch in self.data.channels:
357 375 y = Y[ch]
358 376 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
359 377 plt.legend(bbox_to_anchor=(1.18, 1.0))
360 378 else:
361 379 for ch in self.data.channels:
362 380 y = Y[ch]
363 381 self.axes[0].lines[ch].set_data(x, y)
364 382
365 383
366 384 class PowerProfilePlot(Plot):
367 385
368 386 CODE = 'pow_profile'
369 387 plot_type = 'scatter'
370 388
371 389 def setup(self):
372 390
373 391 self.ncols = 1
374 392 self.nrows = 1
375 393 self.nplots = 1
376 394 self.height = 4
377 395 self.width = 3
378 396 self.ylabel = 'Range [km]'
379 397 self.xlabel = 'Intensity [dB]'
380 398 self.titles = ['Power Profile']
381 399 self.colorbar = False
382 400
383 401 def update(self, dataOut):
384 402
385 403 data = {}
386 404 meta = {}
387 405 data[self.CODE] = dataOut.getPower()
388 406
389 407 return data, meta
390 408
391 409 def plot(self):
392 410
393 411 y = self.data.yrange
394 412 self.y = y
395 413
396 414 x = self.data[-1][self.CODE]
397 415
398 416 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
399 417 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
400 418
401 419 if self.axes[0].firsttime:
402 420 for ch in self.data.channels:
403 421 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
404 422 plt.legend()
405 423 else:
406 424 for ch in self.data.channels:
407 425 self.axes[0].lines[ch].set_data(x[ch], y)
408 426
409 427
410 428 class SpectraCutPlot(Plot):
411 429
412 430 CODE = 'spc_cut'
413 431 plot_type = 'scatter'
414 432 buffering = False
415 433
416 434 def setup(self):
417 435
418 436 self.nplots = len(self.data.channels)
419 437 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
420 438 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
421 439 self.width = 3.4 * self.ncols + 1.5
422 440 self.height = 3 * self.nrows
423 441 self.ylabel = 'Power [dB]'
424 442 self.colorbar = False
425 443 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
426 444
427 445 def update(self, dataOut):
428 446
429 447 data = {}
430 448 meta = {}
431 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
449 spc = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
432 450 data['spc'] = spc
433 451 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
434
452 if self.CODE == 'cut_gaussian_fit':
453 data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
454 data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
435 455 return data, meta
436 456
437 457 def plot(self):
438 458 if self.xaxis == "frequency":
439 459 x = self.data.xrange[0][1:]
440 460 self.xlabel = "Frequency (kHz)"
441 461 elif self.xaxis == "time":
442 462 x = self.data.xrange[1]
443 463 self.xlabel = "Time (ms)"
444 464 else:
445 x = self.data.xrange[2]
465 x = self.data.xrange[2][:-1]
466 self.xlabel = "Velocity (m/s)"
467
468 if self.CODE == 'cut_gaussian_fit':
469 x = self.data.xrange[2][:-1]
446 470 self.xlabel = "Velocity (m/s)"
447 471
448 472 self.titles = []
449 473
450 474 y = self.data.yrange
451 z = self.data[-1]['spc']
475 data = self.data[-1]
476 z = data['spc']
452 477
453 478 if self.height_index:
454 479 index = numpy.array(self.height_index)
455 480 else:
456 481 index = numpy.arange(0, len(y), int((len(y))/9))
457 482
458 483 for n, ax in enumerate(self.axes):
484 if self.CODE == 'cut_gaussian_fit':
485 gau0 = data['gauss_fit0']
486 gau1 = data['gauss_fit1']
459 487 if ax.firsttime:
460 488 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
461 489 self.xmin = self.xmin if self.xmin else -self.xmax
462 490 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
463 491 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
464 ax.plt = ax.plot(x, z[n, :, index].T)
492 ax.plt = ax.plot(x, z[n, :, index].T, lw=0.25)
493 if self.CODE == 'cut_gaussian_fit':
494 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
495 for i, line in enumerate(ax.plt_gau0):
496 line.set_color(ax.plt[i].get_color())
497 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
498 for i, line in enumerate(ax.plt_gau1):
499 line.set_color(ax.plt[i].get_color())
465 500 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
466 501 self.figures[0].legend(ax.plt, labels, loc='center right')
467 502 else:
468 503 for i, line in enumerate(ax.plt):
469 line.set_data(x, z[n, :, i])
504 line.set_data(x, z[n, :, index[i]].T)
505 for i, line in enumerate(ax.plt_gau0):
506 line.set_data(x, gau0[n, :, index[i]].T)
507 line.set_color(ax.plt[i].get_color())
508 for i, line in enumerate(ax.plt_gau1):
509 line.set_data(x, gau1[n, :, index[i]].T)
510 line.set_color(ax.plt[i].get_color())
470 511 self.titles.append('CH {}'.format(n))
471 512
472 513
473 514 class BeaconPhase(Plot):
474 515
475 516 __isConfig = None
476 517 __nsubplots = None
477 518
478 519 PREFIX = 'beacon_phase'
479 520
480 521 def __init__(self):
481 522 Plot.__init__(self)
482 523 self.timerange = 24*60*60
483 524 self.isConfig = False
484 525 self.__nsubplots = 1
485 526 self.counter_imagwr = 0
486 527 self.WIDTH = 800
487 528 self.HEIGHT = 400
488 529 self.WIDTHPROF = 120
489 530 self.HEIGHTPROF = 0
490 531 self.xdata = None
491 532 self.ydata = None
492 533
493 534 self.PLOT_CODE = BEACON_CODE
494 535
495 536 self.FTP_WEI = None
496 537 self.EXP_CODE = None
497 538 self.SUB_EXP_CODE = None
498 539 self.PLOT_POS = None
499 540
500 541 self.filename_phase = None
501 542
502 543 self.figfile = None
503 544
504 545 self.xmin = None
505 546 self.xmax = None
506 547
507 548 def getSubplots(self):
508 549
509 550 ncol = 1
510 551 nrow = 1
511 552
512 553 return nrow, ncol
513 554
514 555 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
515 556
516 557 self.__showprofile = showprofile
517 558 self.nplots = nplots
518 559
519 560 ncolspan = 7
520 561 colspan = 6
521 562 self.__nsubplots = 2
522 563
523 564 self.createFigure(id = id,
524 565 wintitle = wintitle,
525 566 widthplot = self.WIDTH+self.WIDTHPROF,
526 567 heightplot = self.HEIGHT+self.HEIGHTPROF,
527 568 show=show)
528 569
529 570 nrow, ncol = self.getSubplots()
530 571
531 572 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
532 573
533 574 def save_phase(self, filename_phase):
534 575 f = open(filename_phase,'w+')
535 576 f.write('\n\n')
536 577 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
537 578 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
538 579 f.close()
539 580
540 581 def save_data(self, filename_phase, data, data_datetime):
541 582 f=open(filename_phase,'a')
542 583 timetuple_data = data_datetime.timetuple()
543 584 day = str(timetuple_data.tm_mday)
544 585 month = str(timetuple_data.tm_mon)
545 586 year = str(timetuple_data.tm_year)
546 587 hour = str(timetuple_data.tm_hour)
547 588 minute = str(timetuple_data.tm_min)
548 589 second = str(timetuple_data.tm_sec)
549 590 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
550 591 f.close()
551 592
552 593 def plot(self):
553 594 log.warning('TODO: Not yet implemented...')
554 595
555 596 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
556 597 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
557 598 timerange=None,
558 599 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
559 600 server=None, folder=None, username=None, password=None,
560 601 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
561 602
562 603 if dataOut.flagNoData:
563 604 return dataOut
564 605
565 606 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
566 607 return
567 608
568 609 if pairsList == None:
569 610 pairsIndexList = dataOut.pairsIndexList[:10]
570 611 else:
571 612 pairsIndexList = []
572 613 for pair in pairsList:
573 614 if pair not in dataOut.pairsList:
574 615 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
575 616 pairsIndexList.append(dataOut.pairsList.index(pair))
576 617
577 618 if pairsIndexList == []:
578 619 return
579 620
580 621 # if len(pairsIndexList) > 4:
581 622 # pairsIndexList = pairsIndexList[0:4]
582 623
583 624 hmin_index = None
584 625 hmax_index = None
585 626
586 627 if hmin != None and hmax != None:
587 628 indexes = numpy.arange(dataOut.nHeights)
588 629 hmin_list = indexes[dataOut.heightList >= hmin]
589 630 hmax_list = indexes[dataOut.heightList <= hmax]
590 631
591 632 if hmin_list.any():
592 633 hmin_index = hmin_list[0]
593 634
594 635 if hmax_list.any():
595 636 hmax_index = hmax_list[-1]+1
596 637
597 638 x = dataOut.getTimeRange()
598 639
599 640 thisDatetime = dataOut.datatime
600 641
601 642 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
602 643 xlabel = "Local Time"
603 644 ylabel = "Phase (degrees)"
604 645
605 646 update_figfile = False
606 647
607 648 nplots = len(pairsIndexList)
608 649 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
609 650 phase_beacon = numpy.zeros(len(pairsIndexList))
610 651 for i in range(nplots):
611 652 pair = dataOut.pairsList[pairsIndexList[i]]
612 653 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
613 654 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
614 655 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
615 656 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
616 657 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
617 658
618 659 if dataOut.beacon_heiIndexList:
619 660 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
620 661 else:
621 662 phase_beacon[i] = numpy.average(phase)
622 663
623 664 if not self.isConfig:
624 665
625 666 nplots = len(pairsIndexList)
626 667
627 668 self.setup(id=id,
628 669 nplots=nplots,
629 670 wintitle=wintitle,
630 671 showprofile=showprofile,
631 672 show=show)
632 673
633 674 if timerange != None:
634 675 self.timerange = timerange
635 676
636 677 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
637 678
638 679 if ymin == None: ymin = 0
639 680 if ymax == None: ymax = 360
640 681
641 682 self.FTP_WEI = ftp_wei
642 683 self.EXP_CODE = exp_code
643 684 self.SUB_EXP_CODE = sub_exp_code
644 685 self.PLOT_POS = plot_pos
645 686
646 687 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
647 688 self.isConfig = True
648 689 self.figfile = figfile
649 690 self.xdata = numpy.array([])
650 691 self.ydata = numpy.array([])
651 692
652 693 update_figfile = True
653 694
654 695 #open file beacon phase
655 696 path = '%s%03d' %(self.PREFIX, self.id)
656 697 beacon_file = os.path.join(path,'%s.txt'%self.name)
657 698 self.filename_phase = os.path.join(figpath,beacon_file)
658 699 #self.save_phase(self.filename_phase)
659 700
660 701
661 702 #store data beacon phase
662 703 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
663 704
664 705 self.setWinTitle(title)
665 706
666 707
667 708 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
668 709
669 710 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
670 711
671 712 axes = self.axesList[0]
672 713
673 714 self.xdata = numpy.hstack((self.xdata, x[0:1]))
674 715
675 716 if len(self.ydata)==0:
676 717 self.ydata = phase_beacon.reshape(-1,1)
677 718 else:
678 719 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
679 720
680 721
681 722 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
682 723 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
683 724 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
684 725 XAxisAsTime=True, grid='both'
685 726 )
686 727
687 728 self.draw()
688 729
689 730 if dataOut.ltctime >= self.xmax:
690 731 self.counter_imagwr = wr_period
691 732 self.isConfig = False
692 733 update_figfile = True
693 734
694 735 self.save(figpath=figpath,
695 736 figfile=figfile,
696 737 save=save,
697 738 ftp=ftp,
698 739 wr_period=wr_period,
699 740 thisDatetime=thisDatetime,
700 741 update_figfile=update_figfile)
701 742
702 743 return dataOut No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now