##// END OF EJS Templates
Se adiciona script Julia_EEJ.py para test y la operaciones GetSNR
Alexander Valdez -
r1684:1f582a059484
parent child
Show More
@@ -0,0 +1,291
1 #!/usr/bin/env python
2 '''
3 Created on Jul 7, 2014
4
5 @author: roj-idl71
6 '''
7 import os, sys, json
8
9 #path = os.path.dirname(os.getcwd())
10 #path = os.path.join(path, 'source')
11 #sys.path.insert(0, path)
12
13 from schainpy.controller import Project
14
15 if __name__ == '__main__':
16
17 desc = "JULIA raw experiment "
18 filename = "schain.xml"
19
20 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2022/'
21 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2021/'
22 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2017/'
23 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2021/'
24 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2018/'
25 #dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2017/'
26 #dpath = '/home/roberto/Folder_aux/D2019/'
27
28 ##dpath = '/home/roberto/puma/JULIA_EW_IMAGING/JULIA_EW/D2017/'
29 dpath = '/home/soporte/PUMA/JULIA_EW_IMAGING/JULIA_EW/D2017/'
30
31 figpath = '/home/soporte/DATA/Pictures/JULIA/EEJ/Skew_but_dop_is_shift'+'/'+dpath[-5:-1]
32 ppath = "/home/soporte/DATA/MLT/Oblique/2022_03/data_reshape"
33 spcdata = '/home/soporte/DATA/JULIA/EEJ/SPC'
34 path_parameters = '/home/soporte/DATA/JULIA/EEJ/Params'
35 db_range=['25','35']
36 db_range=['10','20']
37 #db_range=['14','20']
38 db_range=['10','23']
39 db_range=['13','20']
40 db_range=['21','30']
41 db_range=['15','30']
42 db_range=['13','28']
43 tiempo=['7','18']
44 altura1=[2,20]
45 altura1=[90,220]
46 velocity=['-80','80']
47 period=60
48 # PROJECT 1
49
50 show_spc = 0
51 save_spc = 0
52 fitting = 1
53 save_params = 1
54 plot_params = 0
55
56 controllerObj = Project()
57 controllerObj.setup(id = '191', name='altura1', description=desc)
58
59 readUnitConfObj1 = controllerObj.addReadUnit(datatype='SpectraReader',
60 path=dpath,
61 startDate='2017/09/01', #Check 21-29 Jun 2021, 18 May 2022
62 endDate='2017/09/30', #05,06,07-01-18
63 #startTime='06:00:00',
64 #endTime='18:00:00',
65 startTime='07:00:00',
66 #startTime='10:00:00',
67 #startTime='08:38:01',
68 #startTime='11:16:00',
69 #startTime='08:13:00',
70 endTime='17:59:59',
71 #startTime='16:30:00',
72 #endTime='17:30:59',
73 online=0,
74 walk=1,
75 expLabel='150EEJ',
76 getByBlock=1,
77 delay=20)
78
79 # opObj00 = readUnitConfObj.addOperation(name='printInfo')
80 # opObj00 = readUnitConfObj.addOperation(name='printNumberOfBlock')
81
82 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=readUnitConfObj1.getId())
83
84 opObj11 = procUnitConfObj1.addOperation(name='selectChannels')
85 opObj11.addParameter(name='channelList', value='2,', format='intlist')
86
87 '''
88 opObj11 = procUnitConfObj1SPC.addOperation(name='removeDC')
89 opObj11.addParameter(name='mode', value='2', format='int')
90 '''
91
92 #opObj11 = procUnitConfObj1.addOperation(name='dopplerFlip') #It fixes the Doppler
93 #opObj11.addParameter(name='chann', value='0', format='int')
94
95 opObj11 = procUnitConfObj1.addOperation(name='removeInterference')
96
97 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
98 #opObj11.addParameter(name='n', value='20', format='int')
99 opObj11.addParameter(name='n', value='1', format='int')
100
101 opObj11 = procUnitConfObj1.addOperation(name='GetSNR', optype='other')
102
103 if save_spc:
104 dataList=['data_spc',
105 'utctime']
106 metadataList=['nFFTPoints','VelRange','normFactor',
107 'heightList','timeZone']
108
109 op221 = procUnitConfObj1.addOperation(name='HDFWriter', optype='external')
110 op221.addParameter(name='path', value=spcdata)
111 #op221.addParameter(name='mode', value=1, format='int')
112 op221.addParameter(name='dataList', value=dataList)
113 op221.addParameter(name='metadataList', value=metadataList)
114 #op221.addParameter(name='blocksPerFile', value=500)
115 op221.addParameter(name='blocksPerFile', value=2000)
116
117 if show_spc:
118 #'''
119 # opObj11 = procUnitConfObj1SPC.addOperation(name='removeInterference')
120
121 opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
122 opObj11.addParameter(name='id', value='1', format='int')
123 opObj11.addParameter(name='wintitle', value='Oblique', format='str')
124 opObj11.addParameter(name='zmin', value=db_range[0], format='int')
125 opObj11.addParameter(name='zmax', value=db_range[1], format='int')
126 opObj11.addParameter(name='xaxis', value='velocity', format='str')
127 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
128 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
129 # opObj11.addParameter(name='xmin', value=velocity[0], format='int')
130 # opObj11.addParameter(name='xmax', value=velocity[1], format='int')
131 opObj11.addParameter(name='showprofile', value='1', format='int')
132 opObj11.addParameter(name='save', value=figpath, format='str')
133 #'''
134 #'''
135 opObj11 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
136 opObj11.addParameter(name='id', value='10', format='int')
137 opObj11.addParameter(name='wintitle', value='JULIA EEJ RTI', format='str')
138 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
139 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
140 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
141 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
142 opObj11.addParameter(name='zmin', value=db_range[0], format='int')
143 opObj11.addParameter(name='zmax', value=db_range[1], format='int')
144 opObj11.addParameter(name='showprofile', value='1', format='int')
145 #opObj11.addParameter(name='save_period', value=40, format='str')
146 opObj11.addParameter(name='save', value=figpath, format='str')
147 #opObj11.addParameter(name='throttle', value=1, format='str')
148 #'''
149 '''
150 opObj11 = procUnitConfObj1.addOperation(name='SnrPlot', optype='other')
151 opObj11.addParameter(name='id', value='10', format='int')
152 opObj11.addParameter(name='wintitle', value='JULIA EEJ SNR', format='str')
153 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
154 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
155 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
156 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
157 opObj11.addParameter(name='zmin', value=0.1, format='int')
158 opObj11.addParameter(name='zmax', value=50, format='int')
159 #opObj11.addParameter(name='showprofile', value='1', format='int')
160 opObj11.addParameter(name='save', value=figpath, format='str')
161 '''
162 if fitting:
163 Dop = 'Max' #Plot and Save the Pos[Max_val] as the Doppler Shift
164 Dop = 'Shift' #Plot and Save the Skew Gaussian Shift as the Doppler Shift
165 opObj11 = procUnitConfObj1.addOperation(name='Oblique_Gauss_Fit', optype='other')
166 opObj11.addParameter(name='mode', value=9, format='int') #Skew
167 opObj11.addParameter(name='Dop', value=Dop)
168 #opObj11.addParameter(name='mode', value=11, format='int')
169
170 if save_params:
171 '''
172 dataList=['powerdB', 'Oblique_params', 'Oblique_param_errors', 'dplr_2_u', 'data_snr',
173 'utctime']
174 metadataList=['VelRange',
175 'heightList','timeZone']
176
177 op221 = procUnitConfObj1.addOperation(name='HDFWriter', optype='external')
178 op221.addParameter(name='path', value=path_parameters)
179 #op221.addParameter(name='mode', value=1, format='int')
180 op221.addParameter(name='dataList', value=dataList)
181 op221.addParameter(name='metadataList', value=metadataList)
182 #op221.addParameter(name='blocksPerFile', value=500)
183 op221.addParameter(name='blocksPerFile', value=2000)
184 '''
185
186 one = {'gdlatr': 'lat', 'gdlonr': 'lon', 'inttms': 'paramInterval'} #reader gdlatr-->lat only 1D
187
188 two = {
189 'snl': 'snl', #DeberΓ­a salir como el original pero mΓ‘s limpio
190 'RANGE': 'heightList', #<----- nmonics
191 'DOPP_T1_EEJ': ('Dop_EEJ_T1', (0)),
192 'DDOPP_T1_EEJ': ('Err_Dop_EEJ_T1', (0)),
193 'SPEC_W_T1_EEJ': ('Spec_W_T1', (0)),
194 'DSPEC_W_T1_EEJ': ('Err_Spec_W_T1', (0)),
195 'DOPP_T2_EEJ': ('Dop_EEJ_T2', (0)),
196 'DDOPP_T2_EEJ': ('Err_Dop_EEJ_T2', (0)),
197 'SPEC_W_T2_EEJ': ('Spec_W_T2', (0)),
198 'DSPEC_W_T2_EEJ': ('Err_Spec_W_T2', (0)),
199 } #writer
200
201
202 ind = ['range']
203
204 meta = {
205 'kinst': 840, #instrumnet code, 840 for JULIA, 14 for JULIA MP CSR
206 'kindat': 1962, #type of data #Este es el nuevo e igual para JULIA y JULIA CSR
207 'catalog': {
208 'principleInvestigator': 'Danny ScipiΓ³n',
209 'expPurpose': 'Equatorial Electrojet Parameters',
210 },
211 'header': {
212 'analyst': 'D. Hysell'
213 }
214 }
215
216
217 op_writer = procUnitConfObj1.addOperation(name='MADWriter', optype='external')
218 #op_writer.addParameter(name='path', value='/home/roberto/DATA/hdf5_outputs/Madrigal/EEJ')
219 #op_writer.addParameter(name='path', value='/home/roberto/DATA/hdf5_outputs/Madrigal/EEJ/Dop_Max_Val')
220 #op_writer.addParameter(name='path', value='/home/roberto/DATA/hdf5_outputs/Madrigal/EEJ/'+Dop+'/'+dpath[-5:-1])
221 #op_writer.addParameter(name='path', value='/home/roberto/DATA/hdf5_outputs/Madrigal/EEJ/CorrectFiles/no_snl_Test/01/'+Dop+'/'+dpath[-5:-1])
222 op_writer.addParameter(name='path', value='/home/soporte/DATA/hdf5_outputs/Madrigal/EEJ/FinalFiles/'+Dop+'/'+dpath[-5:-1])
223 op_writer.addParameter(name='format', value='hdf5', format='str')
224 op_writer.addParameter(name='oneDDict', value=json.dumps(one), format='str')
225 op_writer.addParameter(name='twoDDict', value=json.dumps(two), format='str')
226 op_writer.addParameter(name='ind2DList', value=json.dumps(ind), format='str')
227 op_writer.addParameter(name='metadata', value=json.dumps(meta), format='str')
228
229
230 if plot_params:
231 '''
232 opObj11 = procUnitConfObj1.addOperation(name='DopplerEEJPlot', optype='other')
233 opObj11.addParameter(name='id', value='10', format='int')
234 opObj11.addParameter(name='wintitle', value='Doppler EEJ', format='str')
235 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
236 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
237 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
238 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
239 #opObj11.addParameter(name='zmin', value=-250, format='int')
240 #opObj11.addParameter(name='zmax', value=250, format='int')
241 opObj11.addParameter(name='zlimits', value='(-400,400),(-250,250)')
242 #opObj11.addParameter(name='showprofile', value='1', format='int')
243 opObj11.addParameter(name='save', value=figpath, format='str')
244 #opObj11.addParameter(name='EEJtype', value=1, format='int')
245
246 opObj11 = procUnitConfObj1.addOperation(name='SpcWidthEEJPlot', optype='other')
247 opObj11.addParameter(name='id', value='10', format='int')
248 opObj11.addParameter(name='wintitle', value='Spectral Width EEJ', format='str')
249 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
250 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
251 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
252 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
253 #opObj11.addParameter(name='zmin', value=0., format='int')
254 #opObj11.addParameter(name='zmax', value=250, format='int')
255 opObj11.addParameter(name='zlimits', value='(0.1,100),(0.1,250)')
256 #opObj11.addParameter(name='showprofile', value='1', format='int')
257 opObj11.addParameter(name='save', value=figpath, format='str')
258 #opObj11.addParameter(name='EEJtype', value=1, format='int')
259 '''
260 opObj11 = procUnitConfObj1.addOperation(name='SpectraObliquePlot', optype='other')
261 opObj11.addParameter(name='id', value='1', format='int')
262 opObj11.addParameter(name='wintitle', value='Oblique', format='str')
263 opObj11.addParameter(name='zmin', value=db_range[0], format='int')
264 opObj11.addParameter(name='zmax', value=db_range[1], format='int')
265 opObj11.addParameter(name='xaxis', value='velocity', format='str')
266 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
267 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
268 # opObj11.addParameter(name='xmin', value=velocity[0], format='int')
269 # opObj11.addParameter(name='xmax', value=velocity[1], format='int')
270 opObj11.addParameter(name='showprofile', value='1', format='int')
271 opObj11.addParameter(name='save', value=figpath+'/400', format='str')
272
273 '''
274 opObj11 = procUnitConfObj1.addOperation(name='DopplerEEJPlot', optype='other')
275 opObj11.addParameter(name='id', value='10', format='int')
276 opObj11.addParameter(name='wintitle', value='Doppler EEJ Type II', format='str')
277 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
278 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
279 # opObj11.addParameter(name='ymin', value=altura1[0], format='int')
280 # opObj11.addParameter(name='ymax', value=altura1[1], format='int')
281 opObj11.addParameter(name='zmin', value=-250, format='int')
282 opObj11.addParameter(name='zmax', value=250, format='int')
283 #opObj11.addParameter(name='showprofile', value='1', format='int')
284 opObj11.addParameter(name='save', value=figpath, format='str')
285 opObj11.addParameter(name='EEJtype', value=2, format='int')
286 '''
287
288
289
290
291 controllerObj.start() No newline at end of file
@@ -1,937 +1,967
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13 import numpy
14 14 # repositorio
15 15 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
16 16 from schainpy.model.data.jrodata import Spectra
17 17 from schainpy.model.data.jrodata import hildebrand_sekhon
18 18 from schainpy.utils import log
19 19
20 20
21 21 class SpectraProc(ProcessingUnit):
22 22
23 23 def __init__(self):
24 24
25 25 ProcessingUnit.__init__(self)
26 26
27 27 self.buffer = None
28 28 self.firstdatatime = None
29 29 self.profIndex = 0
30 30 self.dataOut = Spectra()
31 31 self.id_min = None
32 32 self.id_max = None
33 33 self.setupReq = False #Agregar a todas las unidades de proc
34 34
35 35 def __updateSpecFromVoltage(self):
36 36
37 37 self.dataOut.timeZone = self.dataIn.timeZone
38 38 self.dataOut.dstFlag = self.dataIn.dstFlag
39 39 self.dataOut.errorCount = self.dataIn.errorCount
40 40 self.dataOut.useLocalTime = self.dataIn.useLocalTime
41 41 try:
42 42 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
43 43 except:
44 44 pass
45 45 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
46 46 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
47 47 self.dataOut.channelList = self.dataIn.channelList
48 48 self.dataOut.heightList = self.dataIn.heightList
49 49 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
50 50 self.dataOut.nProfiles = self.dataOut.nFFTPoints
51 51 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
52 52 self.dataOut.utctime = self.firstdatatime
53 53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
54 54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
55 55 self.dataOut.flagShiftFFT = False
56 56 self.dataOut.nCohInt = self.dataIn.nCohInt
57 57 self.dataOut.nIncohInt = 1
58 58 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59 self.dataOut.frequency = self.dataIn.frequency
60 60 self.dataOut.realtime = self.dataIn.realtime
61 61 self.dataOut.azimuth = self.dataIn.azimuth
62 62 self.dataOut.zenith = self.dataIn.zenith
63 63 self.dataOut.beam.codeList = self.dataIn.beam.codeList
64 64 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
65 65 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
66 66 self.dataOut.runNextUnit = self.dataIn.runNextUnit
67 67 try:
68 68 self.dataOut.step = self.dataIn.step
69 69 except:
70 70 pass
71 71
72 72 def __getFft(self):
73 73 """
74 74 Convierte valores de Voltaje a Spectra
75 75
76 76 Affected:
77 77 self.dataOut.data_spc
78 78 self.dataOut.data_cspc
79 79 self.dataOut.data_dc
80 80 self.dataOut.heightList
81 81 self.profIndex
82 82 self.buffer
83 83 self.dataOut.flagNoData
84 84 """
85 85 fft_volt = numpy.fft.fft(
86 86 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
87 87 fft_volt = fft_volt.astype(numpy.dtype('complex'))
88 88 dc = fft_volt[:, 0, :]
89 89
90 90 # calculo de self-spectra
91 91 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
92 92 spc = fft_volt * numpy.conjugate(fft_volt)
93 93 spc = spc.real
94 94
95 95 blocksize = 0
96 96 blocksize += dc.size
97 97 blocksize += spc.size
98 98
99 99 cspc = None
100 100 pairIndex = 0
101 101 if self.dataOut.pairsList != None:
102 102 # calculo de cross-spectra
103 103 cspc = numpy.zeros(
104 104 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
105 105 for pair in self.dataOut.pairsList:
106 106 if pair[0] not in self.dataOut.channelList:
107 107 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
108 108 str(pair), str(self.dataOut.channelList)))
109 109 if pair[1] not in self.dataOut.channelList:
110 110 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
111 111 str(pair), str(self.dataOut.channelList)))
112 112
113 113 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
114 114 numpy.conjugate(fft_volt[pair[1], :, :])
115 115 pairIndex += 1
116 116 blocksize += cspc.size
117 117
118 118 self.dataOut.data_spc = spc
119 119 self.dataOut.data_cspc = cspc
120 120 self.dataOut.data_dc = dc
121 121 self.dataOut.blockSize = blocksize
122 122 self.dataOut.flagShiftFFT = False
123 123
124 124 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
125 125
126 126 self.dataIn.runNextUnit = runNextUnit
127 127 if self.dataIn.type == "Spectra":
128 128 self.dataOut.copy(self.dataIn)
129 129 if shift_fft:
130 130 #desplaza a la derecha en el eje 2 determinadas posiciones
131 131 shift = int(self.dataOut.nFFTPoints/2)
132 132 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
133 133
134 134 if self.dataOut.data_cspc is not None:
135 135 #desplaza a la derecha en el eje 2 determinadas posiciones
136 136 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
137 137 if pairsList:
138 138 self.__selectPairs(pairsList)
139 139
140 140 elif self.dataIn.type == "Voltage":
141 141
142 142 self.dataOut.flagNoData = True
143 143
144 144 if nFFTPoints == None:
145 145 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
146 146
147 147 if nProfiles == None:
148 148 nProfiles = nFFTPoints
149 149
150 150 if ippFactor == None:
151 151 self.dataOut.ippFactor = 1
152 152
153 153 self.dataOut.nFFTPoints = nFFTPoints
154 154
155 155 if self.buffer is None:
156 156 self.buffer = numpy.zeros((self.dataIn.nChannels,
157 157 nProfiles,
158 158 self.dataIn.nHeights),
159 159 dtype='complex')
160 160
161 161 if self.dataIn.flagDataAsBlock:
162 162 nVoltProfiles = self.dataIn.data.shape[1]
163 163 if nVoltProfiles == nProfiles:
164 164 self.buffer = self.dataIn.data.copy()
165 165 self.profIndex = nVoltProfiles
166 166
167 167 elif nVoltProfiles < nProfiles:
168 168
169 169 if self.profIndex == 0:
170 170 self.id_min = 0
171 171 self.id_max = nVoltProfiles
172 172
173 173 self.buffer[:, self.id_min:self.id_max,
174 174 :] = self.dataIn.data
175 175 self.profIndex += nVoltProfiles
176 176 self.id_min += nVoltProfiles
177 177 self.id_max += nVoltProfiles
178 178 elif nVoltProfiles > nProfiles:
179 179 self.reader.bypass = True
180 180 if self.profIndex == 0:
181 181 self.id_min = 0
182 182 self.id_max = nProfiles
183 183
184 184 self.buffer = self.dataIn.data[:, self.id_min:self.id_max,:]
185 185 self.profIndex += nProfiles
186 186 self.id_min += nProfiles
187 187 self.id_max += nProfiles
188 188 if self.id_max == nVoltProfiles:
189 189 self.reader.bypass = False
190 190
191 191 else:
192 192 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
193 193 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
194 194 self.dataOut.flagNoData = True
195 195 else:
196 196 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
197 197 self.profIndex += 1
198 198
199 199 if self.firstdatatime == None:
200 200 self.firstdatatime = self.dataIn.utctime
201 201
202 202 if self.profIndex == nProfiles:
203 203 self.__updateSpecFromVoltage()
204 204 if pairsList == None:
205 205 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
206 206 else:
207 207 self.dataOut.pairsList = pairsList
208 208 self.__getFft()
209 209 self.dataOut.flagNoData = False
210 210 self.firstdatatime = None
211 211 #if not self.reader.bypass:
212 212 self.profIndex = 0
213 213 else:
214 214 raise ValueError("The type of input object '%s' is not valid".format(
215 215 self.dataIn.type))
216 216
217 217 def __selectPairs(self, pairsList):
218 218
219 219 if not pairsList:
220 220 return
221 221
222 222 pairs = []
223 223 pairsIndex = []
224 224
225 225 for pair in pairsList:
226 226 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
227 227 continue
228 228 pairs.append(pair)
229 229 pairsIndex.append(pairs.index(pair))
230 230
231 231 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
232 232 self.dataOut.pairsList = pairs
233 233
234 234 return
235 235
236 236 def selectFFTs(self, minFFT, maxFFT ):
237 237 """
238 238 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
239 239 minFFT<= FFT <= maxFFT
240 240 """
241 241
242 242 if (minFFT > maxFFT):
243 243 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
244 244
245 245 if (minFFT < self.dataOut.getFreqRange()[0]):
246 246 minFFT = self.dataOut.getFreqRange()[0]
247 247
248 248 if (maxFFT > self.dataOut.getFreqRange()[-1]):
249 249 maxFFT = self.dataOut.getFreqRange()[-1]
250 250
251 251 minIndex = 0
252 252 maxIndex = 0
253 253 FFTs = self.dataOut.getFreqRange()
254 254
255 255 inda = numpy.where(FFTs >= minFFT)
256 256 indb = numpy.where(FFTs <= maxFFT)
257 257
258 258 try:
259 259 minIndex = inda[0][0]
260 260 except:
261 261 minIndex = 0
262 262
263 263 try:
264 264 maxIndex = indb[0][-1]
265 265 except:
266 266 maxIndex = len(FFTs)
267 267
268 268 self.selectFFTsByIndex(minIndex, maxIndex)
269 269
270 270 return 1
271 271
272 272 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
273 273 newheis = numpy.where(
274 274 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
275 275
276 276 if hei_ref != None:
277 277 newheis = numpy.where(self.dataOut.heightList > hei_ref)
278 278
279 279 minIndex = min(newheis[0])
280 280 maxIndex = max(newheis[0])
281 281 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
282 282 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
283 283
284 284 # determina indices
285 285 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
286 286 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
287 287 avg_dB = 10 * \
288 288 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
289 289 beacon_dB = numpy.sort(avg_dB)[-nheis:]
290 290 beacon_heiIndexList = []
291 291 for val in avg_dB.tolist():
292 292 if val >= beacon_dB[0]:
293 293 beacon_heiIndexList.append(avg_dB.tolist().index(val))
294 294
295 295 data_cspc = None
296 296 if self.dataOut.data_cspc is not None:
297 297 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
298 298
299 299 data_dc = None
300 300 if self.dataOut.data_dc is not None:
301 301 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
302 302
303 303 self.dataOut.data_spc = data_spc
304 304 self.dataOut.data_cspc = data_cspc
305 305 self.dataOut.data_dc = data_dc
306 306 self.dataOut.heightList = heightList
307 307 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
308 308
309 309 return 1
310 310
311 311 def selectFFTsByIndex(self, minIndex, maxIndex):
312 312 """
313 313
314 314 """
315 315
316 316 if (minIndex < 0) or (minIndex > maxIndex):
317 317 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
318 318
319 319 if (maxIndex >= self.dataOut.nProfiles):
320 320 maxIndex = self.dataOut.nProfiles-1
321 321
322 322 #Spectra
323 323 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
324 324
325 325 data_cspc = None
326 326 if self.dataOut.data_cspc is not None:
327 327 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
328 328
329 329 data_dc = None
330 330 if self.dataOut.data_dc is not None:
331 331 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
332 332
333 333 self.dataOut.data_spc = data_spc
334 334 self.dataOut.data_cspc = data_cspc
335 335 self.dataOut.data_dc = data_dc
336 336
337 337 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
338 338 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
339 339 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
340 340
341 341 return 1
342 342
343 343 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
344 344 # validacion de rango
345 345 if minHei == None:
346 346 minHei = self.dataOut.heightList[0]
347 347
348 348 if maxHei == None:
349 349 maxHei = self.dataOut.heightList[-1]
350 350
351 351 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
352 352 print('minHei: %.2f is out of the heights range' % (minHei))
353 353 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
354 354 minHei = self.dataOut.heightList[0]
355 355
356 356 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
357 357 print('maxHei: %.2f is out of the heights range' % (maxHei))
358 358 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
359 359 maxHei = self.dataOut.heightList[-1]
360 360
361 361 # validacion de velocidades
362 362 velrange = self.dataOut.getVelRange(1)
363 363
364 364 if minVel == None:
365 365 minVel = velrange[0]
366 366
367 367 if maxVel == None:
368 368 maxVel = velrange[-1]
369 369
370 370 if (minVel < velrange[0]) or (minVel > maxVel):
371 371 print('minVel: %.2f is out of the velocity range' % (minVel))
372 372 print('minVel is setting to %.2f' % (velrange[0]))
373 373 minVel = velrange[0]
374 374
375 375 if (maxVel > velrange[-1]) or (maxVel < minVel):
376 376 print('maxVel: %.2f is out of the velocity range' % (maxVel))
377 377 print('maxVel is setting to %.2f' % (velrange[-1]))
378 378 maxVel = velrange[-1]
379 379
380 380 # seleccion de indices para rango
381 381 minIndex = 0
382 382 maxIndex = 0
383 383 heights = self.dataOut.heightList
384 384
385 385 inda = numpy.where(heights >= minHei)
386 386 indb = numpy.where(heights <= maxHei)
387 387
388 388 try:
389 389 minIndex = inda[0][0]
390 390 except:
391 391 minIndex = 0
392 392
393 393 try:
394 394 maxIndex = indb[0][-1]
395 395 except:
396 396 maxIndex = len(heights)
397 397
398 398 if (minIndex < 0) or (minIndex > maxIndex):
399 399 raise ValueError("some value in (%d,%d) is not valid" % (
400 400 minIndex, maxIndex))
401 401
402 402 if (maxIndex >= self.dataOut.nHeights):
403 403 maxIndex = self.dataOut.nHeights - 1
404 404
405 405 # seleccion de indices para velocidades
406 406 indminvel = numpy.where(velrange >= minVel)
407 407 indmaxvel = numpy.where(velrange <= maxVel)
408 408 try:
409 409 minIndexVel = indminvel[0][0]
410 410 except:
411 411 minIndexVel = 0
412 412
413 413 try:
414 414 maxIndexVel = indmaxvel[0][-1]
415 415 except:
416 416 maxIndexVel = len(velrange)
417 417
418 418 # seleccion del espectro
419 419 data_spc = self.dataOut.data_spc[:,
420 420 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
421 421 # estimacion de ruido
422 422 noise = numpy.zeros(self.dataOut.nChannels)
423 423
424 424 for channel in range(self.dataOut.nChannels):
425 425 daux = data_spc[channel, :, :]
426 426 sortdata = numpy.sort(daux, axis=None)
427 427 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
428 428
429 429 self.dataOut.noise_estimation = noise.copy()
430 430
431 431 return 1
432 432
433
434 class GetSNR(Operation):
435 '''
436 Written by R. Flores
437 '''
438 """Operation to get SNR.
439
440 Parameters:
441 -----------
442
443 Example
444 --------
445
446 op = proc_unit.addOperation(name='GetSNR', optype='other')
447
448 """
449
450 def __init__(self, **kwargs):
451
452 Operation.__init__(self, **kwargs)
453
454 def run(self,dataOut):
455
456 noise = dataOut.getNoise(ymin_index=-10) #RegiΓ³n superior donde solo deberΓ­a de haber ruido
457 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
458 dataOut.snl = numpy.log10(dataOut.data_snr)
459 dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
460
461 return dataOut
462
433 463 class removeDC(Operation):
434 464
435 465 def run(self, dataOut, mode=2):
436 466 self.dataOut = dataOut
437 467 jspectra = self.dataOut.data_spc
438 468 jcspectra = self.dataOut.data_cspc
439 469
440 470 num_chan = jspectra.shape[0]
441 471 num_hei = jspectra.shape[2]
442 472
443 473 if jcspectra is not None:
444 474 jcspectraExist = True
445 475 num_pairs = jcspectra.shape[0]
446 476 else:
447 477 jcspectraExist = False
448 478
449 479 freq_dc = int(jspectra.shape[1] / 2)
450 480 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
451 481 ind_vel = ind_vel.astype(int)
452 482
453 483 if ind_vel[0] < 0:
454 484 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
455 485
456 486 if mode == 1:
457 487 jspectra[:, freq_dc, :] = (
458 488 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
459 489
460 490 if jcspectraExist:
461 491 jcspectra[:, freq_dc, :] = (
462 492 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
463 493
464 494 if mode == 2:
465 495
466 496 vel = numpy.array([-2, -1, 1, 2])
467 497 xx = numpy.zeros([4, 4])
468 498
469 499 for fil in range(4):
470 500 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
471 501
472 502 xx_inv = numpy.linalg.inv(xx)
473 503 xx_aux = xx_inv[0, :]
474 504
475 505 for ich in range(num_chan):
476 506 yy = jspectra[ich, ind_vel, :]
477 507 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
478 508
479 509 junkid = jspectra[ich, freq_dc, :] <= 0
480 510 cjunkid = sum(junkid)
481 511
482 512 if cjunkid.any():
483 513 jspectra[ich, freq_dc, junkid.nonzero()] = (
484 514 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
485 515
486 516 if jcspectraExist:
487 517 for ip in range(num_pairs):
488 518 yy = jcspectra[ip, ind_vel, :]
489 519 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
490 520
491 521 self.dataOut.data_spc = jspectra
492 522 self.dataOut.data_cspc = jcspectra
493 523
494 524 return self.dataOut
495 525
496 526 class removeInterference(Operation):
497 527
498 528 def removeInterference2(self):
499 529
500 530 cspc = self.dataOut.data_cspc
501 531 spc = self.dataOut.data_spc
502 532 Heights = numpy.arange(cspc.shape[2])
503 533 realCspc = numpy.abs(cspc)
504 534
505 535 for i in range(cspc.shape[0]):
506 536 LinePower= numpy.sum(realCspc[i], axis=0)
507 537 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
508 538 SelectedHeights = Heights[ numpy.where(LinePower < Threshold) ]
509 539 InterferenceSum = numpy.sum(realCspc[i,:,SelectedHeights],axis=0)
510 540 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
511 541 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
512 542
513 543
514 544 InterferenceRange = numpy.where(([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
515 545 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
516 546 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
517 547 cspc[i,InterferenceRange,:] = numpy.NaN
518 548
519 549 self.dataOut.data_cspc = cspc
520 550
521 551 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
522 552
523 553 jspectra = self.dataOut.data_spc
524 554 jcspectra = self.dataOut.data_cspc
525 555 jnoise = self.dataOut.getNoise()
526 556 num_incoh = self.dataOut.nIncohInt
527 557
528 558 num_channel = jspectra.shape[0]
529 559 num_prof = jspectra.shape[1]
530 560 num_hei = jspectra.shape[2]
531 561
532 562 # hei_interf
533 563 if hei_interf is None:
534 564 count_hei = int(num_hei / 2)
535 565 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
536 566 hei_interf = numpy.asarray(hei_interf)[0]
537 567 # nhei_interf
538 568 if (nhei_interf == None):
539 569 nhei_interf = 5
540 570 if (nhei_interf < 1):
541 571 nhei_interf = 1
542 572 if (nhei_interf > count_hei):
543 573 nhei_interf = count_hei
544 574 if (offhei_interf == None):
545 575 offhei_interf = 0
546 576
547 577 ind_hei = list(range(num_hei))
548 578 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
549 579 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
550 580 mask_prof = numpy.asarray(list(range(num_prof)))
551 581 num_mask_prof = mask_prof.size
552 582 comp_mask_prof = [0, num_prof / 2]
553 583
554 584 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
555 585 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
556 586 jnoise = numpy.nan
557 587 noise_exist = jnoise[0] < numpy.Inf
558 588
559 589 # Subrutina de Remocion de la Interferencia
560 590 for ich in range(num_channel):
561 591 # Se ordena los espectros segun su potencia (menor a mayor)
562 592 power = jspectra[ich, mask_prof, :]
563 593 power = power[:, hei_interf]
564 594 power = power.sum(axis=0)
565 595 psort = power.ravel().argsort()
566 596
567 597 # Se estima la interferencia promedio en los Espectros de Potencia empleando
568 598 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
569 599 offhei_interf, nhei_interf + offhei_interf))]]]
570 600
571 601 if noise_exist:
572 602 # tmp_noise = jnoise[ich] / num_prof
573 603 tmp_noise = jnoise[ich]
574 604 junkspc_interf = junkspc_interf - tmp_noise
575 605 #junkspc_interf[:,comp_mask_prof] = 0
576 606
577 607 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
578 608 jspc_interf = jspc_interf.transpose()
579 609 # Calculando el espectro de interferencia promedio
580 610 noiseid = numpy.where(
581 611 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
582 612 noiseid = noiseid[0]
583 613 cnoiseid = noiseid.size
584 614 interfid = numpy.where(
585 615 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
586 616 interfid = interfid[0]
587 617 cinterfid = interfid.size
588 618
589 619 if (cnoiseid > 0):
590 620 jspc_interf[noiseid] = 0
591 621
592 622 # Expandiendo los perfiles a limpiar
593 623 if (cinterfid > 0):
594 624 new_interfid = (
595 625 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
596 626 new_interfid = numpy.asarray(new_interfid)
597 627 new_interfid = {x for x in new_interfid}
598 628 new_interfid = numpy.array(list(new_interfid))
599 629 new_cinterfid = new_interfid.size
600 630 else:
601 631 new_cinterfid = 0
602 632
603 633 for ip in range(new_cinterfid):
604 634 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
605 635 jspc_interf[new_interfid[ip]
606 636 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
607 637
608 638 jspectra[ich, :, ind_hei] = jspectra[ich, :,
609 639 ind_hei] - jspc_interf # Corregir indices
610 640
611 641 # Removiendo la interferencia del punto de mayor interferencia
612 642 ListAux = jspc_interf[mask_prof].tolist()
613 643 maxid = ListAux.index(max(ListAux))
614 644
615 645 if cinterfid > 0:
616 646 for ip in range(cinterfid * (interf == 2) - 1):
617 647 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
618 648 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
619 649 cind = len(ind)
620 650
621 651 if (cind > 0):
622 652 jspectra[ich, interfid[ip], ind] = tmp_noise * \
623 653 (1 + (numpy.random.uniform(cind) - 0.5) /
624 654 numpy.sqrt(num_incoh))
625 655
626 656 ind = numpy.array([-2, -1, 1, 2])
627 657 xx = numpy.zeros([4, 4])
628 658
629 659 for id1 in range(4):
630 660 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
631 661
632 662 xx_inv = numpy.linalg.inv(xx)
633 663 xx = xx_inv[:, 0]
634 664 ind = (ind + maxid + num_mask_prof) % num_mask_prof
635 665 yy = jspectra[ich, mask_prof[ind], :]
636 666 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
637 667 yy.transpose(), xx)
638 668
639 669 indAux = (jspectra[ich, :, :] < tmp_noise *
640 670 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
641 671 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
642 672 (1 - 1 / numpy.sqrt(num_incoh))
643 673
644 674 # Remocion de Interferencia en el Cross Spectra
645 675 if jcspectra is None:
646 676 return jspectra, jcspectra
647 677 num_pairs = int(jcspectra.size / (num_prof * num_hei))
648 678 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
649 679
650 680 for ip in range(num_pairs):
651 681
652 682 #-------------------------------------------
653 683
654 684 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
655 685 cspower = cspower[:, hei_interf]
656 686 cspower = cspower.sum(axis=0)
657 687
658 688 cspsort = cspower.ravel().argsort()
659 689 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
660 690 offhei_interf, nhei_interf + offhei_interf))]]]
661 691 junkcspc_interf = junkcspc_interf.transpose()
662 692 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
663 693
664 694 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
665 695
666 696 median_real = int(numpy.median(numpy.real(
667 697 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
668 698 median_imag = int(numpy.median(numpy.imag(
669 699 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
670 700 comp_mask_prof = [int(e) for e in comp_mask_prof]
671 701 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
672 702 median_real, median_imag)
673 703
674 704 for iprof in range(num_prof):
675 705 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
676 706 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
677 707
678 708 # Removiendo la Interferencia
679 709 jcspectra[ip, :, ind_hei] = jcspectra[ip,
680 710 :, ind_hei] - jcspc_interf
681 711
682 712 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
683 713 maxid = ListAux.index(max(ListAux))
684 714
685 715 ind = numpy.array([-2, -1, 1, 2])
686 716 xx = numpy.zeros([4, 4])
687 717
688 718 for id1 in range(4):
689 719 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
690 720
691 721 xx_inv = numpy.linalg.inv(xx)
692 722 xx = xx_inv[:, 0]
693 723
694 724 ind = (ind + maxid + num_mask_prof) % num_mask_prof
695 725 yy = jcspectra[ip, mask_prof[ind], :]
696 726 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
697 727
698 728 # Guardar Resultados
699 729 self.dataOut.data_spc = jspectra
700 730 self.dataOut.data_cspc = jcspectra
701 731
702 732 return 1
703 733
704 734 def run(self, dataOut, interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None, mode=1):
705 735
706 736 self.dataOut = dataOut
707 737
708 738 if mode == 1:
709 739 self.removeInterference(interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None)
710 740 elif mode == 2:
711 741 self.removeInterference2()
712 742
713 743 return self.dataOut
714 744
715 745
716 746 class deflip(Operation):
717 747
718 748 def run(self, dataOut):
719 749 # arreglo 1: (num_chan, num_profiles, num_heights)
720 750 self.dataOut = dataOut
721 751
722 752 # JULIA-oblicua, indice 2
723 753 # arreglo 2: (num_profiles, num_heights)
724 754 jspectra = self.dataOut.data_spc[2]
725 755 jspectra_tmp=numpy.zeros(jspectra.shape)
726 756 num_profiles=jspectra.shape[0]
727 757 freq_dc = int(num_profiles / 2)
728 758 # Flip con for
729 759 for j in range(num_profiles):
730 760 jspectra_tmp[num_profiles-j-1]= jspectra[j]
731 761 # Intercambio perfil de DC con perfil inmediato anterior
732 762 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
733 763 jspectra_tmp[freq_dc]= jspectra[freq_dc]
734 764 # canal modificado es re-escrito en el arreglo de canales
735 765 self.dataOut.data_spc[2] = jspectra_tmp
736 766
737 767 return self.dataOut
738 768
739 769
740 770 class IncohInt(Operation):
741 771
742 772 __profIndex = 0
743 773 __withOverapping = False
744 774
745 775 __byTime = False
746 776 __initime = None
747 777 __lastdatatime = None
748 778 __integrationtime = None
749 779
750 780 __buffer_spc = None
751 781 __buffer_cspc = None
752 782 __buffer_dc = None
753 783
754 784 __dataReady = False
755 785
756 786 __timeInterval = None
757 787
758 788 n = None
759 789
760 790 def __init__(self):
761 791
762 792 Operation.__init__(self)
763 793
764 794 def setup(self, n=None, timeInterval=None, overlapping=False):
765 795 """
766 796 Set the parameters of the integration class.
767 797
768 798 Inputs:
769 799
770 800 n : Number of coherent integrations
771 801 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
772 802 overlapping :
773 803
774 804 """
775 805
776 806 self.__initime = None
777 807 self.__lastdatatime = 0
778 808
779 809 self.__buffer_spc = 0
780 810 self.__buffer_cspc = 0
781 811 self.__buffer_dc = 0
782 812
783 813 self.__profIndex = 0
784 814 self.__dataReady = False
785 815 self.__byTime = False
786 816
787 817 if n is None and timeInterval is None:
788 818 raise ValueError("n or timeInterval should be specified ...")
789 819
790 820 if n is not None:
791 821 self.n = int(n)
792 822 else:
793 823
794 824 self.__integrationtime = int(timeInterval)
795 825 self.n = None
796 826 self.__byTime = True
797 827
798 828 def putData(self, data_spc, data_cspc, data_dc):
799 829 """
800 830 Add a profile to the __buffer_spc and increase in one the __profileIndex
801 831
802 832 """
803 833
804 834 self.__buffer_spc += data_spc
805 835
806 836 if data_cspc is None:
807 837 self.__buffer_cspc = None
808 838 else:
809 839 self.__buffer_cspc += data_cspc
810 840
811 841 if data_dc is None:
812 842 self.__buffer_dc = None
813 843 else:
814 844 self.__buffer_dc += data_dc
815 845
816 846 self.__profIndex += 1
817 847
818 848 return
819 849
820 850 def pushData(self):
821 851 """
822 852 Return the sum of the last profiles and the profiles used in the sum.
823 853
824 854 Affected:
825 855
826 856 self.__profileIndex
827 857
828 858 """
829 859
830 860 data_spc = self.__buffer_spc
831 861 data_cspc = self.__buffer_cspc
832 862 data_dc = self.__buffer_dc
833 863 n = self.__profIndex
834 864
835 865 self.__buffer_spc = 0
836 866 self.__buffer_cspc = 0
837 867 self.__buffer_dc = 0
838 868 self.__profIndex = 0
839 869
840 870 return data_spc, data_cspc, data_dc, n
841 871
842 872 def byProfiles(self, *args):
843 873
844 874 self.__dataReady = False
845 875 avgdata_spc = None
846 876 avgdata_cspc = None
847 877 avgdata_dc = None
848 878
849 879 self.putData(*args)
850 880
851 881 if self.__profIndex == self.n:
852 882
853 883 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
854 884 self.n = n
855 885 self.__dataReady = True
856 886
857 887 return avgdata_spc, avgdata_cspc, avgdata_dc
858 888
859 889 def byTime(self, datatime, *args):
860 890
861 891 self.__dataReady = False
862 892 avgdata_spc = None
863 893 avgdata_cspc = None
864 894 avgdata_dc = None
865 895
866 896 self.putData(*args)
867 897
868 898 if (datatime - self.__initime) >= self.__integrationtime:
869 899 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
870 900 self.n = n
871 901 self.__dataReady = True
872 902
873 903 return avgdata_spc, avgdata_cspc, avgdata_dc
874 904
875 905 def integrate(self, datatime, *args):
876 906
877 907 if self.__profIndex == 0:
878 908 self.__initime = datatime
879 909
880 910 if self.__byTime:
881 911 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
882 912 datatime, *args)
883 913 else:
884 914 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
885 915
886 916 if not self.__dataReady:
887 917 return None, None, None, None
888 918
889 919 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
890 920
891 921 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
892 922 if n == 1:
893 923 return dataOut
894 924
895 925 dataOut.flagNoData = True
896 926
897 927 if not self.isConfig:
898 928 self.setup(n, timeInterval, overlapping)
899 929 self.isConfig = True
900 930
901 931 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
902 932 dataOut.data_spc,
903 933 dataOut.data_cspc,
904 934 dataOut.data_dc)
905 935
906 936 if self.__dataReady:
907 937
908 938 dataOut.data_spc = avgdata_spc
909 939 dataOut.data_cspc = avgdata_cspc
910 940 dataOut.data_dc = avgdata_dc
911 941 dataOut.nIncohInt *= self.n
912 942 dataOut.utctime = avgdatatime
913 943 dataOut.flagNoData = False
914 944
915 945 return dataOut
916 946
917 947 class dopplerFlip(Operation):
918 948
919 949 def run(self, dataOut):
920 950 # arreglo 1: (num_chan, num_profiles, num_heights)
921 951 self.dataOut = dataOut
922 952 # JULIA-oblicua, indice 2
923 953 # arreglo 2: (num_profiles, num_heights)
924 954 jspectra = self.dataOut.data_spc[2]
925 955 jspectra_tmp = numpy.zeros(jspectra.shape)
926 956 num_profiles = jspectra.shape[0]
927 957 freq_dc = int(num_profiles / 2)
928 958 # Flip con for
929 959 for j in range(num_profiles):
930 960 jspectra_tmp[num_profiles-j-1]= jspectra[j]
931 961 # Intercambio perfil de DC con perfil inmediato anterior
932 962 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
933 963 jspectra_tmp[freq_dc]= jspectra[freq_dc]
934 964 # canal modificado es re-escrito en el arreglo de canales
935 965 self.dataOut.data_spc[2] = jspectra_tmp
936 966
937 967 return self.dataOut No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now