##// END OF EJS Templates
Merge branch 'v3.0-devel' of http://jro-dev.igp.gob.pe/rhodecode/schain into v3.0-devel
Juan C. Espinoza -
r1300:fdbae09ac35b merge
parent child
Show More
@@ -0,0 +1,473
1 import numpy,math,random,time
2 #---------------1 Heredamos JRODatareader
3 from schainpy.model.io.jroIO_base import *
4 #---------------2 Heredamos las propiedades de ProcessingUnit
5 from schainpy.model.proc.jroproc_base import ProcessingUnit,Operation,MPDecorator
6 #---------------3 Importaremos las clases BascicHeader, SystemHeader, RadarControlHeader, ProcessingHeader
7 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader,SystemHeader,RadarControllerHeader, ProcessingHeader
8 #---------------4 Importaremos el objeto Voltge
9 from schainpy.model.data.jrodata import Voltage
10
11 class SimulatorReader(JRODataReader, ProcessingUnit):
12 incIntFactor = 1
13 nFFTPoints = 0
14 FixPP_IncInt = 1
15 FixRCP_IPP = 1000
16 FixPP_CohInt = 1
17 Tau_0 = 250
18 AcqH0_0 = 70
19 H0 = AcqH0_0
20 AcqDH_0 = 1.25
21 DH0 = AcqDH_0
22 Bauds = 32
23 BaudWidth = None
24 FixRCP_TXA = 40
25 FixRCP_TXB = 70
26 fAngle = 2.0*math.pi*(1/16)
27 DC_level = 500
28 stdev = 8
29 Num_Codes = 2
30 #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1])
31 #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0])
32 #Dyn_snCode = numpy.array([Num_Codes,Bauds])
33 Dyn_snCode = None
34 Samples = 200
35 channels = 5
36 pulses = None
37 Reference = None
38 pulse_size = None
39 prof_gen = None
40 Fdoppler = 100
41 Hdoppler = 36
42 Adoppler = 300
43 frequency = 9345
44 nTotalReadFiles = 1000
45
46 def __init__(self):
47 """
48 Inicializador de la clases SimulatorReader para
49 generar datos de voltage simulados.
50 Input:
51 dataOut: Objeto de la clase Voltage.
52 Este Objeto sera utilizado apra almacenar
53 un perfil de datos cada vez qe se haga
54 un requerimiento (getData)
55 """
56 ProcessingUnit.__init__(self)
57 print(" [ START ] init - Metodo Simulator Reader")
58
59 self.isConfig = False
60 self.basicHeaderObj = BasicHeader(LOCALTIME)
61 self.systemHeaderObj = SystemHeader()
62 self.radarControllerHeaderObj = RadarControllerHeader()
63 self.processingHeaderObj = ProcessingHeader()
64 self.profileIndex = 2**32-1
65 self.dataOut = Voltage()
66 #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1])
67 code0 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,1,1,1,-1,1,1,-1,1,-1,-1,-1,1,1,1,-1,1])
68 #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0])
69 code1 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,-1,-1,1,-1,-1,1,-1,1,1,1,-1,-1,-1,1,-1])
70 #self.Dyn_snCode = numpy.array([code0,code1])
71 self.Dyn_snCode = None
72
73 def set_kwargs(self, **kwargs):
74 for key, value in kwargs.items():
75 setattr(self, key, value)
76
77 def __hasNotDataInBuffer(self):
78
79 if self.profileIndex >= self.processingHeaderObj.profilesPerBlock* self.nTxs:
80 if self.nReadBlocks>0:
81 tmp = self.dataOut.utctime
82 tmp_utc = int(self.dataOut.utctime)
83 tmp_milisecond = int((tmp-tmp_utc)*1000)
84 self.basicHeaderObj.utc = tmp_utc
85 self.basicHeaderObj.miliSecond= tmp_milisecond
86 return 1
87 return 0
88
89 def setNextFile(self):
90 """Set the next file to be readed open it and parse de file header"""
91
92 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
93 self.nReadFiles=self.nReadFiles+1
94 if self.nReadFiles ==self.nTotalReadFiles:
95 self.flagNoMoreFiles=1
96 raise schainpy.admin.SchainWarning('No more files to read')
97
98 print('------------------- [Opening file] ------------------------------',self.nReadFiles)
99 self.nReadBlocks = 0
100
101 def __setNewBlock(self):
102 self.setNextFile()
103 if self.flagIsNewFile:
104 return 1
105
106 def readNextBlock(self):
107 while True:
108 self.__setNewBlock()
109 if not(self.readBlock()):
110 return 0
111 self.getBasicHeader()
112 break
113 if self.verbose:
114 print("[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
115 self.processingHeaderObj.dataBlocksPerFile,
116 self.dataOut.datatime.ctime()) )
117 return 1
118
119 def getFirstHeader(self):
120 self.getBasicHeader()
121 self.dataOut.processingHeaderObj = self.processingHeaderObj.copy()
122 self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
123 self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
124
125 self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
126 self.dataOut.heightList = numpy.arange(self.processingHeaderObj.nHeights) * self.processingHeaderObj.deltaHeight + self.processingHeaderObj.firstHeight
127 self.dataOut.channelList = list(range(self.systemHeaderObj.nChannels))
128 self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
129 # asumo q la data no esta decodificada
130 self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode
131 # asumo q la data no esta sin flip
132 self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip
133 self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft
134 self.dataOut.frequency = self.frequency
135
136 def getBasicHeader(self):
137 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
138 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
139
140 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
141 self.dataOut.timeZone = self.basicHeaderObj.timeZone
142 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
143 self.dataOut.errorCount = self.basicHeaderObj.errorCount
144 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
145 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
146
147 def set_RCH(self, expType=2, nTx=1,ipp=None, txA=0, txB=0,
148 nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None,
149 numTaus=0, line6Function=0, line5Function=0, fClock=None,
150 prePulseBefore=0, prePulseAfter=0,
151 codeType=0, nCode=0, nBaud=0, code=None,
152 flip1=0, flip2=0):
153 self.radarControllerHeaderObj.expType = expType
154 self.radarControllerHeaderObj.nTx = nTx
155 self.radarControllerHeaderObj.ipp = float(ipp)
156 self.radarControllerHeaderObj.txA = float(txA)
157 self.radarControllerHeaderObj.txB = float(txB)
158 self.radarControllerHeaderObj.rangeIPP = ipp
159 self.radarControllerHeaderObj.rangeTxA = txA
160 self.radarControllerHeaderObj.rangeTxB = txB
161
162 self.radarControllerHeaderObj.nHeights = int(nHeights)
163 self.radarControllerHeaderObj.firstHeight = numpy.array([firstHeight])
164 self.radarControllerHeaderObj.deltaHeight = numpy.array([deltaHeight])
165 self.radarControllerHeaderObj.samplesWin = numpy.array([nHeights])
166
167
168 self.radarControllerHeaderObj.nWindows = nWindows
169 self.radarControllerHeaderObj.numTaus = numTaus
170 self.radarControllerHeaderObj.codeType = codeType
171 self.radarControllerHeaderObj.line6Function = line6Function
172 self.radarControllerHeaderObj.line5Function = line5Function
173 self.radarControllerHeaderObj.fclock = fClock
174 self.radarControllerHeaderObj.prePulseBefore= prePulseBefore
175 self.radarControllerHeaderObj.prePulseAfter = prePulseAfter
176
177 self.radarControllerHeaderObj.nCode = nCode
178 self.radarControllerHeaderObj.nBaud = nBaud
179 self.radarControllerHeaderObj.code = code
180 self.radarControllerHeaderObj.flip1 = flip1
181 self.radarControllerHeaderObj.flip2 = flip2
182
183 self.radarControllerHeaderObj.code_size = int(numpy.ceil(nBaud / 32.)) * nCode * 4
184
185 if fClock is None and deltaHeight is not None:
186 self.fClock = 0.15 / (deltaHeight * 1e-6)
187
188 def set_PH(self, dtype=0, blockSize=0, profilesPerBlock=0,
189 dataBlocksPerFile=0, nWindows=0, processFlags=0, nCohInt=0,
190 nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0,
191 deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0,
192 code=0, nBaud=None, shif_fft=False, flag_dc=False,
193 flag_cspc=False, flag_decode=False, flag_deflip=False):
194
195 self.processingHeaderObj.profilesPerBlock = profilesPerBlock
196 self.processingHeaderObj.dataBlocksPerFile = dataBlocksPerFile
197 self.processingHeaderObj.nWindows = nWindows
198 self.processingHeaderObj.nCohInt = nCohInt
199 self.processingHeaderObj.nIncohInt = nIncohInt
200 self.processingHeaderObj.totalSpectra = totalSpectra
201 self.processingHeaderObj.nHeights = int(nHeights)
202 self.processingHeaderObj.firstHeight = firstHeight
203 self.processingHeaderObj.deltaHeight = deltaHeight
204 self.processingHeaderObj.samplesWin = nHeights
205
206 def set_BH(self, utc = 0, miliSecond = 0, timeZone = 0):
207 self.basicHeaderObj.utc = utc
208 self.basicHeaderObj.miliSecond = miliSecond
209 self.basicHeaderObj.timeZone = timeZone
210
211 def set_SH(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWidth=0):
212 self.systemHeaderObj.nSamples = nSamples
213 self.systemHeaderObj.nProfiles = nProfiles
214 self.systemHeaderObj.nChannels = nChannels
215 self.systemHeaderObj.adcResolution = adcResolution
216 self.systemHeaderObj.pciDioBusWidth = pciDioBusWidth
217
218 def init_acquisition(self):
219
220 if self.nFFTPoints != 0:
221 self.incIntFactor = m_nProfilesperBlock/self.nFFTPoints
222 if (self.FixPP_IncInt > self.incIntFactor):
223 self.incIntFactor = self.FixPP_IncInt/ self.incIntFactor
224 elif(self.FixPP_IncInt< self.incIntFactor):
225 print("False alert...")
226
227 ProfilesperBlock = self.processingHeaderObj.profilesPerBlock
228
229 self.timeperblock =int(((self.FixRCP_IPP
230 *ProfilesperBlock
231 *self.FixPP_CohInt
232 *self.incIntFactor)
233 /150.0)
234 *0.9
235 +0.5)
236 # para cada canal
237 self.profiles = ProfilesperBlock*self.FixPP_CohInt
238 self.profiles = ProfilesperBlock
239 self.Reference = int((self.Tau_0-self.AcqH0_0)/(self.AcqDH_0)+0.5)
240 self.BaudWidth = int((self.FixRCP_TXA/self.AcqDH_0)/self.Bauds + 0.5 )
241
242 if (self.BaudWidth==0):
243 self.BaudWidth=1
244
245 def init_pulse(self,Num_Codes=Num_Codes,Bauds=Bauds,BaudWidth=BaudWidth,Dyn_snCode=Dyn_snCode):
246
247 Num_Codes = Num_Codes
248 Bauds = Bauds
249 BaudWidth = BaudWidth
250 Dyn_snCode = Dyn_snCode
251
252 if Dyn_snCode:
253 print("EXISTE")
254 else:
255 print("No existe")
256
257 if Dyn_snCode: # if Bauds:
258 pulses = list(range(0,Num_Codes))
259 num_codes = Num_Codes
260 for i in range(num_codes):
261 pulse_size = Bauds*BaudWidth
262 pulses[i] = numpy.zeros(pulse_size)
263 for j in range(Bauds):
264 for k in range(BaudWidth):
265 pulses[i][j*BaudWidth+k] = int(Dyn_snCode[i][j]*600)
266 else:
267 print("sin code")
268 pulses = list(range(1))
269 if self.AcqDH_0>0.149:
270 pulse_size = int(self.FixRCP_TXB/0.15+0.5)
271 else:
272 pulse_size = int((self.FixRCP_TXB/self.AcqDH_0)+0.5) #0.0375
273 pulses[0] = numpy.ones(pulse_size)
274 pulses = 600*pulses[0]
275
276 return pulses,pulse_size
277
278 def jro_GenerateBlockOfData(self,Samples=Samples,DC_level= DC_level,stdev=stdev,
279 Reference= Reference,pulses= pulses,
280 Num_Codes= Num_Codes,pulse_size=pulse_size,
281 prof_gen= prof_gen,H0 = H0,DH0=DH0,
282 Adoppler=Adoppler,Fdoppler= Fdoppler,Hdoppler=Hdoppler):
283 Samples = Samples
284 DC_level = DC_level
285 stdev = stdev
286 m_nR = Reference
287 pulses = pulses
288 num_codes = Num_Codes
289 ps = pulse_size
290 prof_gen = prof_gen
291 channels = self.channels
292 H0 = H0
293 DH0 = DH0
294 ippSec = self.radarControllerHeaderObj.ippSeconds
295 Fdoppler = self.Fdoppler
296 Hdoppler = self.Hdoppler
297 Adoppler = self.Adoppler
298
299 self.datablock = numpy.zeros([channels,prof_gen,Samples],dtype= numpy.complex64)
300 for i in range(channels):
301 for k in range(prof_gen):
302 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·NOISEΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
303 Noise_r = numpy.random.normal(DC_level,stdev,Samples)
304 Noise_i = numpy.random.normal(DC_level,stdev,Samples)
305 Noise = numpy.zeros(Samples,dtype=complex)
306 Noise.real = Noise_r
307 Noise.imag = Noise_i
308 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·PULSOSΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
309 Pulso = numpy.zeros(pulse_size,dtype=complex)
310 Pulso.real = pulses[k%num_codes]
311 Pulso.imag = pulses[k%num_codes]
312 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· PULSES+NOISEΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·
313 InBuffer = numpy.zeros(Samples,dtype=complex)
314 InBuffer[m_nR:m_nR+ps] = Pulso
315 InBuffer = InBuffer+Noise
316 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· ANGLE Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
317 InBuffer.real[m_nR:m_nR+ps] = InBuffer.real[m_nR:m_nR+ps]*(math.cos( self.fAngle)*5)
318 InBuffer.imag[m_nR:m_nR+ps] = InBuffer.imag[m_nR:m_nR+ps]*(math.sin( self.fAngle)*5)
319 InBuffer=InBuffer
320 self.datablock[i][k]= InBuffer
321 #plot_cts(InBuffer,H0=H0,DH0=DH0)
322 #wave_fft(x=InBuffer,plot_show=True)
323 #time.sleep(1)
324 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·DOPPLER SIGNAL...............................................
325 time_vec = numpy.linspace(0,(prof_gen-1)*ippSec,int(prof_gen))+self.nReadBlocks*ippSec*prof_gen+(self.nReadFiles-1)*ippSec*prof_gen
326 fd = Fdoppler #+(600.0/120)*self.nReadBlocks
327 d_signal = Adoppler*numpy.array(numpy.exp(1.0j*2.0*math.pi*fd*time_vec),dtype=numpy.complex64)
328 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· DATABLOCK + DOPPLERΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·...........................
329 HD=int(Hdoppler/self.AcqDH_0)
330 for i in range(12):
331 self.datablock[:,:,HD+i]=self.datablock[:,:,HD+i]+ d_signal # RESULT
332 '''
333 a= numpy.zeros(10)
334 for i in range(10):
335 a[i]=i+self.nReadBlocks+20
336 for i in a:
337 self.datablock[0,:,int(i)]=self.datablock[0,:,int(i)]+ d_signal # RESULT
338 '''
339
340 def readBlock(self):
341
342 self.jro_GenerateBlockOfData(Samples= self.samples,DC_level=self.DC_level,
343 stdev=self.stdev,Reference= self.Reference,
344 pulses = self.pulses,Num_Codes=self.Num_Codes,
345 pulse_size=self.pulse_size,prof_gen=self.profiles,
346 H0=self.H0,DH0=self.DH0)
347
348 self.profileIndex = 0
349 self.flagIsNewFile = 0
350 self.flagIsNewBlock = 1
351 self.nTotalBlocks += 1
352 self.nReadBlocks += 1
353
354 return 1
355
356
357 def getData(self):
358 if self.flagNoMoreFiles:
359 self.dataOut.flagNodata = True
360 return 0
361 self.flagDiscontinuousBlock = 0
362 self.flagIsNewBlock = 0
363 if self.__hasNotDataInBuffer(): # aqui es verdad
364 if not(self.readNextBlock()): # return 1 y por eso el if not salta a getBasic Header
365 return 0
366 self.getFirstHeader() # atributo
367
368 if not self.getByBlock:
369 self.dataOut.flagDataAsBlock = False
370 self.dataOut.data = self.datablock[:, self.profileIndex, :]
371 self.dataOut.profileIndex = self.profileIndex
372 self.profileIndex += 1
373 else:
374 pass
375 self.dataOut.flagNoData = False
376 self.getBasicHeader()
377 self.dataOut.realtime = self.online
378 return self.dataOut.data
379
380
381 def setup(self,frequency=49.92e6,incIntFactor= 1, nFFTPoints = 0, FixPP_IncInt=1,FixRCP_IPP=1000,
382 FixPP_CohInt= 1,Tau_0= 250,AcqH0_0 = 70 ,AcqDH_0=1.25, Bauds= 32,
383 FixRCP_TXA = 40, FixRCP_TXB = 50, fAngle = 2.0*math.pi*(1/16),DC_level= 50,
384 stdev= 8,Num_Codes = 1 , Dyn_snCode = None, samples=200,
385 channels=2,Fdoppler=20,Hdoppler=36,Adoppler=500,nTotalReadFiles=10000,
386 **kwargs):
387
388 self.set_kwargs(**kwargs)
389 self.nReadBlocks = 0
390 self.nReadFiles = 1
391 print('------------------- [Opening file: ] ------------------------------',self.nReadFiles)
392
393 tmp = time.time()
394 tmp_utc = int(tmp)
395 tmp_milisecond = int((tmp-tmp_utc)*1000)
396 print(" SETUP -basicHeaderObj.utc",datetime.datetime.utcfromtimestamp(tmp))
397 if Dyn_snCode is None:
398 Num_Codes=1
399 Bauds =1
400
401
402
403 self.set_BH(utc= tmp_utc,miliSecond= tmp_milisecond,timeZone=300 )
404
405 self.set_RCH( expType=0, nTx=150,ipp=FixRCP_IPP, txA=FixRCP_TXA, txB= FixRCP_TXB,
406 nWindows=1 , nHeights=samples, firstHeight=AcqH0_0, deltaHeight=AcqDH_0,
407 numTaus=1, line6Function=0, line5Function=0, fClock=None,
408 prePulseBefore=0, prePulseAfter=0,
409 codeType=14, nCode=Num_Codes, nBaud=32, code=Dyn_snCode,
410 flip1=0, flip2=0)
411
412 self.set_PH(dtype=0, blockSize=0, profilesPerBlock=300,
413 dataBlocksPerFile=120, nWindows=1, processFlags=0, nCohInt=1,
414 nIncohInt=1, totalSpectra=0, nHeights=samples, firstHeight=AcqH0_0,
415 deltaHeight=AcqDH_0, samplesWin=samples, spectraComb=0, nCode=0,
416 code=0, nBaud=None, shif_fft=False, flag_dc=False,
417 flag_cspc=False, flag_decode=False, flag_deflip=False)
418
419 self.set_SH(nSamples=samples, nProfiles=300, nChannels=channels)
420
421
422 self.frequency = frequency
423 self.incIntFactor = incIntFactor
424 self.nFFTPoints = nFFTPoints
425 self.FixPP_IncInt = FixPP_IncInt
426 self.FixRCP_IPP = FixRCP_IPP
427 self.FixPP_CohInt = FixPP_CohInt
428 self.Tau_0 = Tau_0
429 self.AcqH0_0 = AcqH0_0
430 self.H0 = AcqH0_0
431 self.AcqDH_0 = AcqDH_0
432 self.DH0 = AcqDH_0
433 self.Bauds = Bauds
434 self.FixRCP_TXA = FixRCP_TXA
435 self.FixRCP_TXB = FixRCP_TXB
436 self.fAngle = fAngle
437 self.DC_level = DC_level
438 self.stdev = stdev
439 self.Num_Codes = Num_Codes
440 self.Dyn_snCode = Dyn_snCode
441 self.samples = samples
442 self.channels = channels
443 self.profiles = None
444 self.m_nReference = None
445 self.Baudwidth = None
446 self.Fdoppler = Fdoppler
447 self.Hdoppler = Hdoppler
448 self.Adoppler = Adoppler
449 self.nTotalReadFiles = int(nTotalReadFiles)
450
451 print("IPP ", self.FixRCP_IPP)
452 print("Tau_0 ",self.Tau_0)
453 print("AcqH0_0",self.AcqH0_0)
454 print("samples,window ",self.samples)
455 print("AcqDH_0",AcqDH_0)
456 print("FixRCP_TXA",self.FixRCP_TXA)
457 print("FixRCP_TXB",self.FixRCP_TXB)
458 print("Dyn_snCode",Dyn_snCode)
459 print("Fdoppler", Fdoppler)
460 print("Hdoppler",Hdoppler)
461 print("Vdopplermax",Fdoppler*(3.0e8/self.frequency)/2.0)
462 print("nTotalReadFiles", nTotalReadFiles)
463
464 self.init_acquisition()
465 self.pulses,self.pulse_size=self.init_pulse(Num_Codes=self.Num_Codes,Bauds=self.Bauds,BaudWidth=self.BaudWidth,Dyn_snCode=Dyn_snCode)
466 print(" [ END ] - SETUP metodo")
467 return
468
469 def run(self,**kwargs): # metodo propio
470 if not(self.isConfig):
471 self.setup(**kwargs)
472 self.isConfig = True
473 self.getData()
@@ -0,0 +1,61
1 import os,sys
2 import datetime
3 import time
4 from schainpy.controller import Project
5 path = '/home/alex/Downloads/NEW_WR2'
6 figpath = path
7 desc = "Simulator Test"
8
9 controllerObj = Project()
10
11 controllerObj.setup(id='10',name='Test Simulator',description=desc)
12
13 readUnitConfObj = controllerObj.addReadUnit(datatype='SimulatorReader',
14 frequency=9.345e9,
15 FixRCP_IPP= 60,
16 Tau_0 = 30,
17 AcqH0_0=0,
18 samples=330,
19 AcqDH_0=0.15,
20 FixRCP_TXA=0.15,
21 FixRCP_TXB=0.15,
22 Fdoppler=200.0,
23 Hdoppler=36,
24 Adoppler=300,
25 delay=0,
26 online=0,
27 walk=0)
28
29 #opObj11 = readUnitConfObj.addOperation(name='printInfo')
30
31 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
32
33 #opObj10 = procUnitConfObjA.addOperation(name='selectChannels')
34 #opObj10.addParameter(name='channelList', value=[0,1])
35 #opObj10.addParameter(name='channelList', value='0',format='intlist')
36
37 opObj11 = procUnitConfObjA.addOperation(name='PulsePairVoltage', optype='other')
38 opObj11.addParameter(name='n', value='32', format='int')
39 opObj11.addParameter(name='removeDC', value=1, format='int')
40 '''
41 type="power"
42 opObj10 = procUnitConfObjA.addOperation(name='ScopePlot', optype='external')
43 #opObj10.addParameter(name='id', value='12')
44 opObj10.addParameter(name='wintitle', value=type )
45 opObj10.addParameter(name='type', value=type)
46 '''
47
48 type="WeatherPower"
49 opObj10 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='external')
50 #opObj10.addParameter(name='id', value='12')
51 opObj10.addParameter(name='wintitle', value=type )
52
53 '''
54 type="WeatherVeloity"
55
56 opObj10 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='external')
57 #opObj10.addParameter(name='id', value='12')
58 opObj10.addParameter(name='wintitle', value=type )
59 '''
60
61 controllerObj.start()
@@ -0,0 +1,56
1 import os,sys
2 import datetime
3 import time
4 from schainpy.controller import Project
5 path = '/home/alex/Downloads/NEW_WR2/spc16removeDC'
6 figpath = path
7 desc = "Simulator Test"
8
9 controllerObj = Project()
10
11 controllerObj.setup(id='10',name='Test Simulator',description=desc)
12
13 readUnitConfObj = controllerObj.addReadUnit(datatype='SimulatorReader',
14 frequency=9.345e9,
15 FixRCP_IPP= 60,
16 Tau_0 = 30,
17 AcqH0_0=0,
18 samples=330,
19 AcqDH_0=0.15,
20 FixRCP_TXA=0.15,
21 FixRCP_TXB=0.15,
22 Fdoppler=200.0,
23 Hdoppler=36,
24 Adoppler=300,
25 delay=0,
26 online=0,
27 walk=0,
28 nTotalReadFiles=4)
29
30 opObj11 = readUnitConfObj.addOperation(name='printInfo')
31
32 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
33
34 opObj10 = procUnitConfObjA.addOperation(name='selectChannels')
35 opObj10.addParameter(name='channelList', value=[0,1])
36
37 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
38 procUnitConfObjB.addParameter(name='nFFTPoints', value=200, format='int')
39 procUnitConfObjB.addParameter(name='nProfiles', value=200, format='int')
40
41 opObj11 = procUnitConfObjB.addOperation(name='removeDC')
42 opObj11.addParameter(name='mode', value=2)
43
44 #opObj11 = procUnitConfObjB.addOperation(name='IncohInt', optype='other')
45 #opObj11.addParameter(name='n', value='20', format='int')
46
47 procUnitConfObjC = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjB.getId())
48 procUnitConfObjC.addOperation(name='SpectralMoments')
49
50 opObj11 = procUnitConfObjC.addOperation(name='SpectralMomentsPlot')
51 opObj11.addParameter(name='xmax', value=6)
52 opObj11.addParameter(name='save', value=figpath)
53 opObj11.addParameter(name='showprofile', value=1)
54 opObj11.addParameter(name='save_period', value=10)
55
56 controllerObj.start()
@@ -0,0 +1,40
1 import os,sys
2 import datetime
3 import time
4 from schainpy.controller import Project
5 path = '/home/alex/Downloads/NEW_WR2/spc16removeDC'
6 figpath = path
7 desc = "Simulator Test"
8
9 controllerObj = Project()
10
11 controllerObj.setup(id='10',name='Test Simulator',description=desc)
12
13 readUnitConfObj = controllerObj.addReadUnit(datatype='SimulatorReader',
14 frequency=9.345e9,
15 FixRCP_IPP= 60,
16 Tau_0 = 30,
17 AcqH0_0=0,
18 samples=330,
19 AcqDH_0=0.15,
20 FixRCP_TXA=0.15,
21 FixRCP_TXB=0.15,
22 Fdoppler=200.0,
23 Hdoppler=36,
24 Adoppler=300,
25 delay=0,
26 online=0,
27 walk=0,
28 nTotalReadFiles=4)
29
30 opObj11 = readUnitConfObj.addOperation(name='printInfo')
31 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
32
33 opObj10 = procUnitConfObjA.addOperation(name='selectChannels')
34 opObj10.addParameter(name='channelList', value=[0])
35
36 opObj11 = procUnitConfObjA.addOperation(name='PulsePairVoltage', optype='other')
37 opObj11.addParameter(name='n', value='32', format='int')#10
38 #opObj11.addParameter(name='removeDC', value=1, format='int')
39
40 controllerObj.start()
@@ -1,1387 +1,1394
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10 import json
11 11
12 12 import schainpy.admin
13 13 from schainpy.utils import log
14 14 from .jroheaderIO import SystemHeader, RadarControllerHeader
15 15 from schainpy.model.data import _noise
16 16
17 17
18 18 def getNumpyDtype(dataTypeCode):
19 19
20 20 if dataTypeCode == 0:
21 21 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
22 22 elif dataTypeCode == 1:
23 23 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
24 24 elif dataTypeCode == 2:
25 25 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
26 26 elif dataTypeCode == 3:
27 27 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
28 28 elif dataTypeCode == 4:
29 29 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
30 30 elif dataTypeCode == 5:
31 31 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
32 32 else:
33 33 raise ValueError('dataTypeCode was not defined')
34 34
35 35 return numpyDtype
36 36
37 37
38 38 def getDataTypeCode(numpyDtype):
39 39
40 40 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
41 41 datatype = 0
42 42 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
43 43 datatype = 1
44 44 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
45 45 datatype = 2
46 46 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
47 47 datatype = 3
48 48 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
49 49 datatype = 4
50 50 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
51 51 datatype = 5
52 52 else:
53 53 datatype = None
54 54
55 55 return datatype
56 56
57 57
58 58 def hildebrand_sekhon(data, navg):
59 59 """
60 60 This method is for the objective determination of the noise level in Doppler spectra. This
61 61 implementation technique is based on the fact that the standard deviation of the spectral
62 62 densities is equal to the mean spectral density for white Gaussian noise
63 63
64 64 Inputs:
65 65 Data : heights
66 66 navg : numbers of averages
67 67
68 68 Return:
69 69 mean : noise's level
70 70 """
71 71
72 72 sortdata = numpy.sort(data, axis=None)
73 73 '''
74 74 lenOfData = len(sortdata)
75 75 nums_min = lenOfData*0.2
76 76
77 77 if nums_min <= 5:
78 78
79 79 nums_min = 5
80 80
81 81 sump = 0.
82 82 sumq = 0.
83 83
84 84 j = 0
85 85 cont = 1
86 86
87 87 while((cont == 1)and(j < lenOfData)):
88 88
89 89 sump += sortdata[j]
90 90 sumq += sortdata[j]**2
91 91
92 92 if j > nums_min:
93 93 rtest = float(j)/(j-1) + 1.0/navg
94 94 if ((sumq*j) > (rtest*sump**2)):
95 95 j = j - 1
96 96 sump = sump - sortdata[j]
97 97 sumq = sumq - sortdata[j]**2
98 98 cont = 0
99 99
100 100 j += 1
101 101
102 102 lnoise = sump / j
103 103 '''
104 104 return _noise.hildebrand_sekhon(sortdata, navg)
105 105
106 106
107 107 class Beam:
108 108
109 109 def __init__(self):
110 110 self.codeList = []
111 111 self.azimuthList = []
112 112 self.zenithList = []
113 113
114 114
115 115 class GenericData(object):
116 116
117 117 flagNoData = True
118 118
119 119 def copy(self, inputObj=None):
120 120
121 121 if inputObj == None:
122 122 return copy.deepcopy(self)
123 123
124 124 for key in list(inputObj.__dict__.keys()):
125 125
126 126 attribute = inputObj.__dict__[key]
127 127
128 128 # If this attribute is a tuple or list
129 129 if type(inputObj.__dict__[key]) in (tuple, list):
130 130 self.__dict__[key] = attribute[:]
131 131 continue
132 132
133 133 # If this attribute is another object or instance
134 134 if hasattr(attribute, '__dict__'):
135 135 self.__dict__[key] = attribute.copy()
136 136 continue
137 137
138 138 self.__dict__[key] = inputObj.__dict__[key]
139 139
140 140 def deepcopy(self):
141 141
142 142 return copy.deepcopy(self)
143 143
144 144 def isEmpty(self):
145 145
146 146 return self.flagNoData
147 147
148 148 def isReady(self):
149 149
150 150 return not self.flagNoData
151 151
152 152
153 153 class JROData(GenericData):
154 154
155 155 # m_BasicHeader = BasicHeader()
156 156 # m_ProcessingHeader = ProcessingHeader()
157 157
158 158 systemHeaderObj = SystemHeader()
159 159 radarControllerHeaderObj = RadarControllerHeader()
160 160 # data = None
161 161 type = None
162 162 datatype = None # dtype but in string
163 163 # dtype = None
164 164 # nChannels = None
165 165 # nHeights = None
166 166 nProfiles = None
167 167 heightList = None
168 168 channelList = None
169 169 flagDiscontinuousBlock = False
170 170 useLocalTime = False
171 171 utctime = None
172 172 timeZone = None
173 173 dstFlag = None
174 174 errorCount = None
175 175 blocksize = None
176 176 # nCode = None
177 177 # nBaud = None
178 178 # code = None
179 179 flagDecodeData = False # asumo q la data no esta decodificada
180 180 flagDeflipData = False # asumo q la data no esta sin flip
181 181 flagShiftFFT = False
182 182 # ippSeconds = None
183 183 # timeInterval = None
184 184 nCohInt = None
185 185 # noise = None
186 186 windowOfFilter = 1
187 187 # Speed of ligth
188 188 C = 3e8
189 189 frequency = 49.92e6
190 190 realtime = False
191 191 beacon_heiIndexList = None
192 192 last_block = None
193 193 blocknow = None
194 194 azimuth = None
195 195 zenith = None
196 196 beam = Beam()
197 197 profileIndex = None
198 198 error = None
199 199 data = None
200 200 nmodes = None
201 201
202 202 def __str__(self):
203 203
204 204 return '{} - {}'.format(self.type, self.getDatatime())
205 205
206 206 def getNoise(self):
207 207
208 208 raise NotImplementedError
209 209
210 210 def getNChannels(self):
211 211
212 212 return len(self.channelList)
213 213
214 214 def getChannelIndexList(self):
215 215
216 216 return list(range(self.nChannels))
217 217
218 218 def getNHeights(self):
219 219
220 220 return len(self.heightList)
221 221
222 222 def getHeiRange(self, extrapoints=0):
223 223
224 224 heis = self.heightList
225 225 # deltah = self.heightList[1] - self.heightList[0]
226 226 #
227 227 # heis.append(self.heightList[-1])
228 228
229 229 return heis
230 230
231 231 def getDeltaH(self):
232 232
233 233 delta = self.heightList[1] - self.heightList[0]
234 234
235 235 return delta
236 236
237 237 def getltctime(self):
238 238
239 239 if self.useLocalTime:
240 240 return self.utctime - self.timeZone * 60
241 241
242 242 return self.utctime
243 243
244 244 def getDatatime(self):
245 245
246 246 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
247 247 return datatimeValue
248 248
249 249 def getTimeRange(self):
250 250
251 251 datatime = []
252 252
253 253 datatime.append(self.ltctime)
254 254 datatime.append(self.ltctime + self.timeInterval + 1)
255 255
256 256 datatime = numpy.array(datatime)
257 257
258 258 return datatime
259 259
260 260 def getFmaxTimeResponse(self):
261 261
262 262 period = (10**-6) * self.getDeltaH() / (0.15)
263 263
264 264 PRF = 1. / (period * self.nCohInt)
265 265
266 266 fmax = PRF
267 267
268 268 return fmax
269 269
270 270 def getFmax(self):
271 271 PRF = 1. / (self.ippSeconds * self.nCohInt)
272 272
273 273 fmax = PRF
274 274 return fmax
275 275
276 276 def getVmax(self):
277 277
278 278 _lambda = self.C / self.frequency
279 279
280 280 vmax = self.getFmax() * _lambda / 2
281 281
282 282 return vmax
283 283
284 284 def get_ippSeconds(self):
285 285 '''
286 286 '''
287 287 return self.radarControllerHeaderObj.ippSeconds
288 288
289 289 def set_ippSeconds(self, ippSeconds):
290 290 '''
291 291 '''
292 292
293 293 self.radarControllerHeaderObj.ippSeconds = ippSeconds
294 294
295 295 return
296 296
297 297 def get_dtype(self):
298 298 '''
299 299 '''
300 300 return getNumpyDtype(self.datatype)
301 301
302 302 def set_dtype(self, numpyDtype):
303 303 '''
304 304 '''
305 305
306 306 self.datatype = getDataTypeCode(numpyDtype)
307 307
308 308 def get_code(self):
309 309 '''
310 310 '''
311 311 return self.radarControllerHeaderObj.code
312 312
313 313 def set_code(self, code):
314 314 '''
315 315 '''
316 316 self.radarControllerHeaderObj.code = code
317 317
318 318 return
319 319
320 320 def get_ncode(self):
321 321 '''
322 322 '''
323 323 return self.radarControllerHeaderObj.nCode
324 324
325 325 def set_ncode(self, nCode):
326 326 '''
327 327 '''
328 328 self.radarControllerHeaderObj.nCode = nCode
329 329
330 330 return
331 331
332 332 def get_nbaud(self):
333 333 '''
334 334 '''
335 335 return self.radarControllerHeaderObj.nBaud
336 336
337 337 def set_nbaud(self, nBaud):
338 338 '''
339 339 '''
340 340 self.radarControllerHeaderObj.nBaud = nBaud
341 341
342 342 return
343 343
344 344 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
345 345 channelIndexList = property(
346 346 getChannelIndexList, "I'm the 'channelIndexList' property.")
347 347 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
348 348 #noise = property(getNoise, "I'm the 'nHeights' property.")
349 349 datatime = property(getDatatime, "I'm the 'datatime' property")
350 350 ltctime = property(getltctime, "I'm the 'ltctime' property")
351 351 ippSeconds = property(get_ippSeconds, set_ippSeconds)
352 352 dtype = property(get_dtype, set_dtype)
353 353 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
354 354 code = property(get_code, set_code)
355 355 nCode = property(get_ncode, set_ncode)
356 356 nBaud = property(get_nbaud, set_nbaud)
357 357
358 358
359 359 class Voltage(JROData):
360 360
361 361 # data es un numpy array de 2 dmensiones (canales, alturas)
362 362 data = None
363
363 data_intensity = None
364 data_velocity = None
364 365 def __init__(self):
365 366 '''
366 367 Constructor
367 368 '''
368 369
369 370 self.useLocalTime = True
370 371 self.radarControllerHeaderObj = RadarControllerHeader()
371 372 self.systemHeaderObj = SystemHeader()
372 373 self.type = "Voltage"
373 374 self.data = None
374 375 # self.dtype = None
375 376 # self.nChannels = 0
376 377 # self.nHeights = 0
377 378 self.nProfiles = None
378 379 self.heightList = None
379 380 self.channelList = None
380 381 # self.channelIndexList = None
381 382 self.flagNoData = True
382 383 self.flagDiscontinuousBlock = False
383 384 self.utctime = None
384 385 self.timeZone = None
385 386 self.dstFlag = None
386 387 self.errorCount = None
387 388 self.nCohInt = None
388 389 self.blocksize = None
389 390 self.flagDecodeData = False # asumo q la data no esta decodificada
390 391 self.flagDeflipData = False # asumo q la data no esta sin flip
391 392 self.flagShiftFFT = False
392 393 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
393 394 self.profileIndex = 0
394 395
395 396 def getNoisebyHildebrand(self, channel=None):
396 397 """
397 398 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
398 399
399 400 Return:
400 401 noiselevel
401 402 """
402 403
403 404 if channel != None:
404 405 data = self.data[channel]
405 406 nChannels = 1
406 407 else:
407 408 data = self.data
408 409 nChannels = self.nChannels
409 410
410 411 noise = numpy.zeros(nChannels)
411 412 power = data * numpy.conjugate(data)
412 413
413 414 for thisChannel in range(nChannels):
414 415 if nChannels == 1:
415 416 daux = power[:].real
416 417 else:
417 418 daux = power[thisChannel, :].real
418 419 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
419 420
420 421 return noise
421 422
422 423 def getNoise(self, type=1, channel=None):
423 424
424 425 if type == 1:
425 426 noise = self.getNoisebyHildebrand(channel)
426 427
427 428 return noise
428 429
429 430 def getPower(self, channel=None):
430 431
431 432 if channel != None:
432 433 data = self.data[channel]
433 434 else:
434 435 data = self.data
435 436
436 437 power = data * numpy.conjugate(data)
437 438 powerdB = 10 * numpy.log10(power.real)
438 439 powerdB = numpy.squeeze(powerdB)
439 440
440 441 return powerdB
441 442
442 443 def getTimeInterval(self):
443 444
444 445 timeInterval = self.ippSeconds * self.nCohInt
445 446
446 447 return timeInterval
447 448
448 449 noise = property(getNoise, "I'm the 'nHeights' property.")
449 450 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
450 451
451 452
452 453 class Spectra(JROData):
453 454
454 455 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
455 456 data_spc = None
456 457 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
457 458 data_cspc = None
458 459 # data dc es un numpy array de 2 dmensiones (canales, alturas)
459 460 data_dc = None
460 461 # data power
461 462 data_pwr = None
462 463 nFFTPoints = None
463 464 # nPairs = None
464 465 pairsList = None
465 466 nIncohInt = None
466 467 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
467 468 nCohInt = None # se requiere para determinar el valor de timeInterval
468 469 ippFactor = None
469 470 profileIndex = 0
470 471 plotting = "spectra"
471 472
472 473 def __init__(self):
473 474 '''
474 475 Constructor
475 476 '''
476 477
477 478 self.useLocalTime = True
478 479 self.radarControllerHeaderObj = RadarControllerHeader()
479 480 self.systemHeaderObj = SystemHeader()
480 481 self.type = "Spectra"
481 482 # self.data = None
482 483 # self.dtype = None
483 484 # self.nChannels = 0
484 485 # self.nHeights = 0
485 486 self.nProfiles = None
486 487 self.heightList = None
487 488 self.channelList = None
488 489 # self.channelIndexList = None
489 490 self.pairsList = None
490 491 self.flagNoData = True
491 492 self.flagDiscontinuousBlock = False
492 493 self.utctime = None
493 494 self.nCohInt = None
494 495 self.nIncohInt = None
495 496 self.blocksize = None
496 497 self.nFFTPoints = None
497 498 self.wavelength = None
498 499 self.flagDecodeData = False # asumo q la data no esta decodificada
499 500 self.flagDeflipData = False # asumo q la data no esta sin flip
500 501 self.flagShiftFFT = False
501 502 self.ippFactor = 1
502 503 #self.noise = None
503 504 self.beacon_heiIndexList = []
504 505 self.noise_estimation = None
505 506
506 507 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
507 508 """
508 509 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
509 510
510 511 Return:
511 512 noiselevel
512 513 """
513 514
514 515 noise = numpy.zeros(self.nChannels)
515 516
516 517 for channel in range(self.nChannels):
517 518 daux = self.data_spc[channel,
518 519 xmin_index:xmax_index, ymin_index:ymax_index]
519 520 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
520 521
521 522 return noise
522 523
523 524 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
524 525
525 526 if self.noise_estimation is not None:
526 527 # this was estimated by getNoise Operation defined in jroproc_spectra.py
527 528 return self.noise_estimation
528 529 else:
529 530 noise = self.getNoisebyHildebrand(
530 531 xmin_index, xmax_index, ymin_index, ymax_index)
531 532 return noise
532 533
533 534 def getFreqRangeTimeResponse(self, extrapoints=0):
534 535
535 536 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
536 537 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
537 538
538 539 return freqrange
539 540
540 541 def getAcfRange(self, extrapoints=0):
541 542
542 543 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
543 544 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
544 545
545 546 return freqrange
546 547
547 548 def getFreqRange(self, extrapoints=0):
548 549
549 550 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
550 551 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
551 552
552 553 return freqrange
553 554
554 555 def getVelRange(self, extrapoints=0):
555 556
556 557 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
557 558 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
558 559
559 560 if self.nmodes:
560 561 return velrange/self.nmodes
561 562 else:
562 563 return velrange
563 564
564 565 def getNPairs(self):
565 566
566 567 return len(self.pairsList)
567 568
568 569 def getPairsIndexList(self):
569 570
570 571 return list(range(self.nPairs))
571 572
572 573 def getNormFactor(self):
573 574
574 575 pwcode = 1
575 576
576 577 if self.flagDecodeData:
577 578 pwcode = numpy.sum(self.code[0]**2)
578 579 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
579 580 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
580 581
581 582 return normFactor
582 583
583 584 def getFlagCspc(self):
584 585
585 586 if self.data_cspc is None:
586 587 return True
587 588
588 589 return False
589 590
590 591 def getFlagDc(self):
591 592
592 593 if self.data_dc is None:
593 594 return True
594 595
595 596 return False
596 597
597 598 def getTimeInterval(self):
598 599
599 600 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
600 601 if self.nmodes:
601 602 return self.nmodes*timeInterval
602 603 else:
603 604 return timeInterval
604 605
605 606 def getPower(self):
606 607
607 608 factor = self.normFactor
608 609 z = self.data_spc / factor
609 610 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
610 611 avg = numpy.average(z, axis=1)
611 612
612 613 return 10 * numpy.log10(avg)
613 614
614 615 def getCoherence(self, pairsList=None, phase=False):
615 616
616 617 z = []
617 618 if pairsList is None:
618 619 pairsIndexList = self.pairsIndexList
619 620 else:
620 621 pairsIndexList = []
621 622 for pair in pairsList:
622 623 if pair not in self.pairsList:
623 624 raise ValueError("Pair %s is not in dataOut.pairsList" % (
624 625 pair))
625 626 pairsIndexList.append(self.pairsList.index(pair))
626 627 for i in range(len(pairsIndexList)):
627 628 pair = self.pairsList[pairsIndexList[i]]
628 629 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
629 630 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
630 631 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
631 632 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
632 633 if phase:
633 634 data = numpy.arctan2(avgcoherenceComplex.imag,
634 635 avgcoherenceComplex.real) * 180 / numpy.pi
635 636 else:
636 637 data = numpy.abs(avgcoherenceComplex)
637 638
638 639 z.append(data)
639 640
640 641 return numpy.array(z)
641 642
642 643 def setValue(self, value):
643 644
644 645 print("This property should not be initialized")
645 646
646 647 return
647 648
648 649 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
649 650 pairsIndexList = property(
650 651 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
651 652 normFactor = property(getNormFactor, setValue,
652 653 "I'm the 'getNormFactor' property.")
653 654 flag_cspc = property(getFlagCspc, setValue)
654 655 flag_dc = property(getFlagDc, setValue)
655 656 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
656 657 timeInterval = property(getTimeInterval, setValue,
657 658 "I'm the 'timeInterval' property")
658 659
659 660
660 661 class SpectraHeis(Spectra):
661 662
662 663 data_spc = None
663 664 data_cspc = None
664 665 data_dc = None
665 666 nFFTPoints = None
666 667 # nPairs = None
667 668 pairsList = None
668 669 nCohInt = None
669 670 nIncohInt = None
670 671
671 672 def __init__(self):
672 673
673 674 self.radarControllerHeaderObj = RadarControllerHeader()
674 675
675 676 self.systemHeaderObj = SystemHeader()
676 677
677 678 self.type = "SpectraHeis"
678 679
679 680 # self.dtype = None
680 681
681 682 # self.nChannels = 0
682 683
683 684 # self.nHeights = 0
684 685
685 686 self.nProfiles = None
686 687
687 688 self.heightList = None
688 689
689 690 self.channelList = None
690 691
691 692 # self.channelIndexList = None
692 693
693 694 self.flagNoData = True
694 695
695 696 self.flagDiscontinuousBlock = False
696 697
697 698 # self.nPairs = 0
698 699
699 700 self.utctime = None
700 701
701 702 self.blocksize = None
702 703
703 704 self.profileIndex = 0
704 705
705 706 self.nCohInt = 1
706 707
707 708 self.nIncohInt = 1
708 709
709 710 def getNormFactor(self):
710 711 pwcode = 1
711 712 if self.flagDecodeData:
712 713 pwcode = numpy.sum(self.code[0]**2)
713 714
714 715 normFactor = self.nIncohInt * self.nCohInt * pwcode
715 716
716 717 return normFactor
717 718
718 719 def getTimeInterval(self):
719 720
720 721 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
721 722
722 723 return timeInterval
723 724
724 725 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
725 726 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
726 727
727 728
728 729 class Fits(JROData):
729 730
730 731 heightList = None
731 732 channelList = None
732 733 flagNoData = True
733 734 flagDiscontinuousBlock = False
734 735 useLocalTime = False
735 736 utctime = None
736 737 timeZone = None
737 738 # ippSeconds = None
738 739 # timeInterval = None
739 740 nCohInt = None
740 741 nIncohInt = None
741 742 noise = None
742 743 windowOfFilter = 1
743 744 # Speed of ligth
744 745 C = 3e8
745 746 frequency = 49.92e6
746 747 realtime = False
747 748
748 749 def __init__(self):
749 750
750 751 self.type = "Fits"
751 752
752 753 self.nProfiles = None
753 754
754 755 self.heightList = None
755 756
756 757 self.channelList = None
757 758
758 759 # self.channelIndexList = None
759 760
760 761 self.flagNoData = True
761 762
762 763 self.utctime = None
763 764
764 765 self.nCohInt = 1
765 766
766 767 self.nIncohInt = 1
767 768
768 769 self.useLocalTime = True
769 770
770 771 self.profileIndex = 0
771 772
772 773 # self.utctime = None
773 774 # self.timeZone = None
774 775 # self.ltctime = None
775 776 # self.timeInterval = None
776 777 # self.header = None
777 778 # self.data_header = None
778 779 # self.data = None
779 780 # self.datatime = None
780 781 # self.flagNoData = False
781 782 # self.expName = ''
782 783 # self.nChannels = None
783 784 # self.nSamples = None
784 785 # self.dataBlocksPerFile = None
785 786 # self.comments = ''
786 787 #
787 788
788 789 def getltctime(self):
789 790
790 791 if self.useLocalTime:
791 792 return self.utctime - self.timeZone * 60
792 793
793 794 return self.utctime
794 795
795 796 def getDatatime(self):
796 797
797 798 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
798 799 return datatime
799 800
800 801 def getTimeRange(self):
801 802
802 803 datatime = []
803 804
804 805 datatime.append(self.ltctime)
805 806 datatime.append(self.ltctime + self.timeInterval)
806 807
807 808 datatime = numpy.array(datatime)
808 809
809 810 return datatime
810 811
811 812 def getHeiRange(self):
812 813
813 814 heis = self.heightList
814 815
815 816 return heis
816 817
817 818 def getNHeights(self):
818 819
819 820 return len(self.heightList)
820 821
821 822 def getNChannels(self):
822 823
823 824 return len(self.channelList)
824 825
825 826 def getChannelIndexList(self):
826 827
827 828 return list(range(self.nChannels))
828 829
829 830 def getNoise(self, type=1):
830 831
831 832 #noise = numpy.zeros(self.nChannels)
832 833
833 834 if type == 1:
834 835 noise = self.getNoisebyHildebrand()
835 836
836 837 if type == 2:
837 838 noise = self.getNoisebySort()
838 839
839 840 if type == 3:
840 841 noise = self.getNoisebyWindow()
841 842
842 843 return noise
843 844
844 845 def getTimeInterval(self):
845 846
846 847 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
847 848
848 849 return timeInterval
849 850
850 851 def get_ippSeconds(self):
851 852 '''
852 853 '''
853 854 return self.ipp_sec
854 855
855 856
856 857 datatime = property(getDatatime, "I'm the 'datatime' property")
857 858 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
858 859 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
859 860 channelIndexList = property(
860 861 getChannelIndexList, "I'm the 'channelIndexList' property.")
861 862 noise = property(getNoise, "I'm the 'nHeights' property.")
862 863
863 864 ltctime = property(getltctime, "I'm the 'ltctime' property")
864 865 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
865 866 ippSeconds = property(get_ippSeconds, '')
866 867
867 868 class Correlation(JROData):
868 869
869 870 noise = None
870 871 SNR = None
871 872 #--------------------------------------------------
872 873 mode = None
873 874 split = False
874 875 data_cf = None
875 876 lags = None
876 877 lagRange = None
877 878 pairsList = None
878 879 normFactor = None
879 880 #--------------------------------------------------
880 881 # calculateVelocity = None
881 882 nLags = None
882 883 nPairs = None
883 884 nAvg = None
884 885
885 886 def __init__(self):
886 887 '''
887 888 Constructor
888 889 '''
889 890 self.radarControllerHeaderObj = RadarControllerHeader()
890 891
891 892 self.systemHeaderObj = SystemHeader()
892 893
893 894 self.type = "Correlation"
894 895
895 896 self.data = None
896 897
897 898 self.dtype = None
898 899
899 900 self.nProfiles = None
900 901
901 902 self.heightList = None
902 903
903 904 self.channelList = None
904 905
905 906 self.flagNoData = True
906 907
907 908 self.flagDiscontinuousBlock = False
908 909
909 910 self.utctime = None
910 911
911 912 self.timeZone = None
912 913
913 914 self.dstFlag = None
914 915
915 916 self.errorCount = None
916 917
917 918 self.blocksize = None
918 919
919 920 self.flagDecodeData = False # asumo q la data no esta decodificada
920 921
921 922 self.flagDeflipData = False # asumo q la data no esta sin flip
922 923
923 924 self.pairsList = None
924 925
925 926 self.nPoints = None
926 927
927 928 def getPairsList(self):
928 929
929 930 return self.pairsList
930 931
931 932 def getNoise(self, mode=2):
932 933
933 934 indR = numpy.where(self.lagR == 0)[0][0]
934 935 indT = numpy.where(self.lagT == 0)[0][0]
935 936
936 937 jspectra0 = self.data_corr[:, :, indR, :]
937 938 jspectra = copy.copy(jspectra0)
938 939
939 940 num_chan = jspectra.shape[0]
940 941 num_hei = jspectra.shape[2]
941 942
942 943 freq_dc = jspectra.shape[1] / 2
943 944 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
944 945
945 946 if ind_vel[0] < 0:
946 947 ind_vel[list(range(0, 1))] = ind_vel[list(
947 948 range(0, 1))] + self.num_prof
948 949
949 950 if mode == 1:
950 951 jspectra[:, freq_dc, :] = (
951 952 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
952 953
953 954 if mode == 2:
954 955
955 956 vel = numpy.array([-2, -1, 1, 2])
956 957 xx = numpy.zeros([4, 4])
957 958
958 959 for fil in range(4):
959 960 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
960 961
961 962 xx_inv = numpy.linalg.inv(xx)
962 963 xx_aux = xx_inv[0, :]
963 964
964 965 for ich in range(num_chan):
965 966 yy = jspectra[ich, ind_vel, :]
966 967 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
967 968
968 969 junkid = jspectra[ich, freq_dc, :] <= 0
969 970 cjunkid = sum(junkid)
970 971
971 972 if cjunkid.any():
972 973 jspectra[ich, freq_dc, junkid.nonzero()] = (
973 974 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
974 975
975 976 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
976 977
977 978 return noise
978 979
979 980 def getTimeInterval(self):
980 981
981 982 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
982 983
983 984 return timeInterval
984 985
985 986 def splitFunctions(self):
986 987
987 988 pairsList = self.pairsList
988 989 ccf_pairs = []
989 990 acf_pairs = []
990 991 ccf_ind = []
991 992 acf_ind = []
992 993 for l in range(len(pairsList)):
993 994 chan0 = pairsList[l][0]
994 995 chan1 = pairsList[l][1]
995 996
996 997 # Obteniendo pares de Autocorrelacion
997 998 if chan0 == chan1:
998 999 acf_pairs.append(chan0)
999 1000 acf_ind.append(l)
1000 1001 else:
1001 1002 ccf_pairs.append(pairsList[l])
1002 1003 ccf_ind.append(l)
1003 1004
1004 1005 data_acf = self.data_cf[acf_ind]
1005 1006 data_ccf = self.data_cf[ccf_ind]
1006 1007
1007 1008 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1008 1009
1009 1010 def getNormFactor(self):
1010 1011 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1011 1012 acf_pairs = numpy.array(acf_pairs)
1012 1013 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1013 1014
1014 1015 for p in range(self.nPairs):
1015 1016 pair = self.pairsList[p]
1016 1017
1017 1018 ch0 = pair[0]
1018 1019 ch1 = pair[1]
1019 1020
1020 1021 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1021 1022 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1022 1023 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1023 1024
1024 1025 return normFactor
1025 1026
1026 1027 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1027 1028 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1028 1029
1029 1030
1030 1031 class Parameters(Spectra):
1031 1032
1032 1033 experimentInfo = None # Information about the experiment
1033 1034 # Information from previous data
1034 1035 inputUnit = None # Type of data to be processed
1035 1036 operation = None # Type of operation to parametrize
1036 1037 # normFactor = None #Normalization Factor
1037 1038 groupList = None # List of Pairs, Groups, etc
1038 1039 # Parameters
1039 1040 data_param = None # Parameters obtained
1040 1041 data_pre = None # Data Pre Parametrization
1041 1042 data_SNR = None # Signal to Noise Ratio
1042 1043 # heightRange = None #Heights
1043 1044 abscissaList = None # Abscissa, can be velocities, lags or time
1044 1045 # noise = None #Noise Potency
1045 1046 utctimeInit = None # Initial UTC time
1046 1047 paramInterval = None # Time interval to calculate Parameters in seconds
1047 1048 useLocalTime = True
1048 1049 # Fitting
1049 1050 data_error = None # Error of the estimation
1050 1051 constants = None
1051 1052 library = None
1052 1053 # Output signal
1053 1054 outputInterval = None # Time interval to calculate output signal in seconds
1054 1055 data_output = None # Out signal
1055 1056 nAvg = None
1056 1057 noise_estimation = None
1057 1058 GauSPC = None # Fit gaussian SPC
1058 1059
1059 1060 def __init__(self):
1060 1061 '''
1061 1062 Constructor
1062 1063 '''
1063 1064 self.radarControllerHeaderObj = RadarControllerHeader()
1064 1065
1065 1066 self.systemHeaderObj = SystemHeader()
1066 1067
1067 1068 self.type = "Parameters"
1068 1069
1069 1070 def getTimeRange1(self, interval):
1070 1071
1071 1072 datatime = []
1072 1073
1073 1074 if self.useLocalTime:
1074 1075 time1 = self.utctimeInit - self.timeZone * 60
1075 1076 else:
1076 1077 time1 = self.utctimeInit
1077 1078
1078 1079 datatime.append(time1)
1079 1080 datatime.append(time1 + interval)
1080 1081 datatime = numpy.array(datatime)
1081 1082
1082 1083 return datatime
1083 1084
1084 1085 def getTimeInterval(self):
1085 1086
1086 1087 if hasattr(self, 'timeInterval1'):
1087 1088 return self.timeInterval1
1088 1089 else:
1089 1090 return self.paramInterval
1090 1091
1091 1092 def setValue(self, value):
1092 1093
1093 1094 print("This property should not be initialized")
1094 1095
1095 1096 return
1096 1097
1097 1098 def getNoise(self):
1098 1099
1099 1100 return self.spc_noise
1100 1101
1101 1102 timeInterval = property(getTimeInterval)
1102 1103 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1103 1104
1104 1105
1105 1106 class PlotterData(object):
1106 1107 '''
1107 1108 Object to hold data to be plotted
1108 1109 '''
1109 1110
1110 1111 MAXNUMX = 200
1111 1112 MAXNUMY = 200
1112 1113
1113 1114 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1114 1115
1115 1116 self.key = code
1116 1117 self.throttle = throttle_value
1117 1118 self.exp_code = exp_code
1118 1119 self.buffering = buffering
1119 1120 self.ready = False
1120 1121 self.flagNoData = False
1121 1122 self.localtime = False
1122 1123 self.data = {}
1123 1124 self.meta = {}
1124 1125 self.__times = []
1125 1126 self.__heights = []
1126 1127
1127 1128 if 'snr' in code:
1128 1129 self.plottypes = ['snr']
1129 1130 elif code == 'spc':
1130 1131 self.plottypes = ['spc', 'noise', 'rti']
1131 1132 elif code == 'rti':
1132 1133 self.plottypes = ['noise', 'rti']
1133 1134 else:
1134 1135 self.plottypes = [code]
1135 1136
1136 1137 if 'snr' not in self.plottypes and snr:
1137 1138 self.plottypes.append('snr')
1138 1139
1139 1140 for plot in self.plottypes:
1140 1141 self.data[plot] = {}
1141 1142
1143
1142 1144 def __str__(self):
1143 1145 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1144 1146 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1145 1147
1146 1148 def __len__(self):
1147 1149 return len(self.__times)
1148 1150
1149 1151 def __getitem__(self, key):
1150 1152
1151 1153 if key not in self.data:
1152 1154 raise KeyError(log.error('Missing key: {}'.format(key)))
1153 1155 if 'spc' in key or not self.buffering:
1154 1156 ret = self.data[key]
1155 1157 elif 'scope' in key:
1156 1158 ret = numpy.array(self.data[key][float(self.tm)])
1157 1159 else:
1158 1160 ret = numpy.array([self.data[key][x] for x in self.times])
1159 1161 if ret.ndim > 1:
1160 1162 ret = numpy.swapaxes(ret, 0, 1)
1161 1163 return ret
1162 1164
1163 1165 def __contains__(self, key):
1164 1166 return key in self.data
1165 1167
1166 1168 def setup(self):
1167 1169 '''
1168 1170 Configure object
1169 1171 '''
1170
1171 1172 self.type = ''
1172 1173 self.ready = False
1173 1174 self.data = {}
1174 1175 self.__times = []
1175 1176 self.__heights = []
1176 1177 self.__all_heights = set()
1177 1178 for plot in self.plottypes:
1178 1179 if 'snr' in plot:
1179 1180 plot = 'snr'
1180 1181 elif 'spc_moments' == plot:
1181 1182 plot = 'moments'
1182 1183 self.data[plot] = {}
1183 1184
1184 1185 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1185 1186 self.data['noise'] = {}
1186 1187 self.data['rti'] = {}
1187 1188 if 'noise' not in self.plottypes:
1188 1189 self.plottypes.append('noise')
1189 1190 if 'rti' not in self.plottypes:
1190 1191 self.plottypes.append('rti')
1191 1192
1192 1193 def shape(self, key):
1193 1194 '''
1194 1195 Get the shape of the one-element data for the given key
1195 1196 '''
1196 1197
1197 1198 if len(self.data[key]):
1198 1199 if 'spc' in key or not self.buffering:
1199 1200 return self.data[key].shape
1200 1201 return self.data[key][self.__times[0]].shape
1201 1202 return (0,)
1202 1203
1203 1204 def update(self, dataOut, tm):
1204 1205 '''
1205 1206 Update data object with new dataOut
1206 1207 '''
1207
1208 1208 if tm in self.__times:
1209 1209 return
1210 1210 self.profileIndex = dataOut.profileIndex
1211 1211 self.tm = tm
1212 1212 self.type = dataOut.type
1213 1213 self.parameters = getattr(dataOut, 'parameters', [])
1214 1214
1215 1215 if hasattr(dataOut, 'meta'):
1216 1216 self.meta.update(dataOut.meta)
1217 1217
1218 1218 if hasattr(dataOut, 'pairsList'):
1219 1219 self.pairs = dataOut.pairsList
1220 1220
1221 1221 self.interval = dataOut.getTimeInterval()
1222 1222 self.localtime = dataOut.useLocalTime
1223 1223 if True in ['spc' in ptype for ptype in self.plottypes]:
1224 1224 self.xrange = (dataOut.getFreqRange(1)/1000.,
1225 1225 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1226 1226 self.factor = dataOut.normFactor
1227 1227 self.__heights.append(dataOut.heightList)
1228 1228 self.__all_heights.update(dataOut.heightList)
1229 1229 self.__times.append(tm)
1230
1231 1230 for plot in self.plottypes:
1232 1231 if plot in ('spc', 'spc_moments', 'spc_cut'):
1233 1232 z = dataOut.data_spc/dataOut.normFactor
1234 1233 buffer = 10*numpy.log10(z)
1235 1234 if plot == 'cspc':
1236 1235 z = dataOut.data_spc/dataOut.normFactor
1237 1236 buffer = (dataOut.data_spc, dataOut.data_cspc)
1238 1237 if plot == 'noise':
1239 1238 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1240 1239 if plot in ('rti', 'spcprofile'):
1241 1240 buffer = dataOut.getPower()
1242 1241 if plot == 'snr_db':
1243 1242 buffer = dataOut.data_SNR
1244 1243 if plot == 'snr':
1245 1244 buffer = 10*numpy.log10(dataOut.data_SNR)
1246 1245 if plot == 'dop':
1247 1246 buffer = dataOut.data_DOP
1248 1247 if plot == 'pow':
1249 1248 buffer = 10*numpy.log10(dataOut.data_POW)
1250 1249 if plot == 'width':
1251 1250 buffer = dataOut.data_WIDTH
1252 1251 if plot == 'coh':
1253 1252 buffer = dataOut.getCoherence()
1254 1253 if plot == 'phase':
1255 1254 buffer = dataOut.getCoherence(phase=True)
1256 1255 if plot == 'output':
1257 1256 buffer = dataOut.data_output
1258 1257 if plot == 'param':
1259 1258 buffer = dataOut.data_param
1260 1259 if plot == 'scope':
1261 1260 buffer = dataOut.data
1262 1261 self.flagDataAsBlock = dataOut.flagDataAsBlock
1263 1262 self.nProfiles = dataOut.nProfiles
1263 if plot == 'pp_power':
1264 buffer = dataOut.data_intensity
1265 self.flagDataAsBlock = dataOut.flagDataAsBlock
1266 self.nProfiles = dataOut.nProfiles
1267 if plot == 'pp_velocity':
1268 buffer = dataOut.data_velocity
1269 self.flagDataAsBlock = dataOut.flagDataAsBlock
1270 self.nProfiles = dataOut.nProfiles
1264 1271
1265 1272 if plot == 'spc':
1266 1273 self.data['spc'] = buffer
1267 1274 elif plot == 'cspc':
1268 1275 self.data['spc'] = buffer[0]
1269 1276 self.data['cspc'] = buffer[1]
1270 1277 elif plot == 'spc_moments':
1271 1278 self.data['spc'] = buffer
1272 1279 self.data['moments'][tm] = dataOut.moments
1273 1280 else:
1274 1281 if self.buffering:
1275 1282 self.data[plot][tm] = buffer
1276 1283 else:
1277 1284 self.data[plot] = buffer
1278 1285
1279 1286 if dataOut.channelList is None:
1280 1287 self.channels = range(buffer.shape[0])
1281 1288 else:
1282 1289 self.channels = dataOut.channelList
1283 1290
1284 1291 if buffer is None:
1285 1292 self.flagNoData = True
1286 1293 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1287 1294
1288 1295 def normalize_heights(self):
1289 1296 '''
1290 1297 Ensure same-dimension of the data for different heighList
1291 1298 '''
1292 1299
1293 1300 H = numpy.array(list(self.__all_heights))
1294 1301 H.sort()
1295 1302 for key in self.data:
1296 1303 shape = self.shape(key)[:-1] + H.shape
1297 1304 for tm, obj in list(self.data[key].items()):
1298 1305 h = self.__heights[self.__times.index(tm)]
1299 1306 if H.size == h.size:
1300 1307 continue
1301 1308 index = numpy.where(numpy.in1d(H, h))[0]
1302 1309 dummy = numpy.zeros(shape) + numpy.nan
1303 1310 if len(shape) == 2:
1304 1311 dummy[:, index] = obj
1305 1312 else:
1306 1313 dummy[index] = obj
1307 1314 self.data[key][tm] = dummy
1308 1315
1309 1316 self.__heights = [H for tm in self.__times]
1310 1317
1311 1318 def jsonify(self, plot_name, plot_type, decimate=False):
1312 1319 '''
1313 1320 Convert data to json
1314 1321 '''
1315 1322
1316 1323 tm = self.times[-1]
1317 1324 dy = int(self.heights.size/self.MAXNUMY) + 1
1318 1325 if self.key in ('spc', 'cspc') or not self.buffering:
1319 1326 dx = int(self.data[self.key].shape[1]/self.MAXNUMX) + 1
1320 1327 data = self.roundFloats(
1321 1328 self.data[self.key][::, ::dx, ::dy].tolist())
1322 1329 else:
1323 1330 if self.key is 'noise':
1324 1331 data = [[x] for x in self.roundFloats(self.data[self.key][tm].tolist())]
1325 1332 else:
1326 1333 data = self.roundFloats(self.data[self.key][tm][::, ::dy].tolist())
1327 1334
1328 1335 meta = {}
1329 1336 ret = {
1330 1337 'plot': plot_name,
1331 1338 'code': self.exp_code,
1332 1339 'time': float(tm),
1333 1340 'data': data,
1334 1341 }
1335 1342 meta['type'] = plot_type
1336 1343 meta['interval'] = float(self.interval)
1337 1344 meta['localtime'] = self.localtime
1338 1345 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1339 1346 if 'spc' in self.data or 'cspc' in self.data:
1340 1347 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1341 1348 else:
1342 1349 meta['xrange'] = []
1343 1350
1344 1351 meta.update(self.meta)
1345 1352 ret['metadata'] = meta
1346 1353 return json.dumps(ret)
1347 1354
1348 1355 @property
1349 1356 def times(self):
1350 1357 '''
1351 1358 Return the list of times of the current data
1352 1359 '''
1353 1360
1354 1361 ret = numpy.array(self.__times)
1355 1362 ret.sort()
1356 1363 return ret
1357 1364
1358 1365 @property
1359 1366 def min_time(self):
1360 1367 '''
1361 1368 Return the minimun time value
1362 1369 '''
1363 1370
1364 1371 return self.times[0]
1365 1372
1366 1373 @property
1367 1374 def max_time(self):
1368 1375 '''
1369 1376 Return the maximun time value
1370 1377 '''
1371 1378
1372 1379 return self.times[-1]
1373 1380
1374 1381 @property
1375 1382 def heights(self):
1376 1383 '''
1377 1384 Return the list of heights of the current data
1378 1385 '''
1379 1386
1380 1387 return numpy.array(self.__heights[-1])
1381 1388
1382 1389 @staticmethod
1383 1390 def roundFloats(obj):
1384 1391 if isinstance(obj, list):
1385 1392 return list(map(PlotterData.roundFloats, obj))
1386 1393 elif isinstance(obj, float):
1387 1394 return round(obj, 2)
@@ -1,134 +1,233
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9
10 10 from schainpy.model.graphics.jroplot_base import Plot, plt
11 11
12 12
13 13 class ScopePlot(Plot):
14 14
15 15 '''
16 16 Plot for Scope
17 17 '''
18 18
19 19 CODE = 'scope'
20 20 plot_name = 'Scope'
21 21 plot_type = 'scatter'
22 22
23 23 def setup(self):
24 24
25 25 self.xaxis = 'Range (Km)'
26 26 self.ncols = 1
27 27 self.nrows = 1
28 28 self.nplots = 1
29 29 self.ylabel = 'Intensity [dB]'
30 30 self.titles = ['Scope']
31 31 self.colorbar = False
32 32 self.width = 6
33 33 self.height = 4
34 34
35 35 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
36 36
37 37 yreal = y[channelIndexList,:].real
38 38 yimag = y[channelIndexList,:].imag
39 39 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
40 40 self.xlabel = "Range (Km)"
41 41 self.ylabel = "Intensity - IQ"
42 42
43 43 self.y = yreal
44 44 self.x = x
45 45 self.xmin = min(x)
46 46 self.xmax = max(x)
47 47
48 48
49 49 self.titles[0] = title
50 50
51 51 for i,ax in enumerate(self.axes):
52 52 title = "Channel %d" %(i)
53 53 if ax.firsttime:
54 54 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
55 55 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
56 56 else:
57 57 ax.plt_r.set_data(x, yreal[i,:])
58 58 ax.plt_i.set_data(x, yimag[i,:])
59 59
60 60 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
61 61 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
62 62 yreal = y.real
63 yreal = 10*numpy.log10(yreal)
63 64 self.y = yreal
64 65 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
65 66 self.xlabel = "Range (Km)"
66 67 self.ylabel = "Intensity"
67 68 self.xmin = min(x)
68 69 self.xmax = max(x)
69 70
70 71
71 72 self.titles[0] = title
72 73
73 74 for i,ax in enumerate(self.axes):
74 75 title = "Channel %d" %(i)
75 76
76 77 ychannel = yreal[i,:]
77 78
78 79 if ax.firsttime:
79 80 ax.plt_r = ax.plot(x, ychannel)[0]
80 81 else:
81 82 #pass
82 83 ax.plt_r.set_data(x, ychannel)
83 84
85 def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle):
84 86
85 def plot(self):
86 87
88 y = y[channelIndexList,:]
89 yreal = y.real
90 yreal = 10*numpy.log10(yreal)
91 self.y = yreal
92 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
93 self.xlabel = "Range (Km)"
94 self.ylabel = "Intensity"
95 self.xmin = min(x)
96 self.xmax = max(x)
97
98 self.titles[0] =title
99 for i,ax in enumerate(self.axes):
100 title = "Channel %d" %(i)
101
102 ychannel = yreal[i,:]
103
104 if ax.firsttime:
105 ax.plt_r = ax.plot(x, ychannel)[0]
106 else:
107 #pass
108 ax.plt_r.set_data(x, ychannel)
109
110 def plot_weathervelocity(self, x, y, channelIndexList, thisDatetime, wintitle):
111
112 x = x[channelIndexList,:]
113 yreal = y
114 self.y = yreal
115 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
116 self.xlabel = "Velocity (m/s)"
117 self.ylabel = "Range (Km)"
118 self.xmin = numpy.min(x)
119 self.xmax = numpy.max(x)
120 self.titles[0] =title
121 for i,ax in enumerate(self.axes):
122 title = "Channel %d" %(i)
123 xchannel = x[i,:]
124 if ax.firsttime:
125 ax.plt_r = ax.plot(xchannel, yreal)[0]
126 else:
127 #pass
128 ax.plt_r.set_data(xchannel, yreal)
129
130 def plot(self):
87 131 if self.channels:
88 132 channels = self.channels
89 133 else:
90 134 channels = self.data.channels
91 135
92 136 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
93
94 scope = self.data['scope']
137 if self.CODE == "pp_power":
138 scope = self.data['pp_power']
139 elif self.CODE == "pp_velocity":
140 scope = self.data["pp_velocity"]
141 else:
142 scope =self.data["scope"]
95 143
96 144 if self.data.flagDataAsBlock:
97 145
98 146 for i in range(self.data.nProfiles):
99 147
100 148 wintitle1 = " [Profile = %d] " %i
101
149 if self.CODE =="scope":
102 150 if self.type == "power":
103 151 self.plot_power(self.data.heights,
104 152 scope[:,i,:],
105 153 channels,
106 154 thisDatetime,
107 155 wintitle1
108 156 )
109 157
110 158 if self.type == "iq":
111 159 self.plot_iq(self.data.heights,
112 160 scope[:,i,:],
113 161 channels,
114 162 thisDatetime,
115 163 wintitle1
116 164 )
165 if self.CODE=="pp_power":
166 self.plot_weatherpower(self.data.heights,
167 scope[:,i,:],
168 channels,
169 thisDatetime,
170 wintitle
171 )
172 if self.CODE=="pp_velocity":
173 self.plot_weathervelocity(scope[:,i,:],
174 self.data.heights,
175 channels,
176 thisDatetime,
177 wintitle
178 )
117 179 else:
118 180 wintitle = " [Profile = %d] " %self.data.profileIndex
119
181 if self.CODE== "scope":
120 182 if self.type == "power":
121 183 self.plot_power(self.data.heights,
122 184 scope,
123 185 channels,
124 186 thisDatetime,
125 187 wintitle
126 188 )
127 189
128 190 if self.type == "iq":
129 191 self.plot_iq(self.data.heights,
130 192 scope,
131 193 channels,
132 194 thisDatetime,
133 195 wintitle
134 196 )
197 if self.CODE=="pp_power":
198 self.plot_weatherpower(self.data.heights,
199 scope,
200 channels,
201 thisDatetime,
202 wintitle
203 )
204 if self.CODE=="pp_velocity":
205 self.plot_weathervelocity(scope,
206 self.data.heights,
207 channels,
208 thisDatetime,
209 wintitle
210 )
211
212
213
214 class PulsepairPowerPlot(ScopePlot):
215 '''
216 Plot for
217 '''
218
219 CODE = 'pp_power'
220 plot_name = 'PulsepairPower'
221 plot_type = 'scatter'
222 buffering = False
223
224
225
226 class PulsepairVelocityPlot(ScopePlot):
227 '''
228 Plot for
229 '''
230 CODE = 'pp_velocity'
231 plot_name = 'PulsepairVelocity'
232 plot_type = 'scatter'
233 buffering = False
@@ -1,1409 +1,1557
1 1 import sys
2 import numpy
2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62
63 62 for channel in channelList:
64 63 if channel not in self.dataOut.channelList:
65 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
66 65
67 66 index = self.dataOut.channelList.index(channel)
68 67 channelIndexList.append(index)
69
70 68 self.selectChannelsByIndex(channelIndexList)
71 69 return self.dataOut
72 70
73 71 def selectChannelsByIndex(self, channelIndexList):
74 72 """
75 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
76 74
77 75 Input:
78 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
79 77
80 78 Affected:
81 79 self.dataOut.data
82 80 self.dataOut.channelIndexList
83 81 self.dataOut.nChannels
84 82 self.dataOut.m_ProcessingHeader.totalSpectra
85 83 self.dataOut.systemHeaderObj.numChannels
86 84 self.dataOut.m_ProcessingHeader.blockSize
87 85
88 86 Return:
89 87 None
90 88 """
91 89
92 90 for channelIndex in channelIndexList:
93 91 if channelIndex not in self.dataOut.channelIndexList:
94 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
95 93
96 94 if self.dataOut.type == 'Voltage':
97 95 if self.dataOut.flagDataAsBlock:
98 96 """
99 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
100 98 """
101 99 data = self.dataOut.data[channelIndexList,:,:]
102 100 else:
103 101 data = self.dataOut.data[channelIndexList,:]
104 102
105 103 self.dataOut.data = data
106 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
107 105 self.dataOut.channelList = range(len(channelIndexList))
106
108 107 elif self.dataOut.type == 'Spectra':
109 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
110 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
111 110
112 111 self.dataOut.data_spc = data_spc
113 112 self.dataOut.data_dc = data_dc
114 113
115 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
116 115 self.dataOut.channelList = range(len(channelIndexList))
117 116 self.__selectPairsByChannel(channelIndexList)
118 117
119 118 return 1
120 119
121 120 def __selectPairsByChannel(self, channelList=None):
122 121
123 122 if channelList == None:
124 123 return
125 124
126 125 pairsIndexListSelected = []
127 126 for pairIndex in self.dataOut.pairsIndexList:
128 127 # First pair
129 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
130 129 continue
131 130 # Second pair
132 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
133 132 continue
134 133
135 134 pairsIndexListSelected.append(pairIndex)
136 135
137 136 if not pairsIndexListSelected:
138 137 self.dataOut.data_cspc = None
139 138 self.dataOut.pairsList = []
140 139 return
141 140
142 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
143 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
144 143 for i in pairsIndexListSelected]
145 144
146 145 return
147 146
148 147 class selectHeights(Operation):
149 148
150 149 def run(self, dataOut, minHei=None, maxHei=None):
151 150 """
152 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
153 152 minHei <= height <= maxHei
154 153
155 154 Input:
156 155 minHei : valor minimo de altura a considerar
157 156 maxHei : valor maximo de altura a considerar
158 157
159 158 Affected:
160 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
161 160
162 161 Return:
163 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
164 163 """
165 164
166 165 self.dataOut = dataOut
167 166
168 167 if minHei == None:
169 168 minHei = self.dataOut.heightList[0]
170 169
171 170 if maxHei == None:
172 171 maxHei = self.dataOut.heightList[-1]
173 172
174 173 if (minHei < self.dataOut.heightList[0]):
175 174 minHei = self.dataOut.heightList[0]
176 175
177 176 if (maxHei > self.dataOut.heightList[-1]):
178 177 maxHei = self.dataOut.heightList[-1]
179 178
180 179 minIndex = 0
181 180 maxIndex = 0
182 181 heights = self.dataOut.heightList
183 182
184 183 inda = numpy.where(heights >= minHei)
185 184 indb = numpy.where(heights <= maxHei)
186 185
187 186 try:
188 187 minIndex = inda[0][0]
189 188 except:
190 189 minIndex = 0
191 190
192 191 try:
193 192 maxIndex = indb[0][-1]
194 193 except:
195 194 maxIndex = len(heights)
196 195
197 196 self.selectHeightsByIndex(minIndex, maxIndex)
198 197
199 198 return self.dataOut
200 199
201 200 def selectHeightsByIndex(self, minIndex, maxIndex):
202 201 """
203 202 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
204 203 minIndex <= index <= maxIndex
205 204
206 205 Input:
207 206 minIndex : valor de indice minimo de altura a considerar
208 207 maxIndex : valor de indice maximo de altura a considerar
209 208
210 209 Affected:
211 210 self.dataOut.data
212 211 self.dataOut.heightList
213 212
214 213 Return:
215 214 1 si el metodo se ejecuto con exito caso contrario devuelve 0
216 215 """
217 216
218 217 if self.dataOut.type == 'Voltage':
219 218 if (minIndex < 0) or (minIndex > maxIndex):
220 219 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
221 220
222 221 if (maxIndex >= self.dataOut.nHeights):
223 222 maxIndex = self.dataOut.nHeights
224 223
225 224 #voltage
226 225 if self.dataOut.flagDataAsBlock:
227 226 """
228 227 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
229 228 """
230 229 data = self.dataOut.data[:,:, minIndex:maxIndex]
231 230 else:
232 231 data = self.dataOut.data[:, minIndex:maxIndex]
233 232
234 233 # firstHeight = self.dataOut.heightList[minIndex]
235 234
236 235 self.dataOut.data = data
237 236 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
238 237
239 238 if self.dataOut.nHeights <= 1:
240 239 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
241 240 elif self.dataOut.type == 'Spectra':
242 241 if (minIndex < 0) or (minIndex > maxIndex):
243 242 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
244 243 minIndex, maxIndex))
245 244
246 245 if (maxIndex >= self.dataOut.nHeights):
247 246 maxIndex = self.dataOut.nHeights - 1
248 247
249 248 # Spectra
250 249 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
251 250
252 251 data_cspc = None
253 252 if self.dataOut.data_cspc is not None:
254 253 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
255 254
256 255 data_dc = None
257 256 if self.dataOut.data_dc is not None:
258 257 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
259 258
260 259 self.dataOut.data_spc = data_spc
261 260 self.dataOut.data_cspc = data_cspc
262 261 self.dataOut.data_dc = data_dc
263 262
264 263 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
265 264
266 265 return 1
267 266
268 267
269 268 class filterByHeights(Operation):
270 269
271 270 def run(self, dataOut, window):
272 271
273 272 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
274 273
275 274 if window == None:
276 275 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
277 276
278 277 newdelta = deltaHeight * window
279 278 r = dataOut.nHeights % window
280 279 newheights = (dataOut.nHeights-r)/window
281 280
282 281 if newheights <= 1:
283 282 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
284 283
285 284 if dataOut.flagDataAsBlock:
286 285 """
287 286 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
288 287 """
289 288 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
290 289 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
291 290 buffer = numpy.sum(buffer,3)
292 291
293 292 else:
294 293 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
295 294 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
296 295 buffer = numpy.sum(buffer,2)
297 296
298 297 dataOut.data = buffer
299 298 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
300 299 dataOut.windowOfFilter = window
301 300
302 301 return dataOut
303 302
304 303
305 304 class setH0(Operation):
306 305
307 306 def run(self, dataOut, h0, deltaHeight = None):
308 307
309 308 if not deltaHeight:
310 309 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
311 310
312 311 nHeights = dataOut.nHeights
313 312
314 313 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
315 314
316 315 dataOut.heightList = newHeiRange
317 316
318 317 return dataOut
319 318
320 319
321 320 class deFlip(Operation):
322 321
323 322 def run(self, dataOut, channelList = []):
324 323
325 324 data = dataOut.data.copy()
326 325
327 326 if dataOut.flagDataAsBlock:
328 327 flip = self.flip
329 328 profileList = list(range(dataOut.nProfiles))
330 329
331 330 if not channelList:
332 331 for thisProfile in profileList:
333 332 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
334 333 flip *= -1.0
335 334 else:
336 335 for thisChannel in channelList:
337 336 if thisChannel not in dataOut.channelList:
338 337 continue
339 338
340 339 for thisProfile in profileList:
341 340 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
342 341 flip *= -1.0
343 342
344 343 self.flip = flip
345 344
346 345 else:
347 346 if not channelList:
348 347 data[:,:] = data[:,:]*self.flip
349 348 else:
350 349 for thisChannel in channelList:
351 350 if thisChannel not in dataOut.channelList:
352 351 continue
353 352
354 353 data[thisChannel,:] = data[thisChannel,:]*self.flip
355 354
356 355 self.flip *= -1.
357 356
358 357 dataOut.data = data
359 358
360 359 return dataOut
361 360
362 361
363 362 class setAttribute(Operation):
364 363 '''
365 364 Set an arbitrary attribute(s) to dataOut
366 365 '''
367 366
368 367 def __init__(self):
369 368
370 369 Operation.__init__(self)
371 370 self._ready = False
372 371
373 372 def run(self, dataOut, **kwargs):
374 373
375 374 for key, value in kwargs.items():
376 375 setattr(dataOut, key, value)
377 376
378 377 return dataOut
379 378
380 379
381 380 class interpolateHeights(Operation):
382 381
383 382 def run(self, dataOut, topLim, botLim):
384 383 #69 al 72 para julia
385 384 #82-84 para meteoros
386 385 if len(numpy.shape(dataOut.data))==2:
387 386 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
388 387 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
389 388 #dataOut.data[:,botLim:limSup+1] = sampInterp
390 389 dataOut.data[:,botLim:topLim+1] = sampInterp
391 390 else:
392 391 nHeights = dataOut.data.shape[2]
393 392 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
394 393 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
395 394 f = interpolate.interp1d(x, y, axis = 2)
396 395 xnew = numpy.arange(botLim,topLim+1)
397 396 ynew = f(xnew)
398 397 dataOut.data[:,:,botLim:topLim+1] = ynew
399 398
400 399 return dataOut
401 400
402 401
403 402 class CohInt(Operation):
404 403
405 404 isConfig = False
406 405 __profIndex = 0
407 406 __byTime = False
408 407 __initime = None
409 408 __lastdatatime = None
410 409 __integrationtime = None
411 410 __buffer = None
412 411 __bufferStride = []
413 412 __dataReady = False
414 413 __profIndexStride = 0
415 414 __dataToPutStride = False
416 415 n = None
417 416
418 417 def __init__(self, **kwargs):
419 418
420 419 Operation.__init__(self, **kwargs)
421 420
422 421 # self.isConfig = False
423 422
424 423 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
425 424 """
426 425 Set the parameters of the integration class.
427 426
428 427 Inputs:
429 428
430 429 n : Number of coherent integrations
431 430 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
432 431 overlapping :
433 432 """
434 433
435 434 self.__initime = None
436 435 self.__lastdatatime = 0
437 436 self.__buffer = None
438 437 self.__dataReady = False
439 438 self.byblock = byblock
440 439 self.stride = stride
441 440
442 441 if n == None and timeInterval == None:
443 442 raise ValueError("n or timeInterval should be specified ...")
444 443
445 444 if n != None:
446 445 self.n = n
447 446 self.__byTime = False
448 447 else:
449 448 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
450 449 self.n = 9999
451 450 self.__byTime = True
452 451
453 452 if overlapping:
454 453 self.__withOverlapping = True
455 454 self.__buffer = None
456 455 else:
457 456 self.__withOverlapping = False
458 457 self.__buffer = 0
459 458
460 459 self.__profIndex = 0
461 460
462 461 def putData(self, data):
463 462
464 463 """
465 464 Add a profile to the __buffer and increase in one the __profileIndex
466 465
467 466 """
468 467
469 468 if not self.__withOverlapping:
470 469 self.__buffer += data.copy()
471 470 self.__profIndex += 1
472 471 return
473 472
474 473 #Overlapping data
475 474 nChannels, nHeis = data.shape
476 475 data = numpy.reshape(data, (1, nChannels, nHeis))
477 476
478 477 #If the buffer is empty then it takes the data value
479 478 if self.__buffer is None:
480 479 self.__buffer = data
481 480 self.__profIndex += 1
482 481 return
483 482
484 483 #If the buffer length is lower than n then stakcing the data value
485 484 if self.__profIndex < self.n:
486 485 self.__buffer = numpy.vstack((self.__buffer, data))
487 486 self.__profIndex += 1
488 487 return
489 488
490 489 #If the buffer length is equal to n then replacing the last buffer value with the data value
491 490 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
492 491 self.__buffer[self.n-1] = data
493 492 self.__profIndex = self.n
494 493 return
495 494
496 495
497 496 def pushData(self):
498 497 """
499 498 Return the sum of the last profiles and the profiles used in the sum.
500 499
501 500 Affected:
502 501
503 502 self.__profileIndex
504 503
505 504 """
506 505
507 506 if not self.__withOverlapping:
508 507 data = self.__buffer
509 508 n = self.__profIndex
510 509
511 510 self.__buffer = 0
512 511 self.__profIndex = 0
513 512
514 513 return data, n
515 514
516 515 #Integration with Overlapping
517 516 data = numpy.sum(self.__buffer, axis=0)
518 517 # print data
519 518 # raise
520 519 n = self.__profIndex
521 520
522 521 return data, n
523 522
524 523 def byProfiles(self, data):
525 524
526 525 self.__dataReady = False
527 526 avgdata = None
528 527 # n = None
529 528 # print data
530 529 # raise
531 530 self.putData(data)
532 531
533 532 if self.__profIndex == self.n:
534 533 avgdata, n = self.pushData()
535 534 self.__dataReady = True
536 535
537 536 return avgdata
538 537
539 538 def byTime(self, data, datatime):
540 539
541 540 self.__dataReady = False
542 541 avgdata = None
543 542 n = None
544 543
545 544 self.putData(data)
546 545
547 546 if (datatime - self.__initime) >= self.__integrationtime:
548 547 avgdata, n = self.pushData()
549 548 self.n = n
550 549 self.__dataReady = True
551 550
552 551 return avgdata
553 552
554 553 def integrateByStride(self, data, datatime):
555 554 # print data
556 555 if self.__profIndex == 0:
557 556 self.__buffer = [[data.copy(), datatime]]
558 557 else:
559 558 self.__buffer.append([data.copy(),datatime])
560 559 self.__profIndex += 1
561 560 self.__dataReady = False
562 561
563 562 if self.__profIndex == self.n * self.stride :
564 563 self.__dataToPutStride = True
565 564 self.__profIndexStride = 0
566 565 self.__profIndex = 0
567 566 self.__bufferStride = []
568 567 for i in range(self.stride):
569 568 current = self.__buffer[i::self.stride]
570 569 data = numpy.sum([t[0] for t in current], axis=0)
571 570 avgdatatime = numpy.average([t[1] for t in current])
572 571 # print data
573 572 self.__bufferStride.append((data, avgdatatime))
574 573
575 574 if self.__dataToPutStride:
576 575 self.__dataReady = True
577 576 self.__profIndexStride += 1
578 577 if self.__profIndexStride == self.stride:
579 578 self.__dataToPutStride = False
580 579 # print self.__bufferStride[self.__profIndexStride - 1]
581 580 # raise
582 581 return self.__bufferStride[self.__profIndexStride - 1]
583 582
584 583
585 584 return None, None
586 585
587 586 def integrate(self, data, datatime=None):
588 587
589 588 if self.__initime == None:
590 589 self.__initime = datatime
591 590
592 591 if self.__byTime:
593 592 avgdata = self.byTime(data, datatime)
594 593 else:
595 594 avgdata = self.byProfiles(data)
596 595
597 596
598 597 self.__lastdatatime = datatime
599 598
600 599 if avgdata is None:
601 600 return None, None
602 601
603 602 avgdatatime = self.__initime
604 603
605 604 deltatime = datatime - self.__lastdatatime
606 605
607 606 if not self.__withOverlapping:
608 607 self.__initime = datatime
609 608 else:
610 609 self.__initime += deltatime
611 610
612 611 return avgdata, avgdatatime
613 612
614 613 def integrateByBlock(self, dataOut):
615 614
616 615 times = int(dataOut.data.shape[1]/self.n)
617 616 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
618 617
619 618 id_min = 0
620 619 id_max = self.n
621 620
622 621 for i in range(times):
623 622 junk = dataOut.data[:,id_min:id_max,:]
624 623 avgdata[:,i,:] = junk.sum(axis=1)
625 624 id_min += self.n
626 625 id_max += self.n
627 626
628 627 timeInterval = dataOut.ippSeconds*self.n
629 628 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
630 629 self.__dataReady = True
631 630 return avgdata, avgdatatime
632 631
633 632 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
634 633
635 634 if not self.isConfig:
636 635 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
637 636 self.isConfig = True
638 637
639 638 if dataOut.flagDataAsBlock:
640 639 """
641 640 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
642 641 """
643 642 avgdata, avgdatatime = self.integrateByBlock(dataOut)
644 643 dataOut.nProfiles /= self.n
645 644 else:
646 645 if stride is None:
647 646 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
648 647 else:
649 648 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
650 649
651 650
652 651 # dataOut.timeInterval *= n
653 652 dataOut.flagNoData = True
654 653
655 654 if self.__dataReady:
656 655 dataOut.data = avgdata
657 656 dataOut.nCohInt *= self.n
658 657 dataOut.utctime = avgdatatime
659 658 # print avgdata, avgdatatime
660 659 # raise
661 660 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
662 661 dataOut.flagNoData = False
663 662 return dataOut
664 663
665 664 class Decoder(Operation):
666 665
667 666 isConfig = False
668 667 __profIndex = 0
669 668
670 669 code = None
671 670
672 671 nCode = None
673 672 nBaud = None
674 673
675 674 def __init__(self, **kwargs):
676 675
677 676 Operation.__init__(self, **kwargs)
678 677
679 678 self.times = None
680 679 self.osamp = None
681 680 # self.__setValues = False
682 681 self.isConfig = False
683 682 self.setupReq = False
684 683 def setup(self, code, osamp, dataOut):
685 684
686 685 self.__profIndex = 0
687 686
688 687 self.code = code
689 688
690 689 self.nCode = len(code)
691 690 self.nBaud = len(code[0])
692 691
693 692 if (osamp != None) and (osamp >1):
694 693 self.osamp = osamp
695 694 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
696 695 self.nBaud = self.nBaud*self.osamp
697 696
698 697 self.__nChannels = dataOut.nChannels
699 698 self.__nProfiles = dataOut.nProfiles
700 699 self.__nHeis = dataOut.nHeights
701 700
702 701 if self.__nHeis < self.nBaud:
703 702 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
704 703
705 704 #Frequency
706 705 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
707 706
708 707 __codeBuffer[:,0:self.nBaud] = self.code
709 708
710 709 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
711 710
712 711 if dataOut.flagDataAsBlock:
713 712
714 713 self.ndatadec = self.__nHeis #- self.nBaud + 1
715 714
716 715 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
717 716
718 717 else:
719 718
720 719 #Time
721 720 self.ndatadec = self.__nHeis #- self.nBaud + 1
722 721
723 722 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
724 723
725 724 def __convolutionInFreq(self, data):
726 725
727 726 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
728 727
729 728 fft_data = numpy.fft.fft(data, axis=1)
730 729
731 730 conv = fft_data*fft_code
732 731
733 732 data = numpy.fft.ifft(conv,axis=1)
734 733
735 734 return data
736 735
737 736 def __convolutionInFreqOpt(self, data):
738 737
739 738 raise NotImplementedError
740 739
741 740 def __convolutionInTime(self, data):
742 741
743 742 code = self.code[self.__profIndex]
744 743 for i in range(self.__nChannels):
745 744 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
746 745
747 746 return self.datadecTime
748 747
749 748 def __convolutionByBlockInTime(self, data):
750 749
751 750 repetitions = int(self.__nProfiles / self.nCode)
752 751 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
753 752 junk = junk.flatten()
754 753 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
755 754 profilesList = range(self.__nProfiles)
756 755
757 756 for i in range(self.__nChannels):
758 757 for j in profilesList:
759 758 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
760 759 return self.datadecTime
761 760
762 761 def __convolutionByBlockInFreq(self, data):
763 762
764 763 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
765 764
766 765
767 766 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
768 767
769 768 fft_data = numpy.fft.fft(data, axis=2)
770 769
771 770 conv = fft_data*fft_code
772 771
773 772 data = numpy.fft.ifft(conv,axis=2)
774 773
775 774 return data
776 775
777 776
778 777 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
779 778
780 779 if dataOut.flagDecodeData:
781 780 print("This data is already decoded, recoding again ...")
782 781
783 782 if not self.isConfig:
784 783
785 784 if code is None:
786 785 if dataOut.code is None:
787 786 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
788 787
789 788 code = dataOut.code
790 789 else:
791 790 code = numpy.array(code).reshape(nCode,nBaud)
792 791 self.setup(code, osamp, dataOut)
793 792
794 793 self.isConfig = True
795 794
796 795 if mode == 3:
797 796 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
798 797
799 798 if times != None:
800 799 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
801 800
802 801 if self.code is None:
803 802 print("Fail decoding: Code is not defined.")
804 803 return
805 804
806 805 self.__nProfiles = dataOut.nProfiles
807 806 datadec = None
808 807
809 808 if mode == 3:
810 809 mode = 0
811 810
812 811 if dataOut.flagDataAsBlock:
813 812 """
814 813 Decoding when data have been read as block,
815 814 """
816 815
817 816 if mode == 0:
818 817 datadec = self.__convolutionByBlockInTime(dataOut.data)
819 818 if mode == 1:
820 819 datadec = self.__convolutionByBlockInFreq(dataOut.data)
821 820 else:
822 821 """
823 822 Decoding when data have been read profile by profile
824 823 """
825 824 if mode == 0:
826 825 datadec = self.__convolutionInTime(dataOut.data)
827 826
828 827 if mode == 1:
829 828 datadec = self.__convolutionInFreq(dataOut.data)
830 829
831 830 if mode == 2:
832 831 datadec = self.__convolutionInFreqOpt(dataOut.data)
833 832
834 833 if datadec is None:
835 834 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
836 835
837 836 dataOut.code = self.code
838 837 dataOut.nCode = self.nCode
839 838 dataOut.nBaud = self.nBaud
840 839
841 840 dataOut.data = datadec
842 841
843 842 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
844 843
845 844 dataOut.flagDecodeData = True #asumo q la data esta decodificada
846 845
847 846 if self.__profIndex == self.nCode-1:
848 847 self.__profIndex = 0
849 848 return dataOut
850 849
851 850 self.__profIndex += 1
852 851
853 852 return dataOut
854 853 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
855 854
856 855
857 856 class ProfileConcat(Operation):
858 857
859 858 isConfig = False
860 859 buffer = None
861 860
862 861 def __init__(self, **kwargs):
863 862
864 863 Operation.__init__(self, **kwargs)
865 864 self.profileIndex = 0
866 865
867 866 def reset(self):
868 867 self.buffer = numpy.zeros_like(self.buffer)
869 868 self.start_index = 0
870 869 self.times = 1
871 870
872 871 def setup(self, data, m, n=1):
873 872 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
874 873 self.nHeights = data.shape[1]#.nHeights
875 874 self.start_index = 0
876 875 self.times = 1
877 876
878 877 def concat(self, data):
879 878
880 879 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
881 880 self.start_index = self.start_index + self.nHeights
882 881
883 882 def run(self, dataOut, m):
884 883 dataOut.flagNoData = True
885 884
886 885 if not self.isConfig:
887 886 self.setup(dataOut.data, m, 1)
888 887 self.isConfig = True
889 888
890 889 if dataOut.flagDataAsBlock:
891 890 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
892 891
893 892 else:
894 893 self.concat(dataOut.data)
895 894 self.times += 1
896 895 if self.times > m:
897 896 dataOut.data = self.buffer
898 897 self.reset()
899 898 dataOut.flagNoData = False
900 899 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
901 900 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
902 901 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
903 902 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
904 903 dataOut.ippSeconds *= m
905 904 return dataOut
906 905
907 906 class ProfileSelector(Operation):
908 907
909 908 profileIndex = None
910 909 # Tamanho total de los perfiles
911 910 nProfiles = None
912 911
913 912 def __init__(self, **kwargs):
914 913
915 914 Operation.__init__(self, **kwargs)
916 915 self.profileIndex = 0
917 916
918 917 def incProfileIndex(self):
919 918
920 919 self.profileIndex += 1
921 920
922 921 if self.profileIndex >= self.nProfiles:
923 922 self.profileIndex = 0
924 923
925 924 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
926 925
927 926 if profileIndex < minIndex:
928 927 return False
929 928
930 929 if profileIndex > maxIndex:
931 930 return False
932 931
933 932 return True
934 933
935 934 def isThisProfileInList(self, profileIndex, profileList):
936 935
937 936 if profileIndex not in profileList:
938 937 return False
939 938
940 939 return True
941 940
942 941 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
943 942
944 943 """
945 944 ProfileSelector:
946 945
947 946 Inputs:
948 947 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
949 948
950 949 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
951 950
952 951 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
953 952
954 953 """
955 954
956 955 if rangeList is not None:
957 956 if type(rangeList[0]) not in (tuple, list):
958 957 rangeList = [rangeList]
959 958
960 959 dataOut.flagNoData = True
961 960
962 961 if dataOut.flagDataAsBlock:
963 962 """
964 963 data dimension = [nChannels, nProfiles, nHeis]
965 964 """
966 965 if profileList != None:
967 966 dataOut.data = dataOut.data[:,profileList,:]
968 967
969 968 if profileRangeList != None:
970 969 minIndex = profileRangeList[0]
971 970 maxIndex = profileRangeList[1]
972 971 profileList = list(range(minIndex, maxIndex+1))
973 972
974 973 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
975 974
976 975 if rangeList != None:
977 976
978 977 profileList = []
979 978
980 979 for thisRange in rangeList:
981 980 minIndex = thisRange[0]
982 981 maxIndex = thisRange[1]
983 982
984 983 profileList.extend(list(range(minIndex, maxIndex+1)))
985 984
986 985 dataOut.data = dataOut.data[:,profileList,:]
987 986
988 987 dataOut.nProfiles = len(profileList)
989 988 dataOut.profileIndex = dataOut.nProfiles - 1
990 989 dataOut.flagNoData = False
991 990
992 991 return dataOut
993 992
994 993 """
995 994 data dimension = [nChannels, nHeis]
996 995 """
997 996
998 997 if profileList != None:
999 998
1000 999 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1001 1000
1002 1001 self.nProfiles = len(profileList)
1003 1002 dataOut.nProfiles = self.nProfiles
1004 1003 dataOut.profileIndex = self.profileIndex
1005 1004 dataOut.flagNoData = False
1006 1005
1007 1006 self.incProfileIndex()
1008 1007 return dataOut
1009 1008
1010 1009 if profileRangeList != None:
1011 1010
1012 1011 minIndex = profileRangeList[0]
1013 1012 maxIndex = profileRangeList[1]
1014 1013
1015 1014 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1016 1015
1017 1016 self.nProfiles = maxIndex - minIndex + 1
1018 1017 dataOut.nProfiles = self.nProfiles
1019 1018 dataOut.profileIndex = self.profileIndex
1020 1019 dataOut.flagNoData = False
1021 1020
1022 1021 self.incProfileIndex()
1023 1022 return dataOut
1024 1023
1025 1024 if rangeList != None:
1026 1025
1027 1026 nProfiles = 0
1028 1027
1029 1028 for thisRange in rangeList:
1030 1029 minIndex = thisRange[0]
1031 1030 maxIndex = thisRange[1]
1032 1031
1033 1032 nProfiles += maxIndex - minIndex + 1
1034 1033
1035 1034 for thisRange in rangeList:
1036 1035
1037 1036 minIndex = thisRange[0]
1038 1037 maxIndex = thisRange[1]
1039 1038
1040 1039 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1041 1040
1042 1041 self.nProfiles = nProfiles
1043 1042 dataOut.nProfiles = self.nProfiles
1044 1043 dataOut.profileIndex = self.profileIndex
1045 1044 dataOut.flagNoData = False
1046 1045
1047 1046 self.incProfileIndex()
1048 1047
1049 1048 break
1050 1049
1051 1050 return dataOut
1052 1051
1053 1052
1054 1053 if beam != None: #beam is only for AMISR data
1055 1054 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1056 1055 dataOut.flagNoData = False
1057 1056 dataOut.profileIndex = self.profileIndex
1058 1057
1059 1058 self.incProfileIndex()
1060 1059
1061 1060 return dataOut
1062 1061
1063 1062 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1064 1063
1065 1064
1066 1065 class Reshaper(Operation):
1067 1066
1068 1067 def __init__(self, **kwargs):
1069 1068
1070 1069 Operation.__init__(self, **kwargs)
1071 1070
1072 1071 self.__buffer = None
1073 1072 self.__nitems = 0
1074 1073
1075 1074 def __appendProfile(self, dataOut, nTxs):
1076 1075
1077 1076 if self.__buffer is None:
1078 1077 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1079 1078 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1080 1079
1081 1080 ini = dataOut.nHeights * self.__nitems
1082 1081 end = ini + dataOut.nHeights
1083 1082
1084 1083 self.__buffer[:, ini:end] = dataOut.data
1085 1084
1086 1085 self.__nitems += 1
1087 1086
1088 1087 return int(self.__nitems*nTxs)
1089 1088
1090 1089 def __getBuffer(self):
1091 1090
1092 1091 if self.__nitems == int(1./self.__nTxs):
1093 1092
1094 1093 self.__nitems = 0
1095 1094
1096 1095 return self.__buffer.copy()
1097 1096
1098 1097 return None
1099 1098
1100 1099 def __checkInputs(self, dataOut, shape, nTxs):
1101 1100
1102 1101 if shape is None and nTxs is None:
1103 1102 raise ValueError("Reshaper: shape of factor should be defined")
1104 1103
1105 1104 if nTxs:
1106 1105 if nTxs < 0:
1107 1106 raise ValueError("nTxs should be greater than 0")
1108 1107
1109 1108 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1110 1109 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1111 1110
1112 1111 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1113 1112
1114 1113 return shape, nTxs
1115 1114
1116 1115 if len(shape) != 2 and len(shape) != 3:
1117 1116 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1118 1117
1119 1118 if len(shape) == 2:
1120 1119 shape_tuple = [dataOut.nChannels]
1121 1120 shape_tuple.extend(shape)
1122 1121 else:
1123 1122 shape_tuple = list(shape)
1124 1123
1125 1124 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1126 1125
1127 1126 return shape_tuple, nTxs
1128 1127
1129 1128 def run(self, dataOut, shape=None, nTxs=None):
1130 1129
1131 1130 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1132 1131
1133 1132 dataOut.flagNoData = True
1134 1133 profileIndex = None
1135 1134
1136 1135 if dataOut.flagDataAsBlock:
1137 1136
1138 1137 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1139 1138 dataOut.flagNoData = False
1140 1139
1141 1140 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1142 1141
1143 1142 else:
1144 1143
1145 1144 if self.__nTxs < 1:
1146 1145
1147 1146 self.__appendProfile(dataOut, self.__nTxs)
1148 1147 new_data = self.__getBuffer()
1149 1148
1150 1149 if new_data is not None:
1151 1150 dataOut.data = new_data
1152 1151 dataOut.flagNoData = False
1153 1152
1154 1153 profileIndex = dataOut.profileIndex*nTxs
1155 1154
1156 1155 else:
1157 1156 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1158 1157
1159 1158 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1160 1159
1161 1160 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1162 1161
1163 1162 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1164 1163
1165 1164 dataOut.profileIndex = profileIndex
1166 1165
1167 1166 dataOut.ippSeconds /= self.__nTxs
1168 1167
1169 1168 return dataOut
1170 1169
1171 1170 class SplitProfiles(Operation):
1172 1171
1173 1172 def __init__(self, **kwargs):
1174 1173
1175 1174 Operation.__init__(self, **kwargs)
1176 1175
1177 1176 def run(self, dataOut, n):
1178 1177
1179 1178 dataOut.flagNoData = True
1180 1179 profileIndex = None
1181 1180
1182 1181 if dataOut.flagDataAsBlock:
1183 1182
1184 1183 #nchannels, nprofiles, nsamples
1185 1184 shape = dataOut.data.shape
1186 1185
1187 1186 if shape[2] % n != 0:
1188 1187 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1189 1188
1190 1189 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1191 1190
1192 1191 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1193 1192 dataOut.flagNoData = False
1194 1193
1195 1194 profileIndex = int(dataOut.nProfiles/n) - 1
1196 1195
1197 1196 else:
1198 1197
1199 1198 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1200 1199
1201 1200 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1202 1201
1203 1202 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1204 1203
1205 1204 dataOut.nProfiles = int(dataOut.nProfiles*n)
1206 1205
1207 1206 dataOut.profileIndex = profileIndex
1208 1207
1209 1208 dataOut.ippSeconds /= n
1210 1209
1211 1210 return dataOut
1212 1211
1213 1212 class CombineProfiles(Operation):
1214 1213 def __init__(self, **kwargs):
1215 1214
1216 1215 Operation.__init__(self, **kwargs)
1217 1216
1218 1217 self.__remData = None
1219 1218 self.__profileIndex = 0
1220 1219
1221 1220 def run(self, dataOut, n):
1222 1221
1223 1222 dataOut.flagNoData = True
1224 1223 profileIndex = None
1225 1224
1226 1225 if dataOut.flagDataAsBlock:
1227 1226
1228 1227 #nchannels, nprofiles, nsamples
1229 1228 shape = dataOut.data.shape
1230 1229 new_shape = shape[0], shape[1]/n, shape[2]*n
1231 1230
1232 1231 if shape[1] % n != 0:
1233 1232 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1234 1233
1235 1234 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1236 1235 dataOut.flagNoData = False
1237 1236
1238 1237 profileIndex = int(dataOut.nProfiles*n) - 1
1239 1238
1240 1239 else:
1241 1240
1242 1241 #nchannels, nsamples
1243 1242 if self.__remData is None:
1244 1243 newData = dataOut.data
1245 1244 else:
1246 1245 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1247 1246
1248 1247 self.__profileIndex += 1
1249 1248
1250 1249 if self.__profileIndex < n:
1251 1250 self.__remData = newData
1252 1251 #continue
1253 1252 return
1254 1253
1255 1254 self.__profileIndex = 0
1256 1255 self.__remData = None
1257 1256
1258 1257 dataOut.data = newData
1259 1258 dataOut.flagNoData = False
1260 1259
1261 1260 profileIndex = dataOut.profileIndex/n
1262 1261
1263 1262
1264 1263 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1265 1264
1266 1265 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1267 1266
1268 1267 dataOut.nProfiles = int(dataOut.nProfiles/n)
1269 1268
1270 1269 dataOut.profileIndex = profileIndex
1271 1270
1272 1271 dataOut.ippSeconds *= n
1273 1272
1274 1273 return dataOut
1274
1275 class PulsePairVoltage(Operation):
1276 '''
1277 Function PulsePair(Signal Power, Velocity)
1278 The real component of Lag[0] provides Intensity Information
1279 The imag component of Lag[1] Phase provides Velocity Information
1280
1281 Configuration Parameters:
1282 nPRF = Number of Several PRF
1283 theta = Degree Azimuth angel Boundaries
1284
1285 Input:
1286 self.dataOut
1287 lag[N]
1288 Affected:
1289 self.dataOut.spc
1290 '''
1291 isConfig = False
1292 __profIndex = 0
1293 __initime = None
1294 __lastdatatime = None
1295 __buffer = None
1296 __buffer2 = []
1297 __buffer3 = None
1298 __dataReady = False
1299 n = None
1300 __nch = 0
1301 __nHeis = 0
1302 removeDC = False
1303 ipp = None
1304 lambda_ = 0
1305
1306 def __init__(self,**kwargs):
1307 Operation.__init__(self,**kwargs)
1308
1309 def setup(self, dataOut, n = None, removeDC=False):
1310 '''
1311 n= Numero de PRF's de entrada
1312 '''
1313 self.__initime = None
1314 self.__lastdatatime = 0
1315 self.__dataReady = False
1316 self.__buffer = 0
1317 self.__buffer2 = []
1318 self.__buffer3 = 0
1319 self.__profIndex = 0
1320
1321 self.__nch = dataOut.nChannels
1322 self.__nHeis = dataOut.nHeights
1323 self.removeDC = removeDC
1324 self.lambda_ = 3.0e8/(9345.0e6)
1325 self.ippSec = dataOut.ippSeconds
1326 self.nCohInt = dataOut.nCohInt
1327 print("IPPseconds",dataOut.ippSeconds)
1328
1329 print("ELVALOR DE n es:", n)
1330 if n == None:
1331 raise ValueError("n should be specified.")
1332
1333 if n != None:
1334 if n<2:
1335 raise ValueError("n should be greater than 2")
1336
1337 self.n = n
1338 self.__nProf = n
1339
1340 self.__buffer = numpy.zeros((dataOut.nChannels,
1341 n,
1342 dataOut.nHeights),
1343 dtype='complex')
1344
1345 def putData(self,data):
1346 '''
1347 Add a profile to he __buffer and increase in one the __profiel Index
1348 '''
1349 self.__buffer[:,self.__profIndex,:]= data
1350 self.__profIndex += 1
1351 return
1352
1353 def pushData(self):
1354 '''
1355 Return the PULSEPAIR and the profiles used in the operation
1356 Affected : self.__profileIndex
1357 '''
1358
1359 if self.removeDC==True:
1360 mean = numpy.mean(self.__buffer,1)
1361 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1362 dc= numpy.tile(tmp,[1,self.__nProf,1])
1363 self.__buffer = self.__buffer - dc
1364
1365 data_intensity = numpy.sum(self.__buffer*numpy.conj(self.__buffer),1)/(self.n*self.nCohInt)#*self.nCohInt)
1366 pair1 = self.__buffer[:,1:,:]*numpy.conjugate(self.__buffer[:,:-1,:])
1367 angle = numpy.angle(numpy.sum(pair1,1))*180/(math.pi)
1368 #print(angle.shape)#print("__ANGLE__") #print("angle",angle[:,:10])
1369 data_velocity = (self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(numpy.sum(pair1,1))#self.ippSec*self.nCohInt
1370 n = self.__profIndex
1371
1372 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1373 self.__profIndex = 0
1374 return data_intensity,data_velocity,n
1375
1376 def pulsePairbyProfiles(self,data):
1377
1378 self.__dataReady = False
1379 data_intensity = None
1380 data_velocity = None
1381 self.putData(data)
1382 if self.__profIndex == self.n:
1383 data_intensity, data_velocity, n = self.pushData()
1384 self.__dataReady = True
1385
1386 return data_intensity, data_velocity
1387
1388 def pulsePairOp(self, data, datatime= None):
1389
1390 if self.__initime == None:
1391 self.__initime = datatime
1392
1393 data_intensity, data_velocity = self.pulsePairbyProfiles(data)
1394 self.__lastdatatime = datatime
1395
1396 if data_intensity is None:
1397 return None, None, None
1398
1399 avgdatatime = self.__initime
1400 deltatime = datatime - self.__lastdatatime
1401 self.__initime = datatime
1402
1403 return data_intensity, data_velocity, avgdatatime
1404
1405 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1406
1407 if not self.isConfig:
1408 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1409 self.isConfig = True
1410 data_intensity, data_velocity, avgdatatime = self.pulsePairOp(dataOut.data, dataOut.utctime)
1411 dataOut.flagNoData = True
1412
1413 if self.__dataReady:
1414 dataOut.nCohInt *= self.n
1415 dataOut.data_intensity = data_intensity #valor para intensidad
1416 dataOut.data_velocity = data_velocity #valor para velocidad
1417 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1418 dataOut.utctime = avgdatatime
1419 dataOut.flagNoData = False
1420 return dataOut
1421
1422
1275 1423 # import collections
1276 1424 # from scipy.stats import mode
1277 1425 #
1278 1426 # class Synchronize(Operation):
1279 1427 #
1280 1428 # isConfig = False
1281 1429 # __profIndex = 0
1282 1430 #
1283 1431 # def __init__(self, **kwargs):
1284 1432 #
1285 1433 # Operation.__init__(self, **kwargs)
1286 1434 # # self.isConfig = False
1287 1435 # self.__powBuffer = None
1288 1436 # self.__startIndex = 0
1289 1437 # self.__pulseFound = False
1290 1438 #
1291 1439 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1292 1440 #
1293 1441 # #Read data
1294 1442 #
1295 1443 # powerdB = dataOut.getPower(channel = channel)
1296 1444 # noisedB = dataOut.getNoise(channel = channel)[0]
1297 1445 #
1298 1446 # self.__powBuffer.extend(powerdB.flatten())
1299 1447 #
1300 1448 # dataArray = numpy.array(self.__powBuffer)
1301 1449 #
1302 1450 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1303 1451 #
1304 1452 # maxValue = numpy.nanmax(filteredPower)
1305 1453 #
1306 1454 # if maxValue < noisedB + 10:
1307 1455 # #No se encuentra ningun pulso de transmision
1308 1456 # return None
1309 1457 #
1310 1458 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1311 1459 #
1312 1460 # if len(maxValuesIndex) < 2:
1313 1461 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1314 1462 # return None
1315 1463 #
1316 1464 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1317 1465 #
1318 1466 # #Seleccionar solo valores con un espaciamiento de nSamples
1319 1467 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1320 1468 #
1321 1469 # if len(pulseIndex) < 2:
1322 1470 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1323 1471 # return None
1324 1472 #
1325 1473 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1326 1474 #
1327 1475 # #remover senales que se distancien menos de 10 unidades o muestras
1328 1476 # #(No deberian existir IPP menor a 10 unidades)
1329 1477 #
1330 1478 # realIndex = numpy.where(spacing > 10 )[0]
1331 1479 #
1332 1480 # if len(realIndex) < 2:
1333 1481 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1334 1482 # return None
1335 1483 #
1336 1484 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1337 1485 # realPulseIndex = pulseIndex[realIndex]
1338 1486 #
1339 1487 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1340 1488 #
1341 1489 # print "IPP = %d samples" %period
1342 1490 #
1343 1491 # self.__newNSamples = dataOut.nHeights #int(period)
1344 1492 # self.__startIndex = int(realPulseIndex[0])
1345 1493 #
1346 1494 # return 1
1347 1495 #
1348 1496 #
1349 1497 # def setup(self, nSamples, nChannels, buffer_size = 4):
1350 1498 #
1351 1499 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1352 1500 # maxlen = buffer_size*nSamples)
1353 1501 #
1354 1502 # bufferList = []
1355 1503 #
1356 1504 # for i in range(nChannels):
1357 1505 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1358 1506 # maxlen = buffer_size*nSamples)
1359 1507 #
1360 1508 # bufferList.append(bufferByChannel)
1361 1509 #
1362 1510 # self.__nSamples = nSamples
1363 1511 # self.__nChannels = nChannels
1364 1512 # self.__bufferList = bufferList
1365 1513 #
1366 1514 # def run(self, dataOut, channel = 0):
1367 1515 #
1368 1516 # if not self.isConfig:
1369 1517 # nSamples = dataOut.nHeights
1370 1518 # nChannels = dataOut.nChannels
1371 1519 # self.setup(nSamples, nChannels)
1372 1520 # self.isConfig = True
1373 1521 #
1374 1522 # #Append new data to internal buffer
1375 1523 # for thisChannel in range(self.__nChannels):
1376 1524 # bufferByChannel = self.__bufferList[thisChannel]
1377 1525 # bufferByChannel.extend(dataOut.data[thisChannel])
1378 1526 #
1379 1527 # if self.__pulseFound:
1380 1528 # self.__startIndex -= self.__nSamples
1381 1529 #
1382 1530 # #Finding Tx Pulse
1383 1531 # if not self.__pulseFound:
1384 1532 # indexFound = self.__findTxPulse(dataOut, channel)
1385 1533 #
1386 1534 # if indexFound == None:
1387 1535 # dataOut.flagNoData = True
1388 1536 # return
1389 1537 #
1390 1538 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1391 1539 # self.__pulseFound = True
1392 1540 # self.__startIndex = indexFound
1393 1541 #
1394 1542 # #If pulse was found ...
1395 1543 # for thisChannel in range(self.__nChannels):
1396 1544 # bufferByChannel = self.__bufferList[thisChannel]
1397 1545 # #print self.__startIndex
1398 1546 # x = numpy.array(bufferByChannel)
1399 1547 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1400 1548 #
1401 1549 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1402 1550 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1403 1551 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1404 1552 #
1405 1553 # dataOut.data = self.__arrayBuffer
1406 1554 #
1407 1555 # self.__startIndex += self.__newNSamples
1408 1556 #
1409 1557 # return
General Comments 0
You need to be logged in to leave comments. Login now