##// END OF EJS Templates
schain
RSS
nuevo
joabAM -
r1472:c5c71a942112
parent child
Show More
Side by Side Unified
Context file:
r1472:c5c71a942112 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,82 +1,130
1 1 #include <Python.h>
2 2 #include <numpy/arrayobject.h>
3 3 #include <math.h>
4 4
5 5
6 6 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
7 7 double navg;
8 8 PyObject *data_obj, *data_array;
9 9
10 10 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
11 11 return NULL;
12 12 }
13 13
14 14 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
15 15
16 16 if (data_array == NULL) {
17 17 Py_XDECREF(data_array);
18 18 Py_XDECREF(data_obj);
19 19 return NULL;
20 20 }
21 21 double *sortdata = (double*)PyArray_DATA(data_array);
22 22 int lenOfData = (int)PyArray_SIZE(data_array) ;
23 23 double nums_min = lenOfData*0.2;
24 24 if (nums_min <= 5) nums_min = 5;
25 25 double sump = 0;
26 26 double sumq = 0;
27 27 long j = 0;
28 28 int cont = 1;
29 29 double rtest = 0;
30 30 while ((cont == 1) && (j < lenOfData)) {
31 31 sump = sump + sortdata[j];
32 32 sumq = sumq + pow(sortdata[j], 2);
33 33 if (j > nums_min) {
34 34 rtest = (double)j/(j-1) + 1/navg;
35 35 if ((sumq*j) > (rtest*pow(sump, 2))) {
36 36 j = j - 1;
37 37 sump = sump - sortdata[j];
38 38 sumq = sumq - pow(sortdata[j],2);
39 39 cont = 0;
40 40 }
41 41 }
42 42 j = j + 1;
43 43 }
44 44
45 45 double lnoise = sump / j;
46 46
47 47 Py_DECREF(data_array);
48 48
49 49 // return PyLong_FromLong(lnoise);
50 50 return PyFloat_FromDouble(lnoise);
51 51 }
52 /*
53 static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) {
54 double navg;
55 double th;
56 PyObject *data_obj, *data_array;
57
58 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg, &th)) {
59 return NULL;
60 }
61
62 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
63
64 if (data_array == NULL) {
65 Py_XDECREF(data_array);
66 Py_XDECREF(data_obj);
67 return NULL;
68 }
69 double *sortdata = (double*)PyArray_DATA(data_array);
70 int lenOfData = (int)PyArray_SIZE(data_array) ;
71 double nums_min = lenOfData*th;
72 if (nums_min <= 5) nums_min = 5;
73 double sump = 0;
74 double sumq = 0;
75 long j = 0;
76 int cont = 1;
77 double rtest = 0;
78 while ((cont == 1) && (j < lenOfData)) {
79 sump = sump + sortdata[j];
80 sumq = sumq + pow(sortdata[j], 2);
81 if (j > nums_min) {
82 rtest = (double)j/(j-1) + 1/navg;
83 if ((sumq*j) > (rtest*pow(sump, 2))) {
84 j = j - 1;
85 sump = sump - sortdata[j];
86 sumq = sumq - pow(sortdata[j],2);
87 cont = 0;
88 }
89 }
90 j = j + 1;
91 }
52 92
93 //double lnoise = sump / j;
94
95 Py_DECREF(data_array);
96
97 // return PyLong_FromLong(lnoise);
98 return PyFloat_FromDouble(j,sortID);
99 }
100 */
53 101
54 102 static PyMethodDef noiseMethods[] = {
55 103 { "hildebrand_sekhon", hildebrand_sekhon, METH_VARARGS, "Get noise with hildebrand_sekhon algorithm" },
56 104 { NULL, NULL, 0, NULL }
57 105 };
58 106
59 107 #if PY_MAJOR_VERSION >= 3
60 108
61 109 static struct PyModuleDef noisemodule = {
62 110 PyModuleDef_HEAD_INIT,
63 111 "_noise",
64 112 "Get noise with hildebrand_sekhon algorithm",
65 113 -1,
66 114 noiseMethods
67 115 };
68 116
69 117 #endif
70 118
71 119 #if PY_MAJOR_VERSION >= 3
72 120 PyMODINIT_FUNC PyInit__noise(void) {
73 121 Py_Initialize();
74 122 import_array();
75 123 return PyModule_Create(&noisemodule);
76 124 }
77 125 #else
78 126 PyMODINIT_FUNC init_noise() {
79 127 Py_InitModule("_noise", noiseMethods);
80 128 import_array();
81 129 }
82 130 #endif
Show file before | Show file after | Hide comments
@@ -1,357 +1,359
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 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot
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 38 colormap = 'jet'
39 39 plot_type = 'pcolor'
40 40
41 41
42 42 class SnrPlot(RTIPlot):
43 43 '''
44 44 Plot for SNR Data
45 45 '''
46 46
47 47 CODE = 'snr'
48 48 colormap = 'jet'
49 49
50 50 def update(self, dataOut):
51 51 if len(self.channelList) == 0:
52 52 self.update_list(dataOut)
53 data = {}
53
54 54 meta = {}
55 data['snr'] = 10*numpy.log10(dataOut.data_snr)
56
55 data = {
56 'snr': 10 * numpy.log10(dataOut.data_snr)
57 }
58 #print(data['snr'])
57 59 return data, meta
58 60
59 61 class DopplerPlot(RTIPlot):
60 62 '''
61 63 Plot for DOPPLER Data (1st moment)
62 64 '''
63 65
64 66 CODE = 'dop'
65 67 colormap = 'jet'
66 68
67 69 def update(self, dataOut):
68 70 self.update_list(dataOut)
69 71 data = {
70 72 'dop': 10*numpy.log10(dataOut.data_dop)
71 73 }
72 74
73 75 return data, {}
74 76
75 77 class PowerPlot(RTIPlot):
76 78 '''
77 79 Plot for Power Data (0 moment)
78 80 '''
79 81
80 82 CODE = 'pow'
81 83 colormap = 'jet'
82 84
83 85 def update(self, dataOut):
84 86 self.update_list(dataOut)
85 87 data = {
86 88 'pow': 10*numpy.log10(dataOut.data_pow)
87 89 }
88 90 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
89 91 return data, {}
90 92
91 93 class SpectralWidthPlot(RTIPlot):
92 94 '''
93 95 Plot for Spectral Width Data (2nd moment)
94 96 '''
95 97
96 98 CODE = 'width'
97 99 colormap = 'jet'
98 100
99 101 def update(self, dataOut):
100 102 self.update_list(dataOut)
101 103 data = {
102 104 'width': dataOut.data_width
103 105 }
104 106 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
105 107 return data, {}
106 108
107 109 class SkyMapPlot(Plot):
108 110 '''
109 111 Plot for meteors detection data
110 112 '''
111 113
112 114 CODE = 'param'
113 115
114 116 def setup(self):
115 117
116 118 self.ncols = 1
117 119 self.nrows = 1
118 120 self.width = 7.2
119 121 self.height = 7.2
120 122 self.nplots = 1
121 123 self.xlabel = 'Zonal Zenith Angle (deg)'
122 124 self.ylabel = 'Meridional Zenith Angle (deg)'
123 125 self.polar = True
124 126 self.ymin = -180
125 127 self.ymax = 180
126 128 self.colorbar = False
127 129
128 130 def plot(self):
129 131
130 132 arrayParameters = numpy.concatenate(self.data['param'])
131 133 error = arrayParameters[:, -1]
132 134 indValid = numpy.where(error == 0)[0]
133 135 finalMeteor = arrayParameters[indValid, :]
134 136 finalAzimuth = finalMeteor[:, 3]
135 137 finalZenith = finalMeteor[:, 4]
136 138
137 139 x = finalAzimuth * numpy.pi / 180
138 140 y = finalZenith
139 141
140 142 ax = self.axes[0]
141 143
142 144 if ax.firsttime:
143 145 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
144 146 else:
145 147 ax.plot.set_data(x, y)
146 148
147 149 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
148 150 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
149 151 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
150 152 dt2,
151 153 len(x))
152 154 self.titles[0] = title
153 155
154 156
155 157 class GenericRTIPlot(Plot):
156 158 '''
157 159 Plot for data_xxxx object
158 160 '''
159 161
160 162 CODE = 'param'
161 163 colormap = 'viridis'
162 164 plot_type = 'pcolorbuffer'
163 165
164 166 def setup(self):
165 167 self.xaxis = 'time'
166 168 self.ncols = 1
167 169 self.nrows = self.data.shape('param')[0]
168 170 self.nplots = self.nrows
169 171 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
170 172
171 173 if not self.xlabel:
172 174 self.xlabel = 'Time'
173 175
174 176 self.ylabel = 'Height [km]'
175 177 if not self.titles:
176 178 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
177 179
178 180 def update(self, dataOut):
179 181
180 182 data = {
181 183 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
182 184 }
183 185
184 186 meta = {}
185 187
186 188 return data, meta
187 189
188 190 def plot(self):
189 191 # self.data.normalize_heights()
190 192 self.x = self.data.times
191 193 self.y = self.data.yrange
192 194 self.z = self.data['param']
193 195
194 196 self.z = numpy.ma.masked_invalid(self.z)
195 197
196 198 if self.decimation is None:
197 199 x, y, z = self.fill_gaps(self.x, self.y, self.z)
198 200 else:
199 201 x, y, z = self.fill_gaps(*self.decimate())
200 202
201 203 for n, ax in enumerate(self.axes):
202 204
203 205 self.zmax = self.zmax if self.zmax is not None else numpy.max(
204 206 self.z[n])
205 207 self.zmin = self.zmin if self.zmin is not None else numpy.min(
206 208 self.z[n])
207 209
208 210 if ax.firsttime:
209 211 if self.zlimits is not None:
210 212 self.zmin, self.zmax = self.zlimits[n]
211 213
212 214 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
213 215 vmin=self.zmin,
214 216 vmax=self.zmax,
215 217 cmap=self.cmaps[n]
216 218 )
217 219 else:
218 220 if self.zlimits is not None:
219 221 self.zmin, self.zmax = self.zlimits[n]
220 222 ax.collections.remove(ax.collections[0])
221 223 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
222 224 vmin=self.zmin,
223 225 vmax=self.zmax,
224 226 cmap=self.cmaps[n]
225 227 )
226 228
227 229
228 230 class PolarMapPlot(Plot):
229 231 '''
230 232 Plot for weather radar
231 233 '''
232 234
233 235 CODE = 'param'
234 236 colormap = 'seismic'
235 237
236 238 def setup(self):
237 239 self.ncols = 1
238 240 self.nrows = 1
239 241 self.width = 9
240 242 self.height = 8
241 243 self.mode = self.data.meta['mode']
242 244 if self.channels is not None:
243 245 self.nplots = len(self.channels)
244 246 self.nrows = len(self.channels)
245 247 else:
246 248 self.nplots = self.data.shape(self.CODE)[0]
247 249 self.nrows = self.nplots
248 250 self.channels = list(range(self.nplots))
249 251 if self.mode == 'E':
250 252 self.xlabel = 'Longitude'
251 253 self.ylabel = 'Latitude'
252 254 else:
253 255 self.xlabel = 'Range (km)'
254 256 self.ylabel = 'Height (km)'
255 257 self.bgcolor = 'white'
256 258 self.cb_labels = self.data.meta['units']
257 259 self.lat = self.data.meta['latitude']
258 260 self.lon = self.data.meta['longitude']
259 261 self.xmin, self.xmax = float(
260 262 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
261 263 self.ymin, self.ymax = float(
262 264 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
263 265 # self.polar = True
264 266
265 267 def plot(self):
266 268
267 269 for n, ax in enumerate(self.axes):
268 270 data = self.data['param'][self.channels[n]]
269 271
270 272 zeniths = numpy.linspace(
271 273 0, self.data.meta['max_range'], data.shape[1])
272 274 if self.mode == 'E':
273 275 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
274 276 r, theta = numpy.meshgrid(zeniths, azimuths)
275 277 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
276 278 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
277 279 x = km2deg(x) + self.lon
278 280 y = km2deg(y) + self.lat
279 281 else:
280 282 azimuths = numpy.radians(self.data.yrange)
281 283 r, theta = numpy.meshgrid(zeniths, azimuths)
282 284 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
283 285 self.y = zeniths
284 286
285 287 if ax.firsttime:
286 288 if self.zlimits is not None:
287 289 self.zmin, self.zmax = self.zlimits[n]
288 290 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
289 291 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
290 292 vmin=self.zmin,
291 293 vmax=self.zmax,
292 294 cmap=self.cmaps[n])
293 295 else:
294 296 if self.zlimits is not None:
295 297 self.zmin, self.zmax = self.zlimits[n]
296 298 ax.collections.remove(ax.collections[0])
297 299 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
298 300 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
299 301 vmin=self.zmin,
300 302 vmax=self.zmax,
301 303 cmap=self.cmaps[n])
302 304
303 305 if self.mode == 'A':
304 306 continue
305 307
306 308 # plot district names
307 309 f = open('/data/workspace/schain_scripts/distrito.csv')
308 310 for line in f:
309 311 label, lon, lat = [s.strip() for s in line.split(',') if s]
310 312 lat = float(lat)
311 313 lon = float(lon)
312 314 # ax.plot(lon, lat, '.b', ms=2)
313 315 ax.text(lon, lat, label.decode('utf8'), ha='center',
314 316 va='bottom', size='8', color='black')
315 317
316 318 # plot limites
317 319 limites = []
318 320 tmp = []
319 321 for line in open('/data/workspace/schain_scripts/lima.csv'):
320 322 if '#' in line:
321 323 if tmp:
322 324 limites.append(tmp)
323 325 tmp = []
324 326 continue
325 327 values = line.strip().split(',')
326 328 tmp.append((float(values[0]), float(values[1])))
327 329 for points in limites:
328 330 ax.add_patch(
329 331 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
330 332
331 333 # plot Cuencas
332 334 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
333 335 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
334 336 values = [line.strip().split(',') for line in f]
335 337 points = [(float(s[0]), float(s[1])) for s in values]
336 338 ax.add_patch(Polygon(points, ec='b', fc='none'))
337 339
338 340 # plot grid
339 341 for r in (15, 30, 45, 60):
340 342 ax.add_artist(plt.Circle((self.lon, self.lat),
341 343 km2deg(r), color='0.6', fill=False, lw=0.2))
342 344 ax.text(
343 345 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
344 346 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
345 347 '{}km'.format(r),
346 348 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
347 349
348 350 if self.mode == 'E':
349 351 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
350 352 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
351 353 else:
352 354 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
353 355 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
354 356
355 357 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
356 358 self.titles = ['{} {}'.format(
357 359 self.data.parameters[x], title) for x in self.channels]
Show file before | Show file after | Hide comments
@@ -1,962 +1,966
1 1 # Copyright (c) 2012-2020 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 from itertools import combinations
14 14
15 15
16 16 class SpectraPlot(Plot):
17 17 '''
18 18 Plot for Spectra data
19 19 '''
20 20
21 21 CODE = 'spc'
22 22 colormap = 'jet'
23 23 plot_type = 'pcolor'
24 24 buffering = False
25 25 channelList = []
26 26
27 27 def setup(self):
28 28
29 29 self.nplots = len(self.data.channels)
30 30 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
31 31 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
32 32 self.height = 2.6 * self.nrows
33 33
34 34 self.cb_label = 'dB'
35 35 if self.showprofile:
36 36 self.width = 4 * self.ncols
37 37 else:
38 38 self.width = 3.5 * self.ncols
39 39 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
40 40 self.ylabel = 'Range [km]'
41 41
42 42
43 43 def update_list(self,dataOut):
44 44 if len(self.channelList) == 0:
45 45 self.channelList = dataOut.channelList
46 46
47 47 def update(self, dataOut):
48 48
49 49 self.update_list(dataOut)
50 50 data = {}
51 51 meta = {}
52 52 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
53 53 data['spc'] = spc
54 54 data['rti'] = dataOut.getPower()
55 55 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
56 56 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
57 57 if self.CODE == 'spc_moments':
58 58 data['moments'] = dataOut.moments
59 59
60 60 return data, meta
61 61
62 62 def plot(self):
63 63 if self.xaxis == "frequency":
64 64 x = self.data.xrange[0]
65 65 self.xlabel = "Frequency (kHz)"
66 66 elif self.xaxis == "time":
67 67 x = self.data.xrange[1]
68 68 self.xlabel = "Time (ms)"
69 69 else:
70 70 x = self.data.xrange[2]
71 71 self.xlabel = "Velocity (m/s)"
72 72
73 73 if self.CODE == 'spc_moments':
74 74 x = self.data.xrange[2]
75 75 self.xlabel = "Velocity (m/s)"
76 76
77 77 self.titles = []
78 78 y = self.data.yrange
79 79 self.y = y
80 80
81 81 data = self.data[-1]
82 82 z = data['spc']
83 83
84 84 for n, ax in enumerate(self.axes):
85 85 noise = data['noise'][n]
86 86 if self.CODE == 'spc_moments':
87 87 mean = data['moments'][n, 1]
88 88 if ax.firsttime:
89 89 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
90 90 self.xmin = self.xmin if self.xmin else -self.xmax
91 91 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
92 92 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
93 93 ax.plt = ax.pcolormesh(x, y, z[n].T,
94 94 vmin=self.zmin,
95 95 vmax=self.zmax,
96 96 cmap=plt.get_cmap(self.colormap)
97 97 )
98 98
99 99 if self.showprofile:
100 100 ax.plt_profile = self.pf_axes[n].plot(
101 101 data['rti'][n], y)[0]
102 102 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
103 103 color="k", linestyle="dashed", lw=1)[0]
104 104 if self.CODE == 'spc_moments':
105 105 ax.plt_mean = ax.plot(mean, y, color='k')[0]
106 106 else:
107 107 ax.plt.set_array(z[n].T.ravel())
108 108 if self.showprofile:
109 109 ax.plt_profile.set_data(data['rti'][n], y)
110 110 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
111 111 if self.CODE == 'spc_moments':
112 112 ax.plt_mean.set_data(mean, y)
113 113 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
114 114
115 115
116 116 class CrossSpectraPlot(Plot):
117 117
118 118 CODE = 'cspc'
119 119 colormap = 'jet'
120 120 plot_type = 'pcolor'
121 121 zmin_coh = None
122 122 zmax_coh = None
123 123 zmin_phase = None
124 124 zmax_phase = None
125 125 realChannels = None
126 126 crossPairs = None
127 127
128 128 def setup(self):
129 129
130 130 self.ncols = 4
131 131 self.nplots = len(self.data.pairs) * 2
132 132 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
133 133 self.width = 3.1 * self.ncols
134 134 self.height = 2.6 * self.nrows
135 135 self.ylabel = 'Range [km]'
136 136 self.showprofile = False
137 137 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
138 138
139 139 def update(self, dataOut):
140 140
141 141 data = {}
142 142 meta = {}
143 143
144 144 spc = dataOut.data_spc
145 145 cspc = dataOut.data_cspc
146 146 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
147 147 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
148 148 meta['pairs'] = rawPairs
149 149
150 150 if self.crossPairs == None:
151 151 self.crossPairs = dataOut.pairsList
152 152
153 153 tmp = []
154 154
155 155 for n, pair in enumerate(meta['pairs']):
156 156
157 157 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
158 158 coh = numpy.abs(out)
159 159 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
160 160 tmp.append(coh)
161 161 tmp.append(phase)
162 162
163 163 data['cspc'] = numpy.array(tmp)
164 164
165 165 return data, meta
166 166
167 167 def plot(self):
168 168
169 169 if self.xaxis == "frequency":
170 170 x = self.data.xrange[0]
171 171 self.xlabel = "Frequency (kHz)"
172 172 elif self.xaxis == "time":
173 173 x = self.data.xrange[1]
174 174 self.xlabel = "Time (ms)"
175 175 else:
176 176 x = self.data.xrange[2]
177 177 self.xlabel = "Velocity (m/s)"
178 178
179 179 self.titles = []
180 180
181 181 y = self.data.yrange
182 182 self.y = y
183 183
184 184 data = self.data[-1]
185 185 cspc = data['cspc']
186 186
187 187 for n in range(len(self.data.pairs)):
188 188
189 189 pair = self.crossPairs[n]
190 190
191 191 coh = cspc[n*2]
192 192 phase = cspc[n*2+1]
193 193 ax = self.axes[2 * n]
194 194
195 195 if ax.firsttime:
196 196 ax.plt = ax.pcolormesh(x, y, coh.T,
197 197 vmin=self.zmin_coh,
198 198 vmax=self.zmax_coh,
199 199 cmap=plt.get_cmap(self.colormap_coh)
200 200 )
201 201 else:
202 202 ax.plt.set_array(coh.T.ravel())
203 203 self.titles.append(
204 204 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
205 205
206 206 ax = self.axes[2 * n + 1]
207 207 if ax.firsttime:
208 208 ax.plt = ax.pcolormesh(x, y, phase.T,
209 209 vmin=-180,
210 210 vmax=180,
211 211 cmap=plt.get_cmap(self.colormap_phase)
212 212 )
213 213 else:
214 214 ax.plt.set_array(phase.T.ravel())
215 215
216 216 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
217 217
218 218
219 219 class RTIPlot(Plot):
220 220 '''
221 221 Plot for RTI data
222 222 '''
223 223
224 224 CODE = 'rti'
225 225 colormap = 'jet'
226 226 plot_type = 'pcolorbuffer'
227 227 titles = None
228 228 channelList = []
229 229
230 230 def setup(self):
231 231 self.xaxis = 'time'
232 232 self.ncols = 1
233 233 #print("dataChannels ",self.data.channels)
234 234 self.nrows = len(self.data.channels)
235 235 self.nplots = len(self.data.channels)
236 236 self.ylabel = 'Range [km]'
237 237 self.xlabel = 'Time'
238 238 self.cb_label = 'dB'
239 239 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
240 240 self.titles = ['{} Channel {}'.format(
241 241 self.CODE.upper(), x) for x in range(self.nplots)]
242 242
243 243 def update_list(self,dataOut):
244 244
245 245 self.channelList = dataOut.channelList
246 246
247 247
248 248 def update(self, dataOut):
249 249 if len(self.channelList) == 0:
250 250 self.update_list(dataOut)
251 251 data = {}
252 252 meta = {}
253 253 data['rti'] = dataOut.getPower()
254 254 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
255 255 return data, meta
256 256
257 257 def plot(self):
258 258
259 259 self.x = self.data.times
260 260 self.y = self.data.yrange
261
261 262 self.z = self.data[self.CODE]
262 263 self.z = numpy.array(self.z, dtype=float)
263 264 self.z = numpy.ma.masked_invalid(self.z)
264 265
265 266 try:
266 267 if self.channelList != None:
267 268 self.titles = ['{} Channel {}'.format(
268 269 self.CODE.upper(), x) for x in self.channelList]
269 270 except:
270 271 if self.channelList.any() != None:
271 272 self.titles = ['{} Channel {}'.format(
272 273 self.CODE.upper(), x) for x in self.channelList]
273 274 if self.decimation is None:
274 275 x, y, z = self.fill_gaps(self.x, self.y, self.z)
275 276 else:
276 277 x, y, z = self.fill_gaps(*self.decimate())
277 dummy_var = self.axes #ExtraΓ±amente esto actualiza el valor axes
278
279 #dummy_var = self.axes #ExtraΓ±amente esto actualiza el valor axes
278 280 for n, ax in enumerate(self.axes):
279 281 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
280 282 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
281 283 data = self.data[-1]
284
282 285 if ax.firsttime:
283 286 ax.plt = ax.pcolormesh(x, y, z[n].T,
284 287 vmin=self.zmin,
285 288 vmax=self.zmax,
286 289 cmap=plt.get_cmap(self.colormap)
287 290 )
288 291 if self.showprofile:
289 292 ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0]
290
291 293 if "noise" in self.data:
294
292 295 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
293 296 color="k", linestyle="dashed", lw=1)[0]
294 297 else:
295 298 ax.collections.remove(ax.collections[0])
296 299 ax.plt = ax.pcolormesh(x, y, z[n].T,
297 300 vmin=self.zmin,
298 301 vmax=self.zmax,
299 302 cmap=plt.get_cmap(self.colormap)
300 303 )
301 304 if self.showprofile:
302 305 ax.plot_profile.set_data(data[self.CODE][n], self.y)
303 306 if "noise" in self.data:
307
304 308 ax.plot_noise.set_data(numpy.repeat(
305 309 data['noise'][n], len(self.y)), self.y)
306 310
307 311
308 312 class CoherencePlot(RTIPlot):
309 313 '''
310 314 Plot for Coherence data
311 315 '''
312 316
313 317 CODE = 'coh'
314 318
315 319 def setup(self):
316 320 self.xaxis = 'time'
317 321 self.ncols = 1
318 322 self.nrows = len(self.data.pairs)
319 323 self.nplots = len(self.data.pairs)
320 324 self.ylabel = 'Range [km]'
321 325 self.xlabel = 'Time'
322 326 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
323 327 if self.CODE == 'coh':
324 328 self.cb_label = ''
325 329 self.titles = [
326 330 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
327 331 else:
328 332 self.cb_label = 'Degrees'
329 333 self.titles = [
330 334 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
331 335
332 336 def update(self, dataOut):
333 337 self.update_list(dataOut)
334 338 data = {}
335 339 meta = {}
336 340 data['coh'] = dataOut.getCoherence()
337 341 meta['pairs'] = dataOut.pairsList
338 342
339 343
340 344 return data, meta
341 345
342 346 class PhasePlot(CoherencePlot):
343 347 '''
344 348 Plot for Phase map data
345 349 '''
346 350
347 351 CODE = 'phase'
348 352 colormap = 'seismic'
349 353
350 354 def update(self, dataOut):
351 355
352 356 data = {}
353 357 meta = {}
354 358 data['phase'] = dataOut.getCoherence(phase=True)
355 359 meta['pairs'] = dataOut.pairsList
356 360
357 361 return data, meta
358 362
359 363 class NoisePlot(Plot):
360 364 '''
361 365 Plot for noise
362 366 '''
363 367
364 368 CODE = 'noise'
365 369 plot_type = 'scatterbuffer'
366 370
367 371 def setup(self):
368 372 self.xaxis = 'time'
369 373 self.ncols = 1
370 374 self.nrows = 1
371 375 self.nplots = 1
372 376 self.ylabel = 'Intensity [dB]'
373 377 self.xlabel = 'Time'
374 378 self.titles = ['Noise']
375 379 self.colorbar = False
376 380 self.plots_adjust.update({'right': 0.85 })
377 381
378 382 def update(self, dataOut):
379 383
380 384 data = {}
381 385 meta = {}
382 386 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
383 387 data['noise'] = noise
384 388 meta['yrange'] = numpy.array([])
385 389
386 390 return data, meta
387 391
388 392 def plot(self):
389 393
390 394 x = self.data.times
391 395 xmin = self.data.min_time
392 396 xmax = xmin + self.xrange * 60 * 60
393 397 Y = self.data['noise']
394 398
395 399 if self.axes[0].firsttime:
396 400 if self.ymin is None: self.ymin = numpy.nanmin(Y) - 5
397 401 if self.ymax is None: self.ymax = numpy.nanmax(Y) + 5
398 402 for ch in self.data.channels:
399 403 y = Y[ch]
400 404 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
401 405 plt.legend(bbox_to_anchor=(1.18, 1.0))
402 406 else:
403 407 for ch in self.data.channels:
404 408 y = Y[ch]
405 409 self.axes[0].lines[ch].set_data(x, y)
406 410
407 411
408 412 class PowerProfilePlot(Plot):
409 413
410 414 CODE = 'pow_profile'
411 415 plot_type = 'scatter'
412 416
413 417 def setup(self):
414 418
415 419 self.ncols = 1
416 420 self.nrows = 1
417 421 self.nplots = 1
418 422 self.height = 4
419 423 self.width = 3
420 424 self.ylabel = 'Range [km]'
421 425 self.xlabel = 'Intensity [dB]'
422 426 self.titles = ['Power Profile']
423 427 self.colorbar = False
424 428
425 429 def update(self, dataOut):
426 430
427 431 data = {}
428 432 meta = {}
429 433 data[self.CODE] = dataOut.getPower()
430 434
431 435 return data, meta
432 436
433 437 def plot(self):
434 438
435 439 y = self.data.yrange
436 440 self.y = y
437 441
438 442 x = self.data[-1][self.CODE]
439 443
440 444 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
441 445 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
442 446
443 447 if self.axes[0].firsttime:
444 448 for ch in self.data.channels:
445 449 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
446 450 plt.legend()
447 451 else:
448 452 for ch in self.data.channels:
449 453 self.axes[0].lines[ch].set_data(x[ch], y)
450 454
451 455
452 456 class SpectraCutPlot(Plot):
453 457
454 458 CODE = 'spc_cut'
455 459 plot_type = 'scatter'
456 460 buffering = False
457 461 heights = []
458 462 channelList = []
459 463 maintitle = "Spectra Cuts"
460 464
461 465 def setup(self):
462 466
463 467 self.nplots = len(self.data.channels)
464 468 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
465 469 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
466 470 self.width = 3.6 * self.ncols + 1.5
467 471 self.height = 3 * self.nrows
468 472 self.ylabel = 'Power [dB]'
469 473 self.colorbar = False
470 474 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
471 475 if self.selectedHeight:
472 476 self.maintitle = "Spectra Cut for %d km " %(int(self.selectedHeight))
473 477
474 478 def update(self, dataOut):
475 479 if len(self.channelList) == 0:
476 480 self.channelList = dataOut.channelList
477 481
478 482 self.heights = dataOut.heightList
479 483 if self.selectedHeight:
480 484 index_list = numpy.where(self.heights >= self.selectedHeight)
481 485 self.height_index = index_list[0]
482 486 self.height_index = self.height_index[0]
483 487 #print(self.height_index)
484 488 data = {}
485 489 meta = {}
486 490 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
487 491 if self.selectedHeight:
488 492 data['spc'] = spc[:,:,self.height_index]
489 493 else:
490 494 data['spc'] = spc
491 495 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
492 496
493 497 return data, meta
494 498
495 499 def plot(self):
496 500 if self.xaxis == "frequency":
497 501 x = self.data.xrange[0][1:]
498 502 self.xlabel = "Frequency (kHz)"
499 503 elif self.xaxis == "time":
500 504 x = self.data.xrange[1]
501 505 self.xlabel = "Time (ms)"
502 506 else:
503 507 x = self.data.xrange[2]
504 508 self.xlabel = "Velocity (m/s)"
505 509
506 510 self.titles = []
507 511
508 512 y = self.data.yrange
509 513 z = self.data[-1]['spc']
510 514 #print(z.shape)
511 515 if self.height_index:
512 516 index = numpy.array(self.height_index)
513 517 else:
514 518 index = numpy.arange(0, len(y), int((len(y))/9))
515 519
516 520 for n, ax in enumerate(self.axes):
517 521 if ax.firsttime:
518 522 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
519 523 self.xmin = self.xmin if self.xmin else -self.xmax
520 524 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
521 525 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
522 526 if self.selectedHeight:
523 527 ax.plt = ax.plot(x, z[n,:])
524 528
525 529 else:
526 530 ax.plt = ax.plot(x, z[n, :, index].T)
527 531 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
528 532 self.figures[0].legend(ax.plt, labels, loc='center right')
529 533 else:
530 534 for i, line in enumerate(ax.plt):
531 535 if self.selectedHeight:
532 536 line.set_data(x, z[n, :])
533 537 else:
534 538 line.set_data(x, z[n, :, index[i]])
535 539 self.titles.append('CH {}'.format(self.channelList[n]))
536 540 plt.suptitle(self.maintitle)
537 541
538 542 class BeaconPhase(Plot):
539 543
540 544 __isConfig = None
541 545 __nsubplots = None
542 546
543 547 PREFIX = 'beacon_phase'
544 548
545 549 def __init__(self):
546 550 Plot.__init__(self)
547 551 self.timerange = 24*60*60
548 552 self.isConfig = False
549 553 self.__nsubplots = 1
550 554 self.counter_imagwr = 0
551 555 self.WIDTH = 800
552 556 self.HEIGHT = 400
553 557 self.WIDTHPROF = 120
554 558 self.HEIGHTPROF = 0
555 559 self.xdata = None
556 560 self.ydata = None
557 561
558 562 self.PLOT_CODE = BEACON_CODE
559 563
560 564 self.FTP_WEI = None
561 565 self.EXP_CODE = None
562 566 self.SUB_EXP_CODE = None
563 567 self.PLOT_POS = None
564 568
565 569 self.filename_phase = None
566 570
567 571 self.figfile = None
568 572
569 573 self.xmin = None
570 574 self.xmax = None
571 575
572 576 def getSubplots(self):
573 577
574 578 ncol = 1
575 579 nrow = 1
576 580
577 581 return nrow, ncol
578 582
579 583 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
580 584
581 585 self.__showprofile = showprofile
582 586 self.nplots = nplots
583 587
584 588 ncolspan = 7
585 589 colspan = 6
586 590 self.__nsubplots = 2
587 591
588 592 self.createFigure(id = id,
589 593 wintitle = wintitle,
590 594 widthplot = self.WIDTH+self.WIDTHPROF,
591 595 heightplot = self.HEIGHT+self.HEIGHTPROF,
592 596 show=show)
593 597
594 598 nrow, ncol = self.getSubplots()
595 599
596 600 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
597 601
598 602 def save_phase(self, filename_phase):
599 603 f = open(filename_phase,'w+')
600 604 f.write('\n\n')
601 605 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
602 606 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
603 607 f.close()
604 608
605 609 def save_data(self, filename_phase, data, data_datetime):
606 610 f=open(filename_phase,'a')
607 611 timetuple_data = data_datetime.timetuple()
608 612 day = str(timetuple_data.tm_mday)
609 613 month = str(timetuple_data.tm_mon)
610 614 year = str(timetuple_data.tm_year)
611 615 hour = str(timetuple_data.tm_hour)
612 616 minute = str(timetuple_data.tm_min)
613 617 second = str(timetuple_data.tm_sec)
614 618 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
615 619 f.close()
616 620
617 621 def plot(self):
618 622 log.warning('TODO: Not yet implemented...')
619 623
620 624 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
621 625 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
622 626 timerange=None,
623 627 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
624 628 server=None, folder=None, username=None, password=None,
625 629 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
626 630
627 631 if dataOut.flagNoData:
628 632 return dataOut
629 633
630 634 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
631 635 return
632 636
633 637 if pairsList == None:
634 638 pairsIndexList = dataOut.pairsIndexList[:10]
635 639 else:
636 640 pairsIndexList = []
637 641 for pair in pairsList:
638 642 if pair not in dataOut.pairsList:
639 643 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
640 644 pairsIndexList.append(dataOut.pairsList.index(pair))
641 645
642 646 if pairsIndexList == []:
643 647 return
644 648
645 649 # if len(pairsIndexList) > 4:
646 650 # pairsIndexList = pairsIndexList[0:4]
647 651
648 652 hmin_index = None
649 653 hmax_index = None
650 654
651 655 if hmin != None and hmax != None:
652 656 indexes = numpy.arange(dataOut.nHeights)
653 657 hmin_list = indexes[dataOut.heightList >= hmin]
654 658 hmax_list = indexes[dataOut.heightList <= hmax]
655 659
656 660 if hmin_list.any():
657 661 hmin_index = hmin_list[0]
658 662
659 663 if hmax_list.any():
660 664 hmax_index = hmax_list[-1]+1
661 665
662 666 x = dataOut.getTimeRange()
663 667
664 668 thisDatetime = dataOut.datatime
665 669
666 670 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
667 671 xlabel = "Local Time"
668 672 ylabel = "Phase (degrees)"
669 673
670 674 update_figfile = False
671 675
672 676 nplots = len(pairsIndexList)
673 677 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
674 678 phase_beacon = numpy.zeros(len(pairsIndexList))
675 679 for i in range(nplots):
676 680 pair = dataOut.pairsList[pairsIndexList[i]]
677 681 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
678 682 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
679 683 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
680 684 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
681 685 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
682 686
683 687 if dataOut.beacon_heiIndexList:
684 688 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
685 689 else:
686 690 phase_beacon[i] = numpy.average(phase)
687 691
688 692 if not self.isConfig:
689 693
690 694 nplots = len(pairsIndexList)
691 695
692 696 self.setup(id=id,
693 697 nplots=nplots,
694 698 wintitle=wintitle,
695 699 showprofile=showprofile,
696 700 show=show)
697 701
698 702 if timerange != None:
699 703 self.timerange = timerange
700 704
701 705 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
702 706
703 707 if ymin == None: ymin = 0
704 708 if ymax == None: ymax = 360
705 709
706 710 self.FTP_WEI = ftp_wei
707 711 self.EXP_CODE = exp_code
708 712 self.SUB_EXP_CODE = sub_exp_code
709 713 self.PLOT_POS = plot_pos
710 714
711 715 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
712 716 self.isConfig = True
713 717 self.figfile = figfile
714 718 self.xdata = numpy.array([])
715 719 self.ydata = numpy.array([])
716 720
717 721 update_figfile = True
718 722
719 723 #open file beacon phase
720 724 path = '%s%03d' %(self.PREFIX, self.id)
721 725 beacon_file = os.path.join(path,'%s.txt'%self.name)
722 726 self.filename_phase = os.path.join(figpath,beacon_file)
723 727 #self.save_phase(self.filename_phase)
724 728
725 729
726 730 #store data beacon phase
727 731 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
728 732
729 733 self.setWinTitle(title)
730 734
731 735
732 736 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
733 737
734 738 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
735 739
736 740 axes = self.axesList[0]
737 741
738 742 self.xdata = numpy.hstack((self.xdata, x[0:1]))
739 743
740 744 if len(self.ydata)==0:
741 745 self.ydata = phase_beacon.reshape(-1,1)
742 746 else:
743 747 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
744 748
745 749
746 750 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
747 751 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
748 752 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
749 753 XAxisAsTime=True, grid='both'
750 754 )
751 755
752 756 self.draw()
753 757
754 758 if dataOut.ltctime >= self.xmax:
755 759 self.counter_imagwr = wr_period
756 760 self.isConfig = False
757 761 update_figfile = True
758 762
759 763 self.save(figpath=figpath,
760 764 figfile=figfile,
761 765 save=save,
762 766 ftp=ftp,
763 767 wr_period=wr_period,
764 768 thisDatetime=thisDatetime,
765 769 update_figfile=update_figfile)
766 770
767 771 return dataOut
768 772
769 773 class NoiselessSpectraPlot(Plot):
770 774 '''
771 775 Plot for Spectra data, subtracting
772 776 the noise in all channels, using for
773 777 amisr-14 data
774 778 '''
775 779
776 780 CODE = 'noiseless_spc'
777 781 colormap = 'nipy_spectral'
778 782 plot_type = 'pcolor'
779 783 buffering = False
780 784 channelList = []
781 785
782 786 def setup(self):
783 787
784 788 self.nplots = len(self.data.channels)
785 789 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
786 790 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
787 791 self.height = 2.6 * self.nrows
788 792
789 793 self.cb_label = 'dB'
790 794 if self.showprofile:
791 795 self.width = 4 * self.ncols
792 796 else:
793 797 self.width = 3.5 * self.ncols
794 798 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
795 799 self.ylabel = 'Range [km]'
796 800
797 801
798 802 def update_list(self,dataOut):
799 803 if len(self.channelList) == 0:
800 804 self.channelList = dataOut.channelList
801 805
802 806 def update(self, dataOut):
803 807
804 808 self.update_list(dataOut)
805 809 data = {}
806 810 meta = {}
807 811 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
808 812 (nch, nff, nh) = dataOut.data_spc.shape
809 813 n1 = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
810 814 noise = numpy.repeat(n1,nff, axis=1).reshape((nch,nff,nh))
811 815 #print(noise.shape, "noise", noise)
812 816
813 817 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - noise
814 818
815 819 data['spc'] = spc
816 820 data['rti'] = dataOut.getPower() - n1
817 821
818 822 data['noise'] = n0
819 823 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
820 824
821 825 return data, meta
822 826
823 827 def plot(self):
824 828 if self.xaxis == "frequency":
825 829 x = self.data.xrange[0]
826 830 self.xlabel = "Frequency (kHz)"
827 831 elif self.xaxis == "time":
828 832 x = self.data.xrange[1]
829 833 self.xlabel = "Time (ms)"
830 834 else:
831 835 x = self.data.xrange[2]
832 836 self.xlabel = "Velocity (m/s)"
833 837
834 838 self.titles = []
835 839 y = self.data.yrange
836 840 self.y = y
837 841
838 842 data = self.data[-1]
839 843 z = data['spc']
840 844
841 845 for n, ax in enumerate(self.axes):
842 846 noise = data['noise'][n]
843 847
844 848 if ax.firsttime:
845 849 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
846 850 self.xmin = self.xmin if self.xmin else -self.xmax
847 851 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
848 852 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
849 853 ax.plt = ax.pcolormesh(x, y, z[n].T,
850 854 vmin=self.zmin,
851 855 vmax=self.zmax,
852 856 cmap=plt.get_cmap(self.colormap)
853 857 )
854 858
855 859 if self.showprofile:
856 860 ax.plt_profile = self.pf_axes[n].plot(
857 861 data['rti'][n], y)[0]
858 862 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
859 863 color="k", linestyle="dashed", lw=1)[0]
860 864
861 865 else:
862 866 ax.plt.set_array(z[n].T.ravel())
863 867 if self.showprofile:
864 868 ax.plt_profile.set_data(data['rti'][n], y)
865 869 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
866 870
867 871 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
868 872
869 873
870 874 class NoiselessRTIPlot(Plot):
871 875 '''
872 876 Plot for RTI data
873 877 '''
874 878
875 879 CODE = 'noiseless_rti'
876 880 colormap = 'jet'
877 881 plot_type = 'pcolorbuffer'
878 882 titles = None
879 883 channelList = []
880 884
881 885 def setup(self):
882 886 self.xaxis = 'time'
883 887 self.ncols = 1
884 888 #print("dataChannels ",self.data.channels)
885 889 self.nrows = len(self.data.channels)
886 890 self.nplots = len(self.data.channels)
887 891 self.ylabel = 'Range [km]'
888 892 self.xlabel = 'Time'
889 893 self.cb_label = 'dB'
890 894 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
891 895 self.titles = ['{} Channel {}'.format(
892 896 self.CODE.upper(), x) for x in range(self.nplots)]
893 897
894 898 def update_list(self,dataOut):
895 899
896 900 self.channelList = dataOut.channelList
897 901
898 902
899 903 def update(self, dataOut):
900 904 if len(self.channelList) == 0:
901 905 self.update_list(dataOut)
902 906 data = {}
903 907 meta = {}
904 908
905 909 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
906 910 (nch, nff, nh) = dataOut.data_spc.shape
907 911 noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
908 912
909 913
910 914 data['noiseless_rti'] = dataOut.getPower() - noise
911 915 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
912 916 return data, meta
913 917
914 918 def plot(self):
915 919
916 920 self.x = self.data.times
917 921 self.y = self.data.yrange
918 922 self.z = self.data['noiseless_rti']
919 923 self.z = numpy.array(self.z, dtype=float)
920 924 self.z = numpy.ma.masked_invalid(self.z)
921 925
922 926 try:
923 927 if self.channelList != None:
924 928 self.titles = ['{} Channel {}'.format(
925 929 self.CODE.upper(), x) for x in self.channelList]
926 930 except:
927 931 if self.channelList.any() != None:
928 932 self.titles = ['{} Channel {}'.format(
929 933 self.CODE.upper(), x) for x in self.channelList]
930 934 if self.decimation is None:
931 935 x, y, z = self.fill_gaps(self.x, self.y, self.z)
932 936 else:
933 937 x, y, z = self.fill_gaps(*self.decimate())
934 938 dummy_var = self.axes #ExtraΓ±amente esto actualiza el valor axes
935 939 for n, ax in enumerate(self.axes):
936 940 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
937 941 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
938 942 data = self.data[-1]
939 943 if ax.firsttime:
940 944 ax.plt = ax.pcolormesh(x, y, z[n].T,
941 945 vmin=self.zmin,
942 946 vmax=self.zmax,
943 947 cmap=plt.get_cmap(self.colormap)
944 948 )
945 949 if self.showprofile:
946 950 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
947 951
948 952 if "noise" in self.data:
949 953 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
950 954 color="k", linestyle="dashed", lw=1)[0]
951 955 else:
952 956 ax.collections.remove(ax.collections[0])
953 957 ax.plt = ax.pcolormesh(x, y, z[n].T,
954 958 vmin=self.zmin,
955 959 vmax=self.zmax,
956 960 cmap=plt.get_cmap(self.colormap)
957 961 )
958 962 if self.showprofile:
959 963 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
960 964 if "noise" in self.data:
961 965 ax.plot_noise.set_data(numpy.repeat(
962 966 data['noise'][n], len(self.y)), self.y)
Show file before | Show file after | Hide comments
@@ -1,685 +1,685
1 1 import os
2 2 import time
3 3 import datetime
4 4
5 5 import numpy
6 6 import h5py
7 7
8 8 import schainpy.admin
9 9 from schainpy.model.data.jrodata import *
10 10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 11 from schainpy.model.io.jroIO_base import *
12 12 from schainpy.utils import log
13 13
14 14
15 15 class HDFReader(Reader, ProcessingUnit):
16 16 """Processing unit to read HDF5 format files
17 17
18 18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 20 attributes.
21 21 It is possible to read any HDF5 file by given the structure in the `description`
22 22 parameter, also you can add extra values to metadata with the parameter `extras`.
23 23
24 24 Parameters:
25 25 -----------
26 26 path : str
27 27 Path where files are located.
28 28 startDate : date
29 29 Start date of the files
30 30 endDate : list
31 31 End date of the files
32 32 startTime : time
33 33 Start time of the files
34 34 endTime : time
35 35 End time of the files
36 36 description : dict, optional
37 37 Dictionary with the description of the HDF5 file
38 38 extras : dict, optional
39 39 Dictionary with extra metadata to be be added to `dataOut`
40 40
41 41 Examples
42 42 --------
43 43
44 44 desc = {
45 45 'Data': {
46 46 'data_output': ['u', 'v', 'w'],
47 47 'utctime': 'timestamps',
48 48 } ,
49 49 'Metadata': {
50 50 'heightList': 'heights'
51 51 }
52 52 }
53 53
54 54 desc = {
55 55 'Data': {
56 56 'data_output': 'winds',
57 57 'utctime': 'timestamps'
58 58 },
59 59 'Metadata': {
60 60 'heightList': 'heights'
61 61 }
62 62 }
63 63
64 64 extras = {
65 65 'timeZone': 300
66 66 }
67 67
68 68 reader = project.addReadUnit(
69 69 name='HDFReader',
70 70 path='/path/to/files',
71 71 startDate='2019/01/01',
72 72 endDate='2019/01/31',
73 73 startTime='00:00:00',
74 74 endTime='23:59:59',
75 75 # description=json.dumps(desc),
76 76 # extras=json.dumps(extras),
77 77 )
78 78
79 79 """
80 80
81 81 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
82 82
83 83 def __init__(self):
84 84 ProcessingUnit.__init__(self)
85 85
86 86 self.ext = ".hdf5"
87 87 self.optchar = "D"
88 88 self.meta = {}
89 89 self.data = {}
90 90 self.open_file = h5py.File
91 91 self.open_mode = 'r'
92 92 self.description = {}
93 93 self.extras = {}
94 94 self.filefmt = "*%Y%j***"
95 95 self.folderfmt = "*%Y%j"
96 96 self.utcoffset = 0
97 97
98 98 self.dataOut = Parameters()
99 99 self.dataOut.error=False ## NOTE: Importante definir esto antes inicio
100 100 self.dataOut.flagNoData = True
101 101
102 102 def setup(self, **kwargs):
103 103
104 104 self.set_kwargs(**kwargs)
105 105 if not self.ext.startswith('.'):
106 106 self.ext = '.{}'.format(self.ext)
107 107
108 108 if self.online:
109 109 log.log("Searching files in online mode...", self.name)
110 110
111 111 for nTries in range(self.nTries):
112 112 fullpath = self.searchFilesOnLine(self.path, self.startDate,
113 113 self.endDate, self.expLabel, self.ext, self.walk,
114 114 self.filefmt, self.folderfmt)
115 115 pathname, filename = os.path.split(fullpath)
116 116
117 117 try:
118 118 fullpath = next(fullpath)
119 119
120 120 except:
121 121 fullpath = None
122 122
123 123 if fullpath:
124 124 break
125 125
126 126 log.warning(
127 127 'Waiting {} sec for a valid file in {}: try {} ...'.format(
128 128 self.delay, self.path, nTries + 1),
129 129 self.name)
130 130 time.sleep(self.delay)
131 131
132 132 if not(fullpath):
133 133 raise schainpy.admin.SchainError(
134 134 'There isn\'t any valid file in {}'.format(self.path))
135 135
136 136 pathname, filename = os.path.split(fullpath)
137 137 self.year = int(filename[1:5])
138 138 self.doy = int(filename[5:8])
139 139 self.set = int(filename[8:11]) - 1
140 140 else:
141 141 log.log("Searching files in {}".format(self.path), self.name)
142 142 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
143 143 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
144 144
145 145 self.setNextFile()
146 146
147 147
148 148
149 149
150 150 def readFirstHeader(self):
151 151 '''Read metadata and data'''
152 152
153 153 self.__readMetadata()
154 154 self.__readData()
155 155 self.__setBlockList()
156 156
157 157 for attr in self.meta:
158 158 setattr(self.dataOut, attr, self.meta[attr])
159 159 self.blockIndex = 0
160 160
161 161 return
162 162
163 163 def __setBlockList(self):
164 164 '''
165 165 Selects the data within the times defined
166 166
167 167 self.fp
168 168 self.startTime
169 169 self.endTime
170 170 self.blockList
171 171 self.blocksPerFile
172 172
173 173 '''
174 174
175 175 startTime = self.startTime
176 176 endTime = self.endTime
177 177 thisUtcTime = self.data['utctime'] + self.utcoffset
178 178 self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
179 179 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
180 180 self.startFileDatetime = thisDatetime
181 181 thisDate = thisDatetime.date()
182 182 thisTime = thisDatetime.time()
183 183
184 184 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
185 185 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
186 186
187 187 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
188 188
189 189 self.blockList = ind
190 190 self.blocksPerFile = len(ind)
191 191 self.blocksPerFile = len(thisUtcTime)
192 192 return
193 193
194 194 def __readMetadata(self):
195 195 '''
196 196 Reads Metadata
197 197 '''
198 198
199 199 meta = {}
200 200
201 201 if self.description:
202 202 for key, value in self.description['Metadata'].items():
203 203 meta[key] = self.fp[value][()]
204 204 else:
205 205 grp = self.fp['Metadata']
206 206 for name in grp:
207 207 meta[name] = grp[name][()]
208 208
209 209 if self.extras:
210 210 for key, value in self.extras.items():
211 211 meta[key] = value
212 212 self.meta = meta
213 213
214 214 return
215 215
216 216
217 217
218 218 def checkForRealPath(self, nextFile, nextDay):
219 219
220 220 # print("check FRP")
221 221 # dt = self.startFileDatetime + datetime.timedelta(1)
222 222 # filename = '{}.{}{}'.format(self.path, dt.strftime('%Y%m%d'), self.ext)
223 223 # fullfilename = os.path.join(self.path, filename)
224 224 # print("check Path ",fullfilename,filename)
225 225 # if os.path.exists(fullfilename):
226 226 # return fullfilename, filename
227 227 # return None, filename
228 228 return None,None
229 229
230 230 def __readData(self):
231 231
232 232 data = {}
233 233
234 234 if self.description:
235 235 for key, value in self.description['Data'].items():
236 236 if isinstance(value, str):
237 237 if isinstance(self.fp[value], h5py.Dataset):
238 238 data[key] = self.fp[value][()]
239 239 elif isinstance(self.fp[value], h5py.Group):
240 240 array = []
241 241 for ch in self.fp[value]:
242 242 array.append(self.fp[value][ch][()])
243 243 data[key] = numpy.array(array)
244 244 elif isinstance(value, list):
245 245 array = []
246 246 for ch in value:
247 247 array.append(self.fp[ch][()])
248 248 data[key] = numpy.array(array)
249 249 else:
250 250 grp = self.fp['Data']
251 251 for name in grp:
252 252 if isinstance(grp[name], h5py.Dataset):
253 253 array = grp[name][()]
254 254 elif isinstance(grp[name], h5py.Group):
255 255 array = []
256 256 for ch in grp[name]:
257 257 array.append(grp[name][ch][()])
258 258 array = numpy.array(array)
259 259 else:
260 260 log.warning('Unknown type: {}'.format(name))
261 261
262 262 if name in self.description:
263 263 key = self.description[name]
264 264 else:
265 265 key = name
266 266 data[key] = array
267 267
268 268 self.data = data
269 269 return
270 270
271 271 def getData(self):
272 272 if not self.isDateTimeInRange(self.startFileDatetime, self.startDate, self.endDate, self.startTime, self.endTime):
273 273 self.dataOut.flagNoData = True
274 274 self.blockIndex = self.blocksPerFile
275 275 self.dataOut.error = True # TERMINA EL PROGRAMA
276 276 return
277 277 for attr in self.data:
278 278
279 279 if self.data[attr].ndim == 1:
280 280 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
281 281 else:
282 282 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
283 283
284 284
285 285 self.blockIndex += 1
286 286
287 287 if self.blockIndex == 1:
288 288 log.log("Block No. {}/{} -> {}".format(
289 289 self.blockIndex,
290 290 self.blocksPerFile,
291 291 self.dataOut.datatime.ctime()), self.name)
292 292 else:
293 293 log.log("Block No. {}/{} ".format(
294 294 self.blockIndex,
295 295 self.blocksPerFile),self.name)
296 296
297 297 if self.blockIndex == self.blocksPerFile:
298 298 self.setNextFile()
299 299
300 300 self.dataOut.flagNoData = False
301 301
302 302
303 303 def run(self, **kwargs):
304 304
305 305 if not(self.isConfig):
306 306 self.setup(**kwargs)
307 307 self.isConfig = True
308 308
309 309 self.getData()
310 310
311 #@MPDecorator
312 class HDFWrite(Operation):
311 @MPDecorator
312 class HDFWriter(Operation):
313 313 """Operation to write HDF5 files.
314 314
315 315 The HDF5 file contains by default two groups Data and Metadata where
316 316 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
317 317 parameters, data attributes are normaly time dependent where the metadata
318 318 are not.
319 319 It is possible to customize the structure of the HDF5 file with the
320 320 optional description parameter see the examples.
321 321
322 322 Parameters:
323 323 -----------
324 324 path : str
325 325 Path where files will be saved.
326 326 blocksPerFile : int
327 327 Number of blocks per file
328 328 metadataList : list
329 329 List of the dataOut attributes that will be saved as metadata
330 330 dataList : int
331 331 List of the dataOut attributes that will be saved as data
332 332 setType : bool
333 333 If True the name of the files corresponds to the timestamp of the data
334 334 description : dict, optional
335 335 Dictionary with the desired description of the HDF5 file
336 336
337 337 Examples
338 338 --------
339 339
340 340 desc = {
341 341 'data_output': {'winds': ['z', 'w', 'v']},
342 342 'utctime': 'timestamps',
343 343 'heightList': 'heights'
344 344 }
345 345 desc = {
346 346 'data_output': ['z', 'w', 'v'],
347 347 'utctime': 'timestamps',
348 348 'heightList': 'heights'
349 349 }
350 350 desc = {
351 351 'Data': {
352 352 'data_output': 'winds',
353 353 'utctime': 'timestamps'
354 354 },
355 355 'Metadata': {
356 356 'heightList': 'heights'
357 357 }
358 358 }
359 359
360 360 writer = proc_unit.addOperation(name='HDFWriter')
361 361 writer.addParameter(name='path', value='/path/to/file')
362 362 writer.addParameter(name='blocksPerFile', value='32')
363 363 writer.addParameter(name='metadataList', value='heightList,timeZone')
364 364 writer.addParameter(name='dataList',value='data_output,utctime')
365 365 # writer.addParameter(name='description',value=json.dumps(desc))
366 366
367 367 """
368 368
369 369 ext = ".hdf5"
370 370 optchar = "D"
371 371 filename = None
372 372 path = None
373 373 setFile = None
374 374 fp = None
375 375 firsttime = True
376 376 #Configurations
377 377 blocksPerFile = None
378 378 blockIndex = None
379 379 dataOut = None #eval ??????
380 380 #Data Arrays
381 381 dataList = None
382 382 metadataList = None
383 383 currentDay = None
384 384 lastTime = None
385 385 timeZone = "ut"
386 386 hourLimit = 3
387 387 breakDays = True
388 388
389 389 def __init__(self):
390 390
391 391 Operation.__init__(self)
392 392
393 393
394 394 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None,
395 395 description={},timeZone = "ut",hourLimit = 3, breakDays=True):
396 396 self.path = path
397 397 self.blocksPerFile = blocksPerFile
398 398 self.metadataList = metadataList
399 399 self.dataList = [s.strip() for s in dataList]
400 400 self.setType = setType
401 401 self.description = description
402 402 self.timeZone = timeZone
403 403 self.hourLimit = hourLimit
404 404 self.breakDays = breakDays
405 405
406 406 if self.metadataList is None:
407 407 self.metadataList = self.dataOut.metadata_list
408 408
409 409 tableList = []
410 410 dsList = []
411 411
412 412 for i in range(len(self.dataList)):
413 413 dsDict = {}
414 414 if hasattr(self.dataOut, self.dataList[i]):
415 415 dataAux = getattr(self.dataOut, self.dataList[i])
416 416 dsDict['variable'] = self.dataList[i]
417 417 else:
418 418 log.warning('Attribute {} not found in dataOut'.format(self.dataList[i]),self.name)
419 419 continue
420 420
421 421 if dataAux is None:
422 422 continue
423 423 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
424 424 dsDict['nDim'] = 0
425 425 else:
426 426 dsDict['nDim'] = len(dataAux.shape)
427 427 dsDict['shape'] = dataAux.shape
428 428 dsDict['dsNumber'] = dataAux.shape[0]
429 429 dsDict['dtype'] = dataAux.dtype
430 430
431 431 dsList.append(dsDict)
432 432
433 433 self.blockIndex = 0
434 434 self.dsList = dsList
435 435 self.currentDay = self.dataOut.datatime.date()
436 436
437 437
438 438 def timeFlag(self):
439 439 currentTime = self.dataOut.utctime
440 440 timeTuple = None
441 441 if self.timeZone == "lt":
442 442 timeTuple = time.localtime(currentTime)
443 443 else :
444 444 timeTuple = time.gmtime(currentTime)
445 445
446 446 dataDay = timeTuple.tm_yday
447 447
448 448 if self.lastTime is None:
449 449 self.lastTime = currentTime
450 450 self.currentDay = dataDay
451 451 return False
452 452
453 453 timeDiff = currentTime - self.lastTime
454 454
455 455 #Si el dia es diferente o si la diferencia entre un dato y otro supera self.hourLimit
456 456 if (dataDay != self.currentDay) and self.breakDays:
457 457 self.currentDay = dataDay
458 458 return True
459 459 elif timeDiff > self.hourLimit*60*60:
460 460 self.lastTime = currentTime
461 461 return True
462 462 else:
463 463 self.lastTime = currentTime
464 464 return False
465 465
466 466 def run(self, dataOut,**kwargs):
467 467
468 self.dataOut = dataOut
468 self.dataOut = dataOut.copy()
469 469 if not(self.isConfig):
470 470 self.setup(**kwargs)
471 471
472 472 self.isConfig = True
473 473 self.setNextFile()
474 474
475 475 self.putData()
476 476
477 return self.dataOut
477 #return self.dataOut
478 478
479 479 def setNextFile(self):
480 480
481 481 ext = self.ext
482 482 path = self.path
483 483 setFile = self.setFile
484 484 timeTuple = None
485 485 if self.timeZone == "lt":
486 486 timeTuple = time.localtime(self.dataOut.utctime)
487 487 elif self.timeZone == "ut":
488 488 timeTuple = time.gmtime(self.dataOut.utctime)
489 489 #print("path: ",timeTuple)
490 490 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
491 491 fullpath = os.path.join(path, subfolder)
492 492
493 493 if os.path.exists(fullpath):
494 494 filesList = os.listdir(fullpath)
495 495 filesList = [k for k in filesList if k.startswith(self.optchar)]
496 496 if len( filesList ) > 0:
497 497 filesList = sorted(filesList, key=str.lower)
498 498 filen = filesList[-1]
499 499 # el filename debera tener el siguiente formato
500 500 # 0 1234 567 89A BCDE (hex)
501 501 # x YYYY DDD SSS .ext
502 502 if isNumber(filen[8:11]):
503 503 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
504 504 else:
505 505 setFile = -1
506 506 else:
507 507 setFile = -1 #inicializo mi contador de seteo
508 508 else:
509 509 os.makedirs(fullpath)
510 510 setFile = -1 #inicializo mi contador de seteo
511 511
512 512 if self.setType is None:
513 513 setFile += 1
514 514 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
515 515 timeTuple.tm_year,
516 516 timeTuple.tm_yday,
517 517 setFile,
518 518 ext )
519 519 else:
520 520 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
521 521 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
522 522 timeTuple.tm_year,
523 523 timeTuple.tm_yday,
524 524 setFile,
525 525 ext )
526 526
527 527 self.filename = os.path.join( path, subfolder, file )
528 528
529 529
530 530
531 531 def getLabel(self, name, x=None):
532 532
533 533 if x is None:
534 534 if 'Data' in self.description:
535 535 data = self.description['Data']
536 536 if 'Metadata' in self.description:
537 537 data.update(self.description['Metadata'])
538 538 else:
539 539 data = self.description
540 540 if name in data:
541 541 if isinstance(data[name], str):
542 542 return data[name]
543 543 elif isinstance(data[name], list):
544 544 return None
545 545 elif isinstance(data[name], dict):
546 546 for key, value in data[name].items():
547 547 return key
548 548 return name
549 549 else:
550 550 if 'Metadata' in self.description:
551 551 meta = self.description['Metadata']
552 552 else:
553 553 meta = self.description
554 554 if name in meta:
555 555 if isinstance(meta[name], list):
556 556 return meta[name][x]
557 557 elif isinstance(meta[name], dict):
558 558 for key, value in meta[name].items():
559 559 return value[x]
560 560 if 'cspc' in name:
561 561 return 'pair{:02d}'.format(x)
562 562 else:
563 563 return 'channel{:02d}'.format(x)
564 564
565 565 def writeMetadata(self, fp):
566 566
567 567 if self.description:
568 568 if 'Metadata' in self.description:
569 569 grp = fp.create_group('Metadata')
570 570 else:
571 571 grp = fp
572 572 else:
573 573 grp = fp.create_group('Metadata')
574 574
575 575 for i in range(len(self.metadataList)):
576 576 if not hasattr(self.dataOut, self.metadataList[i]):
577 577 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
578 578 continue
579 579 value = getattr(self.dataOut, self.metadataList[i])
580 580 if isinstance(value, bool):
581 581 if value is True:
582 582 value = 1
583 583 else:
584 584 value = 0
585 585 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
586 586 return
587 587
588 588 def writeData(self, fp):
589 589
590 590 if self.description:
591 591 if 'Data' in self.description:
592 592 grp = fp.create_group('Data')
593 593 else:
594 594 grp = fp
595 595 else:
596 596 grp = fp.create_group('Data')
597 597
598 598 dtsets = []
599 599 data = []
600 600
601 601 for dsInfo in self.dsList:
602 602 if dsInfo['nDim'] == 0:
603 603 ds = grp.create_dataset(
604 604 self.getLabel(dsInfo['variable']),
605 605 (self.blocksPerFile, ),
606 606 chunks=True,
607 607 dtype=numpy.float64)
608 608 dtsets.append(ds)
609 609 data.append((dsInfo['variable'], -1))
610 610 else:
611 611 label = self.getLabel(dsInfo['variable'])
612 612 if label is not None:
613 613 sgrp = grp.create_group(label)
614 614 else:
615 615 sgrp = grp
616 616 for i in range(dsInfo['dsNumber']):
617 617 ds = sgrp.create_dataset(
618 618 self.getLabel(dsInfo['variable'], i),
619 619 (self.blocksPerFile, ) + dsInfo['shape'][1:],
620 620 chunks=True,
621 621 dtype=dsInfo['dtype'])
622 622 dtsets.append(ds)
623 623 data.append((dsInfo['variable'], i))
624 624 fp.flush()
625 625
626 626 log.log('Creating file: {}'.format(fp.filename), self.name)
627 627
628 628 self.ds = dtsets
629 629 self.data = data
630 630 self.firsttime = True
631 631
632 632 return
633 633
634 634 def putData(self):
635 635
636 636 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
637 637 self.closeFile()
638 638 self.setNextFile()
639 639 self.dataOut.flagNoData = False
640 640 self.blockIndex = 0
641 641 return
642 642
643 643
644 644
645 645 if self.blockIndex == 0:
646 646 #Escribir metadata Aqui???
647 647 #Setting HDF5 File
648 648 self.fp = h5py.File(self.filename, 'w')
649 649 #write metadata
650 650 self.writeMetadata(self.fp)
651 651 #Write data
652 652 self.writeData(self.fp)
653 653 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
654 654 elif (self.blockIndex % 10 ==0):
655 655 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
656 656 else:
657 657
658 658 log.log('Block No. {}/{}'.format(self.blockIndex+1, self.blocksPerFile), self.name)
659 659
660 660 for i, ds in enumerate(self.ds):
661 661 attr, ch = self.data[i]
662 662 if ch == -1:
663 663 ds[self.blockIndex] = getattr(self.dataOut, attr)
664 664 else:
665 665 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
666 666
667 667 self.blockIndex += 1
668 668
669 669 self.fp.flush()
670 670 self.dataOut.flagNoData = True
671 671
672 672
673 673 def closeFile(self):
674 674
675 675 if self.blockIndex != self.blocksPerFile:
676 676 for ds in self.ds:
677 677 ds.resize(self.blockIndex, axis=0)
678 678
679 679 if self.fp:
680 680 self.fp.flush()
681 681 self.fp.close()
682 682
683 683 def close(self):
684 684
685 685 self.closeFile()
Show file before | Show file after | Hide comments
@@ -1,211 +1,211
1 1 '''
2 2 Base clases to create Processing units and operations, the MPDecorator
3 3 must be used in plotting and writing operations to allow to run as an
4 4 external process.
5 5 '''
6 6 import os
7 7 import inspect
8 8 import zmq
9 9 import time
10 10 import pickle
11 11 import traceback
12 12 from threading import Thread
13 13 from multiprocessing import Process, Queue
14 14 from schainpy.utils import log
15 15 import copy
16 16 QUEUE_SIZE = int(os.environ.get('QUEUE_MAX_SIZE', '10'))
17 17 class ProcessingUnit(object):
18 18 '''
19 19 Base class to create Signal Chain Units
20 20 '''
21 21
22 22 proc_type = 'processing'
23 23
24 24 def __init__(self):
25 25
26 26 self.dataIn = None
27 27 self.dataOut = None
28 28 self.isConfig = False
29 29 self.operations = []
30 30
31 31 def setInput(self, unit):
32 32
33 33 self.dataIn = unit.dataOut
34 34
35 35
36 36 def getAllowedArgs(self):
37 37 if hasattr(self, '__attrs__'):
38 38 return self.__attrs__
39 39 else:
40 40 return inspect.getargspec(self.run).args
41 41
42 42 def addOperation(self, conf, operation):
43 43 '''
44 44 '''
45 45
46 46 self.operations.append((operation, conf.type, conf.getKwargs()))
47 47
48 48 def getOperationObj(self, objId):
49 49
50 50 if objId not in list(self.operations.keys()):
51 51 return None
52 52
53 53 return self.operations[objId]
54 54
55 55 def call(self, **kwargs):
56 56 '''
57 57 '''
58 58
59 59 try:
60 60
61 61 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
62 62 return self.dataIn.isReady()
63 63 elif self.dataIn is None or not self.dataIn.error:
64 64 self.run(**kwargs)
65 65 elif self.dataIn.error:
66 66 self.dataOut.error = self.dataIn.error
67 67 self.dataOut.flagNoData = True
68 68 except:
69 69
70 70 err = traceback.format_exc()
71 71 if 'SchainWarning' in err:
72 72 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), self.name)
73 73 elif 'SchainError' in err:
74 74 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), self.name)
75 75 else:
76 76 log.error(err, self.name)
77 77 self.dataOut.error = True
78 78
79 79 for op, optype, opkwargs in self.operations:
80 80 if optype == 'other' and self.dataOut.isReady():
81 81 try:
82 82 self.dataOut = op.run(self.dataOut, **opkwargs)
83 83 except Exception as e:
84 84 print(e)
85 85 self.dataOut.error = True
86 86 return 'Error'
87 87 elif optype == 'external' and self.dataOut.isReady() :
88 88 op.queue.put(copy.deepcopy(self.dataOut))
89 89 elif optype == 'external' and self.dataOut.error:
90 90 op.queue.put(copy.deepcopy(self.dataOut))
91 91
92 return 'Error' if self.dataOut.error else True#self.dataOut.isReady()
92 return 'Error' if self.dataOut.error else self.dataOut.isReady()
93 93
94 94 def setup(self):
95 95
96 96 raise NotImplementedError
97 97
98 98 def run(self):
99 99
100 100 raise NotImplementedError
101 101
102 102 def close(self):
103 103
104 104 return
105 105
106 106
107 107 class Operation(object):
108 108
109 109 '''
110 110 '''
111 111
112 112 proc_type = 'operation'
113 113
114 114 def __init__(self):
115 115
116 116 self.id = None
117 117 self.isConfig = False
118 118
119 119 if not hasattr(self, 'name'):
120 120 self.name = self.__class__.__name__
121 121
122 122 def getAllowedArgs(self):
123 123 if hasattr(self, '__attrs__'):
124 124 return self.__attrs__
125 125 else:
126 126 return inspect.getargspec(self.run).args
127 127
128 128 def setup(self):
129 129
130 130 self.isConfig = True
131 131
132 132 raise NotImplementedError
133 133
134 134 def run(self, dataIn, **kwargs):
135 135 """
136 136 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
137 137 atributos del objeto dataIn.
138 138
139 139 Input:
140 140
141 141 dataIn : objeto del tipo JROData
142 142
143 143 Return:
144 144
145 145 None
146 146
147 147 Affected:
148 148 __buffer : buffer de recepcion de datos.
149 149
150 150 """
151 151 if not self.isConfig:
152 152 self.setup(**kwargs)
153 153
154 154 raise NotImplementedError
155 155
156 156 def close(self):
157 157
158 158 return
159 159
160 160
161 161 def MPDecorator(BaseClass):
162 162 """
163 163 Multiprocessing class decorator
164 164
165 165 This function add multiprocessing features to a BaseClass.
166 166 """
167 167
168 168 class MPClass(BaseClass, Process):
169 169
170 170 def __init__(self, *args, **kwargs):
171 171 super(MPClass, self).__init__()
172 172 Process.__init__(self)
173 173
174 174 self.args = args
175 175 self.kwargs = kwargs
176 176 self.t = time.time()
177 177 self.op_type = 'external'
178 178 self.name = BaseClass.__name__
179 179 self.__doc__ = BaseClass.__doc__
180 180
181 181 if 'plot' in self.name.lower() and not self.name.endswith('_'):
182 182 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
183 183
184 184 self.start_time = time.time()
185 185 self.err_queue = args[3]
186 186 self.queue = Queue(maxsize=QUEUE_SIZE)
187 187 self.myrun = BaseClass.run
188 188
189 189 def run(self):
190 190
191 191 while True:
192 192
193 193 dataOut = self.queue.get()
194 194
195 195 if not dataOut.error:
196 196 try:
197 197 BaseClass.run(self, dataOut, **self.kwargs)
198 198 except:
199 199 err = traceback.format_exc()
200 200 log.error(err, self.name)
201 201 else:
202 202 break
203 203
204 204 self.close()
205 205
206 206 def close(self):
207 207
208 208 BaseClass.close(self)
209 209 log.success('Done...(Time:{:4.2f} secs)'.format(time.time()-self.start_time), self.name)
210 210
211 211 return MPClass
Show file before | Show file after | Hide comments
@@ -1,1814 +1,1823
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15 import math
16 16
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 18 from schainpy.model.data.jrodata import Spectra
19 19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 20 from schainpy.utils import log
21 21
22 22 from scipy.optimize import curve_fit
23 23
24 24 class SpectraProc(ProcessingUnit):
25 25
26 26 def __init__(self):
27 27
28 28 ProcessingUnit.__init__(self)
29 29
30 30 self.buffer = None
31 31 self.firstdatatime = None
32 32 self.profIndex = 0
33 33 self.dataOut = Spectra()
34 34 self.id_min = None
35 35 self.id_max = None
36 36 self.setupReq = False #Agregar a todas las unidades de proc
37 37
38 38 def __updateSpecFromVoltage(self):
39 39
40 40 self.dataOut.timeZone = self.dataIn.timeZone
41 41 self.dataOut.dstFlag = self.dataIn.dstFlag
42 42 self.dataOut.errorCount = self.dataIn.errorCount
43 43 self.dataOut.useLocalTime = self.dataIn.useLocalTime
44 44 try:
45 45 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
46 46 except:
47 47 pass
48 48 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
49 49 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
50 50 self.dataOut.channelList = self.dataIn.channelList
51 51 self.dataOut.heightList = self.dataIn.heightList
52 52 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
53 53 self.dataOut.nProfiles = self.dataOut.nFFTPoints
54 54 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
55 55 self.dataOut.utctime = self.firstdatatime
56 56 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
57 57 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
58 58 self.dataOut.flagShiftFFT = False
59 59 self.dataOut.nCohInt = self.dataIn.nCohInt
60 60 self.dataOut.nIncohInt = 1
61 61 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
62 62 self.dataOut.frequency = self.dataIn.frequency
63 63 self.dataOut.realtime = self.dataIn.realtime
64 64 self.dataOut.azimuth = self.dataIn.azimuth
65 65 self.dataOut.zenith = self.dataIn.zenith
66 66 self.dataOut.codeList = self.dataIn.codeList
67 67 self.dataOut.azimuthList = self.dataIn.azimuthList
68 68 self.dataOut.elevationList = self.dataIn.elevationList
69 69
70 70
71 71 def __getFft(self):
72 72 """
73 73 Convierte valores de Voltaje a Spectra
74 74
75 75 Affected:
76 76 self.dataOut.data_spc
77 77 self.dataOut.data_cspc
78 78 self.dataOut.data_dc
79 79 self.dataOut.heightList
80 80 self.profIndex
81 81 self.buffer
82 82 self.dataOut.flagNoData
83 83 """
84 84 fft_volt = numpy.fft.fft(
85 85 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
86 86 fft_volt = fft_volt.astype(numpy.dtype('complex'))
87 87 dc = fft_volt[:, 0, :]
88 88
89 89 # calculo de self-spectra
90 90 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
91 91 spc = fft_volt * numpy.conjugate(fft_volt)
92 92 spc = spc.real
93 93
94 94 blocksize = 0
95 95 blocksize += dc.size
96 96 blocksize += spc.size
97 97
98 98 cspc = None
99 99 pairIndex = 0
100 100 if self.dataOut.pairsList != None:
101 101 # calculo de cross-spectra
102 102 cspc = numpy.zeros(
103 103 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
104 104 for pair in self.dataOut.pairsList:
105 105 if pair[0] not in self.dataOut.channelList:
106 106 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
107 107 str(pair), str(self.dataOut.channelList)))
108 108 if pair[1] not in self.dataOut.channelList:
109 109 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
110 110 str(pair), str(self.dataOut.channelList)))
111 111
112 112 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
113 113 numpy.conjugate(fft_volt[pair[1], :, :])
114 114 pairIndex += 1
115 115 blocksize += cspc.size
116 116
117 117 self.dataOut.data_spc = spc
118 118 self.dataOut.data_cspc = cspc
119 119 self.dataOut.data_dc = dc
120 120 self.dataOut.blockSize = blocksize
121 121 self.dataOut.flagShiftFFT = False
122 122
123 123 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
124 124
125 125 if self.dataIn.type == "Spectra":
126 126
127 127 try:
128 128 self.dataOut.copy(self.dataIn)
129 129
130 130 except Exception as e:
131 131 print(e)
132 132
133 133 if shift_fft:
134 134 #desplaza a la derecha en el eje 2 determinadas posiciones
135 135 shift = int(self.dataOut.nFFTPoints/2)
136 136 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
137 137
138 138 if self.dataOut.data_cspc is not None:
139 139 #desplaza a la derecha en el eje 2 determinadas posiciones
140 140 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
141 141 if pairsList:
142 142 self.__selectPairs(pairsList)
143 143
144 144
145 145 elif self.dataIn.type == "Voltage":
146 146
147 147 self.dataOut.flagNoData = True
148 148
149 149 if nFFTPoints == None:
150 150 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
151 151
152 152 if nProfiles == None:
153 153 nProfiles = nFFTPoints
154 154
155 155 if ippFactor == None:
156 156 self.dataOut.ippFactor = 1
157 157
158 158 self.dataOut.nFFTPoints = nFFTPoints
159 159
160 160 if self.buffer is None:
161 161 self.buffer = numpy.zeros((self.dataIn.nChannels,
162 162 nProfiles,
163 163 self.dataIn.nHeights),
164 164 dtype='complex')
165 165
166 166 if self.dataIn.flagDataAsBlock:
167 167 nVoltProfiles = self.dataIn.data.shape[1]
168 168
169 169 if nVoltProfiles == nProfiles:
170 170 self.buffer = self.dataIn.data.copy()
171 171 self.profIndex = nVoltProfiles
172 172
173 173 elif nVoltProfiles < nProfiles:
174 174
175 175 if self.profIndex == 0:
176 176 self.id_min = 0
177 177 self.id_max = nVoltProfiles
178 178
179 179 self.buffer[:, self.id_min:self.id_max,
180 180 :] = self.dataIn.data
181 181 self.profIndex += nVoltProfiles
182 182 self.id_min += nVoltProfiles
183 183 self.id_max += nVoltProfiles
184 184 else:
185 185 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
186 186 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
187 187 self.dataOut.flagNoData = True
188 188 else:
189 189 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
190 190 self.profIndex += 1
191 191
192 192 if self.firstdatatime == None:
193 193 self.firstdatatime = self.dataIn.utctime
194 194
195 195 if self.profIndex == nProfiles:
196 196 self.__updateSpecFromVoltage()
197 197 if pairsList == None:
198 198 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
199 199 else:
200 200 self.dataOut.pairsList = pairsList
201 201 self.__getFft()
202 202 self.dataOut.flagNoData = False
203 203 self.firstdatatime = None
204 204 self.profIndex = 0
205 self.dataOut.noise_estimation = None
205
206 206 else:
207 207 raise ValueError("The type of input object '%s' is not valid".format(
208 208 self.dataIn.type))
209 209
210 210 def __selectPairs(self, pairsList):
211 211
212 212 if not pairsList:
213 213 return
214 214
215 215 pairs = []
216 216 pairsIndex = []
217 217
218 218 for pair in pairsList:
219 219 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
220 220 continue
221 221 pairs.append(pair)
222 222 pairsIndex.append(pairs.index(pair))
223 223
224 224 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
225 225 self.dataOut.pairsList = pairs
226 226
227 227 return
228 228
229 229 def selectFFTs(self, minFFT, maxFFT ):
230 230 """
231 231 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
232 232 minFFT<= FFT <= maxFFT
233 233 """
234 234
235 235 if (minFFT > maxFFT):
236 236 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
237 237
238 238 if (minFFT < self.dataOut.getFreqRange()[0]):
239 239 minFFT = self.dataOut.getFreqRange()[0]
240 240
241 241 if (maxFFT > self.dataOut.getFreqRange()[-1]):
242 242 maxFFT = self.dataOut.getFreqRange()[-1]
243 243
244 244 minIndex = 0
245 245 maxIndex = 0
246 246 FFTs = self.dataOut.getFreqRange()
247 247
248 248 inda = numpy.where(FFTs >= minFFT)
249 249 indb = numpy.where(FFTs <= maxFFT)
250 250
251 251 try:
252 252 minIndex = inda[0][0]
253 253 except:
254 254 minIndex = 0
255 255
256 256 try:
257 257 maxIndex = indb[0][-1]
258 258 except:
259 259 maxIndex = len(FFTs)
260 260
261 261 self.selectFFTsByIndex(minIndex, maxIndex)
262 262
263 263 return 1
264 264
265 265 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
266 266 newheis = numpy.where(
267 267 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
268 268
269 269 if hei_ref != None:
270 270 newheis = numpy.where(self.dataOut.heightList > hei_ref)
271 271
272 272 minIndex = min(newheis[0])
273 273 maxIndex = max(newheis[0])
274 274 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
275 275 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
276 276
277 277 # determina indices
278 278 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
279 279 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
280 280 avg_dB = 10 * \
281 281 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
282 282 beacon_dB = numpy.sort(avg_dB)[-nheis:]
283 283 beacon_heiIndexList = []
284 284 for val in avg_dB.tolist():
285 285 if val >= beacon_dB[0]:
286 286 beacon_heiIndexList.append(avg_dB.tolist().index(val))
287 287
288 288 #data_spc = data_spc[:,:,beacon_heiIndexList]
289 289 data_cspc = None
290 290 if self.dataOut.data_cspc is not None:
291 291 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
292 292 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
293 293
294 294 data_dc = None
295 295 if self.dataOut.data_dc is not None:
296 296 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
297 297 #data_dc = data_dc[:,beacon_heiIndexList]
298 298
299 299 self.dataOut.data_spc = data_spc
300 300 self.dataOut.data_cspc = data_cspc
301 301 self.dataOut.data_dc = data_dc
302 302 self.dataOut.heightList = heightList
303 303 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
304 304
305 305 return 1
306 306
307 307 def selectFFTsByIndex(self, minIndex, maxIndex):
308 308 """
309 309
310 310 """
311 311
312 312 if (minIndex < 0) or (minIndex > maxIndex):
313 313 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
314 314
315 315 if (maxIndex >= self.dataOut.nProfiles):
316 316 maxIndex = self.dataOut.nProfiles-1
317 317
318 318 #Spectra
319 319 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
320 320
321 321 data_cspc = None
322 322 if self.dataOut.data_cspc is not None:
323 323 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
324 324
325 325 data_dc = None
326 326 if self.dataOut.data_dc is not None:
327 327 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
328 328
329 329 self.dataOut.data_spc = data_spc
330 330 self.dataOut.data_cspc = data_cspc
331 331 self.dataOut.data_dc = data_dc
332 332
333 333 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
334 334 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
335 335 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
336 336
337 337 return 1
338 338
339 339 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
340 340 # validacion de rango
341 341 if minHei == None:
342 342 minHei = self.dataOut.heightList[0]
343 343
344 344 if maxHei == None:
345 345 maxHei = self.dataOut.heightList[-1]
346 346
347 347 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
348 348 print('minHei: %.2f is out of the heights range' % (minHei))
349 349 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
350 350 minHei = self.dataOut.heightList[0]
351 351
352 352 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
353 353 print('maxHei: %.2f is out of the heights range' % (maxHei))
354 354 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
355 355 maxHei = self.dataOut.heightList[-1]
356 356
357 357 # validacion de velocidades
358 358 velrange = self.dataOut.getVelRange(1)
359 359
360 360 if minVel == None:
361 361 minVel = velrange[0]
362 362
363 363 if maxVel == None:
364 364 maxVel = velrange[-1]
365 365
366 366 if (minVel < velrange[0]) or (minVel > maxVel):
367 367 print('minVel: %.2f is out of the velocity range' % (minVel))
368 368 print('minVel is setting to %.2f' % (velrange[0]))
369 369 minVel = velrange[0]
370 370
371 371 if (maxVel > velrange[-1]) or (maxVel < minVel):
372 372 print('maxVel: %.2f is out of the velocity range' % (maxVel))
373 373 print('maxVel is setting to %.2f' % (velrange[-1]))
374 374 maxVel = velrange[-1]
375 375
376 376 # seleccion de indices para rango
377 377 minIndex = 0
378 378 maxIndex = 0
379 379 heights = self.dataOut.heightList
380 380
381 381 inda = numpy.where(heights >= minHei)
382 382 indb = numpy.where(heights <= maxHei)
383 383
384 384 try:
385 385 minIndex = inda[0][0]
386 386 except:
387 387 minIndex = 0
388 388
389 389 try:
390 390 maxIndex = indb[0][-1]
391 391 except:
392 392 maxIndex = len(heights)
393 393
394 394 if (minIndex < 0) or (minIndex > maxIndex):
395 395 raise ValueError("some value in (%d,%d) is not valid" % (
396 396 minIndex, maxIndex))
397 397
398 398 if (maxIndex >= self.dataOut.nHeights):
399 399 maxIndex = self.dataOut.nHeights - 1
400 400
401 401 # seleccion de indices para velocidades
402 402 indminvel = numpy.where(velrange >= minVel)
403 403 indmaxvel = numpy.where(velrange <= maxVel)
404 404 try:
405 405 minIndexVel = indminvel[0][0]
406 406 except:
407 407 minIndexVel = 0
408 408
409 409 try:
410 410 maxIndexVel = indmaxvel[0][-1]
411 411 except:
412 412 maxIndexVel = len(velrange)
413 413
414 414 # seleccion del espectro
415 415 data_spc = self.dataOut.data_spc[:,
416 416 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
417 417 # estimacion de ruido
418 418 noise = numpy.zeros(self.dataOut.nChannels)
419 419
420 420 for channel in range(self.dataOut.nChannels):
421 421 daux = data_spc[channel, :, :]
422 422 sortdata = numpy.sort(daux, axis=None)
423 423 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
424 424
425 425 self.dataOut.noise_estimation = noise.copy()
426 426
427 427 return 1
428 428
429 429 class removeDC(Operation):
430 430
431 431 def run(self, dataOut, mode=2):
432 432 self.dataOut = dataOut
433 433 jspectra = self.dataOut.data_spc
434 434 jcspectra = self.dataOut.data_cspc
435 435
436 436 num_chan = jspectra.shape[0]
437 437 num_hei = jspectra.shape[2]
438 438
439 439 if jcspectra is not None:
440 440 jcspectraExist = True
441 441 num_pairs = jcspectra.shape[0]
442 442 else:
443 443 jcspectraExist = False
444 444
445 445 freq_dc = int(jspectra.shape[1] / 2)
446 446 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
447 447 ind_vel = ind_vel.astype(int)
448 448
449 449 if ind_vel[0] < 0:
450 450 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
451 451
452 452 if mode == 1:
453 453 jspectra[:, freq_dc, :] = (
454 454 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
455 455
456 456 if jcspectraExist:
457 457 jcspectra[:, freq_dc, :] = (
458 458 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
459 459
460 460 if mode == 2:
461 461
462 462 vel = numpy.array([-2, -1, 1, 2])
463 463 xx = numpy.zeros([4, 4])
464 464
465 465 for fil in range(4):
466 466 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
467 467
468 468 xx_inv = numpy.linalg.inv(xx)
469 469 xx_aux = xx_inv[0, :]
470 470
471 471 for ich in range(num_chan):
472 472 yy = jspectra[ich, ind_vel, :]
473 473 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
474 474
475 475 junkid = jspectra[ich, freq_dc, :] <= 0
476 476 cjunkid = sum(junkid)
477 477
478 478 if cjunkid.any():
479 479 jspectra[ich, freq_dc, junkid.nonzero()] = (
480 480 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
481 481
482 482 if jcspectraExist:
483 483 for ip in range(num_pairs):
484 484 yy = jcspectra[ip, ind_vel, :]
485 485 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
486 486
487 487 self.dataOut.data_spc = jspectra
488 488 self.dataOut.data_cspc = jcspectra
489 489
490 490 return self.dataOut
491 491
492 492 class getNoise(Operation):
493
493 494 def __init__(self):
494 495
495 496 Operation.__init__(self)
496 497
497 def run(self, dataOut, minHei=None, maxHei=None, minVel=None, maxVel=None, minFreq= None, maxFreq=None,):
498 self.dataOut = dataOut.copy()
499 print("1: ",dataOut.noise_estimation, dataOut.normFactor)
498 def run(self, dataOut, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None):
499 self.dataOut = dataOut
500 500
501 501 if minHei == None:
502 502 minHei = self.dataOut.heightList[0]
503 503
504 504 if maxHei == None:
505 505 maxHei = self.dataOut.heightList[-1]
506 506
507 507 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
508 508 print('minHei: %.2f is out of the heights range' % (minHei))
509 509 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
510 510 minHei = self.dataOut.heightList[0]
511 511
512 512 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
513 513 print('maxHei: %.2f is out of the heights range' % (maxHei))
514 514 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
515 515 maxHei = self.dataOut.heightList[-1]
516 516
517 517
518 518 #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
519 519 minIndexFFT = 0
520 520 maxIndexFFT = 0
521 521 # validacion de velocidades
522 522 indminPoint = None
523 523 indmaxPoint = None
524 524
525 525 if minVel == None and maxVel == None:
526 526
527 527 freqrange = self.dataOut.getFreqRange(1)
528 528
529 529 if minFreq == None:
530 530 minFreq = freqrange[0]
531 531
532 532 if maxFreq == None:
533 533 maxFreq = freqrange[-1]
534 534
535 535 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
536 536 print('minFreq: %.2f is out of the frequency range' % (minFreq))
537 537 print('minFreq is setting to %.2f' % (freqrange[0]))
538 538 minFreq = freqrange[0]
539 539
540 540 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
541 541 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
542 542 print('maxFreq is setting to %.2f' % (freqrange[-1]))
543 543 maxFreq = freqrange[-1]
544 544
545 545 indminPoint = numpy.where(freqrange >= minFreq)
546 546 indmaxPoint = numpy.where(freqrange <= maxFreq)
547 547
548 548 else:
549 549 velrange = self.dataOut.getVelRange(1)
550 550
551 551 if minVel == None:
552 552 minVel = velrange[0]
553 553
554 554 if maxVel == None:
555 555 maxVel = velrange[-1]
556 556
557 557 if (minVel < velrange[0]) or (minVel > maxVel):
558 558 print('minVel: %.2f is out of the velocity range' % (minVel))
559 559 print('minVel is setting to %.2f' % (velrange[0]))
560 560 minVel = velrange[0]
561 561
562 562 if (maxVel > velrange[-1]) or (maxVel < minVel):
563 563 print('maxVel: %.2f is out of the velocity range' % (maxVel))
564 564 print('maxVel is setting to %.2f' % (velrange[-1]))
565 565 maxVel = velrange[-1]
566 566
567 567 indminPoint = numpy.where(velrange >= minVel)
568 568 indmaxPoint = numpy.where(velrange <= maxVel)
569 569
570 570
571 571 # seleccion de indices para rango
572 572 minIndex = 0
573 573 maxIndex = 0
574 574 heights = self.dataOut.heightList
575 575
576 576 inda = numpy.where(heights >= minHei)
577 577 indb = numpy.where(heights <= maxHei)
578 578
579 579 try:
580 580 minIndex = inda[0][0]
581 581 except:
582 582 minIndex = 0
583 583
584 584 try:
585 585 maxIndex = indb[0][-1]
586 586 except:
587 587 maxIndex = len(heights)
588 588
589 589 if (minIndex < 0) or (minIndex > maxIndex):
590 590 raise ValueError("some value in (%d,%d) is not valid" % (
591 591 minIndex, maxIndex))
592 592
593 593 if (maxIndex >= self.dataOut.nHeights):
594 594 maxIndex = self.dataOut.nHeights - 1
595 595 #############################################################3
596 596 # seleccion de indices para velocidades
597 597
598 598 try:
599 599 minIndexFFT = indminPoint[0][0]
600 600 except:
601 601 minIndexFFT = 0
602 602
603 603 try:
604 604 maxIndexFFT = indmaxPoint[0][-1]
605 605 except:
606 606 maxIndexFFT = len( self.dataOut.getFreqRange(1))
607 607
608 608 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
609 self.dataOut.noise_estimation = None
609 610 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
610 611
611 self.dataOut.noise_estimation = noise.copy()
612 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
612 613 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
613 614 return self.dataOut
614 615
615 616
616 617
617 618 # import matplotlib.pyplot as plt
618 619
619 620 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
620 621 z = (x - a1) / a2
621 622 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
622 623 return y
623 624
624 625
625 626 class CleanRayleigh(Operation):
626 627
627 628 def __init__(self):
628 629
629 630 Operation.__init__(self)
630 631 self.i=0
631 632 self.isConfig = False
632 633 self.__dataReady = False
633 634 self.__profIndex = 0
634 635 self.byTime = False
635 636 self.byProfiles = False
636 637
637 638 self.bloques = None
638 639 self.bloque0 = None
639 640
640 641 self.index = 0
641 642
642 643 self.buffer = 0
643 644 self.buffer2 = 0
644 645 self.buffer3 = 0
645 646
646 647
647 648 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
648 649
649 650 self.nChannels = dataOut.nChannels
650 651 self.nProf = dataOut.nProfiles
651 652 self.nPairs = dataOut.data_cspc.shape[0]
652 653 self.pairsArray = numpy.array(dataOut.pairsList)
653 654 self.spectra = dataOut.data_spc
654 655 self.cspectra = dataOut.data_cspc
655 656 self.heights = dataOut.heightList #alturas totales
656 657 self.nHeights = len(self.heights)
657 658 self.min_hei = min_hei
658 659 self.max_hei = max_hei
659 660 if (self.min_hei == None):
660 661 self.min_hei = 0
661 662 if (self.max_hei == None):
662 663 self.max_hei = dataOut.heightList[-1]
663 664 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
664 665 self.heightsClean = self.heights[self.hval] #alturas filtradas
665 666 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
666 667 self.nHeightsClean = len(self.heightsClean)
667 668 self.channels = dataOut.channelList
668 669 self.nChan = len(self.channels)
669 670 self.nIncohInt = dataOut.nIncohInt
670 671 self.__initime = dataOut.utctime
671 672 self.maxAltInd = self.hval[-1]+1
672 673 self.minAltInd = self.hval[0]
673 674
674 675 self.crosspairs = dataOut.pairsList
675 676 self.nPairs = len(self.crosspairs)
676 677 self.normFactor = dataOut.normFactor
677 678 self.nFFTPoints = dataOut.nFFTPoints
678 679 self.ippSeconds = dataOut.ippSeconds
679 680 self.currentTime = self.__initime
680 681 self.pairsArray = numpy.array(dataOut.pairsList)
681 682 self.factor_stdv = factor_stdv
682 683
683 684 if n != None :
684 685 self.byProfiles = True
685 686 self.nIntProfiles = n
686 687 else:
687 688 self.__integrationtime = timeInterval
688 689
689 690 self.__dataReady = False
690 691 self.isConfig = True
691 692
692 693
693 694
694 695 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
695 696
696 697 if not self.isConfig :
697 698
698 699 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
699 700
700 701 tini=dataOut.utctime
701 702
702 703 if self.byProfiles:
703 704 if self.__profIndex == self.nIntProfiles:
704 705 self.__dataReady = True
705 706 else:
706 707 if (tini - self.__initime) >= self.__integrationtime:
707 708
708 709 self.__dataReady = True
709 710 self.__initime = tini
710 711
711 712 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
712 713
713 714 if self.__dataReady:
714 715
715 716 self.__profIndex = 0
716 717 jspc = self.buffer
717 718 jcspc = self.buffer2
718 719 #jnoise = self.buffer3
719 720 self.buffer = dataOut.data_spc
720 721 self.buffer2 = dataOut.data_cspc
721 722 #self.buffer3 = dataOut.noise
722 723 self.currentTime = dataOut.utctime
723 724 if numpy.any(jspc) :
724 725 #print( jspc.shape, jcspc.shape)
725 726 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
726 727 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
727 728 self.__dataReady = False
728 729 #print( jspc.shape, jcspc.shape)
729 730 dataOut.flagNoData = False
730 731 else:
731 732 dataOut.flagNoData = True
732 733 self.__dataReady = False
733 734 return dataOut
734 735 else:
735 736 #print( len(self.buffer))
736 737 if numpy.any(self.buffer):
737 738 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
738 739 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
739 740 self.buffer3 += dataOut.data_dc
740 741 else:
741 742 self.buffer = dataOut.data_spc
742 743 self.buffer2 = dataOut.data_cspc
743 744 self.buffer3 = dataOut.data_dc
744 745 #print self.index, self.fint
745 746 #print self.buffer2.shape
746 747 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
747 748 self.__profIndex += 1
748 749 return dataOut ## NOTE: REV
749 750
750 751
751 752 #index = tini.tm_hour*12+tini.tm_min/5
752 753 '''REVISAR'''
753 754 # jspc = jspc/self.nFFTPoints/self.normFactor
754 755 # jcspc = jcspc/self.nFFTPoints/self.normFactor
755 756
756 757
757 758
758 759 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
759 760 dataOut.data_spc = tmp_spectra
760 761 dataOut.data_cspc = tmp_cspectra
761 762
762 763 #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
763 764
764 765 dataOut.data_dc = self.buffer3
765 766 dataOut.nIncohInt *= self.nIntProfiles
766 767 dataOut.utctime = self.currentTime #tiempo promediado
767 768 #print("Time: ",time.localtime(dataOut.utctime))
768 769 # dataOut.data_spc = sat_spectra
769 770 # dataOut.data_cspc = sat_cspectra
770 771 self.buffer = 0
771 772 self.buffer2 = 0
772 773 self.buffer3 = 0
773 774
774 775 return dataOut
775 776
776 777 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
777 778 #print("OP cleanRayleigh")
778 779 #import matplotlib.pyplot as plt
779 780 #for k in range(149):
780 781 #channelsProcssd = []
781 782 #channelA_ok = False
782 783 #rfunc = cspectra.copy() #self.bloques
783 784 rfunc = spectra.copy()
784 785 #rfunc = cspectra
785 786 #val_spc = spectra*0.0 #self.bloque0*0.0
786 787 #val_cspc = cspectra*0.0 #self.bloques*0.0
787 788 #in_sat_spectra = spectra.copy() #self.bloque0
788 789 #in_sat_cspectra = cspectra.copy() #self.bloques
789 790
790 791
791 792 ###ONLY FOR TEST:
792 793 raxs = math.ceil(math.sqrt(self.nPairs))
793 794 caxs = math.ceil(self.nPairs/raxs)
794 795 if self.nPairs <4:
795 796 raxs = 2
796 797 caxs = 2
797 798 #print(raxs, caxs)
798 799 fft_rev = 14 #nFFT to plot
799 800 hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
800 801 hei_rev = hei_rev[0]
801 802 #print(hei_rev)
802 803
803 804 #print numpy.absolute(rfunc[:,0,0,14])
804 805
805 806 gauss_fit, covariance = None, None
806 807 for ih in range(self.minAltInd,self.maxAltInd):
807 808 for ifreq in range(self.nFFTPoints):
808 809 '''
809 810 ###ONLY FOR TEST:
810 811 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
811 812 fig, axs = plt.subplots(raxs, caxs)
812 813 fig2, axs2 = plt.subplots(raxs, caxs)
813 814 col_ax = 0
814 815 row_ax = 0
815 816 '''
816 817 #print(self.nPairs)
817 818 for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
818 819 # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
819 820 # continue
820 821 # if not self.crosspairs[ii][0] in channelsProcssd:
821 822 # channelA_ok = True
822 823 #print("pair: ",self.crosspairs[ii])
823 824 '''
824 825 ###ONLY FOR TEST:
825 826 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
826 827 col_ax = 0
827 828 row_ax += 1
828 829 '''
829 830 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
830 831 #print(func2clean.shape)
831 832 val = (numpy.isfinite(func2clean)==True).nonzero()
832 833
833 834 if len(val)>0: #limitador
834 835 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
835 836 if min_val <= -40 :
836 837 min_val = -40
837 838 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
838 839 if max_val >= 200 :
839 840 max_val = 200
840 841 #print min_val, max_val
841 842 step = 1
842 843 #print("Getting bins and the histogram")
843 844 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
844 845 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
845 846 #print(len(y_dist),len(binstep[:-1]))
846 847 #print(row_ax,col_ax, " ..")
847 848 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
848 849 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
849 850 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
850 851 parg = [numpy.amax(y_dist),mean,sigma]
851 852
852 853 newY = None
853 854
854 855 try :
855 856 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
856 857 mode = gauss_fit[1]
857 858 stdv = gauss_fit[2]
858 859 #print(" FIT OK",gauss_fit)
859 860 '''
860 861 ###ONLY FOR TEST:
861 862 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
862 863 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
863 864 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
864 865 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
865 866 axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
866 867 '''
867 868 except:
868 869 mode = mean
869 870 stdv = sigma
870 871 #print("FIT FAIL")
871 872 #continue
872 873
873 874
874 875 #print(mode,stdv)
875 876 #Removing echoes greater than mode + std_factor*stdv
876 877 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
877 878 #noval tiene los indices que se van a remover
878 879 #print("Chan ",ii," novals: ",len(noval[0]))
879 880 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
880 881 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
881 882 #print(novall)
882 883 #print(" ",self.pairsArray[ii])
883 884 #cross_pairs = self.pairsArray[ii]
884 885 #Getting coherent echoes which are removed.
885 886 # if len(novall[0]) > 0:
886 887 #
887 888 # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
888 889 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
889 890 # val_cspc[novall[0],ii,ifreq,ih] = 1
890 891 #print("OUT NOVALL 1")
891 892 try:
892 893 pair = (self.channels[ii],self.channels[ii + 1])
893 894 except:
894 895 pair = (99,99)
895 896 #print("par ", pair)
896 897 if ( pair in self.crosspairs):
897 898 q = self.crosspairs.index(pair)
898 899 #print("estΓ‘ aqui: ", q, (ii,ii + 1))
899 900 new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
900 901 cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
901 902
902 903 #if channelA_ok:
903 904 #chA = self.channels.index(cross_pairs[0])
904 905 new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
905 906 spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
906 907 #channelA_ok = False
907 908
908 909 # chB = self.channels.index(cross_pairs[1])
909 910 # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
910 911 # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
911 912 #
912 913 # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
913 914 # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
914 915 '''
915 916 ###ONLY FOR TEST:
916 917 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
917 918 func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
918 919 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
919 920 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
920 921 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
921 922 axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
922 923 '''
923 924 '''
924 925 ###ONLY FOR TEST:
925 926 col_ax += 1 #contador de ploteo columnas
926 927 ##print(col_ax)
927 928 ###ONLY FOR TEST:
928 929 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
929 930 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
930 931 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
931 932 fig.suptitle(title)
932 933 fig2.suptitle(title2)
933 934 plt.show()
934 935 '''
935 936 ##################################################################################################
936 937
937 938 #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
938 939 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
939 940 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
940 941 for ih in range(self.nHeights):
941 942 for ifreq in range(self.nFFTPoints):
942 943 for ich in range(self.nChan):
943 944 tmp = spectra[:,ich,ifreq,ih]
944 945 valid = (numpy.isfinite(tmp[:])==True).nonzero()
945 946
946 947 if len(valid[0]) >0 :
947 948 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
948 949
949 950 for icr in range(self.nPairs):
950 951 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
951 952 valid = (numpy.isfinite(tmp)==True).nonzero()
952 953 if len(valid[0]) > 0:
953 954 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
954 955
955 956 return out_spectra, out_cspectra
956 957
957 958 def REM_ISOLATED_POINTS(self,array,rth):
958 959 # import matplotlib.pyplot as plt
959 960 if rth == None :
960 961 rth = 4
961 962 #print("REM ISO")
962 963 num_prof = len(array[0,:,0])
963 964 num_hei = len(array[0,0,:])
964 965 n2d = len(array[:,0,0])
965 966
966 967 for ii in range(n2d) :
967 968 #print ii,n2d
968 969 tmp = array[ii,:,:]
969 970 #print tmp.shape, array[ii,101,:],array[ii,102,:]
970 971
971 972 # fig = plt.figure(figsize=(6,5))
972 973 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
973 974 # ax = fig.add_axes([left, bottom, width, height])
974 975 # x = range(num_prof)
975 976 # y = range(num_hei)
976 977 # cp = ax.contour(y,x,tmp)
977 978 # ax.clabel(cp, inline=True,fontsize=10)
978 979 # plt.show()
979 980
980 981 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
981 982 tmp = numpy.reshape(tmp,num_prof*num_hei)
982 983 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
983 984 indxs2 = (tmp > 0).nonzero()
984 985
985 986 indxs1 = (indxs1[0])
986 987 indxs2 = indxs2[0]
987 988 #indxs1 = numpy.array(indxs1[0])
988 989 #indxs2 = numpy.array(indxs2[0])
989 990 indxs = None
990 991 #print indxs1 , indxs2
991 992 for iv in range(len(indxs2)):
992 993 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
993 994 #print len(indxs2), indv
994 995 if len(indv[0]) > 0 :
995 996 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
996 997 # print indxs
997 998 indxs = indxs[1:]
998 999 #print(indxs, len(indxs))
999 1000 if len(indxs) < 4 :
1000 1001 array[ii,:,:] = 0.
1001 1002 return
1002 1003
1003 1004 xpos = numpy.mod(indxs ,num_hei)
1004 1005 ypos = (indxs / num_hei)
1005 1006 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
1006 1007 #print sx
1007 1008 xpos = xpos[sx]
1008 1009 ypos = ypos[sx]
1009 1010
1010 1011 # *********************************** Cleaning isolated points **********************************
1011 1012 ic = 0
1012 1013 while True :
1013 1014 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
1014 1015 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
1015 1016 #plt.plot(r)
1016 1017 #plt.show()
1017 1018 no_coh1 = (numpy.isfinite(r)==True).nonzero()
1018 1019 no_coh2 = (r <= rth).nonzero()
1019 1020 #print r, no_coh1, no_coh2
1020 1021 no_coh1 = numpy.array(no_coh1[0])
1021 1022 no_coh2 = numpy.array(no_coh2[0])
1022 1023 no_coh = None
1023 1024 #print valid1 , valid2
1024 1025 for iv in range(len(no_coh2)):
1025 1026 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
1026 1027 if len(indv[0]) > 0 :
1027 1028 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
1028 1029 no_coh = no_coh[1:]
1029 1030 #print len(no_coh), no_coh
1030 1031 if len(no_coh) < 4 :
1031 1032 #print xpos[ic], ypos[ic], ic
1032 1033 # plt.plot(r)
1033 1034 # plt.show()
1034 1035 xpos[ic] = numpy.nan
1035 1036 ypos[ic] = numpy.nan
1036 1037
1037 1038 ic = ic + 1
1038 1039 if (ic == len(indxs)) :
1039 1040 break
1040 1041 #print( xpos, ypos)
1041 1042
1042 1043 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
1043 1044 #print indxs[0]
1044 1045 if len(indxs[0]) < 4 :
1045 1046 array[ii,:,:] = 0.
1046 1047 return
1047 1048
1048 1049 xpos = xpos[indxs[0]]
1049 1050 ypos = ypos[indxs[0]]
1050 1051 for i in range(0,len(ypos)):
1051 1052 ypos[i]=int(ypos[i])
1052 1053 junk = tmp
1053 1054 tmp = junk*0.0
1054 1055
1055 1056 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
1056 1057 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
1057 1058
1058 1059 #print array.shape
1059 1060 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
1060 1061 #print tmp.shape
1061 1062
1062 1063 # fig = plt.figure(figsize=(6,5))
1063 1064 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1064 1065 # ax = fig.add_axes([left, bottom, width, height])
1065 1066 # x = range(num_prof)
1066 1067 # y = range(num_hei)
1067 1068 # cp = ax.contour(y,x,array[ii,:,:])
1068 1069 # ax.clabel(cp, inline=True,fontsize=10)
1069 1070 # plt.show()
1070 1071 return array
1071 1072
1072 1073
1073 1074 class IntegrationFaradaySpectra(Operation):
1074 1075
1075 1076 __profIndex = 0
1076 1077 __withOverapping = False
1077 1078
1078 1079 __byTime = False
1079 1080 __initime = None
1080 1081 __lastdatatime = None
1081 1082 __integrationtime = None
1082 1083
1083 1084 __buffer_spc = None
1084 1085 __buffer_cspc = None
1085 1086 __buffer_dc = None
1086 1087
1087 1088 __dataReady = False
1088 1089
1089 1090 __timeInterval = None
1090 1091
1091 1092 n = None
1092 1093
1093 1094 def __init__(self):
1094 1095
1095 1096 Operation.__init__(self)
1096 1097
1097 1098 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None):
1098 1099 """
1099 1100 Set the parameters of the integration class.
1100 1101
1101 1102 Inputs:
1102 1103
1103 1104 n : Number of coherent integrations
1104 1105 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1105 1106 overlapping :
1106 1107
1107 1108 """
1108 1109
1109 1110 self.__initime = None
1110 1111 self.__lastdatatime = 0
1111 1112
1112 1113 self.__buffer_spc = []
1113 1114 self.__buffer_cspc = []
1114 1115 self.__buffer_dc = 0
1115 1116
1116 1117 self.__profIndex = 0
1117 1118 self.__dataReady = False
1118 1119 self.__byTime = False
1119 1120
1120 1121 #self.ByLags = dataOut.ByLags ###REDEFINIR
1121 1122 self.ByLags = False
1122 1123
1123 1124 if DPL != None:
1124 1125 self.DPL=DPL
1125 1126 else:
1126 1127 #self.DPL=dataOut.DPL ###REDEFINIR
1127 1128 self.DPL=0
1128 1129
1129 1130 if n is None and timeInterval is None:
1130 1131 raise ValueError("n or timeInterval should be specified ...")
1131 1132
1132 1133 if n is not None:
1133 1134 self.n = int(n)
1134 1135 else:
1135 1136
1136 1137 self.__integrationtime = int(timeInterval)
1137 1138 self.n = None
1138 1139 self.__byTime = True
1139 1140
1140 1141 def putData(self, data_spc, data_cspc, data_dc):
1141 1142 """
1142 1143 Add a profile to the __buffer_spc and increase in one the __profileIndex
1143 1144
1144 1145 """
1145 1146
1146 1147 self.__buffer_spc.append(data_spc)
1147 1148
1148 1149 if data_cspc is None:
1149 1150 self.__buffer_cspc = None
1150 1151 else:
1151 1152 self.__buffer_cspc.append(data_cspc)
1152 1153
1153 1154 if data_dc is None:
1154 1155 self.__buffer_dc = None
1155 1156 else:
1156 1157 self.__buffer_dc += data_dc
1157 1158
1158 1159 self.__profIndex += 1
1159 1160
1160 1161 return
1161 1162
1162 1163 def hildebrand_sekhon_Integration(self,data,navg):
1163 1164
1164 1165 sortdata = numpy.sort(data, axis=None)
1165 1166 sortID=data.argsort()
1166 1167 lenOfData = len(sortdata)
1167 1168 nums_min = lenOfData*0.75
1168 1169 if nums_min <= 5:
1169 1170 nums_min = 5
1170 1171 sump = 0.
1171 1172 sumq = 0.
1172 1173 j = 0
1173 1174 cont = 1
1174 1175 while((cont == 1)and(j < lenOfData)):
1175 1176 sump += sortdata[j]
1176 1177 sumq += sortdata[j]**2
1177 1178 if j > nums_min:
1178 1179 rtest = float(j)/(j-1) + 1.0/navg
1179 1180 if ((sumq*j) > (rtest*sump**2)):
1180 1181 j = j - 1
1181 1182 sump = sump - sortdata[j]
1182 1183 sumq = sumq - sortdata[j]**2
1183 1184 cont = 0
1184 1185 j += 1
1185 1186 #lnoise = sump / j
1186 1187
1187 1188 return j,sortID
1188 1189
1189 1190 def pushData(self):
1190 1191 """
1191 1192 Return the sum of the last profiles and the profiles used in the sum.
1192 1193
1193 1194 Affected:
1194 1195
1195 1196 self.__profileIndex
1196 1197
1197 1198 """
1198 1199 bufferH=None
1199 1200 buffer=None
1200 1201 buffer1=None
1201 1202 buffer_cspc=None
1202 1203 self.__buffer_spc=numpy.array(self.__buffer_spc)
1203 1204 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1205
1204 1206 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1205 1207 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1206 1208 for k in range(7,self.nHeights):
1207 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1209 try:
1210 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1211 except Exception as e:
1212 print(e)
1208 1213 outliers_IDs_cspc=[]
1209 1214 cspc_outliers_exist=False
1210 1215 for i in range(self.nChannels):#dataOut.nChannels):
1211 1216
1212 1217 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1213 1218 indexes=[]
1214 1219 #sortIDs=[]
1215 1220 outliers_IDs=[]
1216 1221
1217 1222 for j in range(self.nProfiles):
1218 1223 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1219 1224 # continue
1220 1225 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1221 1226 # continue
1222 1227 buffer=buffer1[:,j]
1223 1228 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1224 1229
1225 1230 indexes.append(index)
1226 1231 #sortIDs.append(sortID)
1227 1232 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1228 1233
1229 1234 outliers_IDs=numpy.array(outliers_IDs)
1230 1235 outliers_IDs=outliers_IDs.ravel()
1231 1236 outliers_IDs=numpy.unique(outliers_IDs)
1232 1237 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1233 1238 indexes=numpy.array(indexes)
1234 1239 indexmin=numpy.min(indexes)
1235 1240
1236 1241 if indexmin != buffer1.shape[0]:
1237 1242 cspc_outliers_exist=True
1238 1243 ###sortdata=numpy.sort(buffer1,axis=0)
1239 1244 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
1240 1245 lt=outliers_IDs
1241 1246 avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1242 1247
1243 1248 for p in list(outliers_IDs):
1244 1249 buffer1[p,:]=avg
1245 1250
1246 1251 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1247 1252 ###cspc IDs
1248 1253 #indexmin_cspc+=indexmin_cspc
1249 1254 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1250 1255
1251 1256 #if not breakFlag:
1252 1257 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1253 1258 if cspc_outliers_exist:
1254 1259 #sortdata=numpy.sort(buffer_cspc,axis=0)
1255 1260 #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0)
1256 1261 lt=outliers_IDs_cspc
1257 1262
1258 1263 avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1259 1264 for p in list(outliers_IDs_cspc):
1260 1265 buffer_cspc[p,:]=avg
1261 1266
1262 1267 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1263 1268 #else:
1264 1269 #break
1265 1270
1266 1271
1267 1272
1268 1273
1269 1274 buffer=None
1270 1275 bufferH=None
1271 1276 buffer1=None
1272 1277 buffer_cspc=None
1273 1278
1274 1279 #print("cpsc",self.__buffer_cspc[:,0,0,0,0])
1275 1280 #print(self.__profIndex)
1276 1281 #exit()
1277 1282
1278 1283 buffer=None
1279 1284 #print(self.__buffer_spc[:,1,3,20,0])
1280 1285 #print(self.__buffer_spc[:,1,5,37,0])
1281 1286 data_spc = numpy.sum(self.__buffer_spc,axis=0)
1282 1287 data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1283 1288
1284 1289 #print(numpy.shape(data_spc))
1285 1290 #data_spc[1,4,20,0]=numpy.nan
1286 1291
1287 1292 #data_cspc = self.__buffer_cspc
1288 1293 data_dc = self.__buffer_dc
1289 1294 n = self.__profIndex
1290 1295
1291 1296 self.__buffer_spc = []
1292 1297 self.__buffer_cspc = []
1293 1298 self.__buffer_dc = 0
1294 1299 self.__profIndex = 0
1295
1300 #print("pushData Done")
1296 1301 return data_spc, data_cspc, data_dc, n
1297 1302
1298 1303 def byProfiles(self, *args):
1299 1304
1300 1305 self.__dataReady = False
1301 1306 avgdata_spc = None
1302 1307 avgdata_cspc = None
1303 1308 avgdata_dc = None
1304 1309
1305 1310 self.putData(*args)
1306 1311
1307 1312 if self.__profIndex == self.n:
1308 1313
1309 1314 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1310 1315 self.n = n
1311 1316 self.__dataReady = True
1312 1317
1313 1318 return avgdata_spc, avgdata_cspc, avgdata_dc
1314 1319
1315 1320 def byTime(self, datatime, *args):
1316 1321
1317 1322 self.__dataReady = False
1318 1323 avgdata_spc = None
1319 1324 avgdata_cspc = None
1320 1325 avgdata_dc = None
1321 1326
1322 1327 self.putData(*args)
1323 1328
1324 1329 if (datatime - self.__initime) >= self.__integrationtime:
1325 1330 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1326 1331 self.n = n
1327 1332 self.__dataReady = True
1328 1333
1329 1334 return avgdata_spc, avgdata_cspc, avgdata_dc
1330 1335
1331 1336 def integrate(self, datatime, *args):
1332 1337
1333 1338 if self.__profIndex == 0:
1334 1339 self.__initime = datatime
1335 1340
1336 1341 if self.__byTime:
1337 1342 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1338 1343 datatime, *args)
1339 1344 else:
1340 1345 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1341 1346
1342 1347 if not self.__dataReady:
1343 1348 return None, None, None, None
1344 1349
1345 1350 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1346 1351
1347 1352 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False):
1348 1353 if n == 1:
1349 1354 return dataOut
1350 1355
1351 1356 dataOut.flagNoData = True
1352 1357
1353 1358 if not self.isConfig:
1354 1359 self.setup(dataOut, n, timeInterval, overlapping,DPL )
1355 1360 self.isConfig = True
1356 1361
1357 1362 if not self.ByLags:
1358 1363 self.nProfiles=dataOut.nProfiles
1359 1364 self.nChannels=dataOut.nChannels
1360 1365 self.nHeights=dataOut.nHeights
1361 1366 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1362 1367 dataOut.data_spc,
1363 1368 dataOut.data_cspc,
1364 1369 dataOut.data_dc)
1365 1370 else:
1366 1371 self.nProfiles=dataOut.nProfiles
1367 1372 self.nChannels=dataOut.nChannels
1368 1373 self.nHeights=dataOut.nHeights
1369 1374 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1370 1375 dataOut.dataLag_spc,
1371 1376 dataOut.dataLag_cspc,
1372 1377 dataOut.dataLag_dc)
1373 1378
1374 1379 if self.__dataReady:
1375 1380
1376 1381 if not self.ByLags:
1377
1378 dataOut.data_spc = numpy.squeeze(avgdata_spc)
1382 if self.nChannels == 1:
1383 dataOut.data_spc = avgdata_spc
1384 else:
1385 dataOut.data_spc = numpy.squeeze(avgdata_spc)
1379 1386 dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1380 1387 dataOut.data_dc = avgdata_dc
1388
1389
1381 1390 else:
1382 1391 dataOut.dataLag_spc = avgdata_spc
1383 1392 dataOut.dataLag_cspc = avgdata_cspc
1384 1393 dataOut.dataLag_dc = avgdata_dc
1385 1394
1386 1395 dataOut.data_spc=dataOut.dataLag_spc[:,:,:,dataOut.LagPlot]
1387 1396 dataOut.data_cspc=dataOut.dataLag_cspc[:,:,:,dataOut.LagPlot]
1388 1397 dataOut.data_dc=dataOut.dataLag_dc[:,:,dataOut.LagPlot]
1389 1398
1390 1399
1391 1400 dataOut.nIncohInt *= self.n
1392 1401 dataOut.utctime = avgdatatime
1393 1402 dataOut.flagNoData = False
1394 1403
1395 1404 return dataOut
1396 1405
1397 1406 class removeInterference(Operation):
1398 1407
1399 1408 def removeInterference2(self):
1400 1409
1401 1410 cspc = self.dataOut.data_cspc
1402 1411 spc = self.dataOut.data_spc
1403 1412 Heights = numpy.arange(cspc.shape[2])
1404 1413 realCspc = numpy.abs(cspc)
1405 1414
1406 1415 for i in range(cspc.shape[0]):
1407 1416 LinePower= numpy.sum(realCspc[i], axis=0)
1408 1417 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1409 1418 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1410 1419 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1411 1420 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1412 1421 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1413 1422
1414 1423
1415 1424 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1416 1425 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1417 1426 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1418 1427 cspc[i,InterferenceRange,:] = numpy.NaN
1419 1428
1420 1429 self.dataOut.data_cspc = cspc
1421 1430
1422 1431 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1423 1432
1424 1433 jspectra = self.dataOut.data_spc
1425 1434 jcspectra = self.dataOut.data_cspc
1426 1435 jnoise = self.dataOut.getNoise()
1427 1436 num_incoh = self.dataOut.nIncohInt
1428 1437
1429 1438 num_channel = jspectra.shape[0]
1430 1439 num_prof = jspectra.shape[1]
1431 1440 num_hei = jspectra.shape[2]
1432 1441
1433 1442 # hei_interf
1434 1443 if hei_interf is None:
1435 1444 count_hei = int(num_hei / 2)
1436 1445 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1437 1446 hei_interf = numpy.asarray(hei_interf)[0]
1438 1447 # nhei_interf
1439 1448 if (nhei_interf == None):
1440 1449 nhei_interf = 5
1441 1450 if (nhei_interf < 1):
1442 1451 nhei_interf = 1
1443 1452 if (nhei_interf > count_hei):
1444 1453 nhei_interf = count_hei
1445 1454 if (offhei_interf == None):
1446 1455 offhei_interf = 0
1447 1456
1448 1457 ind_hei = list(range(num_hei))
1449 1458 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1450 1459 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1451 1460 mask_prof = numpy.asarray(list(range(num_prof)))
1452 1461 num_mask_prof = mask_prof.size
1453 1462 comp_mask_prof = [0, num_prof / 2]
1454 1463
1455 1464 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1456 1465 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1457 1466 jnoise = numpy.nan
1458 1467 noise_exist = jnoise[0] < numpy.Inf
1459 1468
1460 1469 # Subrutina de Remocion de la Interferencia
1461 1470 for ich in range(num_channel):
1462 1471 # Se ordena los espectros segun su potencia (menor a mayor)
1463 1472 power = jspectra[ich, mask_prof, :]
1464 1473 power = power[:, hei_interf]
1465 1474 power = power.sum(axis=0)
1466 1475 psort = power.ravel().argsort()
1467 1476
1468 1477 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1469 1478 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1470 1479 offhei_interf, nhei_interf + offhei_interf))]]]
1471 1480
1472 1481 if noise_exist:
1473 1482 # tmp_noise = jnoise[ich] / num_prof
1474 1483 tmp_noise = jnoise[ich]
1475 1484 junkspc_interf = junkspc_interf - tmp_noise
1476 1485 #junkspc_interf[:,comp_mask_prof] = 0
1477 1486
1478 1487 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1479 1488 jspc_interf = jspc_interf.transpose()
1480 1489 # Calculando el espectro de interferencia promedio
1481 1490 noiseid = numpy.where(
1482 1491 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1483 1492 noiseid = noiseid[0]
1484 1493 cnoiseid = noiseid.size
1485 1494 interfid = numpy.where(
1486 1495 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1487 1496 interfid = interfid[0]
1488 1497 cinterfid = interfid.size
1489 1498
1490 1499 if (cnoiseid > 0):
1491 1500 jspc_interf[noiseid] = 0
1492 1501
1493 1502 # Expandiendo los perfiles a limpiar
1494 1503 if (cinterfid > 0):
1495 1504 new_interfid = (
1496 1505 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1497 1506 new_interfid = numpy.asarray(new_interfid)
1498 1507 new_interfid = {x for x in new_interfid}
1499 1508 new_interfid = numpy.array(list(new_interfid))
1500 1509 new_cinterfid = new_interfid.size
1501 1510 else:
1502 1511 new_cinterfid = 0
1503 1512
1504 1513 for ip in range(new_cinterfid):
1505 1514 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1506 1515 jspc_interf[new_interfid[ip]
1507 1516 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1508 1517
1509 1518 jspectra[ich, :, ind_hei] = jspectra[ich, :,
1510 1519 ind_hei] - jspc_interf # Corregir indices
1511 1520
1512 1521 # Removiendo la interferencia del punto de mayor interferencia
1513 1522 ListAux = jspc_interf[mask_prof].tolist()
1514 1523 maxid = ListAux.index(max(ListAux))
1515 1524
1516 1525 if cinterfid > 0:
1517 1526 for ip in range(cinterfid * (interf == 2) - 1):
1518 1527 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1519 1528 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1520 1529 cind = len(ind)
1521 1530
1522 1531 if (cind > 0):
1523 1532 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1524 1533 (1 + (numpy.random.uniform(cind) - 0.5) /
1525 1534 numpy.sqrt(num_incoh))
1526 1535
1527 1536 ind = numpy.array([-2, -1, 1, 2])
1528 1537 xx = numpy.zeros([4, 4])
1529 1538
1530 1539 for id1 in range(4):
1531 1540 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1532 1541
1533 1542 xx_inv = numpy.linalg.inv(xx)
1534 1543 xx = xx_inv[:, 0]
1535 1544 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1536 1545 yy = jspectra[ich, mask_prof[ind], :]
1537 1546 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
1538 1547 yy.transpose(), xx)
1539 1548
1540 1549 indAux = (jspectra[ich, :, :] < tmp_noise *
1541 1550 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1542 1551 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1543 1552 (1 - 1 / numpy.sqrt(num_incoh))
1544 1553
1545 1554 # Remocion de Interferencia en el Cross Spectra
1546 1555 if jcspectra is None:
1547 1556 return jspectra, jcspectra
1548 1557 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1549 1558 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1550 1559
1551 1560 for ip in range(num_pairs):
1552 1561
1553 1562 #-------------------------------------------
1554 1563
1555 1564 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1556 1565 cspower = cspower[:, hei_interf]
1557 1566 cspower = cspower.sum(axis=0)
1558 1567
1559 1568 cspsort = cspower.ravel().argsort()
1560 1569 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1561 1570 offhei_interf, nhei_interf + offhei_interf))]]]
1562 1571 junkcspc_interf = junkcspc_interf.transpose()
1563 1572 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1564 1573
1565 1574 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1566 1575
1567 1576 median_real = int(numpy.median(numpy.real(
1568 1577 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1569 1578 median_imag = int(numpy.median(numpy.imag(
1570 1579 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1571 1580 comp_mask_prof = [int(e) for e in comp_mask_prof]
1572 1581 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1573 1582 median_real, median_imag)
1574 1583
1575 1584 for iprof in range(num_prof):
1576 1585 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1577 1586 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1578 1587
1579 1588 # Removiendo la Interferencia
1580 1589 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1581 1590 :, ind_hei] - jcspc_interf
1582 1591
1583 1592 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1584 1593 maxid = ListAux.index(max(ListAux))
1585 1594
1586 1595 ind = numpy.array([-2, -1, 1, 2])
1587 1596 xx = numpy.zeros([4, 4])
1588 1597
1589 1598 for id1 in range(4):
1590 1599 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1591 1600
1592 1601 xx_inv = numpy.linalg.inv(xx)
1593 1602 xx = xx_inv[:, 0]
1594 1603
1595 1604 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1596 1605 yy = jcspectra[ip, mask_prof[ind], :]
1597 1606 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1598 1607
1599 1608 # Guardar Resultados
1600 1609 self.dataOut.data_spc = jspectra
1601 1610 self.dataOut.data_cspc = jcspectra
1602 1611
1603 1612 return 1
1604 1613
1605 1614 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
1606 1615
1607 1616 self.dataOut = dataOut
1608 1617
1609 1618 if mode == 1:
1610 1619 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1611 1620 elif mode == 2:
1612 1621 self.removeInterference2()
1613 1622
1614 1623 return self.dataOut
1615 1624
1616 1625
1617 1626 class IncohInt(Operation):
1618 1627
1619 1628 __profIndex = 0
1620 1629 __withOverapping = False
1621 1630
1622 1631 __byTime = False
1623 1632 __initime = None
1624 1633 __lastdatatime = None
1625 1634 __integrationtime = None
1626 1635
1627 1636 __buffer_spc = None
1628 1637 __buffer_cspc = None
1629 1638 __buffer_dc = None
1630 1639
1631 1640 __dataReady = False
1632 1641
1633 1642 __timeInterval = None
1634 1643
1635 1644 n = None
1636 1645
1637 1646 def __init__(self):
1638 1647
1639 1648 Operation.__init__(self)
1640 1649
1641 1650 def setup(self, n=None, timeInterval=None, overlapping=False):
1642 1651 """
1643 1652 Set the parameters of the integration class.
1644 1653
1645 1654 Inputs:
1646 1655
1647 1656 n : Number of coherent integrations
1648 1657 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1649 1658 overlapping :
1650 1659
1651 1660 """
1652 1661
1653 1662 self.__initime = None
1654 1663 self.__lastdatatime = 0
1655 1664
1656 1665 self.__buffer_spc = 0
1657 1666 self.__buffer_cspc = 0
1658 1667 self.__buffer_dc = 0
1659 1668
1660 1669 self.__profIndex = 0
1661 1670 self.__dataReady = False
1662 1671 self.__byTime = False
1663 1672
1664 1673 if n is None and timeInterval is None:
1665 1674 raise ValueError("n or timeInterval should be specified ...")
1666 1675
1667 1676 if n is not None:
1668 1677 self.n = int(n)
1669 1678 else:
1670 1679
1671 1680 self.__integrationtime = int(timeInterval)
1672 1681 self.n = None
1673 1682 self.__byTime = True
1674 1683
1675 1684 def putData(self, data_spc, data_cspc, data_dc):
1676 1685 """
1677 1686 Add a profile to the __buffer_spc and increase in one the __profileIndex
1678 1687
1679 1688 """
1680 1689
1681 1690 self.__buffer_spc += data_spc
1682 1691
1683 1692 if data_cspc is None:
1684 1693 self.__buffer_cspc = None
1685 1694 else:
1686 1695 self.__buffer_cspc += data_cspc
1687 1696
1688 1697 if data_dc is None:
1689 1698 self.__buffer_dc = None
1690 1699 else:
1691 1700 self.__buffer_dc += data_dc
1692 1701
1693 1702 self.__profIndex += 1
1694 1703
1695 1704 return
1696 1705
1697 1706 def pushData(self):
1698 1707 """
1699 1708 Return the sum of the last profiles and the profiles used in the sum.
1700 1709
1701 1710 Affected:
1702 1711
1703 1712 self.__profileIndex
1704 1713
1705 1714 """
1706 1715
1707 1716 data_spc = self.__buffer_spc
1708 1717 data_cspc = self.__buffer_cspc
1709 1718 data_dc = self.__buffer_dc
1710 1719 n = self.__profIndex
1711 1720
1712 1721 self.__buffer_spc = 0
1713 1722 self.__buffer_cspc = 0
1714 1723 self.__buffer_dc = 0
1715 1724 self.__profIndex = 0
1716 1725
1717 1726 return data_spc, data_cspc, data_dc, n
1718 1727
1719 1728 def byProfiles(self, *args):
1720 1729
1721 1730 self.__dataReady = False
1722 1731 avgdata_spc = None
1723 1732 avgdata_cspc = None
1724 1733 avgdata_dc = None
1725 1734
1726 1735 self.putData(*args)
1727 1736
1728 1737 if self.__profIndex == self.n:
1729 1738
1730 1739 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1731 1740 self.n = n
1732 1741 self.__dataReady = True
1733 1742
1734 1743 return avgdata_spc, avgdata_cspc, avgdata_dc
1735 1744
1736 1745 def byTime(self, datatime, *args):
1737 1746
1738 1747 self.__dataReady = False
1739 1748 avgdata_spc = None
1740 1749 avgdata_cspc = None
1741 1750 avgdata_dc = None
1742 1751
1743 1752 self.putData(*args)
1744 1753
1745 1754 if (datatime - self.__initime) >= self.__integrationtime:
1746 1755 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1747 1756 self.n = n
1748 1757 self.__dataReady = True
1749 1758
1750 1759 return avgdata_spc, avgdata_cspc, avgdata_dc
1751 1760
1752 1761 def integrate(self, datatime, *args):
1753 1762
1754 1763 if self.__profIndex == 0:
1755 1764 self.__initime = datatime
1756 1765
1757 1766 if self.__byTime:
1758 1767 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1759 1768 datatime, *args)
1760 1769 else:
1761 1770 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1762 1771
1763 1772 if not self.__dataReady:
1764 1773 return None, None, None, None
1765 1774
1766 1775 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1767 1776
1768 1777 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1769 1778 if n == 1:
1770 1779 return dataOut
1771 1780
1772 1781 dataOut.flagNoData = True
1773 1782
1774 1783 if not self.isConfig:
1775 1784 self.setup(n, timeInterval, overlapping)
1776 1785 self.isConfig = True
1777 1786
1778 1787 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1779 1788 dataOut.data_spc,
1780 1789 dataOut.data_cspc,
1781 1790 dataOut.data_dc)
1782 1791
1783 1792 if self.__dataReady:
1784 1793
1785 1794 dataOut.data_spc = avgdata_spc
1786 1795 dataOut.data_cspc = avgdata_cspc
1787 1796 dataOut.data_dc = avgdata_dc
1788 1797 dataOut.nIncohInt *= self.n
1789 1798 dataOut.utctime = avgdatatime
1790 1799 dataOut.flagNoData = False
1791 1800
1792 1801 return dataOut
1793 1802
1794 1803 class dopplerFlip(Operation):
1795 1804
1796 1805 def run(self, dataOut):
1797 1806 # arreglo 1: (num_chan, num_profiles, num_heights)
1798 1807 self.dataOut = dataOut
1799 1808 # JULIA-oblicua, indice 2
1800 1809 # arreglo 2: (num_profiles, num_heights)
1801 1810 jspectra = self.dataOut.data_spc[2]
1802 1811 jspectra_tmp = numpy.zeros(jspectra.shape)
1803 1812 num_profiles = jspectra.shape[0]
1804 1813 freq_dc = int(num_profiles / 2)
1805 1814 # Flip con for
1806 1815 for j in range(num_profiles):
1807 1816 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1808 1817 # Intercambio perfil de DC con perfil inmediato anterior
1809 1818 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1810 1819 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1811 1820 # canal modificado es re-escrito en el arreglo de canales
1812 1821 self.dataOut.data_spc[2] = jspectra_tmp
1813 1822
1814 1823 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now