##// END OF EJS Templates
schain
RSS
fix diferent blocks per file size at the end,-> resize arrays
joabAM -
r1555:73660288fffc
parent child
Show More
Side by Side Unified
Context file:
r1555:73660288fffc Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,495 +1,518
1 1 """
2 2 Utilities for IO modules
3 3 @modified: Joab Apaza
4 4 @email: roj-op01@igp.gob.pe, joab.apaza32@gmail.com
5 5
6 6 """
7 7 ################################################################################
8 8 ################################################################################
9 9 import os
10 10 from datetime import datetime
11 11 import numpy
12 12
13 13 from schainpy.model.data.jrodata import Parameters
14 14 import itertools
15 15 import numpy
16 16 import h5py
17 17 import re
18 18 import time
19 19 ################################################################################
20 20 ################################################################################
21 21 ################################################################################
22 22 def folder_in_range(folder, start_date, end_date, pattern):
23 23 """
24 24 Check whether folder is bettwen start_date and end_date
25 25
26 26 Args:
27 27 folder (str): Folder to check
28 28 start_date (date): Initial date
29 29 end_date (date): Final date
30 30 pattern (str): Datetime format of the folder
31 31 Returns:
32 32 bool: True for success, False otherwise
33 33 """
34 34 try:
35 35 dt = datetime.strptime(folder, pattern)
36 36 except:
37 37 raise ValueError('Folder {} does not match {} format'.format(folder, pattern))
38 38 return start_date <= dt.date() <= end_date
39 39
40 40 ################################################################################
41 41 ################################################################################
42 42 ################################################################################
43 43 def getHei_index( minHei, maxHei, heightList):
44 44 if (minHei < heightList[0]):
45 45 minHei = heightList[0]
46 46
47 47 if (maxHei > heightList[-1]):
48 48 maxHei = heightList[-1]
49 49
50 50 minIndex = 0
51 51 maxIndex = 0
52 52 heights = numpy.asarray(heightList)
53 53
54 54 inda = numpy.where(heights >= minHei)
55 55 indb = numpy.where(heights <= maxHei)
56 56
57 57 try:
58 58 minIndex = inda[0][0]
59 59 except:
60 60 minIndex = 0
61 61
62 62 try:
63 63 maxIndex = indb[0][-1]
64 64 except:
65 65 maxIndex = len(heightList)
66 66 return minIndex,maxIndex
67 67
68 68 ################################################################################
69 69 ################################################################################
70 70 ################################################################################
71 71 class MergeH5(object):
72 72 """Processing unit to read HDF5 format files
73 73
74 74 This unit reads HDF5 files created with `HDFWriter` operation when channels area
75 75 processed by separated. Then merge all channels in a single files.
76 76
77 77 "example"
78 78 nChannels = 4
79 79 pathOut = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/merged"
80 80 p0 = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/d2022240_Ch0"
81 81 p1 = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/d2022240_Ch1"
82 82 p2 = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/d2022240_Ch2"
83 83 p3 = "/home/soporte/Data/OutTest/clean2D/perpAndObliq/byChannels/d2022240_Ch3"
84 84 list = ['data_spc','data_cspc','nIncohInt','utctime']
85 85 merger = MergeH5(nChannels,pathOut,list, p0, p1,p2,p3)
86 86 merger.run()
87 87
88 88 """
89 89
90 90 # #__attrs__ = ['paths', 'nChannels']
91 91 isConfig = False
92 92 inPaths = None
93 93 nChannels = None
94 94 ch_dataIn = []
95 95
96 96 channelList = []
97 97
98 98 def __init__(self,nChannels, pOut, dataList, *args):
99 99
100 100 self.inPaths = [p for p in args]
101 101 #print(self.inPaths)
102 102 if len(self.inPaths) != nChannels:
103 103 print("ERROR, number of channels different from iput paths {} != {}".format(nChannels, len(args)))
104 104 return
105 105
106 106 self.pathOut = pOut
107 107 self.dataList = dataList
108 108 self.nChannels = len(self.inPaths)
109 109 self.ch_dataIn = [Parameters() for p in args]
110 110 self.dataOut = Parameters()
111 111 self.channelList = [n for n in range(nChannels)]
112 112 self.blocksPerFile = None
113 113 self.date = None
114 114 self.ext = ".hdf5$"
115 115 self.dataList = dataList
116 116 self.optchar = "D"
117 117 self.meta = {}
118 118 self.data = {}
119 119 self.open_file = h5py.File
120 120 self.open_mode = 'r'
121 121 self.description = {}
122 122 self.extras = {}
123 123 self.filefmt = "*%Y%j***"
124 124 self.folderfmt = "*%Y%j"
125 125
126 126
127 127 def setup(self):
128 128
129 129
130 130 # if not self.ext.startswith('.'):
131 131 # self.ext = '.{}'.format(self.ext)
132 132
133 133 self.filenameList = self.searchFiles(self.inPaths, None)
134 134 self.nfiles = len(self.filenameList[0])
135 135
136 136
137 137
138 138 def searchFiles(self, paths, date, walk=True):
139 139 # self.paths = path
140 140 #self.date = startDate
141 141 #self.walk = walk
142 142 filenameList = [[] for n in range(self.nChannels)]
143 143 ch = 0
144 144 for path in paths:
145 145 if os.path.exists(path):
146 146 print("Searching files in {}".format(path))
147 147 filenameList[ch] = self.getH5files(path, walk)
148 148 print("Found: ")
149 149 for f in filenameList[ch]:
150 150 print(f)
151 151 else:
152 152 self.status = 0
153 153 print('Path:%s does not exists'%path)
154 154 return 0
155 155 ch+=1
156 156 return filenameList
157 157
158 158 def getH5files(self, path, walk):
159 159
160 160 dirnameList = []
161 161 pat = '(\d)+.'+self.ext
162 162 if walk:
163 163 for root, dirs, files in os.walk(path):
164 164 for dir in dirs:
165 165 #print(os.path.join(root,dir))
166 166 files = [re.search(pat,x) for x in os.listdir(os.path.join(root,dir))]
167 167 #print(files)
168 168 files = [x for x in files if x!=None]
169 169 files = [x.string for x in files]
170 170 files = [os.path.join(root,dir,x) for x in files if x!=None]
171 171 files.sort()
172 172
173 173 dirnameList += files
174 174 return dirnameList
175 175 else:
176 176
177 177 dirnameList = [re.search(pat,x) for x in os.listdir(path)]
178 178 dirnameList = [x for x in dirnameList if x!=None]
179 179 dirnameList = [x.string for x in dirnameList]
180 180 dirnameList = [x for x in dirnameList if x!=None]
181 181 dirnameList.sort()
182 182
183 183 return dirnameList
184 184
185 185
186 186 def readFile(self,fp,ch):
187 187
188 188 '''Read metadata and data'''
189 189
190 190 self.readMetadata(fp)
191 191 data = self.readData(fp)
192 192
193 193
194 194 for attr in self.meta:
195 195 setattr(self.ch_dataIn[ch], attr, self.meta[attr])
196 196
197 197 self.fill_dataIn(data, self.ch_dataIn[ch])
198 198
199
199 200 return
200 201
201 202
202 203 def readMetadata(self, fp):
203 204 '''
204 205 Reads Metadata
205 206 '''
206 207 meta = {}
207 208 self.metadataList = []
208 209 grp = fp['Metadata']
209 210 for name in grp:
210 211 meta[name] = grp[name][()]
211 212 self.metadataList.append(name)
212 213
213 214 for k in meta:
214 215 if ('List' in k):
215 216 meta[k] = meta[k].tolist()
216 217 if not self.meta:
217 218 self.meta = meta
218 219 self.meta["channelList"] =[n for n in range(self.nChannels)]
219 220 return 1
220 221 else:
221 222 if len(self.meta) == len(meta):
222 223 for k in meta:
223 224 if 'List' in k and 'channel' not in k and "height" not in k:
224 225 self.meta[k] += meta[k]
225 226
226 227 return 1
227 228 else:
228 229 return 0
229 230
230 231
231 232
232 233 def fill_dataIn(self,data, dataIn):
233 234
234 235 for attr in data:
235 236 if data[attr].ndim == 1:
236 237 setattr(dataIn, attr, data[attr][:])
237 238 else:
238 setattr(dataIn, attr, data[attr][:,:])
239 setattr(dataIn, attr, numpy.squeeze(data[attr][:,:]))
239 240 #print("shape in", dataIn.data_spc.shape, len(dataIn.data_spc))
240 241 if dataIn.data_spc.ndim > 3:
241 242 dataIn.data_spc = dataIn.data_spc[0]
242 243 #print("shape in", dataIn.data_spc.shape)
243 if self.blocksPerFile == None:
244 self.blocksPerFile = len(dataIn.data_spc) #blocks, ch, fft, hei
245 print("blocks per file: ", self.blocksPerFile)
244
245
246
247 def getBlocksPerFile(self):
248 b = numpy.zeros(self.nChannels)
249 for i in range(self.nChannels):
250 b[i] = self.ch_dataIn[i].data_spc.shape[0] #number of blocks
251
252 self.blocksPerFile = int(b.min())
253 iresh_ch = numpy.where(b > self.blocksPerFile)[0]
254 if len(iresh_ch) > 0:
255 for ich in iresh_ch:
256 for i in range(len(self.dataList)):
257 if hasattr(self.ch_dataIn[ich], self.dataList[i]):
258 # print("reshaping ", self.dataList[i])
259 # print(getattr(self.ch_dataIn[ich], self.dataList[i]).shape)
260 dataAux = getattr(self.ch_dataIn[ich], self.dataList[i])
261 setattr(self.ch_dataIn[ich], self.dataList[i], None)
262 setattr(self.ch_dataIn[ich], self.dataList[i], dataAux[0:self.blocksPerFile])
263 # print(getattr(self.ch_dataIn[ich], self.dataList[i]).shape)
264 else:
265 return
246 266
247 267
248 268 def getLabel(self, name, x=None):
249 269
250 270 if x is None:
251 271 if 'Data' in self.description:
252 272 data = self.description['Data']
253 273 if 'Metadata' in self.description:
254 274 data.update(self.description['Metadata'])
255 275 else:
256 276 data = self.description
257 277 if name in data:
258 278 if isinstance(data[name], str):
259 279 return data[name]
260 280 elif isinstance(data[name], list):
261 281 return None
262 282 elif isinstance(data[name], dict):
263 283 for key, value in data[name].items():
264 284 return key
265 285 return name
266 286 else:
267 287 if 'Metadata' in self.description:
268 288 meta = self.description['Metadata']
269 289 else:
270 290 meta = self.description
271 291 if name in meta:
272 292 if isinstance(meta[name], list):
273 293 return meta[name][x]
274 294 elif isinstance(meta[name], dict):
275 295 for key, value in meta[name].items():
276 296 return value[x]
277 297
278 298 if 'cspc' in name:
279 299 return 'pair{:02d}'.format(x)
280 300 else:
281 301 return 'channel{:02d}'.format(x)
282 302
283 303 def readData(self, fp):
284 304 #print("read fp: ", fp)
285 305 data = {}
286 306
287 307 grp = fp['Data']
288 308
289 309 for name in grp:
290 310 if isinstance(grp[name], h5py.Dataset):
291 311 array = grp[name][()]
292 312 elif isinstance(grp[name], h5py.Group):
293 313 array = []
294 314 for ch in grp[name]:
295 315 array.append(grp[name][ch][()])
296 316 array = numpy.array(array)
297 317 else:
298 318 print('Unknown type: {}'.format(name))
299 319 data[name] = array
300 320
301 321 return data
302 322
303 323 def getDataOut(self):
304 324
305 325 self.dataOut = self.ch_dataIn[0].copy() #dataIn #blocks, fft, hei
306 326 self.dataOut.data_spc = None
307 327 self.dataOut.utctim = None
308 328 self.dataOut.nIncohInt = None
309 329 #print(self.ch_dataIn[0].data_spc.shape)
310 330 spc = [data.data_spc for data in self.ch_dataIn]
311 331 self.dataOut.data_spc = numpy.stack(spc, axis=1) #blocks, ch, fft, hei
312 332 time = [data.utctime for data in self.ch_dataIn]
313 333 time = numpy.asarray(time).mean(axis=0)
314 334 #time = numpy.reshape(time, (len(time),1))
315 335 time = numpy.squeeze(time)
316 336 self.dataOut.utctime = time
317 337 ints = [data.nIncohInt for data in self.ch_dataIn]
318 338 self.dataOut.nIncohInt = numpy.stack(ints, axis=1)
319 339
320 print("nIncohInt 1: ",self.dataOut.nIncohInt.shape)
321 340
322 341 if self.dataOut.nIncohInt.ndim > 3:
323 342 aux = self.dataOut.nIncohInt
324 343 self.dataOut.nIncohInt = None
325 344 self.dataOut.nIncohInt = aux[0]
326 345
327 346
328 347 if self.dataOut.nIncohInt.ndim < 3:
329 348 nIncohInt = numpy.repeat(self.dataOut.nIncohInt, self.dataOut.nHeights).reshape(self.blocksPerFile,self.nChannels, self.dataOut.nHeights)
330 349 #nIncohInt = numpy.reshape(nIncohInt, (self.blocksPerFile,self.nChannels, self.dataOut.nHeights))
331 350 self.dataOut.nIncohInt = None
332 351 self.dataOut.nIncohInt = nIncohInt
333 352
334 353 if (self.dataOut.nIncohInt.shape)[0]==self.nChannels: ## ch,blocks, hei
335 354 self.dataOut.nIncohInt = numpy.swapaxes(self.dataOut.nIncohInt, 0, 1) ## blocks,ch, hei
336 355
337 print("nIncohInt 2: ", self.dataOut.nIncohInt.shape)
338 356 #print("utcTime: ", time.shape)
339 357 #print("data_spc ",self.dataOut.data_spc.shape)
340 358 pairsList = [pair for pair in itertools.combinations(self.channelList, 2)]
341 359 #print("PairsList: ", pairsList)
342 360 self.dataOut.pairsList = pairsList
343 361 cspc = []
344 362
345 363 for i, j in pairsList:
346 364 cspc.append(self.ch_dataIn[i].data_spc*numpy.conjugate(self.ch_dataIn[j].data_spc)) #blocks, fft, hei
347 365
348 366 cspc = numpy.asarray(cspc) # # pairs, blocks, fft, hei
349 367 #print("cspc: ",cspc.shape)
350 368 self.dataOut.data_cspc = numpy.swapaxes(cspc, 0, 1) ## blocks, pairs, fft, hei
351 369 #print("dataOut.data_cspc: ",self.dataOut.data_cspc.shape)
352 370
353 371 def writeMetadata(self, fp):
354 372
355 373
356 374 grp = fp.create_group('Metadata')
357 375
358 376 for i in range(len(self.metadataList)):
359 377 if not hasattr(self.dataOut, self.metadataList[i]):
360 378 print('Metadata: `{}` not found'.format(self.metadataList[i]))
361 379 continue
362 380 value = getattr(self.dataOut, self.metadataList[i])
363 381 if isinstance(value, bool):
364 382 if value is True:
365 383 value = 1
366 384 else:
367 385 value = 0
368 386 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
369 387 return
370 388
371 389 def getDsList(self):
372 390
373 391 dsList =[]
374 392 for i in range(len(self.dataList)):
375 393 dsDict = {}
376 394 if hasattr(self.dataOut, self.dataList[i]):
377 395 dataAux = getattr(self.dataOut, self.dataList[i])
378 396 dsDict['variable'] = self.dataList[i]
379 397 else:
380 398 print('Attribute {} not found in dataOut'.format(self.dataList[i]))
381 399 continue
382 400
383 401 if dataAux is None:
384 402 continue
385 403 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
386 404 dsDict['nDim'] = 0
387 405 else:
388 406
389 407 dsDict['nDim'] = len(dataAux.shape) -1
390 408 dsDict['shape'] = dataAux.shape
391 409
392 410 if len(dsDict['shape'])>=2:
393 411 dsDict['dsNumber'] = dataAux.shape[1]
394 412 else:
395 413 dsDict['dsNumber'] = 1
396 414 dsDict['dtype'] = dataAux.dtype
397 415 # if len(dataAux.shape) == 4:
398 416 # dsDict['nDim'] = len(dataAux.shape) -1
399 417 # dsDict['shape'] = dataAux.shape
400 418 # dsDict['dsNumber'] = dataAux.shape[1]
401 419 # dsDict['dtype'] = dataAux.dtype
402 420 # else:
403 421 # dsDict['nDim'] = len(dataAux.shape)
404 422 # dsDict['shape'] = dataAux.shape
405 423 # dsDict['dsNumber'] = dataAux.shape[0]
406 424 # dsDict['dtype'] = dataAux.dtype
407 425
408 426 dsList.append(dsDict)
409 427 #print(dsList)
410 428 self.dsList = dsList
411 429
412 430 def clean_dataIn(self):
413 431 for ch in range(self.nChannels):
414 432 self.ch_dataIn[ch].data_spc = None
415 433 self.ch_dataIn[ch].utctime = None
416 434 self.ch_dataIn[ch].nIncohInt = None
417 435 self.meta ={}
418 436 self.blocksPerFile = None
419 437
420 438 def writeData(self, outFilename):
421 439
422 440 self.getDsList()
423 441
424 442 fp = h5py.File(outFilename, 'w')
425 443 self.writeMetadata(fp)
426 444 grp = fp.create_group('Data')
427 445
428 446 dtsets = []
429 447 data = []
430 448 for dsInfo in self.dsList:
431 449 if dsInfo['nDim'] == 0:
432 450 ds = grp.create_dataset(
433 451 self.getLabel(dsInfo['variable']),(self.blocksPerFile, ),chunks=True,dtype=numpy.float64)
434 452 dtsets.append(ds)
435 453 data.append((dsInfo['variable'], -1))
436 454 else:
437 455 label = self.getLabel(dsInfo['variable'])
438 456 if label is not None:
439 457 sgrp = grp.create_group(label)
440 458 else:
441 459 sgrp = grp
442 460 k = -1*(dsInfo['nDim'] - 1)
443 461 #print(k, dsInfo['shape'], dsInfo['shape'][k:])
444 462 for i in range(dsInfo['dsNumber']):
445 463 ds = sgrp.create_dataset(
446 464 self.getLabel(dsInfo['variable'], i),(self.blocksPerFile, ) + dsInfo['shape'][k:],
447 465 chunks=True,
448 466 dtype=dsInfo['dtype'])
449 467 dtsets.append(ds)
450 468 data.append((dsInfo['variable'], i))
451 469
452 470 #print("\n",dtsets)
453 471
454 472 print('Creating merged file: {}'.format(fp.filename))
455 473
456 474 for i, ds in enumerate(dtsets):
457 475 attr, ch = data[i]
458 476 if ch == -1:
459 477 ds[:] = getattr(self.dataOut, attr)
460 478 else:
461 479 #print(ds, getattr(self.dataOut, attr)[ch].shape)
462 480 aux = getattr(self.dataOut, attr)# block, ch, ...
463 481 aux = numpy.swapaxes(aux,0,1) # ch, blocks, ...
464 482 #print(ds.shape, aux.shape)
465 483 #ds[:] = getattr(self.dataOut, attr)[ch]
466 484 ds[:] = aux[ch]
467 485
468 486 fp.flush()
469 487 fp.close()
470 488 self.clean_dataIn()
471 489 return
472 490
473 491
474 492
475 493 def run(self):
476 494
477 495 if not(self.isConfig):
478 496 self.setup()
479 497 self.isConfig = True
480 498
481 499 for nf in range(self.nfiles):
482 500 name = None
483 501 for ch in range(self.nChannels):
484 502 name = self.filenameList[ch][nf]
485 503 filename = os.path.join(self.inPaths[ch], name)
486 504 fp = h5py.File(filename, 'r')
487 505 #print("Opening file: ",filename)
488 506 self.readFile(fp,ch)
489 507 fp.close()
508
509 if self.blocksPerFile == None:
510 self.getBlocksPerFile()
511 print("blocks per file: ", self.blocksPerFile)
512
490 513 self.getDataOut()
491 514 name = name[-16:]
492 515 #print("Final name out: ", name)
493 516 outFile = os.path.join(self.pathOut, name)
494 517 #print("Outfile: ", outFile)
495 518 self.writeData(outFile)
Show file before | Show file after | Hide comments
@@ -1,2160 +1,2164
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15 import math
16 16
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 18 from schainpy.model.data.jrodata import Spectra
19 19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 20 from schainpy.model.data import _noise
21 21
22 22 from schainpy.utils import log
23 23 import matplotlib.pyplot as plt
24 24 #from scipy.optimize import curve_fit
25 25 from schainpy.model.io.utilsIO import getHei_index
26 26
27 27 class SpectraProc(ProcessingUnit):
28 28
29 29 def __init__(self):
30 30
31 31 ProcessingUnit.__init__(self)
32 32
33 33 self.buffer = None
34 34 self.firstdatatime = None
35 35 self.profIndex = 0
36 36 self.dataOut = Spectra()
37 37 self.id_min = None
38 38 self.id_max = None
39 39 self.setupReq = False #Agregar a todas las unidades de proc
40 40
41 41 def __updateSpecFromVoltage(self):
42 42
43 43
44 44
45 45 self.dataOut.timeZone = self.dataIn.timeZone
46 46 self.dataOut.dstFlag = self.dataIn.dstFlag
47 47 self.dataOut.errorCount = self.dataIn.errorCount
48 48 self.dataOut.useLocalTime = self.dataIn.useLocalTime
49 49 try:
50 50 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
51 51 except:
52 52 pass
53 53 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
54 54 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
55 55 self.dataOut.channelList = self.dataIn.channelList
56 56 self.dataOut.heightList = self.dataIn.heightList
57 57 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
58 58 self.dataOut.nProfiles = self.dataOut.nFFTPoints
59 59 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
60 60 self.dataOut.utctime = self.firstdatatime
61 61 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
62 62 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
63 63 self.dataOut.flagShiftFFT = False
64 64 self.dataOut.nCohInt = self.dataIn.nCohInt
65 65 self.dataOut.nIncohInt = 1
66 66 self.dataOut.radar_ipp = self.dataIn.radar_ipp
67 67 self.dataOut.sampled_heightsFFT = self.dataIn.sampled_heightsFFT
68 68 self.dataOut.pulseLength_TxA = self.dataIn.pulseLength_TxA
69 69 self.dataOut.deltaHeight = self.dataIn.deltaHeight
70 70 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
71 71 self.dataOut.frequency = self.dataIn.frequency
72 72 self.dataOut.realtime = self.dataIn.realtime
73 73 self.dataOut.azimuth = self.dataIn.azimuth
74 74 self.dataOut.zenith = self.dataIn.zenith
75 75 self.dataOut.codeList = self.dataIn.codeList
76 76 self.dataOut.azimuthList = self.dataIn.azimuthList
77 77 self.dataOut.elevationList = self.dataIn.elevationList
78 78
79 79
80 80 def __getFft(self):
81 81 # print("fft donw")
82 82 """
83 83 Convierte valores de Voltaje a Spectra
84 84
85 85 Affected:
86 86 self.dataOut.data_spc
87 87 self.dataOut.data_cspc
88 88 self.dataOut.data_dc
89 89 self.dataOut.heightList
90 90 self.profIndex
91 91 self.buffer
92 92 self.dataOut.flagNoData
93 93 """
94 94 fft_volt = numpy.fft.fft(
95 95 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
96 96 fft_volt = fft_volt.astype(numpy.dtype('complex'))
97 97 dc = fft_volt[:, 0, :]
98 98
99 99 # calculo de self-spectra
100 100 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
101 101 spc = fft_volt * numpy.conjugate(fft_volt)
102 102 spc = spc.real
103 103
104 104 blocksize = 0
105 105 blocksize += dc.size
106 106 blocksize += spc.size
107 107
108 108 cspc = None
109 109 pairIndex = 0
110 110 if self.dataOut.pairsList != None:
111 111 # calculo de cross-spectra
112 112 cspc = numpy.zeros(
113 113 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
114 114 for pair in self.dataOut.pairsList:
115 115 if pair[0] not in self.dataOut.channelList:
116 116 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
117 117 str(pair), str(self.dataOut.channelList)))
118 118 if pair[1] not in self.dataOut.channelList:
119 119 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
120 120 str(pair), str(self.dataOut.channelList)))
121 121
122 122 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
123 123 numpy.conjugate(fft_volt[pair[1], :, :])
124 124 pairIndex += 1
125 125 blocksize += cspc.size
126 126
127 127 self.dataOut.data_spc = spc
128 128 self.dataOut.data_cspc = cspc
129 129 self.dataOut.data_dc = dc
130 130 self.dataOut.blockSize = blocksize
131 131 self.dataOut.flagShiftFFT = False
132 132
133 133 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, zeroPad=False):
134 134 #print("run spc proc")
135 135
136 136 try:
137 137 type = self.dataIn.type.decode("utf-8")
138 138 self.dataIn.type = type
139 139 except:
140 140 pass
141 141 if self.dataIn.type == "Spectra":
142 142
143 143 try:
144 144 self.dataOut.copy(self.dataIn)
145 145 self.dataOut.nProfiles = self.dataOut.nFFTPoints
146 146 #self.dataOut.nHeights = len(self.dataOut.heightList)
147 147 except Exception as e:
148 148 print("Error dataIn ",e)
149 149
150 150 if shift_fft:
151 151 #desplaza a la derecha en el eje 2 determinadas posiciones
152 152 shift = int(self.dataOut.nFFTPoints/2)
153 153 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
154 154
155 155 if self.dataOut.data_cspc is not None:
156 156 #desplaza a la derecha en el eje 2 determinadas posiciones
157 157 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
158 158 if pairsList:
159 159 self.__selectPairs(pairsList)
160 160
161 161
162 162 elif self.dataIn.type == "Voltage":
163 163
164 164 self.dataOut.flagNoData = True
165 165
166 166 if nFFTPoints == None:
167 167 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
168 168
169 169 if nProfiles == None:
170 170 nProfiles = nFFTPoints
171 171
172 172 if ippFactor == None:
173 173 self.dataOut.ippFactor = 1
174 174
175 175 self.dataOut.nFFTPoints = nFFTPoints
176 176 #print(" volts ch,prof, h: ", self.dataIn.data.shape)
177 177 if self.buffer is None:
178 178 if not zeroPad:
179 179 self.buffer = numpy.zeros((self.dataIn.nChannels,
180 180 nProfiles,
181 181 self.dataIn.nHeights),
182 182 dtype='complex')
183 183 else:
184 184 self.buffer = numpy.zeros((self.dataIn.nChannels,
185 185 nFFTPoints,
186 186 self.dataIn.nHeights),
187 187 dtype='complex')
188 188
189 189 if self.dataIn.flagDataAsBlock:
190 190 nVoltProfiles = self.dataIn.data.shape[1]
191 191
192 192 if nVoltProfiles == nProfiles or zeroPad:
193 193 self.buffer = self.dataIn.data.copy()
194 194 self.profIndex = nVoltProfiles
195 195
196 196 elif nVoltProfiles < nProfiles:
197 197
198 198 if self.profIndex == 0:
199 199 self.id_min = 0
200 200 self.id_max = nVoltProfiles
201 201
202 202 self.buffer[:, self.id_min:self.id_max,
203 203 :] = self.dataIn.data
204 204 self.profIndex += nVoltProfiles
205 205 self.id_min += nVoltProfiles
206 206 self.id_max += nVoltProfiles
207 207 else:
208 208 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
209 209 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
210 210 self.dataOut.flagNoData = True
211 211 else:
212 212 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
213 213 self.profIndex += 1
214 214
215 215 if self.firstdatatime == None:
216 216 self.firstdatatime = self.dataIn.utctime
217 217
218 218 if self.profIndex == nProfiles or zeroPad:
219 219
220 220 self.__updateSpecFromVoltage()
221 221
222 222 if pairsList == None:
223 223 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
224 224 else:
225 225 self.dataOut.pairsList = pairsList
226 226 self.__getFft()
227 227 self.dataOut.flagNoData = False
228 228 self.firstdatatime = None
229 229 self.profIndex = 0
230 230
231 231 elif self.dataIn.type == "Parameters":
232 232
233 233 self.dataOut.data_spc = self.dataIn.data_spc
234 234 self.dataOut.data_cspc = self.dataIn.data_cspc
235 235 self.dataOut.data_outlier = self.dataIn.data_outlier
236 236 self.dataOut.nProfiles = self.dataIn.nProfiles
237 237 self.dataOut.nIncohInt = self.dataIn.nIncohInt
238 238 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
239 239 self.dataOut.ippFactor = self.dataIn.ippFactor
240 240 self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
241 241 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
242 242 self.dataOut.ipp = self.dataIn.ipp
243 243 #self.dataOut.abscissaList = self.dataIn.getVelRange(1)
244 244 #self.dataOut.spc_noise = self.dataIn.getNoise()
245 245 #self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
246 246 # self.dataOut.normFactor = self.dataIn.normFactor
247 247 if hasattr(self.dataIn, 'channelList'):
248 248 self.dataOut.channelList = self.dataIn.channelList
249 249 if hasattr(self.dataIn, 'pairsList'):
250 250 self.dataOut.pairsList = self.dataIn.pairsList
251 251 self.dataOut.groupList = self.dataIn.pairsList
252 252
253 253 self.dataOut.flagNoData = False
254 254
255 255 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
256 256 self.dataOut.ChanDist = self.dataIn.ChanDist
257 257 else: self.dataOut.ChanDist = None
258 258
259 259 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
260 260 # self.dataOut.VelRange = self.dataIn.VelRange
261 261 #else: self.dataOut.VelRange = None
262 262
263 263
264 264
265 265 else:
266 266 raise ValueError("The type of input object {} is not valid".format(
267 267 self.dataIn.type))
268 268 #print("spc proc Done", self.dataOut.data_spc.shape)
269 269
270 270 def __selectPairs(self, pairsList):
271 271
272 272 if not pairsList:
273 273 return
274 274
275 275 pairs = []
276 276 pairsIndex = []
277 277
278 278 for pair in pairsList:
279 279 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
280 280 continue
281 281 pairs.append(pair)
282 282 pairsIndex.append(pairs.index(pair))
283 283
284 284 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
285 285 self.dataOut.pairsList = pairs
286 286
287 287 return
288 288
289 289 def selectFFTs(self, minFFT, maxFFT ):
290 290 """
291 291 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
292 292 minFFT<= FFT <= maxFFT
293 293 """
294 294
295 295 if (minFFT > maxFFT):
296 296 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
297 297
298 298 if (minFFT < self.dataOut.getFreqRange()[0]):
299 299 minFFT = self.dataOut.getFreqRange()[0]
300 300
301 301 if (maxFFT > self.dataOut.getFreqRange()[-1]):
302 302 maxFFT = self.dataOut.getFreqRange()[-1]
303 303
304 304 minIndex = 0
305 305 maxIndex = 0
306 306 FFTs = self.dataOut.getFreqRange()
307 307
308 308 inda = numpy.where(FFTs >= minFFT)
309 309 indb = numpy.where(FFTs <= maxFFT)
310 310
311 311 try:
312 312 minIndex = inda[0][0]
313 313 except:
314 314 minIndex = 0
315 315
316 316 try:
317 317 maxIndex = indb[0][-1]
318 318 except:
319 319 maxIndex = len(FFTs)
320 320
321 321 self.selectFFTsByIndex(minIndex, maxIndex)
322 322
323 323 return 1
324 324
325 325 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
326 326 newheis = numpy.where(
327 327 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
328 328
329 329 if hei_ref != None:
330 330 newheis = numpy.where(self.dataOut.heightList > hei_ref)
331 331
332 332 minIndex = min(newheis[0])
333 333 maxIndex = max(newheis[0])
334 334 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
335 335 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
336 336
337 337 # determina indices
338 338 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
339 339 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
340 340 avg_dB = 10 * \
341 341 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
342 342 beacon_dB = numpy.sort(avg_dB)[-nheis:]
343 343 beacon_heiIndexList = []
344 344 for val in avg_dB.tolist():
345 345 if val >= beacon_dB[0]:
346 346 beacon_heiIndexList.append(avg_dB.tolist().index(val))
347 347
348 348 #data_spc = data_spc[:,:,beacon_heiIndexList]
349 349 data_cspc = None
350 350 if self.dataOut.data_cspc is not None:
351 351 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
352 352 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
353 353
354 354 data_dc = None
355 355 if self.dataOut.data_dc is not None:
356 356 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
357 357 #data_dc = data_dc[:,beacon_heiIndexList]
358 358
359 359 self.dataOut.data_spc = data_spc
360 360 self.dataOut.data_cspc = data_cspc
361 361 self.dataOut.data_dc = data_dc
362 362 self.dataOut.heightList = heightList
363 363 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
364 364
365 365 return 1
366 366
367 367 def selectFFTsByIndex(self, minIndex, maxIndex):
368 368 """
369 369
370 370 """
371 371
372 372 if (minIndex < 0) or (minIndex > maxIndex):
373 373 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
374 374
375 375 if (maxIndex >= self.dataOut.nProfiles):
376 376 maxIndex = self.dataOut.nProfiles-1
377 377
378 378 #Spectra
379 379 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
380 380
381 381 data_cspc = None
382 382 if self.dataOut.data_cspc is not None:
383 383 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
384 384
385 385 data_dc = None
386 386 if self.dataOut.data_dc is not None:
387 387 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
388 388
389 389 self.dataOut.data_spc = data_spc
390 390 self.dataOut.data_cspc = data_cspc
391 391 self.dataOut.data_dc = data_dc
392 392
393 393 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
394 394 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
395 395 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
396 396
397 397 return 1
398 398
399 399 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
400 400 # validacion de rango
401 401 if minHei == None:
402 402 minHei = self.dataOut.heightList[0]
403 403
404 404 if maxHei == None:
405 405 maxHei = self.dataOut.heightList[-1]
406 406
407 407 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
408 408 print('minHei: %.2f is out of the heights range' % (minHei))
409 409 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
410 410 minHei = self.dataOut.heightList[0]
411 411
412 412 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
413 413 print('maxHei: %.2f is out of the heights range' % (maxHei))
414 414 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
415 415 maxHei = self.dataOut.heightList[-1]
416 416
417 417 # validacion de velocidades
418 418 velrange = self.dataOut.getVelRange(1)
419 419
420 420 if minVel == None:
421 421 minVel = velrange[0]
422 422
423 423 if maxVel == None:
424 424 maxVel = velrange[-1]
425 425
426 426 if (minVel < velrange[0]) or (minVel > maxVel):
427 427 print('minVel: %.2f is out of the velocity range' % (minVel))
428 428 print('minVel is setting to %.2f' % (velrange[0]))
429 429 minVel = velrange[0]
430 430
431 431 if (maxVel > velrange[-1]) or (maxVel < minVel):
432 432 print('maxVel: %.2f is out of the velocity range' % (maxVel))
433 433 print('maxVel is setting to %.2f' % (velrange[-1]))
434 434 maxVel = velrange[-1]
435 435
436 436 # seleccion de indices para rango
437 437 minIndex = 0
438 438 maxIndex = 0
439 439 heights = self.dataOut.heightList
440 440
441 441 inda = numpy.where(heights >= minHei)
442 442 indb = numpy.where(heights <= maxHei)
443 443
444 444 try:
445 445 minIndex = inda[0][0]
446 446 except:
447 447 minIndex = 0
448 448
449 449 try:
450 450 maxIndex = indb[0][-1]
451 451 except:
452 452 maxIndex = len(heights)
453 453
454 454 if (minIndex < 0) or (minIndex > maxIndex):
455 455 raise ValueError("some value in (%d,%d) is not valid" % (
456 456 minIndex, maxIndex))
457 457
458 458 if (maxIndex >= self.dataOut.nHeights):
459 459 maxIndex = self.dataOut.nHeights - 1
460 460
461 461 # seleccion de indices para velocidades
462 462 indminvel = numpy.where(velrange >= minVel)
463 463 indmaxvel = numpy.where(velrange <= maxVel)
464 464 try:
465 465 minIndexVel = indminvel[0][0]
466 466 except:
467 467 minIndexVel = 0
468 468
469 469 try:
470 470 maxIndexVel = indmaxvel[0][-1]
471 471 except:
472 472 maxIndexVel = len(velrange)
473 473
474 474 # seleccion del espectro
475 475 data_spc = self.dataOut.data_spc[:,
476 476 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
477 477 # estimacion de ruido
478 478 noise = numpy.zeros(self.dataOut.nChannels)
479 479
480 480 for channel in range(self.dataOut.nChannels):
481 481 daux = data_spc[channel, :, :]
482 482 sortdata = numpy.sort(daux, axis=None)
483 483 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
484 484
485 485 self.dataOut.noise_estimation = noise.copy()
486 486
487 487 return 1
488 488
489 489 class removeDC(Operation):
490 490
491 491 def run(self, dataOut, mode=2):
492 492 self.dataOut = dataOut
493 493 jspectra = self.dataOut.data_spc
494 494 jcspectra = self.dataOut.data_cspc
495 495
496 496 num_chan = jspectra.shape[0]
497 497 num_hei = jspectra.shape[2]
498 498
499 499 if jcspectra is not None:
500 500 jcspectraExist = True
501 501 num_pairs = jcspectra.shape[0]
502 502 else:
503 503 jcspectraExist = False
504 504
505 505 freq_dc = int(jspectra.shape[1] / 2)
506 506 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
507 507 ind_vel = ind_vel.astype(int)
508 508
509 509 if ind_vel[0] < 0:
510 510 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
511 511
512 512 if mode == 1:
513 513 jspectra[:, freq_dc, :] = (
514 514 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
515 515
516 516 if jcspectraExist:
517 517 jcspectra[:, freq_dc, :] = (
518 518 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
519 519
520 520 if mode == 2:
521 521
522 522 vel = numpy.array([-2, -1, 1, 2])
523 523 xx = numpy.zeros([4, 4])
524 524
525 525 for fil in range(4):
526 526 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
527 527
528 528 xx_inv = numpy.linalg.inv(xx)
529 529 xx_aux = xx_inv[0, :]
530 530
531 531 for ich in range(num_chan):
532 532 yy = jspectra[ich, ind_vel, :]
533 533 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
534 534
535 535 junkid = jspectra[ich, freq_dc, :] <= 0
536 536 cjunkid = sum(junkid)
537 537
538 538 if cjunkid.any():
539 539 jspectra[ich, freq_dc, junkid.nonzero()] = (
540 540 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
541 541
542 542 if jcspectraExist:
543 543 for ip in range(num_pairs):
544 544 yy = jcspectra[ip, ind_vel, :]
545 545 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
546 546
547 547 self.dataOut.data_spc = jspectra
548 548 self.dataOut.data_cspc = jcspectra
549 549
550 550 return self.dataOut
551 551
552 552 class getNoiseB(Operation):
553 553
554 554 __slots__ =('offset','warnings', 'isConfig', 'minIndex','maxIndex','minIndexFFT','maxIndexFFT')
555 555 def __init__(self):
556 556
557 557 Operation.__init__(self)
558 558 self.isConfig = False
559 559
560 560 def setup(self, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
561 561
562 562 self.warnings = warnings
563 563 if minHei == None:
564 564 minHei = self.dataOut.heightList[0]
565 565
566 566 if maxHei == None:
567 567 maxHei = self.dataOut.heightList[-1]
568 568
569 569 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
570 570 if self.warnings:
571 571 print('minHei: %.2f is out of the heights range' % (minHei))
572 572 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
573 573 minHei = self.dataOut.heightList[0]
574 574
575 575 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
576 576 if self.warnings:
577 577 print('maxHei: %.2f is out of the heights range' % (maxHei))
578 578 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
579 579 maxHei = self.dataOut.heightList[-1]
580 580
581 581
582 582 #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
583 583 minIndexFFT = 0
584 584 maxIndexFFT = 0
585 585 # validacion de velocidades
586 586 indminPoint = None
587 587 indmaxPoint = None
588 588 if self.dataOut.type == 'Spectra':
589 589 if minVel == None and maxVel == None :
590 590
591 591 freqrange = self.dataOut.getFreqRange(1)
592 592
593 593 if minFreq == None:
594 594 minFreq = freqrange[0]
595 595
596 596 if maxFreq == None:
597 597 maxFreq = freqrange[-1]
598 598
599 599 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
600 600 if self.warnings:
601 601 print('minFreq: %.2f is out of the frequency range' % (minFreq))
602 602 print('minFreq is setting to %.2f' % (freqrange[0]))
603 603 minFreq = freqrange[0]
604 604
605 605 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
606 606 if self.warnings:
607 607 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
608 608 print('maxFreq is setting to %.2f' % (freqrange[-1]))
609 609 maxFreq = freqrange[-1]
610 610
611 611 indminPoint = numpy.where(freqrange >= minFreq)
612 612 indmaxPoint = numpy.where(freqrange <= maxFreq)
613 613
614 614 else:
615 615
616 616 velrange = self.dataOut.getVelRange(1)
617 617
618 618 if minVel == None:
619 619 minVel = velrange[0]
620 620
621 621 if maxVel == None:
622 622 maxVel = velrange[-1]
623 623
624 624 if (minVel < velrange[0]) or (minVel > maxVel):
625 625 if self.warnings:
626 626 print('minVel: %.2f is out of the velocity range' % (minVel))
627 627 print('minVel is setting to %.2f' % (velrange[0]))
628 628 minVel = velrange[0]
629 629
630 630 if (maxVel > velrange[-1]) or (maxVel < minVel):
631 631 if self.warnings:
632 632 print('maxVel: %.2f is out of the velocity range' % (maxVel))
633 633 print('maxVel is setting to %.2f' % (velrange[-1]))
634 634 maxVel = velrange[-1]
635 635
636 636 indminPoint = numpy.where(velrange >= minVel)
637 637 indmaxPoint = numpy.where(velrange <= maxVel)
638 638
639 639
640 640 # seleccion de indices para rango
641 641 minIndex = 0
642 642 maxIndex = 0
643 643 heights = self.dataOut.heightList
644 644
645 645 inda = numpy.where(heights >= minHei)
646 646 indb = numpy.where(heights <= maxHei)
647 647
648 648 try:
649 649 minIndex = inda[0][0]
650 650 except:
651 651 minIndex = 0
652 652
653 653 try:
654 654 maxIndex = indb[0][-1]
655 655 except:
656 656 maxIndex = len(heights)
657 657
658 658 if (minIndex < 0) or (minIndex > maxIndex):
659 659 raise ValueError("some value in (%d,%d) is not valid" % (
660 660 minIndex, maxIndex))
661 661
662 662 if (maxIndex >= self.dataOut.nHeights):
663 663 maxIndex = self.dataOut.nHeights - 1
664 664 #############################################################3
665 665 # seleccion de indices para velocidades
666 666 if self.dataOut.type == 'Spectra':
667 667 try:
668 668 minIndexFFT = indminPoint[0][0]
669 669 except:
670 670 minIndexFFT = 0
671 671
672 672 try:
673 673 maxIndexFFT = indmaxPoint[0][-1]
674 674 except:
675 675 maxIndexFFT = len( self.dataOut.getFreqRange(1))
676 676
677 677 self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT
678 678 self.isConfig = True
679 679 self.offset = 1
680 680 if offset!=None:
681 681 self.offset = 10**(offset/10)
682 682 #print("config getNoiseB Done")
683 683
684 684 def run(self, dataOut, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
685 685 self.dataOut = dataOut
686 686
687 687 if not self.isConfig:
688 688 self.setup(offset, minHei, maxHei,minVel, maxVel, minFreq, maxFreq, warnings)
689 689
690 690 self.dataOut.noise_estimation = None
691 691 noise = None
692 692 #print("data type: ",self.dataOut.type, self.dataOut.nIncohInt, self.dataOut.max_nIncohInt)
693 693 if self.dataOut.type == 'Voltage':
694 694 noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
695 695 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
696 696 elif self.dataOut.type == 'Spectra':
697 697 #print(self.dataOut.nChannels, self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.max_nIncohInt, self.dataOut.nIncohInt)
698 698 noise = numpy.zeros( self.dataOut.nChannels)
699 699 norm = 1
700 700
701 701 for channel in range( self.dataOut.nChannels):
702 702 if not hasattr(self.dataOut.nIncohInt,'__len__'):
703 703 norm = 1
704 704 else:
705 705 norm = self.dataOut.max_nIncohInt[channel]/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
706 706 #print("norm nIncoh: ", norm ,self.dataOut.data_spc.shape, self.dataOut.max_nIncohInt)
707 707 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
708 708 daux = numpy.multiply(daux, norm)
709 709 #print("offset: ", self.offset, 10*numpy.log10(self.offset))
710 710 # noise[channel] = self.getNoiseByMean(daux)/self.offset
711 711 #print(daux.shape, daux)
712 712 #noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
713 713 sortdata = numpy.sort(daux, axis=None)
714 714
715 715 noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt[channel])/self.offset
716 716
717 717
718 718 #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
719 719 else:
720 720 noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
721 721 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
722 722 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
723 723 #print("2: ",self.dataOut.noise_estimation)
724 724 #print(self.dataOut.flagNoData)
725 725 #print("getNoise Done", noise, self.dataOut.nProfiles ,self.dataOut.ippFactor)
726 726 return self.dataOut
727 727
728 728 def getNoiseByMean(self,data):
729 729 #data debe estar ordenado
730 730 data = numpy.mean(data,axis=1)
731 731 sortdata = numpy.sort(data, axis=None)
732 732 #sortID=data.argsort()
733 733 #print(data.shape)
734 734
735 735 pnoise = None
736 736 j = 0
737 737
738 738 mean = numpy.mean(sortdata)
739 739 min = numpy.min(sortdata)
740 740 delta = mean - min
741 741 indexes = numpy.where(sortdata > (mean+delta))[0] #only array of indexes
742 742 #print(len(indexes))
743 743 if len(indexes)==0:
744 744 pnoise = numpy.mean(sortdata)
745 745 else:
746 746 j = indexes[0]
747 747 pnoise = numpy.mean(sortdata[0:j])
748 748
749 749 # from matplotlib import pyplot as plt
750 750 # plt.plot(sortdata)
751 751 # plt.vlines(j,(pnoise-delta),(pnoise+delta), color='r')
752 752 # plt.show()
753 753 #print("noise: ", 10*numpy.log10(pnoise))
754 754 return pnoise
755 755
756 756 def getNoiseByHS(self,data, navg):
757 757 #data debe estar ordenado
758 758 #data = numpy.mean(data,axis=1)
759 759 sortdata = numpy.sort(data, axis=None)
760 760
761 761 lenOfData = len(sortdata)
762 762 nums_min = lenOfData*0.2
763 763
764 764 if nums_min <= 5:
765 765
766 766 nums_min = 5
767 767
768 768 sump = 0.
769 769 sumq = 0.
770 770
771 771 j = 0
772 772 cont = 1
773 773
774 774 while((cont == 1)and(j < lenOfData)):
775 775
776 776 sump += sortdata[j]
777 777 sumq += sortdata[j]**2
778 778 #sumq -= sump**2
779 779 if j > nums_min:
780 780 rtest = float(j)/(j-1) + 1.0/navg
781 781 #if ((sumq*j) > (sump**2)):
782 782 if ((sumq*j) > (rtest*sump**2)):
783 783 j = j - 1
784 784 sump = sump - sortdata[j]
785 785 sumq = sumq - sortdata[j]**2
786 786 cont = 0
787 787
788 788 j += 1
789 789
790 790 lnoise = sump / j
791 791
792 792 return lnoise
793 793
794 794
795 795
796 796 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
797 797 z = (x - a1) / a2
798 798 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
799 799 return y
800 800
801 801
802 802 # class CleanRayleigh(Operation):
803 803 #
804 804 # def __init__(self):
805 805 #
806 806 # Operation.__init__(self)
807 807 # self.i=0
808 808 # self.isConfig = False
809 809 # self.__dataReady = False
810 810 # self.__profIndex = 0
811 811 # self.byTime = False
812 812 # self.byProfiles = False
813 813 #
814 814 # self.bloques = None
815 815 # self.bloque0 = None
816 816 #
817 817 # self.index = 0
818 818 #
819 819 # self.buffer = 0
820 820 # self.buffer2 = 0
821 821 # self.buffer3 = 0
822 822 #
823 823 #
824 824 # def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
825 825 #
826 826 # self.nChannels = dataOut.nChannels
827 827 # self.nProf = dataOut.nProfiles
828 828 # self.nPairs = dataOut.data_cspc.shape[0]
829 829 # self.pairsArray = numpy.array(dataOut.pairsList)
830 830 # self.spectra = dataOut.data_spc
831 831 # self.cspectra = dataOut.data_cspc
832 832 # self.heights = dataOut.heightList #alturas totales
833 833 # self.nHeights = len(self.heights)
834 834 # self.min_hei = min_hei
835 835 # self.max_hei = max_hei
836 836 # if (self.min_hei == None):
837 837 # self.min_hei = 0
838 838 # if (self.max_hei == None):
839 839 # self.max_hei = dataOut.heightList[-1]
840 840 # self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
841 841 # self.heightsClean = self.heights[self.hval] #alturas filtradas
842 842 # self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
843 843 # self.nHeightsClean = len(self.heightsClean)
844 844 # self.channels = dataOut.channelList
845 845 # self.nChan = len(self.channels)
846 846 # self.nIncohInt = dataOut.nIncohInt
847 847 # self.__initime = dataOut.utctime
848 848 # self.maxAltInd = self.hval[-1]+1
849 849 # self.minAltInd = self.hval[0]
850 850 #
851 851 # self.crosspairs = dataOut.pairsList
852 852 # self.nPairs = len(self.crosspairs)
853 853 # self.normFactor = dataOut.normFactor
854 854 # self.nFFTPoints = dataOut.nFFTPoints
855 855 # self.ippSeconds = dataOut.ippSeconds
856 856 # self.currentTime = self.__initime
857 857 # self.pairsArray = numpy.array(dataOut.pairsList)
858 858 # self.factor_stdv = factor_stdv
859 859 #
860 860 # if n != None :
861 861 # self.byProfiles = True
862 862 # self.nIntProfiles = n
863 863 # else:
864 864 # self.__integrationtime = timeInterval
865 865 #
866 866 # self.__dataReady = False
867 867 # self.isConfig = True
868 868 #
869 869 #
870 870 #
871 871 # def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
872 872 # #print("runing cleanRayleigh")
873 873 # if not self.isConfig :
874 874 #
875 875 # self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
876 876 #
877 877 # tini=dataOut.utctime
878 878 #
879 879 # if self.byProfiles:
880 880 # if self.__profIndex == self.nIntProfiles:
881 881 # self.__dataReady = True
882 882 # else:
883 883 # if (tini - self.__initime) >= self.__integrationtime:
884 884 #
885 885 # self.__dataReady = True
886 886 # self.__initime = tini
887 887 #
888 888 # #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
889 889 #
890 890 # if self.__dataReady:
891 891 #
892 892 # self.__profIndex = 0
893 893 # jspc = self.buffer
894 894 # jcspc = self.buffer2
895 895 # #jnoise = self.buffer3
896 896 # self.buffer = dataOut.data_spc
897 897 # self.buffer2 = dataOut.data_cspc
898 898 # #self.buffer3 = dataOut.noise
899 899 # self.currentTime = dataOut.utctime
900 900 # if numpy.any(jspc) :
901 901 # #print( jspc.shape, jcspc.shape)
902 902 # jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
903 903 # try:
904 904 # jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
905 905 # except:
906 906 # print("no cspc")
907 907 # self.__dataReady = False
908 908 # #print( jspc.shape, jcspc.shape)
909 909 # dataOut.flagNoData = False
910 910 # else:
911 911 # dataOut.flagNoData = True
912 912 # self.__dataReady = False
913 913 # return dataOut
914 914 # else:
915 915 # #print( len(self.buffer))
916 916 # if numpy.any(self.buffer):
917 917 # self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
918 918 # try:
919 919 # self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
920 920 # self.buffer3 += dataOut.data_dc
921 921 # except:
922 922 # pass
923 923 # else:
924 924 # self.buffer = dataOut.data_spc
925 925 # self.buffer2 = dataOut.data_cspc
926 926 # self.buffer3 = dataOut.data_dc
927 927 # #print self.index, self.fint
928 928 # #print self.buffer2.shape
929 929 # dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
930 930 # self.__profIndex += 1
931 931 # return dataOut ## NOTE: REV
932 932 #
933 933 #
934 934 # #index = tini.tm_hour*12+tini.tm_min/5
935 935 # '''
936 936 # #REVISAR
937 937 # '''
938 938 # # jspc = jspc/self.nFFTPoints/self.normFactor
939 939 # # jcspc = jcspc/self.nFFTPoints/self.normFactor
940 940 #
941 941 #
942 942 #
943 943 # tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
944 944 # dataOut.data_spc = tmp_spectra
945 945 # dataOut.data_cspc = tmp_cspectra
946 946 #
947 947 # #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
948 948 #
949 949 # dataOut.data_dc = self.buffer3
950 950 # dataOut.nIncohInt *= self.nIntProfiles
951 951 # dataOut.max_nIncohInt = self.nIntProfiles
952 952 # dataOut.utctime = self.currentTime #tiempo promediado
953 953 # #print("Time: ",time.localtime(dataOut.utctime))
954 954 # # dataOut.data_spc = sat_spectra
955 955 # # dataOut.data_cspc = sat_cspectra
956 956 # self.buffer = 0
957 957 # self.buffer2 = 0
958 958 # self.buffer3 = 0
959 959 #
960 960 # return dataOut
961 961 #
962 962 # def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
963 963 # print("OP cleanRayleigh")
964 964 # #import matplotlib.pyplot as plt
965 965 # #for k in range(149):
966 966 # #channelsProcssd = []
967 967 # #channelA_ok = False
968 968 # #rfunc = cspectra.copy() #self.bloques
969 969 # rfunc = spectra.copy()
970 970 # #rfunc = cspectra
971 971 # #val_spc = spectra*0.0 #self.bloque0*0.0
972 972 # #val_cspc = cspectra*0.0 #self.bloques*0.0
973 973 # #in_sat_spectra = spectra.copy() #self.bloque0
974 974 # #in_sat_cspectra = cspectra.copy() #self.bloques
975 975 #
976 976 #
977 977 # ###ONLY FOR TEST:
978 978 # raxs = math.ceil(math.sqrt(self.nPairs))
979 979 # if raxs == 0:
980 980 # raxs = 1
981 981 # caxs = math.ceil(self.nPairs/raxs)
982 982 # if self.nPairs <4:
983 983 # raxs = 2
984 984 # caxs = 2
985 985 # #print(raxs, caxs)
986 986 # fft_rev = 14 #nFFT to plot
987 987 # hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
988 988 # hei_rev = hei_rev[0]
989 989 # #print(hei_rev)
990 990 #
991 991 # #print numpy.absolute(rfunc[:,0,0,14])
992 992 #
993 993 # gauss_fit, covariance = None, None
994 994 # for ih in range(self.minAltInd,self.maxAltInd):
995 995 # for ifreq in range(self.nFFTPoints):
996 996 # '''
997 997 # ###ONLY FOR TEST:
998 998 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
999 999 # fig, axs = plt.subplots(raxs, caxs)
1000 1000 # fig2, axs2 = plt.subplots(raxs, caxs)
1001 1001 # col_ax = 0
1002 1002 # row_ax = 0
1003 1003 # '''
1004 1004 # #print(self.nPairs)
1005 1005 # for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
1006 1006 # # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
1007 1007 # # continue
1008 1008 # # if not self.crosspairs[ii][0] in channelsProcssd:
1009 1009 # # channelA_ok = True
1010 1010 # #print("pair: ",self.crosspairs[ii])
1011 1011 # '''
1012 1012 # ###ONLY FOR TEST:
1013 1013 # if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
1014 1014 # col_ax = 0
1015 1015 # row_ax += 1
1016 1016 # '''
1017 1017 # func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
1018 1018 # #print(func2clean.shape)
1019 1019 # val = (numpy.isfinite(func2clean)==True).nonzero()
1020 1020 #
1021 1021 # if len(val)>0: #limitador
1022 1022 # min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
1023 1023 # if min_val <= -40 :
1024 1024 # min_val = -40
1025 1025 # max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
1026 1026 # if max_val >= 200 :
1027 1027 # max_val = 200
1028 1028 # #print min_val, max_val
1029 1029 # step = 1
1030 1030 # #print("Getting bins and the histogram")
1031 1031 # x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
1032 1032 # y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1033 1033 # #print(len(y_dist),len(binstep[:-1]))
1034 1034 # #print(row_ax,col_ax, " ..")
1035 1035 # #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
1036 1036 # mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
1037 1037 # sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
1038 1038 # parg = [numpy.amax(y_dist),mean,sigma]
1039 1039 #
1040 1040 # newY = None
1041 1041 #
1042 1042 # try :
1043 1043 # gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
1044 1044 # mode = gauss_fit[1]
1045 1045 # stdv = gauss_fit[2]
1046 1046 # #print(" FIT OK",gauss_fit)
1047 1047 # '''
1048 1048 # ###ONLY FOR TEST:
1049 1049 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1050 1050 # newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
1051 1051 # axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1052 1052 # axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1053 1053 # axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1054 1054 # '''
1055 1055 # except:
1056 1056 # mode = mean
1057 1057 # stdv = sigma
1058 1058 # #print("FIT FAIL")
1059 1059 # #continue
1060 1060 #
1061 1061 #
1062 1062 # #print(mode,stdv)
1063 1063 # #Removing echoes greater than mode + std_factor*stdv
1064 1064 # noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
1065 1065 # #noval tiene los indices que se van a remover
1066 1066 # #print("Chan ",ii," novals: ",len(noval[0]))
1067 1067 # if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
1068 1068 # novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
1069 1069 # #print(novall)
1070 1070 # #print(" ",self.pairsArray[ii])
1071 1071 # #cross_pairs = self.pairsArray[ii]
1072 1072 # #Getting coherent echoes which are removed.
1073 1073 # # if len(novall[0]) > 0:
1074 1074 # #
1075 1075 # # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
1076 1076 # # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
1077 1077 # # val_cspc[novall[0],ii,ifreq,ih] = 1
1078 1078 # #print("OUT NOVALL 1")
1079 1079 # try:
1080 1080 # pair = (self.channels[ii],self.channels[ii + 1])
1081 1081 # except:
1082 1082 # pair = (99,99)
1083 1083 # #print("par ", pair)
1084 1084 # if ( pair in self.crosspairs):
1085 1085 # q = self.crosspairs.index(pair)
1086 1086 # #print("estΓ‘ aqui: ", q, (ii,ii + 1))
1087 1087 # new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
1088 1088 # cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
1089 1089 #
1090 1090 # #if channelA_ok:
1091 1091 # #chA = self.channels.index(cross_pairs[0])
1092 1092 # new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
1093 1093 # spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
1094 1094 # #channelA_ok = False
1095 1095 #
1096 1096 # # chB = self.channels.index(cross_pairs[1])
1097 1097 # # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
1098 1098 # # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
1099 1099 # #
1100 1100 # # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
1101 1101 # # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
1102 1102 # '''
1103 1103 # ###ONLY FOR TEST:
1104 1104 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1105 1105 # func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
1106 1106 # y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1107 1107 # axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1108 1108 # axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1109 1109 # axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1110 1110 # '''
1111 1111 # '''
1112 1112 # ###ONLY FOR TEST:
1113 1113 # col_ax += 1 #contador de ploteo columnas
1114 1114 # ##print(col_ax)
1115 1115 # ###ONLY FOR TEST:
1116 1116 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1117 1117 # title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
1118 1118 # title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
1119 1119 # fig.suptitle(title)
1120 1120 # fig2.suptitle(title2)
1121 1121 # plt.show()
1122 1122 # '''
1123 1123 # ##################################################################################################
1124 1124 #
1125 1125 # #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
1126 1126 # out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
1127 1127 # out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
1128 1128 # for ih in range(self.nHeights):
1129 1129 # for ifreq in range(self.nFFTPoints):
1130 1130 # for ich in range(self.nChan):
1131 1131 # tmp = spectra[:,ich,ifreq,ih]
1132 1132 # valid = (numpy.isfinite(tmp[:])==True).nonzero()
1133 1133 #
1134 1134 # if len(valid[0]) >0 :
1135 1135 # out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1136 1136 #
1137 1137 # for icr in range(self.nPairs):
1138 1138 # tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
1139 1139 # valid = (numpy.isfinite(tmp)==True).nonzero()
1140 1140 # if len(valid[0]) > 0:
1141 1141 # out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1142 1142 #
1143 1143 # return out_spectra, out_cspectra
1144 1144 #
1145 1145 # def REM_ISOLATED_POINTS(self,array,rth):
1146 1146 # # import matplotlib.pyplot as plt
1147 1147 # if rth == None :
1148 1148 # rth = 4
1149 1149 # #print("REM ISO")
1150 1150 # num_prof = len(array[0,:,0])
1151 1151 # num_hei = len(array[0,0,:])
1152 1152 # n2d = len(array[:,0,0])
1153 1153 #
1154 1154 # for ii in range(n2d) :
1155 1155 # #print ii,n2d
1156 1156 # tmp = array[ii,:,:]
1157 1157 # #print tmp.shape, array[ii,101,:],array[ii,102,:]
1158 1158 #
1159 1159 # # fig = plt.figure(figsize=(6,5))
1160 1160 # # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1161 1161 # # ax = fig.add_axes([left, bottom, width, height])
1162 1162 # # x = range(num_prof)
1163 1163 # # y = range(num_hei)
1164 1164 # # cp = ax.contour(y,x,tmp)
1165 1165 # # ax.clabel(cp, inline=True,fontsize=10)
1166 1166 # # plt.show()
1167 1167 #
1168 1168 # #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
1169 1169 # tmp = numpy.reshape(tmp,num_prof*num_hei)
1170 1170 # indxs1 = (numpy.isfinite(tmp)==True).nonzero()
1171 1171 # indxs2 = (tmp > 0).nonzero()
1172 1172 #
1173 1173 # indxs1 = (indxs1[0])
1174 1174 # indxs2 = indxs2[0]
1175 1175 # #indxs1 = numpy.array(indxs1[0])
1176 1176 # #indxs2 = numpy.array(indxs2[0])
1177 1177 # indxs = None
1178 1178 # #print indxs1 , indxs2
1179 1179 # for iv in range(len(indxs2)):
1180 1180 # indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
1181 1181 # #print len(indxs2), indv
1182 1182 # if len(indv[0]) > 0 :
1183 1183 # indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
1184 1184 # # print indxs
1185 1185 # indxs = indxs[1:]
1186 1186 # #print(indxs, len(indxs))
1187 1187 # if len(indxs) < 4 :
1188 1188 # array[ii,:,:] = 0.
1189 1189 # return
1190 1190 #
1191 1191 # xpos = numpy.mod(indxs ,num_hei)
1192 1192 # ypos = (indxs / num_hei)
1193 1193 # sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
1194 1194 # #print sx
1195 1195 # xpos = xpos[sx]
1196 1196 # ypos = ypos[sx]
1197 1197 #
1198 1198 # # *********************************** Cleaning isolated points **********************************
1199 1199 # ic = 0
1200 1200 # while True :
1201 1201 # r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
1202 1202 # #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
1203 1203 # #plt.plot(r)
1204 1204 # #plt.show()
1205 1205 # no_coh1 = (numpy.isfinite(r)==True).nonzero()
1206 1206 # no_coh2 = (r <= rth).nonzero()
1207 1207 # #print r, no_coh1, no_coh2
1208 1208 # no_coh1 = numpy.array(no_coh1[0])
1209 1209 # no_coh2 = numpy.array(no_coh2[0])
1210 1210 # no_coh = None
1211 1211 # #print valid1 , valid2
1212 1212 # for iv in range(len(no_coh2)):
1213 1213 # indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
1214 1214 # if len(indv[0]) > 0 :
1215 1215 # no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
1216 1216 # no_coh = no_coh[1:]
1217 1217 # #print len(no_coh), no_coh
1218 1218 # if len(no_coh) < 4 :
1219 1219 # #print xpos[ic], ypos[ic], ic
1220 1220 # # plt.plot(r)
1221 1221 # # plt.show()
1222 1222 # xpos[ic] = numpy.nan
1223 1223 # ypos[ic] = numpy.nan
1224 1224 #
1225 1225 # ic = ic + 1
1226 1226 # if (ic == len(indxs)) :
1227 1227 # break
1228 1228 # #print( xpos, ypos)
1229 1229 #
1230 1230 # indxs = (numpy.isfinite(list(xpos))==True).nonzero()
1231 1231 # #print indxs[0]
1232 1232 # if len(indxs[0]) < 4 :
1233 1233 # array[ii,:,:] = 0.
1234 1234 # return
1235 1235 #
1236 1236 # xpos = xpos[indxs[0]]
1237 1237 # ypos = ypos[indxs[0]]
1238 1238 # for i in range(0,len(ypos)):
1239 1239 # ypos[i]=int(ypos[i])
1240 1240 # junk = tmp
1241 1241 # tmp = junk*0.0
1242 1242 #
1243 1243 # tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
1244 1244 # array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
1245 1245 #
1246 1246 # #print array.shape
1247 1247 # #tmp = numpy.reshape(tmp,(num_prof,num_hei))
1248 1248 # #print tmp.shape
1249 1249 #
1250 1250 # # fig = plt.figure(figsize=(6,5))
1251 1251 # # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1252 1252 # # ax = fig.add_axes([left, bottom, width, height])
1253 1253 # # x = range(num_prof)
1254 1254 # # y = range(num_hei)
1255 1255 # # cp = ax.contour(y,x,array[ii,:,:])
1256 1256 # # ax.clabel(cp, inline=True,fontsize=10)
1257 1257 # # plt.show()
1258 1258 # return array
1259 1259 #
1260 1260
1261 1261 class IntegrationFaradaySpectra(Operation):
1262 1262
1263 1263 __profIndex = 0
1264 1264 __withOverapping = False
1265 1265
1266 1266 __byTime = False
1267 1267 __initime = None
1268 1268 __lastdatatime = None
1269 1269 __integrationtime = None
1270 1270
1271 1271 __buffer_spc = None
1272 1272 __buffer_cspc = None
1273 1273 __buffer_dc = None
1274 1274
1275 1275 __dataReady = False
1276 1276
1277 1277 __timeInterval = None
1278 1278 n_ints = None #matriz de numero de integracions (CH,HEI)
1279 1279 n = None
1280 1280 minHei_ind = None
1281 1281 maxHei_ind = None
1282 1282 navg = 1.0
1283 1283 factor = 0.0
1284 1284 dataoutliers = None # (CHANNELS, HEIGHTS)
1285 1285
1286 1286 def __init__(self):
1287 1287
1288 1288 Operation.__init__(self)
1289 1289
1290 1290 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, avg=1,factor=0.75):
1291 1291 """
1292 1292 Set the parameters of the integration class.
1293 1293
1294 1294 Inputs:
1295 1295
1296 1296 n : Number of coherent integrations
1297 1297 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1298 1298 overlapping :
1299 1299
1300 1300 """
1301 1301
1302 1302 self.__initime = None
1303 1303 self.__lastdatatime = 0
1304 1304
1305 1305 self.__buffer_spc = []
1306 1306 self.__buffer_cspc = []
1307 1307 self.__buffer_dc = 0
1308 1308
1309 1309 self.__profIndex = 0
1310 1310 self.__dataReady = False
1311 1311 self.__byTime = False
1312 1312
1313 1313 self.factor = factor
1314 1314 self.navg = avg
1315 1315 #self.ByLags = dataOut.ByLags ###REDEFINIR
1316 1316 self.ByLags = False
1317 1317 self.maxProfilesInt = 0
1318 1318 self.__nChannels = dataOut.nChannels
1319 1319 if DPL != None:
1320 1320 self.DPL=DPL
1321 1321 else:
1322 1322 #self.DPL=dataOut.DPL ###REDEFINIR
1323 1323 self.DPL=0
1324 1324
1325 1325 if n is None and timeInterval is None:
1326 1326 raise ValueError("n or timeInterval should be specified ...")
1327 1327
1328 1328 if n is not None:
1329 1329 self.n = int(n)
1330 1330 else:
1331 1331 self.__integrationtime = int(timeInterval)
1332 1332 self.n = None
1333 1333 self.__byTime = True
1334 1334
1335 1335 if minHei == None:
1336 1336 minHei = self.dataOut.heightList[0]
1337 1337
1338 1338 if maxHei == None:
1339 1339 maxHei = self.dataOut.heightList[-1]
1340 1340
1341 1341 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1342 1342 print('minHei: %.2f is out of the heights range' % (minHei))
1343 1343 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
1344 1344 minHei = self.dataOut.heightList[0]
1345 1345
1346 1346 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
1347 1347 print('maxHei: %.2f is out of the heights range' % (maxHei))
1348 1348 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
1349 1349 maxHei = self.dataOut.heightList[-1]
1350 1350
1351 1351 ind_list1 = numpy.where(self.dataOut.heightList >= minHei)
1352 1352 ind_list2 = numpy.where(self.dataOut.heightList <= maxHei)
1353 1353 self.minHei_ind = ind_list1[0][0]
1354 1354 self.maxHei_ind = ind_list2[0][-1]
1355 1355 #print("setup rem sats done")
1356 1356
1357 1357 def putData(self, data_spc, data_cspc, data_dc):
1358 1358 """
1359 1359 Add a profile to the __buffer_spc and increase in one the __profileIndex
1360 1360
1361 1361 """
1362 1362
1363 1363 self.__buffer_spc.append(data_spc)
1364 1364
1365 1365 if self.__nChannels < 2:
1366 1366 self.__buffer_cspc = None
1367 1367 else:
1368 1368 self.__buffer_cspc.append(data_cspc)
1369 1369
1370 1370 if data_dc is None:
1371 1371 self.__buffer_dc = None
1372 1372 else:
1373 1373 self.__buffer_dc += data_dc
1374 1374
1375 1375 self.__profIndex += 1
1376 1376
1377 1377 return
1378 1378
1379 1379 def hildebrand_sekhon_Integration(self,sortdata,navg, factor):
1380 1380 #data debe estar ordenado
1381 1381 #sortdata = numpy.sort(data, axis=None)
1382 1382 #sortID=data.argsort()
1383 1383 lenOfData = len(sortdata)
1384 1384 nums_min = lenOfData*factor
1385 1385 if nums_min <= 5:
1386 1386 nums_min = 5
1387 1387 sump = 0.
1388 1388 sumq = 0.
1389 1389 j = 0
1390 1390 cont = 1
1391 1391 while((cont == 1)and(j < lenOfData)):
1392 1392 sump += sortdata[j]
1393 1393 sumq += sortdata[j]**2
1394 1394 if j > nums_min:
1395 1395 rtest = float(j)/(j-1) + 1.0/navg
1396 1396 if ((sumq*j) > (rtest*sump**2)):
1397 1397 j = j - 1
1398 1398 sump = sump - sortdata[j]
1399 1399 sumq = sumq - sortdata[j]**2
1400 1400 cont = 0
1401 1401 j += 1
1402 1402 #lnoise = sump / j
1403 1403 #print("H S done")
1404 1404 #return j,sortID
1405 1405 return j
1406 1406
1407 1407
1408 1408 def pushData(self):
1409 1409 """
1410 1410 Return the sum of the last profiles and the profiles used in the sum.
1411 1411
1412 1412 Affected:
1413 1413
1414 1414 self.__profileIndex
1415 1415
1416 1416 """
1417 1417 bufferH=None
1418 1418 buffer=None
1419 1419 buffer1=None
1420 1420 buffer_cspc=None
1421 1421 #print("aes: ", self.__buffer_cspc)
1422 1422 self.__buffer_spc=numpy.array(self.__buffer_spc)
1423 1423 if self.__nChannels > 1 :
1424 1424 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1425 1425
1426 1426 #print("FREQ_DC",self.__buffer_spc.shape,self.__buffer_cspc.shape)
1427 1427
1428 1428 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1429 1429 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1430 1430
1431 1431 self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
1432 1432
1433 1433 for k in range(self.minHei_ind,self.maxHei_ind):
1434 1434 if self.__nChannels > 1:
1435 1435 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1436 1436
1437 1437 outliers_IDs_cspc=[]
1438 1438 cspc_outliers_exist=False
1439 1439 for i in range(self.nChannels):#dataOut.nChannels):
1440 1440
1441 1441 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1442 1442 indexes=[]
1443 1443 #sortIDs=[]
1444 1444 outliers_IDs=[]
1445 1445
1446 1446 for j in range(self.nProfiles): #frecuencias en el tiempo
1447 1447 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1448 1448 # continue
1449 1449 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1450 1450 # continue
1451 1451 buffer=buffer1[:,j]
1452 1452 sortdata = numpy.sort(buffer, axis=None)
1453 1453
1454 1454 sortID=buffer.argsort()
1455 1455 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1456 1456
1457 1457 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1458 1458
1459 1459 # fig,ax = plt.subplots()
1460 1460 # ax.set_title(str(k)+" "+str(j))
1461 1461 # x=range(len(sortdata))
1462 1462 # ax.scatter(x,sortdata)
1463 1463 # ax.axvline(index)
1464 1464 # plt.show()
1465 1465
1466 1466 indexes.append(index)
1467 1467 #sortIDs.append(sortID)
1468 1468 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1469 1469
1470 1470 #print("Outliers: ",outliers_IDs)
1471 1471 outliers_IDs=numpy.array(outliers_IDs)
1472 1472 outliers_IDs=outliers_IDs.ravel()
1473 1473 outliers_IDs=numpy.unique(outliers_IDs)
1474 1474 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1475 1475 indexes=numpy.array(indexes)
1476 1476 indexmin=numpy.min(indexes)
1477 1477
1478 1478
1479 1479 #print(indexmin,buffer1.shape[0], k)
1480 1480
1481 1481 # fig,ax = plt.subplots()
1482 1482 # ax.plot(sortdata)
1483 1483 # ax2 = ax.twinx()
1484 1484 # x=range(len(indexes))
1485 1485 # #plt.scatter(x,indexes)
1486 1486 # ax2.scatter(x,indexes)
1487 1487 # plt.show()
1488 1488
1489 1489 if indexmin != buffer1.shape[0]:
1490 1490 if self.__nChannels > 1:
1491 1491 cspc_outliers_exist= True
1492 1492
1493 1493 lt=outliers_IDs
1494 1494 #avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1495 1495
1496 1496 for p in list(outliers_IDs):
1497 1497 #buffer1[p,:]=avg
1498 1498 buffer1[p,:] = numpy.NaN
1499 1499
1500 1500 self.dataOutliers[i,k] = len(outliers_IDs)
1501 1501
1502 1502
1503 1503 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1504 1504
1505 1505
1506 1506 if self.__nChannels > 1:
1507 1507 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1508 1508
1509 1509
1510 1510 if self.__nChannels > 1:
1511 1511 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1512 1512 if cspc_outliers_exist:
1513 1513
1514 1514 lt=outliers_IDs_cspc
1515 1515
1516 1516 #avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1517 1517 for p in list(outliers_IDs_cspc):
1518 1518 #buffer_cspc[p,:]=avg
1519 1519 buffer_cspc[p,:] = numpy.NaN
1520 1520
1521 1521 if self.__nChannels > 1:
1522 1522 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1523 1523
1524 1524
1525 1525
1526 1526
1527 1527 nOutliers = len(outliers_IDs)
1528 1528 #print("Outliers n: ",self.dataOutliers,nOutliers)
1529 1529 buffer=None
1530 1530 bufferH=None
1531 1531 buffer1=None
1532 1532 buffer_cspc=None
1533 1533
1534 1534
1535 1535 buffer=None
1536 1536
1537 1537 #data_spc = numpy.sum(self.__buffer_spc,axis=0)
1538 1538 data_spc = numpy.nansum(self.__buffer_spc,axis=0)
1539 1539 if self.__nChannels > 1:
1540 1540 #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1541 1541 data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
1542 1542 else:
1543 1543 data_cspc = None
1544 1544 data_dc = self.__buffer_dc
1545 1545 #(CH, HEIGH)
1546 1546 self.maxProfilesInt = self.__profIndex - 1
1547 1547 n = self.__profIndex - self.dataOutliers # n becomes a matrix
1548 1548
1549 1549 self.__buffer_spc = []
1550 1550 self.__buffer_cspc = []
1551 1551 self.__buffer_dc = 0
1552 1552 self.__profIndex = 0
1553 1553 #print("cleaned ",data_cspc)
1554 1554 return data_spc, data_cspc, data_dc, n
1555 1555
1556 1556 def byProfiles(self, *args):
1557 1557
1558 1558 self.__dataReady = False
1559 1559 avgdata_spc = None
1560 1560 avgdata_cspc = None
1561 1561 avgdata_dc = None
1562 1562
1563 1563 self.putData(*args)
1564 1564
1565 1565 if self.__profIndex == self.n:
1566 1566
1567 1567 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1568 1568 self.n_ints = n
1569 1569 self.__dataReady = True
1570 1570
1571 1571 return avgdata_spc, avgdata_cspc, avgdata_dc
1572 1572
1573 1573 def byTime(self, datatime, *args):
1574 1574
1575 1575 self.__dataReady = False
1576 1576 avgdata_spc = None
1577 1577 avgdata_cspc = None
1578 1578 avgdata_dc = None
1579 1579
1580 1580 self.putData(*args)
1581 1581
1582 1582 if (datatime - self.__initime) >= self.__integrationtime:
1583 1583 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1584 1584 self.n_ints = n
1585 1585 self.__dataReady = True
1586 1586
1587 1587 return avgdata_spc, avgdata_cspc, avgdata_dc
1588 1588
1589 1589 def integrate(self, datatime, *args):
1590 1590
1591 1591 if self.__profIndex == 0:
1592 1592 self.__initime = datatime
1593 1593
1594 1594 if self.__byTime:
1595 1595 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1596 1596 datatime, *args)
1597 1597 else:
1598 1598 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1599 1599
1600 1600 if not self.__dataReady:
1601 1601 return None, None, None, None
1602 1602
1603 1603 #print("integrate", avgdata_cspc)
1604 1604 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1605 1605
1606 1606 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
1607 1607 self.dataOut = dataOut
1608 1608 if n == 1:
1609 1609 return self.dataOut
1610 1610
1611 1611 #print("nchannels", self.dataOut.nChannels)
1612 1612 if self.dataOut.nChannels == 1:
1613 1613 self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
1614 1614 #print("IN spc:", self.dataOut.data_spc.shape, self.dataOut.data_cspc)
1615 1615 if not self.isConfig:
1616 1616 self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor)
1617 1617 self.isConfig = True
1618 1618
1619 1619 if not self.ByLags:
1620 1620 self.nProfiles=self.dataOut.nProfiles
1621 1621 self.nChannels=self.dataOut.nChannels
1622 1622 self.nHeights=self.dataOut.nHeights
1623 1623 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1624 1624 self.dataOut.data_spc,
1625 1625 self.dataOut.data_cspc,
1626 1626 self.dataOut.data_dc)
1627 1627 else:
1628 1628 self.nProfiles=self.dataOut.nProfiles
1629 1629 self.nChannels=self.dataOut.nChannels
1630 1630 self.nHeights=self.dataOut.nHeights
1631 1631 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1632 1632 self.dataOut.dataLag_spc,
1633 1633 self.dataOut.dataLag_cspc,
1634 1634 self.dataOut.dataLag_dc)
1635 1635 self.dataOut.flagNoData = True
1636 1636 if self.__dataReady:
1637 1637
1638 1638 if not self.ByLags:
1639 1639 if self.nChannels == 1:
1640 1640 #print("f int", avgdata_spc.shape)
1641 1641 self.dataOut.data_spc = avgdata_spc
1642 1642 self.dataOut.data_cspc = None
1643 1643 else:
1644 1644 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
1645 1645 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1646 1646 self.dataOut.data_dc = avgdata_dc
1647 1647 self.dataOut.data_outlier = self.dataOutliers
1648 1648
1649 1649 else:
1650 1650 self.dataOut.dataLag_spc = avgdata_spc
1651 1651 self.dataOut.dataLag_cspc = avgdata_cspc
1652 1652 self.dataOut.dataLag_dc = avgdata_dc
1653 1653
1654 1654 self.dataOut.data_spc=self.dataOut.dataLag_spc[:,:,:,self.dataOut.LagPlot]
1655 1655 self.dataOut.data_cspc=self.dataOut.dataLag_cspc[:,:,:,self.dataOut.LagPlot]
1656 1656 self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot]
1657 1657
1658 1658
1659 1659 self.dataOut.nIncohInt *= self.n_ints
1660 1660 #print("maxProfilesInt: ",self.maxProfilesInt)
1661 1661
1662 1662 self.dataOut.utctime = avgdatatime
1663 1663 self.dataOut.flagNoData = False
1664 1664 #print("Faraday Integration DONE...", self.dataOut.data_cspc)
1665 1665 #print(self.dataOut.flagNoData)
1666 1666 return self.dataOut
1667 1667
1668 1668
1669 1669
1670 1670 class removeInterference(Operation):
1671 1671
1672 1672 def removeInterference3(self, min_hei = None, max_hei = None):
1673 1673
1674 1674 jspectra = self.dataOut.data_spc
1675 1675 #jcspectra = self.dataOut.data_cspc
1676 1676 jnoise = self.dataOut.getNoise()
1677 1677 num_incoh = self.dataOut.max_nIncohInt
1678 1678 #print(jspectra.shape)
1679 1679 num_channel, num_prof, num_hei = jspectra.shape
1680 1680 minHei = min_hei
1681 1681 maxHei = max_hei
1682 1682 ########################################################################
1683 1683 if minHei == None or (minHei < self.dataOut.heightList[0]):
1684 1684 minHei = self.dataOut.heightList[0]
1685 1685
1686 1686 if maxHei == None or (maxHei > self.dataOut.heightList[-1]):
1687 1687 maxHei = self.dataOut.heightList[-1]
1688 1688 minIndex = 0
1689 1689 maxIndex = 0
1690 1690 heights = self.dataOut.heightList
1691 1691
1692 1692 inda = numpy.where(heights >= minHei)
1693 1693 indb = numpy.where(heights <= maxHei)
1694 1694
1695 1695 try:
1696 1696 minIndex = inda[0][0]
1697 1697 except:
1698 1698 minIndex = 0
1699 1699 try:
1700 1700 maxIndex = indb[0][-1]
1701 1701 except:
1702 1702 maxIndex = len(heights)
1703 1703
1704 1704 if (minIndex < 0) or (minIndex > maxIndex):
1705 1705 raise ValueError("some value in (%d,%d) is not valid" % (
1706 1706 minIndex, maxIndex))
1707 1707 if (maxIndex >= self.dataOut.nHeights):
1708 1708 maxIndex = self.dataOut.nHeights - 1
1709 1709
1710 1710 ########################################################################
1711 1711
1712 1712
1713 1713 #dataOut.max_nIncohInt * dataOut.nCohInt
1714 1714 norm = self.dataOut.nIncohInt /self.dataOut.max_nIncohInt
1715 1715 #print(norm.shape)
1716 1716 # Subrutina de Remocion de la Interferencia
1717 1717 for ich in range(num_channel):
1718 1718 # Se ordena los espectros segun su potencia (menor a mayor)
1719 1719 #power = jspectra[ich, mask_prof, :]
1720 1720 interf = jspectra[ich, :, minIndex:maxIndex]
1721 1721 #print(interf.shape)
1722 1722 inttef = interf.mean(axis=1)
1723 1723
1724 1724 for hei in range(num_hei):
1725 1725 temp = jspectra[ich,:, hei]
1726 1726 temp -= inttef
1727 1727 temp += jnoise[ich]*norm[ich,hei]
1728 1728 jspectra[ich,:, hei] = temp
1729 1729
1730 1730 # Guardar Resultados
1731 1731 self.dataOut.data_spc = jspectra
1732 1732 #self.dataOut.data_cspc = jcspectra
1733 1733
1734 1734 return 1
1735 1735
1736 1736 def removeInterference2(self):
1737 1737
1738 1738 cspc = self.dataOut.data_cspc
1739 1739 spc = self.dataOut.data_spc
1740 1740 Heights = numpy.arange(cspc.shape[2])
1741 1741 realCspc = numpy.abs(cspc)
1742 1742
1743 1743 for i in range(cspc.shape[0]):
1744 1744 LinePower= numpy.sum(realCspc[i], axis=0)
1745 1745 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1746 1746 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1747 1747 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1748 1748 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1749 1749 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1750 1750
1751 1751
1752 1752 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1753 1753 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1754 1754 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1755 1755 cspc[i,InterferenceRange,:] = numpy.NaN
1756 1756
1757 1757 self.dataOut.data_cspc = cspc
1758 1758
1759 1759 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1760 1760
1761 1761 jspectra = self.dataOut.data_spc
1762 1762 jcspectra = self.dataOut.data_cspc
1763 1763 jnoise = self.dataOut.getNoise()
1764 1764 #num_incoh = self.dataOut.nIncohInt
1765 1765 num_incoh = self.dataOut.max_nIncohInt
1766 1766 #print("spc: ", jspectra.shape, jcspectra)
1767 1767 num_channel = jspectra.shape[0]
1768 1768 num_prof = jspectra.shape[1]
1769 1769 num_hei = jspectra.shape[2]
1770 1770
1771 1771 # hei_interf
1772 1772 if hei_interf is None:
1773 1773 count_hei = int(num_hei / 2) # a half of total ranges
1774 1774 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1775 1775 hei_interf = numpy.asarray(hei_interf)[0]
1776 1776 #print(hei_interf)
1777 1777 # nhei_interf
1778 1778 if (nhei_interf == None):
1779 1779 nhei_interf = 5
1780 1780 if (nhei_interf < 1):
1781 1781 nhei_interf = 1
1782 1782 if (nhei_interf > count_hei):
1783 1783 nhei_interf = count_hei
1784 1784 if (offhei_interf == None):
1785 1785 offhei_interf = 0
1786 1786
1787 1787 ind_hei = list(range(num_hei))
1788 1788 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1789 1789 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1790 1790 mask_prof = numpy.asarray(list(range(num_prof)))
1791 1791 num_mask_prof = mask_prof.size
1792 1792 comp_mask_prof = [0, num_prof / 2]
1793 1793
1794 1794 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1795 1795 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1796 1796 jnoise = numpy.nan
1797 1797 noise_exist = jnoise[0] < numpy.Inf
1798 1798
1799 1799 # Subrutina de Remocion de la Interferencia
1800 1800 for ich in range(num_channel):
1801 1801 # Se ordena los espectros segun su potencia (menor a mayor)
1802 1802 power = jspectra[ich, mask_prof, :]
1803 1803 power = power[:, hei_interf]
1804 1804 power = power.sum(axis=0)
1805 1805 psort = power.ravel().argsort()
1806 1806 print(hei_interf[psort[list(range(offhei_interf, nhei_interf + offhei_interf))]])
1807 1807 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1808 1808 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1809 1809 offhei_interf, nhei_interf + offhei_interf))]]]
1810 1810
1811 1811 if noise_exist:
1812 1812 # tmp_noise = jnoise[ich] / num_prof
1813 1813 tmp_noise = jnoise[ich]
1814 1814 junkspc_interf = junkspc_interf - tmp_noise
1815 1815 #junkspc_interf[:,comp_mask_prof] = 0
1816 1816 print(junkspc_interf.shape)
1817 1817 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1818 1818 jspc_interf = jspc_interf.transpose()
1819 1819 # Calculando el espectro de interferencia promedio
1820 1820 noiseid = numpy.where(jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1821 1821 noiseid = noiseid[0]
1822 1822 cnoiseid = noiseid.size
1823 1823 interfid = numpy.where(jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1824 1824 interfid = interfid[0]
1825 1825 cinterfid = interfid.size
1826 1826
1827 1827 if (cnoiseid > 0):
1828 1828 jspc_interf[noiseid] = 0
1829 1829 # Expandiendo los perfiles a limpiar
1830 1830 if (cinterfid > 0):
1831 1831 new_interfid = (
1832 1832 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1833 1833 new_interfid = numpy.asarray(new_interfid)
1834 1834 new_interfid = {x for x in new_interfid}
1835 1835 new_interfid = numpy.array(list(new_interfid))
1836 1836 new_cinterfid = new_interfid.size
1837 1837 else:
1838 1838 new_cinterfid = 0
1839 1839
1840 1840 for ip in range(new_cinterfid):
1841 1841 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1842 1842 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1843 1843
1844 1844 jspectra[ich, :, ind_hei] = jspectra[ich, :,ind_hei] - jspc_interf # Corregir indices
1845 1845
1846 1846 # Removiendo la interferencia del punto de mayor interferencia
1847 1847 ListAux = jspc_interf[mask_prof].tolist()
1848 1848 maxid = ListAux.index(max(ListAux))
1849 1849 print(cinterfid)
1850 1850 if cinterfid > 0:
1851 1851 for ip in range(cinterfid * (interf == 2) - 1):
1852 1852 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1853 1853 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1854 1854 cind = len(ind)
1855 1855
1856 1856 if (cind > 0):
1857 1857 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1858 1858 (1 + (numpy.random.uniform(cind) - 0.5) /
1859 1859 numpy.sqrt(num_incoh))
1860 1860
1861 1861 ind = numpy.array([-2, -1, 1, 2])
1862 1862 xx = numpy.zeros([4, 4])
1863 1863
1864 1864 for id1 in range(4):
1865 1865 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1866 1866 xx_inv = numpy.linalg.inv(xx)
1867 1867 xx = xx_inv[:, 0]
1868 1868 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1869 1869 yy = jspectra[ich, mask_prof[ind], :]
1870 1870 jspectra[ich, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1871 1871
1872 1872 indAux = (jspectra[ich, :, :] < tmp_noise *
1873 1873 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1874 1874 print(indAux)
1875 1875 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1876 1876 (1 - 1 / numpy.sqrt(num_incoh))
1877 1877
1878 1878 # Remocion de Interferencia en el Cross Spectra
1879 1879 if jcspectra is None:
1880 1880 return jspectra, jcspectra
1881 1881 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1882 1882 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1883 1883
1884 1884 for ip in range(num_pairs):
1885 1885
1886 1886 #-------------------------------------------
1887 1887
1888 1888 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1889 1889 cspower = cspower[:, hei_interf]
1890 1890 cspower = cspower.sum(axis=0)
1891 1891
1892 1892 cspsort = cspower.ravel().argsort()
1893 1893 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1894 1894 offhei_interf, nhei_interf + offhei_interf))]]]
1895 1895 junkcspc_interf = junkcspc_interf.transpose()
1896 1896 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1897 1897
1898 1898 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1899 1899
1900 1900 median_real = int(numpy.median(numpy.real(
1901 1901 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1902 1902 median_imag = int(numpy.median(numpy.imag(
1903 1903 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1904 1904 comp_mask_prof = [int(e) for e in comp_mask_prof]
1905 1905 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1906 1906 median_real, median_imag)
1907 1907
1908 1908 for iprof in range(num_prof):
1909 1909 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1910 1910 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1911 1911
1912 1912 # Removiendo la Interferencia
1913 1913 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1914 1914 :, ind_hei] - jcspc_interf
1915 1915
1916 1916 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1917 1917 maxid = ListAux.index(max(ListAux))
1918 1918
1919 1919 ind = numpy.array([-2, -1, 1, 2])
1920 1920 xx = numpy.zeros([4, 4])
1921 1921
1922 1922 for id1 in range(4):
1923 1923 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1924 1924
1925 1925 xx_inv = numpy.linalg.inv(xx)
1926 1926 xx = xx_inv[:, 0]
1927 1927
1928 1928 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1929 1929 yy = jcspectra[ip, mask_prof[ind], :]
1930 1930 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1931 1931
1932 1932 # Guardar Resultados
1933 1933 self.dataOut.data_spc = jspectra
1934 1934 self.dataOut.data_cspc = jcspectra
1935 1935
1936 1936 return 1
1937 1937
1938 1938 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1, minHei=None, maxHei=None):
1939 1939
1940 1940 self.dataOut = dataOut
1941 1941
1942 1942 if mode == 1:
1943 1943 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1944 1944 elif mode == 2:
1945 1945 self.removeInterference2()
1946 1946 elif mode == 3:
1947 1947 self.removeInterference3(min_hei=minHei, max_hei=maxHei)
1948 1948 return self.dataOut
1949 1949
1950 1950
1951 1951 class IncohInt(Operation):
1952 1952
1953 1953 __profIndex = 0
1954 1954 __withOverapping = False
1955 1955
1956 1956 __byTime = False
1957 1957 __initime = None
1958 1958 __lastdatatime = None
1959 1959 __integrationtime = None
1960 1960
1961 1961 __buffer_spc = None
1962 1962 __buffer_cspc = None
1963 1963 __buffer_dc = None
1964 1964
1965 1965 __dataReady = False
1966 1966
1967 1967 __timeInterval = None
1968 1968 incohInt = 0
1969 1969 nOutliers = 0
1970 1970 n = None
1971 1971
1972 1972 def __init__(self):
1973 1973
1974 1974 Operation.__init__(self)
1975 1975
1976 1976 def setup(self, n=None, timeInterval=None, overlapping=False):
1977 1977 """
1978 1978 Set the parameters of the integration class.
1979 1979
1980 1980 Inputs:
1981 1981
1982 1982 n : Number of coherent integrations
1983 1983 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1984 1984 overlapping :
1985 1985
1986 1986 """
1987 1987
1988 1988 self.__initime = None
1989 1989 self.__lastdatatime = 0
1990 1990
1991 1991 self.__buffer_spc = 0
1992 1992 self.__buffer_cspc = 0
1993 1993 self.__buffer_dc = 0
1994 1994
1995 1995 self.__profIndex = 0
1996 1996 self.__dataReady = False
1997 1997 self.__byTime = False
1998 1998 self.incohInt = 0
1999 1999 self.nOutliers = 0
2000 2000 if n is None and timeInterval is None:
2001 2001 raise ValueError("n or timeInterval should be specified ...")
2002 2002
2003 2003 if n is not None:
2004 2004 self.n = int(n)
2005 2005 else:
2006 2006
2007 2007 self.__integrationtime = int(timeInterval)
2008 2008 self.n = None
2009 2009 self.__byTime = True
2010 2010
2011 2011 def putData(self, data_spc, data_cspc, data_dc):
2012 2012 """
2013 2013 Add a profile to the __buffer_spc and increase in one the __profileIndex
2014 2014
2015 2015 """
2016 2016 if data_spc.all() == numpy.nan :
2017 2017 print("nan ")
2018 2018 return
2019 2019 self.__buffer_spc += data_spc
2020 2020
2021 2021 if data_cspc is None:
2022 2022 self.__buffer_cspc = None
2023 2023 else:
2024 2024 self.__buffer_cspc += data_cspc
2025 2025
2026 2026 if data_dc is None:
2027 2027 self.__buffer_dc = None
2028 2028 else:
2029 2029 self.__buffer_dc += data_dc
2030 2030
2031 2031 self.__profIndex += 1
2032 2032
2033 2033 return
2034 2034
2035 2035 def pushData(self):
2036 2036 """
2037 2037 Return the sum of the last profiles and the profiles used in the sum.
2038 2038
2039 2039 Affected:
2040 2040
2041 2041 self.__profileIndex
2042 2042
2043 2043 """
2044 2044
2045 2045 data_spc = self.__buffer_spc
2046 2046 data_cspc = self.__buffer_cspc
2047 2047 data_dc = self.__buffer_dc
2048 2048 n = self.__profIndex
2049 2049
2050 2050 self.__buffer_spc = 0
2051 2051 self.__buffer_cspc = 0
2052 2052 self.__buffer_dc = 0
2053 2053
2054 2054
2055 2055 return data_spc, data_cspc, data_dc, n
2056 2056
2057 2057 def byProfiles(self, *args):
2058 2058
2059 2059 self.__dataReady = False
2060 2060 avgdata_spc = None
2061 2061 avgdata_cspc = None
2062 2062 avgdata_dc = None
2063 2063
2064 2064 self.putData(*args)
2065 2065
2066 2066 if self.__profIndex == self.n:
2067 2067
2068 2068 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
2069 2069 self.n = n
2070 2070 self.__dataReady = True
2071 2071
2072 2072 return avgdata_spc, avgdata_cspc, avgdata_dc
2073 2073
2074 2074 def byTime(self, datatime, *args):
2075 2075
2076 2076 self.__dataReady = False
2077 2077 avgdata_spc = None
2078 2078 avgdata_cspc = None
2079 2079 avgdata_dc = None
2080 2080
2081 2081 self.putData(*args)
2082 2082
2083 2083 if (datatime - self.__initime) >= self.__integrationtime:
2084 2084 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
2085 2085 self.n = n
2086 2086 self.__dataReady = True
2087 2087
2088 2088 return avgdata_spc, avgdata_cspc, avgdata_dc
2089 2089
2090 2090 def integrate(self, datatime, *args):
2091 2091
2092 2092 if self.__profIndex == 0:
2093 2093 self.__initime = datatime
2094 2094
2095 2095 if self.__byTime:
2096 2096 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
2097 2097 datatime, *args)
2098 2098 else:
2099 2099 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
2100 2100
2101 2101 if not self.__dataReady:
2102 2102 return None, None, None, None
2103 2103
2104 2104 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
2105 2105
2106 2106 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
2107 2107 if n == 1:
2108 2108 return dataOut
2109 2109
2110 2110 if dataOut.flagNoData == True:
2111 2111 return dataOut
2112 2112
2113 2113 dataOut.flagNoData = True
2114 2114
2115 2115 if not self.isConfig:
2116 2116 self.setup(n, timeInterval, overlapping)
2117 2117 self.isConfig = True
2118 2118
2119 2119 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
2120 2120 dataOut.data_spc,
2121 2121 dataOut.data_cspc,
2122 2122 dataOut.data_dc)
2123
2123 2124 self.incohInt += dataOut.nIncohInt
2124 self.nOutliers += dataOut.data_outlier
2125
2126 if isinstance(dataOut.data_outlier,numpy.ndarray) or isinstance(dataOut.data_outlier,int) or isinstance(dataOut.data_outlier, float):
2127 self.nOutliers += dataOut.data_outlier
2128
2125 2129 if self.__dataReady:
2126 2130 #print("prof: ",dataOut.max_nIncohInt,self.__profIndex)
2127 2131 dataOut.data_spc = avgdata_spc
2128 2132 dataOut.data_cspc = avgdata_cspc
2129 2133 dataOut.data_dc = avgdata_dc
2130 2134 dataOut.nIncohInt = self.incohInt
2131 2135 dataOut.data_outlier = self.nOutliers
2132 2136 dataOut.utctime = avgdatatime
2133 2137 dataOut.flagNoData = False
2134 2138 self.incohInt = 0
2135 2139 self.nOutliers = 0
2136 2140 self.__profIndex = 0
2137 2141 #print("IncohInt Done")
2138 2142 return dataOut
2139 2143
2140 2144 class dopplerFlip(Operation):
2141 2145
2142 2146 def run(self, dataOut):
2143 2147 # arreglo 1: (num_chan, num_profiles, num_heights)
2144 2148 self.dataOut = dataOut
2145 2149 # JULIA-oblicua, indice 2
2146 2150 # arreglo 2: (num_profiles, num_heights)
2147 2151 jspectra = self.dataOut.data_spc[2]
2148 2152 jspectra_tmp = numpy.zeros(jspectra.shape)
2149 2153 num_profiles = jspectra.shape[0]
2150 2154 freq_dc = int(num_profiles / 2)
2151 2155 # Flip con for
2152 2156 for j in range(num_profiles):
2153 2157 jspectra_tmp[num_profiles-j-1]= jspectra[j]
2154 2158 # Intercambio perfil de DC con perfil inmediato anterior
2155 2159 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
2156 2160 jspectra_tmp[freq_dc]= jspectra[freq_dc]
2157 2161 # canal modificado es re-escrito en el arreglo de canales
2158 2162 self.dataOut.data_spc[2] = jspectra_tmp
2159 2163
2160 2164 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now