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