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