##// END OF EJS Templates
schain
RSS
jrproc_parametes.py update y script EWDriftsMP_01feb2023.py
sebastianVP -
r1743:6b1a256569db
parent child
Show More
Side by Side Unified
Context file:
r1743:6b1a256569db Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,262
1
2 import os, sys
3 import json
4
5 #from controller import *
6 from schainpy.controller import Project
7
8 desc = "EW DRIFTS MP Experiment"
9 filename = "EWDrifts.xml"
10
11 controllerObj = Project()
12
13 controllerObj.setup(id = '191', name='test01', description=desc)
14
15 #Experimentos
16
17 #path = "/data/dia"
18 #path = '/home/pcondor/data'
19 #path = '/media/pcondor/DATA1/Database/ewdriftsschain2023prue/data'
20 #path = '/data/2024_01/MP_ISR/main_radar/rawdata/d2024023'
21 path = '/data/ISR_JULIA/d2024092'
22 #pathFigure = '/media/pcondor/DATA1/Database/ewdriftsschain2023wh5'
23 pathFile = '/media/pcondor/DATA1/Database/ewdriftsabr2024sch/EW_Drifts_01abr'
24 pathFigure = pathFile
25 pathFileavg = pathFile+'/avg'
26 pathFiledata = pathFile+'/Drifts-data'
27 #pathFileavg = '/media/pcondor/DATA1/Database/ewdriftsschain2023wh5/avg'
28 #pathFiledata = '/media/pcondor/DATA1/Database/ewdriftsschain2023wh5/Drifts-data'
29
30 xmin = 0
31 xmax = 24
32 #------------------------------------------------------------------------------------------------
33 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
34 path=path,
35 startDate='2024/04/01',
36 endDate='2024/04/01',
37 startTime='00:00:00',
38 endTime='23:59:59',
39 online=0,
40 getByBlock=1,
41 walk=0)
42
43 #--------------------------------------------------------------------------------------------------
44
45 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
46
47 #opObj11 = procUnitConfObj0.addOperation(name='selectHeights')
48 # # opObj11.addParameter(name='minHei', value='320.0', format='float')
49 # # opObj11.addParameter(name='maxHei', value='350.0', format='float')
50 #opObj11.addParameter(name='minHei', value='0.01', format='float')
51 #opObj11.addParameter(name='maxHei', value='960.0', format='float')
52
53 opObj11 = procUnitConfObj0.addOperation(name='selectChannels')
54 opObj11.addParameter(name='channelList', value='0,0,1,1', format='intlist')
55
56 #opObj11 = procUnitConfObj0.addOperation(name='Reshaper')
57 #opObj11.addParameter(name='shape', value='(500,980)', format='intlist')
58
59 opObj11 = procUnitConfObj0.addOperation(name='ProfileSelector', optype='other')
60 opObj11.addParameter(name='profileRangeList', value='0,127', format='intlist')
61
62 opObj11 = procUnitConfObj0.addOperation(name='filterByHeights')
63 opObj11.addParameter(name='window', value='10', format='int')
64
65 code=[[-1,-1,1],[1,1,-1]]
66 #code = [[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,-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,-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],[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,-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,-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],[-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,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,-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],[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,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,-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],[-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,-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,-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],[-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,-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,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],[1,1,-1],[-1,-1,1],[-1,-1,1],[1,1,-1],[-1,-1,1],[-1,-1,1],[1,1,-1],[1,1,-1]]
67 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
68 opObj11.addParameter(name='code', value=code, format='floatlist')
69 opObj11.addParameter(name='nCode', value='2', format='int')
70 opObj11.addParameter(name='nBaud', value='3', format='int')
71
72 opObj11 = procUnitConfObj0.addOperation(name='selectHeights')
73 opObj11.addParameter(name='minHei', value='0.0', format='float')
74 opObj11.addParameter(name='maxHei', value='960', format='float')
75
76 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
77 procUnitConfObj1.addParameter(name='nFFTPoints', value='128', format='int')
78 procUnitConfObj1.addParameter(name='nProfiles', value='128', format='int')
79 #procUnitConfObj1.addParameter(name='pairsList', value='(2,3),(4,5)', format='pairsList')#,(2,3)
80 procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(2,3)', format='pairsList')
81
82 #opObj11 = procUnitConfObj1.addOperation(name='selectHeights')
83 # # opObj11.addParameter(name='minHei', value='320.0', format='float')
84 # # opObj11.addParameter(name='maxHei', value='350.0', format='float')
85 #opObj11.addParameter(name='minHei', value='0.0', format='float')
86 #opObj11.addParameter(name='maxHei', value='960.0', format='float')
87
88 #opObj11 = procUnitConfObj1.addOperation(name='selectChannels')
89 #opObj11.addParameter(name='channelList', value='2,3,4,5', format='intlist')
90
91 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
92 opObj11.addParameter(name='n', value='1', format='float')
93 #opObj11.addParameter(name='timeInterval', value='300.0', format='float')
94
95 #opObj13 = procUnitConfObj1.addOperation(name='removeDC')
96
97 #opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
98 #opObj14.addParameter(name='id', value='65', format='int')
99 ## # opObj14.addParameter(name='wintitle', value='Con interf', format='str')
100 #opObj14.addParameter(name='save', value=pathFigure, format='str')
101 ##opObj14.addParameter(name='save_period', value=1, format='int')
102 #opObj14.addParameter(name='zmin', value='10', format='int')
103 #opObj14.addParameter(name='zmax', value='26', format='int')
104 #
105
106 #opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
107 #opObj12.addParameter(name='id', value='63', format='int')
108 #opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
109 #opObj12.addParameter(name='save', value=pathFigure, format='str')
110 #opObj12.addParameter(name='save_period', value=10, format='int')
111 ##opObj12.addParameter(name='figpath', value = pathFigure, format='str')
112 #opObj12.addParameter(name='xmin', value=xmin, format='float')
113 #opObj12.addParameter(name='xmax', value=xmax, format='float')
114 #opObj12.addParameter(name='zmin', value='20', format='int')
115 #opObj12.addParameter(name='zmax', value='36', format='int')
116
117 #--------------------------------------------------------------------------------------------------
118
119 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
120 opObj20 = procUnitConfObj2.addOperation(name='SpectralFitting', optype='other')
121 opObj20.addParameter(name='path', value='/home/pcondor/DIR_MADRIGAL/schain/schainpy/model/proc', format='str')
122 opObj20.addParameter(name='file', value='modelSpectralFitting', format='str')
123 opObj20.addParameter(name='groupList', value='(0,1),(2,3)',format='multiList')
124 opObj20.addParameter(name='taver', value='5')
125 opObj20.addParameter(name='coh_th', value='[1]',format='multiList')
126 opObj20.addParameter(name='hei_th', value='[2000]',format='multiList')
127 #opObj20.addParameter(name='filec', value='weightfit', format='str')
128
129 opObj22 = procUnitConfObj2.addOperation(name='HDFWriter', optype='other')
130 opObj22.addParameter(name='path', value=pathFiledata)
131 opObj22.addParameter(name='blocksPerFile', value='1')
132 opObj22.addParameter(name='metadataList',value='heightList,timeZone')
133 opObj22.addParameter(name='dataList',value='tmp_spectra_i,tmp_cspectra_i,tmp_spectra_c,tmp_cspectra_c,clean_num_aver,coh_num_aver,sat_spectra,sat_cspectra,index,utctime,nIncohInt,nCohInt,nProfiles,nFFTPoints,ippFactor,ippSeconds,paramInterval')
134 ##opObj22.addParameter(name='dataList',value='tmp_spectra_i,tmp_cspectra_i,tmp_spectra_c,tmp_cspectra_c,clean_num_aver,coh_num_aver,index,utctime,nIncohInt,nCohInt,nProfiles,nFFTPoints,normFactor,channelList,ippFactor,ippSeconds')
135
136 #angles :-2.41116 3.01082
137 opObj21 = procUnitConfObj2.addOperation(name='EWDriftsEstimation', optype='other')
138 opObj21.addParameter(name='zenith', value='-2.41116, 3.01082', format='floatlist')
139 opObj21.addParameter(name='zenithCorrection', value='0.0', format='float')
140 opObj21.addParameter(name='fileDrifts', value=pathFile)
141
142 # Drifts en h5
143 one = {'gdlatr': 'lat', 'gdlonr': 'lon', 'spcst':'spcst','pl':'pl','cbadn':'cbadn','inttms': 'inttms','azdir7':'azw','eldir7':'elw','azdir8':'aze','eldir8':'ele','jro14':'jro14','jro15':'jro15','jro16':'jro16','nwlos':'nwlos'}
144 two = {
145 'range': ('params', 0),
146 'gdalt': ('params', 1),
147 'VIPN': ('params', 2),
148 'dvipn': ('params', 3),
149 'vipe': ('params', 4),
150 'dvipe': ('params', 5),
151 'vi7': ('params', 6),
152 'dvi7': ('params', 7),
153 'vi8': ('params', 8),
154 'dvi8': ('params', 9),
155 'PAIWL': ('params', 10),
156 'pacwl': ('params', 11),
157 'pbiwl': ('params', 12),
158 'pbcwl': ('params', 13),
159 'pciel': ('params', 14),
160 'pccel': ('params', 15),
161 'pdiel': ('params', 16),
162 'pdcel': ('params', 17),
163 'jro10': ('params', 18),
164 'jro11': ('params', 19)
165 } #writer
166 ind = ['gdalt']
167
168 meta = {
169 'kinst': 10, #instrument code
170 'kindat': 1910, #type of data
171 'catalog': {
172 'principleInvestigator': 'Danny ScipiΓ³n',
173 'expPurpose': 'Drifts'#,
174 #'sciRemarks': file_contents
175 },
176 'header': {
177 'analyst': 'Danny ScipiΓ³n'
178 }
179 }
180
181 op_writer = procUnitConfObj2.addOperation(name='MADWriter')
182 op_writer.addParameter(name='path', value=pathFile)
183 op_writer.addParameter(name='format', value='hdf5')
184 op_writer.addParameter(name='oneDDict', value=json.dumps(one))
185 op_writer.addParameter(name='twoDDict', value=json.dumps(two))
186 op_writer.addParameter(name='ind2DList', value=json.dumps(ind))
187 op_writer.addParameter(name='metadata', value=json.dumps(meta))
188
189 op_writer = procUnitConfObj2.addOperation(name='setHeightDriftsavg')
190
191 # Avg Drifts
192 one_avg = {'gdlatr': 'lat', 'gdlonr': 'lon', 'spcst':'spcst','pl':'pl','cbadn':'cbadn','inttms': 'inttms'}
193 two_avg = {
194 'range': ('params_avg', 4),
195 'gdalt': ('params_avg', 5),
196 'altav': ('params_avg', 6),
197 'VIPN': ('params_avg', 0),
198 'dvipn': ('params_avg', 1),
199 'vipe': ('params_avg', 2),
200 'dvipe': ('params_avg', 3)
201 }
202 ind_avg = ['gdalt']
203 meta = {
204 'kinst': 10, #instrument code
205 'kindat': 1911, #type of data
206 'catalog': {
207 'principleInvestigator': 'Danny ScipiΓ³n',
208 'expPurpose': 'Drifts'#,
209 #'sciRemarks': file_contents
210 },
211 'header': {
212 'analyst': 'Danny ScipiΓ³n'
213 }
214 }
215 #dataOut.heightList = dataOut.params_avg[4]
216 op_writer = procUnitConfObj2.addOperation(name='MADWriter')
217 op_writer.addParameter(name='path', value=pathFileavg)
218 op_writer.addParameter(name='format', value='hdf5')
219 op_writer.addParameter(name='oneDDict', value=json.dumps(one_avg))
220 op_writer.addParameter(name='twoDDict', value=json.dumps(two_avg))
221 op_writer.addParameter(name='ind2DList', value=json.dumps(ind_avg))
222 op_writer.addParameter(name='metadata', value=json.dumps(meta))
223
224 op_writer = procUnitConfObj2.addOperation(name='setHeightDrifts')
225
226 opObj24 = procUnitConfObj2.addOperation(name='SpectralMomentsPlot', optype='other')
227 opObj24.addParameter(name='id', value='1', format='int')
228 ### # opObj14.addParameter(name='wintitle', value='Spectral Averaged', format='str')
229 opObj24.addParameter(name='save', value=pathFigure, format='str')
230 ###opObj24.addParameter(name='save_period', value=1, format='int')
231 opObj24.addParameter(name='zmin', value='-8', format='int')
232 opObj24.addParameter(name='zmax', value='16', format='int')
233 opObj24.addParameter(name='xaxis', value='Velocity', format='str')
234
235 #
236 titles=('SNR,Vertical Drifts,Zonal Drifts')
237 #titles=('Zonal Drifts,Vertical Drifts')
238 opObj23 = procUnitConfObj2.addOperation(name='GenericRTIPlot')
239 opObj23.addParameter(name='colormaps', value='jet,RdBu_r,RdBu_r')
240 opObj23.addParameter(name='attr_data', value='data_snr1,data_output')
241 #opObj23.addParameter(name='colormaps', value='RdBu,RdBu')
242 #opObj23.addParameter(name='attr_data', value='data_output')
243 opObj23.addParameter(name='wintitle', value='EW Drifts')
244 opObj23.addParameter(name='save', value=pathFigure)
245 opObj23.addParameter(name='titles', value=titles)
246 opObj23.addParameter(name='zfactors', value='1,1,1')
247 opObj23.addParameter(name='zlimits', value='(-5,20),(-50,50),(-150,150)')
248 opObj23.addParameter(name='cb_labels', value='dB,m/s,m/s')
249 #opObj23.addParameter(name='titles', value=titles)
250 #opObj23.addParameter(name='zfactors', value='1,1')
251 #opObj23.addParameter(name='zlimits', value='(-150,150),(-40,40)')
252 #opObj23.addParameter(name='cb_labels', value='m/s,m/s')
253 opObj23.addParameter(name='throttle', value='1')
254 opObj23.addParameter(name='xmin', value=xmin)
255 opObj23.addParameter(name='xmax', value=xmax)
256 #opObj23.addParameter(name='exp_code', value='110', format='int')
257 #opObj23.addParameter(name='server', value='10.10.110.243:4444', format='int')
258 #opObj23.addParameter(name='tag', value= 'jicamarca', format='str')
259
260 #--------------------------------------------------------------------------------------------------
261
262 controllerObj.start()
Show file before | Show file after | Hide comments
@@ -0,0 +1,186
1
2 import os, sys
3 import json
4 #from controller import *
5 from schainpy.controller import Project
6
7 desc = "EW DRIFTS MP Experiment"
8 filename = "EWDrifts.xml"
9
10 controllerObj = Project()
11
12 controllerObj.setup(id = '191', name='test01', description=desc)
13
14 #Experimentos
15
16 #path = '/media/pcondor/DATA1/Database/ewdriftsene2024sch/EW_Drifts_01ene/Drifts-data'
17 path = '/media/soporte/DATA/PERCY_SCHAIN_UPDATE/driftsschain'
18 #pathFigure = '/media/pcondor/DATA1/Database/ewdriftsschain2023proc'
19 pathFile ='/media/soporte/DATA/PERCY_SCHAIN_UPDATE/driftsschain/tmp'
20 #pathFile = '/media/pcondor/DATA1/Database/ewdriftsene2024sch/EW_Drifts_01enetmp'
21 pathFigure = pathFile
22 pathFileavg = pathFile+'/avg'
23 pathFiledata = pathFile+'/Drifts-data'
24
25 xmin = 0
26 xmax = 24
27 #------------------------------------------------------------------------------------------------
28 readUnitConfObj = controllerObj.addReadUnit(datatype='HDFReader',
29 path=path,
30 startDate='2024/01/23',
31 endDate='2024/01/23',
32 startTime='00:00:00',
33 endTime='23:59:59',
34 #online=0,
35 #getByBlock=1,
36 walk=1,
37 utcoffset='-18000')
38
39 #--------------------------------------------------------------------------------------------------
40
41 #--------------------------------------------------------------------------------------------------
42
43 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=readUnitConfObj.getId())
44
45 opObj20 = procUnitConfObj2.addOperation(name='SpectralFitting', optype='other')
46 opObj20.addParameter(name='path', value='/home/pcondor/DIR_MADRIGAL/schain/schainpy/model/proc', format='str')
47 opObj20.addParameter(name='file', value='modelSpectralFitting', format='str')
48 opObj20.addParameter(name='groupList', value='(0,1),(2,3)',format='multiList')
49 opObj20.addParameter(name='taver', value='5')
50 opObj20.addParameter(name='coh_th', value='[1]',format='multiList')
51 opObj20.addParameter(name='hei_th', value='[2000]',format='multiList')
52 opObj20.addParameter(name='proc', value='1')
53 opObj20.addParameter(name='channelList', value='0,0,1,1')
54 opObj20.addParameter(name='filec', value='weightfit', format='str')
55
56 #opObj22 = procUnitConfObj2.addOperation(name='HDFWriter', optype='other')
57 #opObj22.addParameter(name='path', value=pathFile)
58 #opObj22.addParameter(name='blocksPerFile', value='1')
59 #opObj22.addParameter(name='metadataList',value='heightList,timeZone')
60 #opObj22.addParameter(name='dataList',value='tmp_spectra_i,tmp_cspectra_i,tmp_spectra_c,tmp_cspectra_c,clean_num_aver,coh_num_aver,index,utctime')
61 #angles :-2.41116 3.01082
62 opObj21 = procUnitConfObj2.addOperation(name='EWDriftsEstimation', optype='other')
63 opObj21.addParameter(name='zenith', value='-2.41116, 3.01082', format='floatlist')
64 opObj21.addParameter(name='zenithCorrection', value='0.0', format='float')
65 opObj21.addParameter(name='fileDrifts', value=pathFile)
66
67 # Drifts en h5
68 one = {'gdlatr': 'lat', 'gdlonr': 'lon', 'spcst':'spcst','pl':'pl','cbadn':'cbadn','inttms': 'inttms','azdir7':'azw','eldir7':'elw','azdir8':'aze','eldir8':'ele','jro14':'jro14','jro15':'jro15','jro16':'jro16','nwlos':'nwlos'}
69 two = {
70 'range': ('params', 0),
71 'gdalt': ('params', 1),
72 'VIPN': ('params', 2),
73 'dvipn': ('params', 3),
74 'vipe': ('params', 4),
75 'dvipe': ('params', 5),
76 'vi7': ('params', 6),
77 'dvi7': ('params', 7),
78 'vi8': ('params', 8),
79 'dvi8': ('params', 9),
80 'PAIWL': ('params', 10),
81 'pacwl': ('params', 11),
82 'pbiwl': ('params', 12),
83 'pbcwl': ('params', 13),
84 'pciel': ('params', 14),
85 'pccel': ('params', 15),
86 'pdiel': ('params', 16),
87 'pdcel': ('params', 17),
88 'jro10': ('params', 18),
89 'jro11': ('params', 19)
90 } #writer
91 ind = ['gdalt']
92
93 #f=open('/home/roberto/moder_test.txt','r')
94 #file_contents=f.read()
95
96 meta = {
97 'kinst': 10, #instrument code
98 'kindat': 1910, #type of data
99 'catalog': {
100 'principleInvestigator': 'Danny ScipiΓ³n',
101 'expPurpose': 'Drifts'#,
102 #'sciRemarks': file_contents
103 },
104 'header': {
105 'analyst': 'Danny ScipiΓ³n'
106 }
107 }
108 #f.close()
109
110 op_writer = procUnitConfObj2.addOperation(name='MADWriter')
111 op_writer.addParameter(name='path', value=pathFile)
112 op_writer.addParameter(name='format', value='hdf5')
113 op_writer.addParameter(name='oneDDict', value=json.dumps(one))
114 op_writer.addParameter(name='twoDDict', value=json.dumps(two))
115 op_writer.addParameter(name='ind2DList', value=json.dumps(ind))
116 op_writer.addParameter(name='metadata', value=json.dumps(meta))
117
118 op_writer = procUnitConfObj2.addOperation(name='setHeightDriftsavg')
119
120 # Avg Drifts
121 one_avg = {'gdlatr': 'lat', 'gdlonr': 'lon', 'spcst':'spcst','pl':'pl','cbadn':'cbadn','inttms': 'inttms'}
122 two_avg = {
123 'range': ('params_avg', 4),
124 'gdalt': ('params_avg', 5),
125 'altav': ('params_avg', 6),
126 'VIPN': ('params_avg', 0),
127 'dvipn': ('params_avg', 1),
128 'vipe': ('params_avg', 2),
129 'dvipe': ('params_avg', 3)
130 }
131 ind_avg = ['gdalt']
132 meta = {
133 'kinst': 10, #instrument code
134 'kindat': 1911, #type of data
135 'catalog': {
136 'principleInvestigator': 'Danny ScipiΓ³n',
137 'expPurpose': 'Drifts'#,
138 #'sciRemarks': file_contents
139 },
140 'header': {
141 'analyst': 'Danny ScipiΓ³n'
142 }
143 }
144
145 op_writer = procUnitConfObj2.addOperation(name='MADWriter')
146 op_writer.addParameter(name='path', value=pathFileavg)
147 op_writer.addParameter(name='format', value='hdf5')
148 op_writer.addParameter(name='oneDDict', value=json.dumps(one_avg))
149 op_writer.addParameter(name='twoDDict', value=json.dumps(two_avg))
150 op_writer.addParameter(name='ind2DList', value=json.dumps(ind_avg))
151 op_writer.addParameter(name='metadata', value=json.dumps(meta))
152
153 op_writer = procUnitConfObj2.addOperation(name='setHeightDrifts')
154
155 opObj24 = procUnitConfObj2.addOperation(name='SpectralMomentsPlot', optype='other')
156 opObj24.addParameter(name='id', value='1', format='int')
157 ### # opObj14.addParameter(name='wintitle', value='Spectral Averaged', format='str')
158 opObj24.addParameter(name='save', value=pathFigure, format='str')
159 ###opObj24.addParameter(name='save_period', value=1, format='int')
160 opObj24.addParameter(name='zmin', value='-8', format='int')
161 opObj24.addParameter(name='zmax', value='16', format='int')
162 opObj24.addParameter(name='xaxis', value='Velocity', format='str')
163
164 #
165 titles=('SNR,Vertical Drifts,Zonal Drifts')
166 opObj23 = procUnitConfObj2.addOperation(name='GenericRTIPlot')
167 #opObj23.addParameter(name='colormaps', value='jet,RdBu_r,RdBu_r')
168 opObj23.addParameter(name='colormaps', value='jro,seismic,seismic')
169 #opObj23.addParameter(name='colormaps', value='jro,bwr,bwr')
170 opObj23.addParameter(name='attr_data', value='data_snr1,data_output')
171 opObj23.addParameter(name='wintitle', value='EW Drifts')
172 opObj23.addParameter(name='save', value=pathFigure)
173 opObj23.addParameter(name='titles', value=titles)
174 opObj23.addParameter(name='zfactors', value='1,1,1')
175 opObj23.addParameter(name='zlimits', value='(0,13),(-50,50),(-150,150)')
176 opObj23.addParameter(name='cb_labels', value='dB,m/s,m/s')
177 opObj23.addParameter(name='throttle', value='1')
178 opObj23.addParameter(name='xmin', value=xmin)
179 opObj23.addParameter(name='xmax', value=xmax)
180 #opObj23.addParameter(name='exp_code', value='110', format='int')
181 #opObj23.addParameter(name='server', value='10.10.110.243:4444', format='int')
182 #opObj23.addParameter(name='tag', value= 'jicamarca', format='str')
183
184 #--------------------------------------------------------------------------------------------------
185
186 controllerObj.start()
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (2380 lines changed) Show them Hide them
@@ -1,3 +1,4
1 # MASTER
1 import numpy
2 import numpy
2 import math
3 import math
3 from scipy import optimize, interpolate, signal, stats, ndimage
4 from scipy import optimize, interpolate, signal, stats, ndimage
@@ -17,7 +18,7 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papamete
17 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
18 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
18 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
19 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
19 from schainpy.model.data.jrodata import Spectra
20 from schainpy.model.data.jrodata import Spectra
20 from scipy import asarray as ar,exp
21 #from scipy import asarray as ar,exp
21 from scipy.optimize import fmin, curve_fit
22 from scipy.optimize import fmin, curve_fit
22 from schainpy.utils import log
23 from schainpy.utils import log
23 import warnings
24 import warnings
@@ -180,7 +181,6 class ParametersProc(ProcessingUnit):
180 self.__updateObjFromInput()
181 self.__updateObjFromInput()
181 self.dataOut.utctimeInit = self.dataIn.utctime
182 self.dataOut.utctimeInit = self.dataIn.utctime
182 self.dataOut.paramInterval = self.dataIn.timeInterval
183 self.dataOut.paramInterval = self.dataIn.timeInterval
183
184 return
184 return
185
185
186
186
@@ -242,7 +242,6 class RemoveWideGC(Operation):
242 continue
242 continue
243 junk3 = numpy.squeeze(numpy.diff(j1index))
243 junk3 = numpy.squeeze(numpy.diff(j1index))
244 junk4 = numpy.squeeze(numpy.diff(j2index))
244 junk4 = numpy.squeeze(numpy.diff(j2index))
245
246 valleyindex = j2index[numpy.where(junk4>1)]
245 valleyindex = j2index[numpy.where(junk4>1)]
247 peakindex = j1index[numpy.where(junk3>1)]
246 peakindex = j1index[numpy.where(junk3>1)]
248
247
@@ -252,7 +251,6 class RemoveWideGC(Operation):
252 if numpy.size(isvalid) >1 :
251 if numpy.size(isvalid) >1 :
253 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
252 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
254 isvalid = isvalid[vindex]
253 isvalid = isvalid[vindex]
255
256 # clutter peak
254 # clutter peak
257 gcpeak = peakindex[isvalid]
255 gcpeak = peakindex[isvalid]
258 vl = numpy.where(valleyindex < gcpeak)
256 vl = numpy.where(valleyindex < gcpeak)
@@ -302,8 +300,6 class SpectralFilters(Operation):
302 VelRange = dataOut.spc_range[2]
300 VelRange = dataOut.spc_range[2]
303
301
304 # novalid corresponds to data within the Negative and PositiveLimit
302 # novalid corresponds to data within the Negative and PositiveLimit
305
306
307 # Removing novalid data from the spectra
303 # Removing novalid data from the spectra
308 for i in range(self.Num_Chn):
304 for i in range(self.Num_Chn):
309 self.spc[i,novalid,:] = dataOut.noise[i]
305 self.spc[i,novalid,:] = dataOut.noise[i]
@@ -435,7 +431,6 class GaussianFit(Operation):
435 # print ('stop 2.1')
431 # print ('stop 2.1')
436 fatspectra=1.0
432 fatspectra=1.0
437 # noise per channel.... we might want to use the noise at each range
433 # noise per channel.... we might want to use the noise at each range
438
439 # wnoise = noise_ #/ spc_norm_max #commented by D. ScipiΓ³n 19.03.2021
434 # wnoise = noise_ #/ spc_norm_max #commented by D. ScipiΓ³n 19.03.2021
440 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
435 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
441 #if wnoise>1.1*pnoise: # to be tested later
436 #if wnoise>1.1*pnoise: # to be tested later
@@ -617,7 +612,7 class GaussianFit(Operation):
617 if Amplitude1<0.05:
612 if Amplitude1<0.05:
618 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
613 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
619
614
620 # print ('stop 16 ')
615 # print ('stop 16 ')
621 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0)/width0)**p0)
616 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0)/width0)**p0)
622 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1)/width1)**p1)
617 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1)/width1)**p1)
623 # SPCparam = (SPC_ch1,SPC_ch2)
618 # SPCparam = (SPC_ch1,SPC_ch2)
@@ -714,7 +709,6 class Oblique_Gauss_Fit(Operation):
714
709
715 return gaussian(x, popt[0], popt[1], popt[2], popt[3]), popt[0], popt[1], popt[2], popt[3]
710 return gaussian(x, popt[0], popt[1], popt[2], popt[3]), popt[0], popt[1], popt[2], popt[3]
716
711
717
718 def Gauss_fit_2(self,spc,x,nGauss):
712 def Gauss_fit_2(self,spc,x,nGauss):
719
713
720
714
@@ -1099,7 +1093,6 class Oblique_Gauss_Fit(Operation):
1099
1093
1100 # fit
1094 # fit
1101 bounds=([0,-numpy.inf,0,0,-400,0,0,0],[numpy.inf,-340,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1095 bounds=([0,-numpy.inf,0,0,-400,0,0,0],[numpy.inf,-340,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1102
1103 params_scale = [spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,1,spc_max]
1096 params_scale = [spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,1,spc_max]
1104 x0_value = numpy.array([spc_max,-400,30,spc_max/4,-200,150,1,1.0e7])
1097 x0_value = numpy.array([spc_max,-400,30,spc_max/4,-200,150,1,1.0e7])
1105 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1098 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
@@ -1881,8 +1874,7 class PrecipitationProc(Operation):
1881 self.i=0
1874 self.i=0
1882
1875
1883 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
1876 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
1884 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350,SNRdBlimit=-30,
1877 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350,SNRdBlimit=-30,channel=None):
1885 channel=None):
1886
1878
1887 # print ('Entering PrecepitationProc ... ')
1879 # print ('Entering PrecepitationProc ... ')
1888
1880
@@ -1915,7 +1907,7 class PrecipitationProc(Operation):
1915 self.Lambda = Lambda
1907 self.Lambda = Lambda
1916 self.aL = aL
1908 self.aL = aL
1917 self.tauW = tauW
1909 self.tauW = tauW
1918 self.ThetaT = ThetaT
1910 self.ThetaT = ThetaT
1919 self.ThetaR = ThetaR
1911 self.ThetaR = ThetaR
1920 self.GSys = 10**(36.63/10) # Ganancia de los LNA 36.63 dB
1912 self.GSys = 10**(36.63/10) # Ganancia de los LNA 36.63 dB
1921 self.lt = 10**(1.67/10) # Perdida en cables Tx 1.67 dB
1913 self.lt = 10**(1.67/10) # Perdida en cables Tx 1.67 dB
@@ -1958,7 +1950,7 class PrecipitationProc(Operation):
1958 ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
1950 ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
1959 ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
1951 ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
1960 # Radar Cross Section
1952 # Radar Cross Section
1961 sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
1953 sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
1962 # Drop Size Distribution
1954 # Drop Size Distribution
1963 DSD = ETAn / sigmaD
1955 DSD = ETAn / sigmaD
1964 # Equivalente Reflectivy
1956 # Equivalente Reflectivy
@@ -1979,7 +1971,7 class PrecipitationProc(Operation):
1979 dataOut.data_output = RR[8]
1971 dataOut.data_output = RR[8]
1980 dataOut.data_param = numpy.ones([3,self.Num_Hei])
1972 dataOut.data_param = numpy.ones([3,self.Num_Hei])
1981 dataOut.channelList = [0,1,2]
1973 dataOut.channelList = [0,1,2]
1982
1974
1983 dataOut.data_param[0]=10*numpy.log10(Ze_org)
1975 dataOut.data_param[0]=10*numpy.log10(Ze_org)
1984 dataOut.data_param[1]=-W
1976 dataOut.data_param[1]=-W
1985 dataOut.data_param[2]=RR
1977 dataOut.data_param[2]=RR
@@ -2054,8 +2046,7 class FullSpectralAnalysis(Operation):
2054 Parameters affected: Winds, height range, SNR
2046 Parameters affected: Winds, height range, SNR
2055
2047
2056 """
2048 """
2057 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRdBlimit=-30,
2049 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRdBlimit=-30,minheight=None, maxheight=None, NegativeLimit=None, PositiveLimit=None):
2058 minheight=None, maxheight=None, NegativeLimit=None, PositiveLimit=None):
2059
2050
2060 spc = dataOut.data_pre[0].copy()
2051 spc = dataOut.data_pre[0].copy()
2061 cspc = dataOut.data_pre[1]
2052 cspc = dataOut.data_pre[1]
@@ -2098,14 +2089,14 class FullSpectralAnalysis(Operation):
2098
2089
2099 if Height >= range_min and Height < range_max:
2090 if Height >= range_min and Height < range_max:
2100 # error_code will be useful in future analysis
2091 # error_code will be useful in future analysis
2101 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList,
2092 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList,
2102 ChanDist, Height, dataOut.noise, dataOut.spc_range, dbSNR[Height], SNRdBlimit, NegativeLimit, PositiveLimit,dataOut.frequency)
2093 ChanDist, Height, dataOut.noise, dataOut.spc_range, dbSNR[Height], SNRdBlimit, NegativeLimit, PositiveLimit,dataOut.frequency)
2103
2094
2104 if abs(Vzon) < 100. and abs(Vmer) < 100.:
2095 if abs(Vzon) < 100. and abs(Vmer) < 100.:
2105 velocityX[Height] = Vzon
2096 velocityX[Height] = Vzon
2106 velocityY[Height] = -Vmer
2097 velocityY[Height] = -Vmer
2107 velocityZ[Height] = Vver
2098 velocityZ[Height] = Vver
2108
2099
2109 # Censoring data with SNR threshold
2100 # Censoring data with SNR threshold
2110 dbSNR [dbSNR < SNRdBlimit] = numpy.NaN
2101 dbSNR [dbSNR < SNRdBlimit] = numpy.NaN
2111
2102
@@ -2205,7 +2196,7 class FullSpectralAnalysis(Operation):
2205 xSamples = xFrec # the frequency range is taken
2196 xSamples = xFrec # the frequency range is taken
2206 delta_x = xSamples[1] - xSamples[0] # delta_f or delta_x
2197 delta_x = xSamples[1] - xSamples[0] # delta_f or delta_x
2207
2198
2208 # only consider velocities with in NegativeLimit and PositiveLimit
2199 # only consider velocities with in NegativeLimit and PositiveLimit
2209 if (NegativeLimit is None):
2200 if (NegativeLimit is None):
2210 NegativeLimit = numpy.min(xVel)
2201 NegativeLimit = numpy.min(xVel)
2211 if (PositiveLimit is None):
2202 if (PositiveLimit is None):
@@ -2622,7 +2613,7 class SpectralMoments(Operation):
2622
2613
2623 if (type1 is None): type1 = 0
2614 if (type1 is None): type1 = 0
2624 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
2615 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
2625 if (snrth is None): snrth = -20.0
2616 if (snrth is None): snrth = -3 #-20.0
2626 if (dc is None): dc = 0
2617 if (dc is None): dc = 0
2627 if (aliasing is None): aliasing = 0
2618 if (aliasing is None): aliasing = 0
2628 if (oldfd is None): oldfd = 0
2619 if (oldfd is None): oldfd = 0
@@ -3146,76 +3137,20 class SpectralFitting(Operation):
3146 __dataReady = False
3137 __dataReady = False
3147 bloques = None
3138 bloques = None
3148 bloque0 = None
3139 bloque0 = None
3149 index = 0
3150 fint = 0
3151 buffer = 0
3152 buffer2 = 0
3153 buffer3 = 0
3154
3140
3155 def __init__(self):
3141 def __init__(self):
3156 Operation.__init__(self)
3142 Operation.__init__(self)
3157 self.i=0
3143 self.i=0
3158 self.isConfig = False
3144 self.isConfig = False
3159
3145
3160
3146 def setup(self,nChan,nProf,nHei,nBlocks):
3161 def setup(self,dataOut,groupList,path,file,filec):
3162 self.__dataReady = False
3147 self.__dataReady = False
3163 # new
3148 self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex)
3164 self.nChannels = dataOut.nChannels
3149 self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks])
3165 self.channels = dataOut.channelList
3166 self.nHeights = dataOut.heightList.size
3167 self.heights = dataOut.heightList
3168 self.nProf = dataOut.nProfiles
3169 self.nIncohInt = dataOut.nIncohInt
3170 self.absc = dataOut.abscissaList[:-1]
3171
3172
3173 #To be inserted as a parameter
3174 try:
3175 self.groupArray = numpy.array(groupList)#groupArray = numpy.array([[0,1],[2,3]])
3176 except:
3177 print("Please insert groupList. Example (0,1),(2,3) format multilist")
3178 dataOut.groupList = self.groupArray
3179 self.crosspairs = dataOut.groupList
3180 self.nPairs = len(self.crosspairs)
3181 self.nGroups = self.groupArray.shape[0]
3182
3183 #List of possible combinations
3184
3185 self.listComb = itertools.combinations(numpy.arange(self.groupArray.shape[1]),2)
3186 self.indCross = numpy.zeros(len(list(self.listComb)), dtype = 'int')
3187
3188 #Parameters Array
3189 dataOut.data_param = None
3190 dataOut.data_paramC = None
3191 dataOut.clean_num_aver = None
3192 dataOut.coh_num_aver = None
3193 dataOut.tmp_spectra_i = None
3194 dataOut.tmp_cspectra_i = None
3195 dataOut.tmp_spectra_c = None
3196 dataOut.tmp_cspectra_c = None
3197 dataOut.index = None
3198
3199 if path != None:
3200 sys.path.append(path)
3201 self.library = importlib.import_module(file)
3202 if filec != None:
3203 self.weightf = importlib.import_module(filec)
3204
3205 #Set constants
3206 self.constants = self.library.setConstants(dataOut)
3207 dataOut.constants = self.constants
3208 self.M = dataOut.normFactor
3209 self.N = dataOut.nFFTPoints
3210 self.ippSeconds = dataOut.ippSeconds
3211 self.K = dataOut.nIncohInt
3212 self.pairsArray = numpy.array(dataOut.pairsList)
3213 self.snrth = 20
3214
3150
3215 def __calculateMoments(self,oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
3151 def __calculateMoments(self,oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
3216
3217 if (nicoh is None): nicoh = 1
3152 if (nicoh is None): nicoh = 1
3218 if (graph is None): graph = 0
3153 if (graph is None): graph = 0
3219 if (smooth is None): smooth = 0
3154 if (smooth is None): smooth = 0
3220 elif (self.smooth < 3): smooth = 0
3155 elif (self.smooth < 3): smooth = 0
3221
3156
@@ -3226,65 +3161,61 class SpectralFitting(Operation):
3226 if (aliasing is None): aliasing = 0
3161 if (aliasing is None): aliasing = 0
3227 if (oldfd is None): oldfd = 0
3162 if (oldfd is None): oldfd = 0
3228 if (wwauto is None): wwauto = 0
3163 if (wwauto is None): wwauto = 0
3229
3164
3230 if (n0 < 1.e-20): n0 = 1.e-20
3165 if (n0 < 1.e-20): n0 = 1.e-20
3231
3232 freq = oldfreq
3166 freq = oldfreq
3233 vec_power = numpy.zeros(oldspec.shape[1])
3167 vec_power = numpy.zeros(oldspec.shape[1])
3234 vec_fd = numpy.zeros(oldspec.shape[1])
3168 vec_fd = numpy.zeros(oldspec.shape[1])
3235 vec_w = numpy.zeros(oldspec.shape[1])
3169 vec_w = numpy.zeros(oldspec.shape[1])
3236 vec_snr = numpy.zeros(oldspec.shape[1])
3170 vec_snr = numpy.zeros(oldspec.shape[1])
3237
3238 oldspec = numpy.ma.masked_invalid(oldspec)
3171 oldspec = numpy.ma.masked_invalid(oldspec)
3239
3172
3240 for ind in range(oldspec.shape[1]):
3173 for ind in range(oldspec.shape[1]):
3241
3242 spec = oldspec[:,ind]
3174 spec = oldspec[:,ind]
3243 aux = spec*fwindow
3175 aux = spec*fwindow
3244 max_spec = aux.max()
3176 max_spec = aux.max()
3245 m = list(aux).index(max_spec)
3177 m = list(aux).index(max_spec)
3246
3178 #Smooth
3247 #Smooth
3248 if (smooth == 0): spec2 = spec
3179 if (smooth == 0): spec2 = spec
3249 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
3180 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
3250
3181
3251 # Calculo de Momentos
3182 # Calculo de Momentos
3252 bb = spec2[list(range(m,spec2.size))]
3183 bb = spec2[list(range(m,spec2.size))]
3253 bb = (bb<n0).nonzero()
3184 bb = (bb<n0).nonzero()
3254 bb = bb[0]
3185 bb = bb[0]
3255
3186
3256 ss = spec2[list(range(0,m + 1))]
3187 ss = spec2[list(range(0,m + 1))]
3257 ss = (ss<n0).nonzero()
3188 ss = (ss<n0).nonzero()
3258 ss = ss[0]
3189 ss = ss[0]
3259
3190
3260 if (bb.size == 0):
3191 if (bb.size == 0):
3261 bb0 = spec.size - 1 - m
3192 bb0 = spec.size - 1 - m
3262 else:
3193 else:
3263 bb0 = bb[0] - 1
3194 bb0 = bb[0] - 1
3264 if (bb0 < 0):
3195 if (bb0 < 0):
3265 bb0 = 0
3196 bb0 = 0
3266
3197
3267 if (ss.size == 0): ss1 = 1
3198 if (ss.size == 0): ss1 = 1
3268 else: ss1 = max(ss) + 1
3199 else: ss1 = max(ss) + 1
3269
3200
3270 if (ss1 > m): ss1 = m
3201 if (ss1 > m): ss1 = m
3271
3202
3272 valid = numpy.asarray(list(range(int(m + bb0 - ss1 + 1)))) + ss1
3203 valid = numpy.asarray(list(range(int(m + bb0 - ss1 + 1)))) + ss1
3273 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
3204 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
3274 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
3205 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
3275 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
3206 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
3276 snr = (spec2.mean()-n0)/n0
3207 snr = (spec2.mean()-n0)/n0
3277
3208
3278 if (snr < 1.e-20) :
3209 if (snr < 1.e-20) :
3279 snr = 1.e-20
3210 snr = 1.e-20
3280
3211
3281 vec_power[ind] = power
3212 vec_power[ind] = power
3282 vec_fd[ind] = fd
3213 vec_fd[ind] = fd
3283 vec_w[ind] = w
3214 vec_w[ind] = w
3284 vec_snr[ind] = snr
3215 vec_snr[ind] = snr
3285
3216
3286 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
3217 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
3287 return moments
3218 return moments
3288
3219
3289 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
3220 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
3290
3221
@@ -3305,12 +3236,12 class SpectralFitting(Operation):
3305 coh_spectra = numpy.zeros([nChan, nProf, nHei], dtype='float')
3236 coh_spectra = numpy.zeros([nChan, nProf, nHei], dtype='float')
3306 coh_cspectra = numpy.zeros([nPairs, nProf, nHei], dtype='complex')
3237 coh_cspectra = numpy.zeros([nPairs, nProf, nHei], dtype='complex')
3307 coh_aver = numpy.zeros([nChan, nHei])
3238 coh_aver = numpy.zeros([nChan, nHei])
3308
3239
3309 incoh_spectra = numpy.zeros([nChan, nProf, nHei], dtype='float')
3240 incoh_spectra = numpy.zeros([nChan, nProf, nHei], dtype='float')
3310 incoh_cspectra = numpy.zeros([nPairs, nProf, nHei], dtype='complex')
3241 incoh_cspectra = numpy.zeros([nPairs, nProf, nHei], dtype='complex')
3311 incoh_aver = numpy.zeros([nChan, nHei])
3242 incoh_aver = numpy.zeros([nChan, nHei])
3312 power = numpy.sum(spectra, axis=1)
3243 power = numpy.sum(spectra, axis=1)
3313
3244
3314 if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65
3245 if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65
3315 if hei_th == None : hei_th = numpy.array([60,300,650])
3246 if hei_th == None : hei_th = numpy.array([60,300,650])
3316 for ic in range(nPairs):
3247 for ic in range(nPairs):
@@ -3328,7 +3259,7 class SpectralFitting(Operation):
3328 if len(indv[0]) == 0 :
3259 if len(indv[0]) == 0 :
3329 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
3260 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
3330 if len(valid)>0:
3261 if len(valid)>0:
3331 my_coh_aver[pair[0],valid]=1
3262 my_coh_aver[pair[0],valid]=1
3332 my_coh_aver[pair[1],valid]=1
3263 my_coh_aver[pair[1],valid]=1
3333 # si la coherencia es mayor a la coherencia threshold los datos se toman
3264 # si la coherencia es mayor a la coherencia threshold los datos se toman
3334 coh = numpy.squeeze(numpy.nansum(cspectra[ic,:,:], axis=0)/numpy.sqrt(numpy.nansum(spectra[pair[0],:,:], axis=0)*numpy.nansum(spectra[pair[1],:,:], axis=0)))
3265 coh = numpy.squeeze(numpy.nansum(cspectra[ic,:,:], axis=0)/numpy.sqrt(numpy.nansum(spectra[pair[0],:,:], axis=0)*numpy.nansum(spectra[pair[1],:,:], axis=0)))
@@ -3341,7 +3272,7 class SpectralFitting(Operation):
3341 if len(valid)>0:
3272 if len(valid)>0:
3342 my_coh_aver[pair[0],hvalid[valid]] =1
3273 my_coh_aver[pair[0],hvalid[valid]] =1
3343 my_coh_aver[pair[1],hvalid[valid]] =1
3274 my_coh_aver[pair[1],hvalid[valid]] =1
3344
3275
3345 coh_echoes = (my_coh_aver[pair[0],:] == 1).nonzero()
3276 coh_echoes = (my_coh_aver[pair[0],:] == 1).nonzero()
3346 incoh_echoes = (my_coh_aver[pair[0],:] != 1).nonzero()
3277 incoh_echoes = (my_coh_aver[pair[0],:] != 1).nonzero()
3347 incoh_echoes = incoh_echoes[0]
3278 incoh_echoes = incoh_echoes[0]
@@ -3352,7 +3283,7 class SpectralFitting(Operation):
3352 my_incoh_aver[pair[0],incoh_echoes] = 1
3283 my_incoh_aver[pair[0],incoh_echoes] = 1
3353 my_incoh_aver[pair[1],incoh_echoes] = 1
3284 my_incoh_aver[pair[1],incoh_echoes] = 1
3354
3285
3355
3286
3356 for ic in range(nPairs):
3287 for ic in range(nPairs):
3357 pair = crosspairs[ic]
3288 pair = crosspairs[ic]
3358
3289
@@ -3390,8 +3321,8 class SpectralFitting(Operation):
3390 incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
3321 incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
3391 incoh_aver[pair[0],incoh_echoes]=1
3322 incoh_aver[pair[0],incoh_echoes]=1
3392 incoh_aver[pair[1],incoh_echoes]=1
3323 incoh_aver[pair[1],incoh_echoes]=1
3393 return my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver
3394
3324
3325 return my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver
3395
3326
3396 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
3327 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
3397
3328
@@ -3402,9 +3333,10 class SpectralFitting(Operation):
3402 nChan = len(channels)
3333 nChan = len(channels)
3403 crosspairs = dataOut.groupList
3334 crosspairs = dataOut.groupList
3404 nPairs = len(crosspairs)
3335 nPairs = len(crosspairs)
3405
3336
3406 absc = dataOut.abscissaList[:-1]
3337 absc = dataOut.abscissaList[:-1]
3407 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
3338 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
3339
3408 clean_coh_spectra = spectra.copy()
3340 clean_coh_spectra = spectra.copy()
3409 clean_coh_cspectra = cspectra.copy()
3341 clean_coh_cspectra = cspectra.copy()
3410 clean_coh_aver = coh_aver.copy()
3342 clean_coh_aver = coh_aver.copy()
@@ -3433,14 +3365,14 class SpectralFitting(Operation):
3433 # Checking spectral widths
3365 # Checking spectral widths
3434 if (spwd[pair[0],ih] > spwd_th[0]) or (spwd[pair[1],ih] > spwd_th[0]) :
3366 if (spwd[pair[0],ih] > spwd_th[0]) or (spwd[pair[1],ih] > spwd_th[0]) :
3435 # satelite
3367 # satelite
3436 clean_coh_spectra[pair,ih,:] = 0.0
3368 clean_coh_spectra[pair,:,ih] = 0.0
3437 clean_coh_cspectra[ic,ih,:] = 0.0
3369 clean_coh_cspectra[ic,:,ih] = 0.0
3438 clean_coh_aver[pair,ih] = 0
3370 clean_coh_aver[pair,ih] = 0
3439 else :
3371 else :
3440 if ((spwd[pair[0],ih] < spwd_th[1]) or (spwd[pair[1],ih] < spwd_th[1])) :
3372 if ((spwd[pair[0],ih] < spwd_th[1]) or (spwd[pair[1],ih] < spwd_th[1])) :
3441 # Especial event like sun.
3373 # Especial event like sun.
3442 clean_coh_spectra[pair,ih,:] = 0.0
3374 clean_coh_spectra[pair,:,ih] = 0.0
3443 clean_coh_cspectra[ic,ih,:] = 0.0
3375 clean_coh_cspectra[ic,:,ih] = 0.0
3444 clean_coh_aver[pair,ih] = 0
3376 clean_coh_aver[pair,ih] = 0
3445
3377
3446 return clean_coh_spectra, clean_coh_cspectra, clean_coh_aver
3378 return clean_coh_spectra, clean_coh_cspectra, clean_coh_aver
@@ -3449,9 +3381,9 class SpectralFitting(Operation):
3449
3381
3450 rfunc = cspectra.copy()
3382 rfunc = cspectra.copy()
3451 n_funct = len(rfunc[0,:,0,0])
3383 n_funct = len(rfunc[0,:,0,0])
3452 val_spc = spectra*0.0
3384 val_spc = spectra*0.0
3453 val_cspc = cspectra*0.0
3385 val_cspc = cspectra*0.0
3454 in_sat_spectra = spectra.copy()
3386 in_sat_spectra = spectra.copy()
3455 in_sat_cspectra = cspectra.copy()
3387 in_sat_cspectra = cspectra.copy()
3456
3388
3457 min_hei = 200
3389 min_hei = 200
@@ -3464,13 +3396,14 class SpectralFitting(Operation):
3464 nPairs = len(crosspairs)
3396 nPairs = len(crosspairs)
3465 hval=(heights >= min_hei).nonzero()
3397 hval=(heights >= min_hei).nonzero()
3466 ih=hval[0]
3398 ih=hval[0]
3399
3467 for ih in range(hval[0][0],nHei):
3400 for ih in range(hval[0][0],nHei):
3468 for ifreq in range(nProf):
3401 for ifreq in range(nProf):
3469 for ii in range(n_funct):
3402 for ii in range(n_funct):
3470
3403
3471 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
3404 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
3472 val = (numpy.isfinite(func2clean)==True).nonzero()
3405 val = (numpy.isfinite(func2clean)==True).nonzero()
3473 if len(val)>0:
3406 if len(val)>0:
3474 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
3407 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
3475 if min_val <= -40 : min_val = -40
3408 if min_val <= -40 : min_val = -40
3476 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
3409 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
@@ -3478,22 +3411,22 class SpectralFitting(Operation):
3478 step = 1
3411 step = 1
3479 #Getting bins and the histogram
3412 #Getting bins and the histogram
3480 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
3413 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
3481 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
3414 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
3482 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
3415 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
3483 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
3416 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
3484 parg = [numpy.amax(y_dist),mean,sigma]
3417 parg = [numpy.amax(y_dist),mean,sigma]
3485 try :
3418 try :
3486 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
3419 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
3487 mode = gauss_fit[1]
3420 mode = gauss_fit[1]
3488 stdv = gauss_fit[2]
3421 stdv = gauss_fit[2]
3489 except:
3422 except:
3490 mode = mean
3423 mode = mean
3491 stdv = sigma
3424 stdv = sigma
3492
3425
3493 #Removing echoes greater than mode + 3*stdv
3426 #Removing echoes greater than mode + 3*stdv
3494 factor_stdv = 2.5
3427 factor_stdv = 2.5
3495 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
3428 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
3496
3429
3497 if len(noval[0]) > 0:
3430 if len(noval[0]) > 0:
3498 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
3431 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
3499 cross_pairs = crosspairs[ii]
3432 cross_pairs = crosspairs[ii]
@@ -3512,11 +3445,12 class SpectralFitting(Operation):
3512 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
3445 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
3513 for ih in range(nHei):
3446 for ih in range(nHei):
3514 for ifreq in range(nProf):
3447 for ifreq in range(nProf):
3515 for ich in range(nChan):
3448 for ich in range(nChan):
3516 tmp = spectra[:,ich,ifreq,ih]
3449 tmp = spectra[:,ich,ifreq,ih]
3517 valid = (numpy.isfinite(tmp[:])==True).nonzero()
3450 valid = (numpy.isfinite(tmp[:])==True).nonzero()
3518 if len(valid[0]) >0 :
3451 if len(valid[0]) >0 :
3519 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3452 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3453
3520 for icr in range(nPairs):
3454 for icr in range(nPairs):
3521 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
3455 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
3522 valid = (numpy.isfinite(tmp)==True).nonzero()
3456 valid = (numpy.isfinite(tmp)==True).nonzero()
@@ -3525,10 +3459,10 class SpectralFitting(Operation):
3525 #Removing fake coherent echoes (at least 4 points around the point)
3459 #Removing fake coherent echoes (at least 4 points around the point)
3526 val_spectra = numpy.sum(val_spc,0)
3460 val_spectra = numpy.sum(val_spc,0)
3527 val_cspectra = numpy.sum(val_cspc,0)
3461 val_cspectra = numpy.sum(val_cspc,0)
3528
3462
3529 val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
3463 val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
3530 val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
3464 val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
3531
3465
3532 for i in range(nChan):
3466 for i in range(nChan):
3533 for j in range(nProf):
3467 for j in range(nProf):
3534 for k in range(nHei):
3468 for k in range(nHei):
@@ -3544,10 +3478,11 class SpectralFitting(Operation):
3544 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
3478 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
3545 tmp_sat_cspectra = cspectra.copy()
3479 tmp_sat_cspectra = cspectra.copy()
3546 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
3480 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
3481
3547 val = (val_spc > 0).nonzero()
3482 val = (val_spc > 0).nonzero()
3548 if len(val[0]) > 0:
3483 if len(val[0]) > 0:
3549 tmp_sat_spectra[val] = in_sat_spectra[val]
3484 tmp_sat_spectra[val] = in_sat_spectra[val]
3550
3485
3551 val = (val_cspc > 0).nonzero()
3486 val = (val_cspc > 0).nonzero()
3552 if len(val[0]) > 0:
3487 if len(val[0]) > 0:
3553 tmp_sat_cspectra[val] = in_sat_cspectra[val]
3488 tmp_sat_cspectra[val] = in_sat_cspectra[val]
@@ -3559,7 +3494,7 class SpectralFitting(Operation):
3559 for ifreq in range(nProf):
3494 for ifreq in range(nProf):
3560 for ich in range(nChan):
3495 for ich in range(nChan):
3561 tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
3496 tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
3562 valid = (numpy.isfinite(tmp)).nonzero()
3497 valid = (numpy.isfinite(tmp)).nonzero()
3563 if len(valid[0]) > 0:
3498 if len(valid[0]) > 0:
3564 sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3499 sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3565
3500
@@ -3568,40 +3503,45 class SpectralFitting(Operation):
3568 valid = (numpy.isfinite(tmp)).nonzero()
3503 valid = (numpy.isfinite(tmp)).nonzero()
3569 if len(valid[0]) > 0:
3504 if len(valid[0]) > 0:
3570 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3505 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3506
3571 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
3507 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
3572
3508
3573 def REM_ISOLATED_POINTS(self,array,rth):
3509 def REM_ISOLATED_POINTS(self,array,rth):
3574 if rth == None : rth = 4
3510 if rth == None : rth = 4
3575 num_prof = len(array[0,:,0])
3511 num_prof = len(array[0,:,0])
3576 num_hei = len(array[0,0,:])
3512 num_hei = len(array[0,0,:])
3577 n2d = len(array[:,0,0])
3513 n2d = len(array[:,0,0])
3578
3514
3579 for ii in range(n2d) :
3515 for ii in range(n2d) :
3580 tmp = array[ii,:,:]
3516 tmp = array[ii,:,:]
3581 tmp = numpy.reshape(tmp,num_prof*num_hei)
3517 tmp = numpy.reshape(tmp,num_prof*num_hei)
3582 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
3518 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
3583 indxs2 = (tmp > 0).nonzero()
3519 indxs2 = (tmp > 0).nonzero()
3584 indxs1 = (indxs1[0])
3520 indxs1 = (indxs1[0])
3585 indxs2 = indxs2[0]
3521 indxs2 = indxs2[0]
3586 indxs = None
3522 indxs = None
3523
3587 for iv in range(len(indxs2)):
3524 for iv in range(len(indxs2)):
3588 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
3525 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
3589 if len(indv[0]) > 0 :
3526 if len(indv[0]) > 0 :
3590 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
3527 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
3528
3591 indxs = indxs[1:]
3529 indxs = indxs[1:]
3592 if len(indxs) < 4 :
3530 if len(indxs) < 4 :
3593 array[ii,:,:] = 0.
3531 array[ii,:,:] = 0.
3594 return
3532 return
3595
3533
3596 xpos = numpy.mod(indxs ,num_hei)
3534 xpos = numpy.mod(indxs ,num_prof)
3597 ypos = (indxs / num_hei)
3535 ypos = (indxs / num_prof)
3598 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
3536 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
3599 xpos = xpos[sx]
3537 xpos = xpos[sx]
3600 ypos = ypos[sx]
3538 ypos = ypos[sx]
3601 # *********************************** Cleaning isolated points **********************************
3539
3540 # *********************************** Cleaning isolated points **********************************
3602 ic = 0
3541 ic = 0
3603 while True :
3542 while True :
3604 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
3543 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
3544
3605 no_coh1 = (numpy.isfinite(r)==True).nonzero()
3545 no_coh1 = (numpy.isfinite(r)==True).nonzero()
3606 no_coh2 = (r <= rth).nonzero()
3546 no_coh2 = (r <= rth).nonzero()
3607 no_coh1 = numpy.array(no_coh1[0])
3547 no_coh1 = numpy.array(no_coh1[0])
@@ -3615,32 +3555,33 class SpectralFitting(Operation):
3615 if len(no_coh) < 4 :
3555 if len(no_coh) < 4 :
3616 xpos[ic] = numpy.nan
3556 xpos[ic] = numpy.nan
3617 ypos[ic] = numpy.nan
3557 ypos[ic] = numpy.nan
3618
3558
3619 ic = ic + 1
3559 ic = ic + 1
3620 if (ic == len(indxs)) :
3560 if (ic == len(indxs)) :
3621 break
3561 break
3622 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
3562 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
3623 if len(indxs[0]) < 4 :
3563 if len(indxs[0]) < 4 :
3624 array[ii,:,:] = 0.
3564 array[ii,:,:] = 0.
3625 return
3565 return
3626
3566
3627 xpos = xpos[indxs[0]]
3567 xpos = xpos[indxs[0]]
3628 ypos = ypos[indxs[0]]
3568 ypos = ypos[indxs[0]]
3629 for i in range(0,len(ypos)):
3569 for i in range(0,len(ypos)):
3630 ypos[i]=int(ypos[i])
3570 ypos[i]=int(ypos[i])
3631 junk = tmp
3571 junk = tmp
3632 tmp = junk*0.0
3572 tmp = junk*0.0
3633
3573
3634 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
3574 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
3635 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
3575 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
3576
3636 return array
3577 return array
3637
3578
3638 def moments(self,doppler,yarray,npoints):
3579 def moments(self,doppler,yarray,npoints):
3639 ytemp = yarray
3580 ytemp = yarray
3640 val = (ytemp > 0).nonzero()
3581 val = (ytemp > 0).nonzero()
3641 val = val[0]
3582 val = val[0]
3642 if len(val) == 0 : val = range(npoints-1)
3583 if len(val) == 0 : val = range(npoints-1)
3643
3584
3644 ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]])
3585 ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]])
3645 ytemp[len(ytemp):] = [ynew]
3586 ytemp[len(ytemp):] = [ynew]
3646
3587
@@ -3653,317 +3594,1492 class SpectralFitting(Operation):
3653 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3594 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3654 return [fmom,numpy.sqrt(smom)]
3595 return [fmom,numpy.sqrt(smom)]
3655
3596
3656 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None, filec=None,coh_th=None, hei_th=None,taver=None,proc=None,nhei=None,nprofs=None,ipp=None,channelList=None):
3597 def windowing_single_old(self,spc,x,A,B,C,D,nFFTPoints):
3657 if not self.isConfig:
3598 '''
3658 self.setup(dataOut = dataOut,groupList=groupList,path=path,file=file,filec=filec)
3599 Written by R. Flores
3659 self.isConfig = True
3600 '''
3660
3601 from scipy.optimize import curve_fit,fmin
3661 if not numpy.any(proc):
3662 if numpy.any(taver):
3663 taver = int(taver)
3664 else :
3665 taver = 5
3666 tini = time.localtime(dataOut.utctime)
3667 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3668 self.index = 0
3669 jspc = self.buffer
3670 jcspc = self.buffer2
3671 jnoise = self.buffer3
3672 self.buffer = dataOut.data_spc
3673 self.buffer2 = dataOut.data_cspc
3674 self.buffer3 = dataOut.noise
3675 self.fint = 1
3676 if numpy.any(jspc) :
3677 jspc = numpy.reshape(jspc ,(int(len(jspc) / self.nChannels) , self.nChannels ,self.nProf,self.nHeights ))
3678 jcspc = numpy.reshape(jcspc ,(int(len(jcspc) /int(self.nChannels/2)),int(self.nChannels/2),self.nProf,self.nHeights ))
3679 jnoise = numpy.reshape(jnoise,(int(len(jnoise)/ self.nChannels) , self.nChannels))
3680 else:
3681 dataOut.flagNoData = True
3682 return dataOut
3683 else :
3684 if (tini.tm_min % taver) == 0 :
3685 self.fint = 1
3686 else :
3687 self.fint = 0
3688
3689 self.index += 1
3690 if numpy.any(self.buffer):
3691 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
3692 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
3693 self.buffer3 = numpy.concatenate((self.buffer3,dataOut.noise), axis=0)
3694 else:
3695 self.buffer = dataOut.data_spc
3696 self.buffer2 = dataOut.data_cspc
3697 self.buffer3 = dataOut.noise
3698 dataOut.flagNoData = True
3699 return dataOut
3700
3602
3701 jnoise = jnoise/self.N# creo que falta dividirlo entre N
3603 def gaussian(x, a, b, c, d):
3702 noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
3604 val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d
3703 index = tini.tm_hour*12+tini.tm_min/taver
3605 return val
3704 dataOut.index = index
3705 jspc = jspc/self.N/self.N
3706 jcspc = jcspc/self.N/self.N
3707
3606
3708 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3607 def R_gaussian(x, a, b, c):
3709 jspectra = tmp_spectra * len(jspc[:,0,0,0])
3608 N = int(numpy.shape(x)[0])
3710 jcspectra = tmp_cspectra * len(jspc[:,0,0,0])
3609 val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi)
3610 return val
3711
3611
3712 my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, self.snrth,coh_th, hei_th)
3612 def T(x,N):
3713 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(self.snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
3613 T = 1-abs(x)/N
3614 return T
3714
3615
3715 dataOut.data_spc = incoh_spectra
3616 def R_T_spc_fun(x, a, b, c, d, nFFTPoints):
3716 dataOut.data_cspc = incoh_cspectra
3717 clean_num_aver = incoh_aver * len(jspc[:,0,0,0])
3718 coh_num_aver = clean_coh_aver* len(jspc[:,0,0,0])
3719 dataOut.clean_num_aver = clean_num_aver
3720 dataOut.coh_num_aver = coh_num_aver
3721
3617
3722 else:
3618 N = int(numpy.shape(x)[0])
3723 clean_num_aver = dataOut.clean_num_aver
3724 coh_num_aver = dataOut.coh_num_aver
3725 dataOut.data_spc = dataOut.tmp_spectra_i
3726 dataOut.data_cspc = dataOut.tmp_cspectra_i
3727 clean_coh_spectra = dataOut.tmp_spectra_c
3728 clean_coh_cspectra = dataOut.tmp_cspectra_c
3729 jspectra = dataOut.data_spc+clean_coh_spectra
3730 nHeights = len(dataOut.heightList) # nhei
3731 nProf = int(dataOut.nProfiles)
3732 dataOut.nProfiles = nProf
3733 dataOut.data_param = None
3734 dataOut.data_paramC = None
3735 dataOut.code = numpy.array([[-1.,-1.,1.],[1.,1.,-1.]])
3736 #M=600
3737 #N=200
3738 dataOut.flagDecodeData=True
3739 M = int(dataOut.normFactor)
3740 N = int(dataOut.nFFTPoints)
3741 dataOut.nFFTPoints = N
3742 dataOut.nIncohInt= int(dataOut.nIncohInt)
3743 dataOut.nProfiles = int(dataOut.nProfiles)
3744 dataOut.nCohInt = int(dataOut.nCohInt)
3745 #dataOut.nFFTPoints=nprofs
3746 #dataOut.normFactor = nprofs
3747 dataOut.channelList = channelList
3748 #dataOut.ippFactor=1
3749 #ipp = ipp/150*1.e-3
3750 vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
3751 #dataOut.ippSeconds=ipp
3752 absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf
3753 if path != None:
3754 sys.path.append(path)
3755 self.library = importlib.import_module(file)
3756 constants = self.library.setConstants(dataOut)
3757 constants['M'] = M
3758 dataOut.constants = constants
3759
3760 #List of possible combinations
3761 listComb = itertools.combinations(numpy.arange(self.groupArray.shape[1]),2)
3762 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3763 if dataOut.data_paramC is None:
3764 dataOut.data_paramC = numpy.zeros((self.nGroups*4, self.nHeights ,2))*numpy.nan
3765 for i in range(self.nGroups):
3766 coord = self.groupArray[i,:]
3767 #Input data array
3768 data = dataOut.data_spc[coord,:,:]/(self.M*self.N)
3769 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3770
3619
3771 #Cross Spectra data array for Covariance Matrixes
3620 x_max = x[-1]
3772 ind = 0
3773 for pairs in listComb:
3774 pairsSel = numpy.array([coord[x],coord[y]])
3775 indCross[ind] = int(numpy.where(numpy.all(self.pairsArray == pairsSel, axis = 1))[0][0])
3776 ind += 1
3777 dataCross = dataOut.data_cspc[indCross,:,:]/(self.M*self.N)
3778 dataCross = dataCross**2
3779 nhei = self.nHeights
3780 poweri = numpy.sum(dataOut.data_spc[:,1:self.nProf-0,:],axis=1)/clean_num_aver[:,:]
3781
3621
3782 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
3622 x_pos = x[nFFTPoints:]
3783 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(self.nProf-1)
3623 x_neg = x[:nFFTPoints]
3784 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(self.nProf-1)
3785 n0 = n0i
3786 n1= n1i
3787 my_noises[2*i+0] = n0
3788 my_noises[2*i+1] = n1
3789 snrth = -15.0 # -4 -16 -25
3790 snrth = 10**(snrth/10.0)
3791 jvelr = numpy.zeros(self.nHeights, dtype = 'float')
3792 hvalid = [0]
3793 coh2 = abs(dataOut.data_cspc[i,1:self.nProf,:])**2/(dataOut.data_spc[0+i*2,1:self.nProf-0,:]*dataOut.data_spc[1+i*2,1:self.nProf-0,:])
3794
3624
3795 for h in range(self.nHeights):
3625 R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0])
3796 smooth = clean_num_aver[i+1,h]
3626 R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1])
3797 signalpn0 = (dataOut.data_spc[i*2,1:(self.nProf-0),h])/smooth
3627 #print(T(x_pos,x[-1]),x_pos,x[-1])
3798 signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth
3628 #print(R_T_neg_1.shape,R_T_pos_1.shape)
3799 signal0 = signalpn0-n0
3629 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
3800 signal1 = signalpn1-n1
3630 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
3801 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3631 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
3802 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3632 max_val_1 = numpy.max(R_T_spc_1)
3803 gamma = coh2[:,h]
3633 R_T_spc_1 = R_T_spc_1*a/max_val_1
3804 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
3634 print("R_T_spc_1: ", R_T_spc_1)
3805 if len(indxs) >0:
3806 if numpy.nanmean(gamma) > 0.07:
3807 maxp0 = numpy.argmax(signal0*gamma)
3808 maxp1 = numpy.argmax(signal1*gamma)
3809 #print('usa gamma',numpy.nanmean(gamma))
3810 else:
3811 maxp0 = numpy.argmax(signal0)
3812 maxp1 = numpy.argmax(signal1)
3813 jvelr[h] = (self.absc[maxp0]+self.absc[maxp1])/2.
3814 else: jvelr[h] = self.absc[0]
3815 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
3816 #print(maxp0,absc[maxp0],snr0,jvelr[h])
3817
3818 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
3819 else: fd0 = numpy.nan
3820 for h in range(self.nHeights):
3821 d = data[:,h]
3822 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
3823 signalpn0 = (dataOut.data_spc[i*2,1:(self.nProf-0),h])/smooth
3824 signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth
3825 signal0 = signalpn0-n0
3826 signal1 = signalpn1-n1
3827 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3828 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3829 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
3830 #Covariance Matrix
3831 D = numpy.diag(d**2)
3832 ind = 0
3833 for pairs in listComb:
3834 #Coordinates in Covariance Matrix
3835 x = pairs[0]
3836 y = pairs[1]
3837 #Channel Index
3838 S12 = dataCross[ind,:,h]
3839 D12 = numpy.diag(S12)
3840 #Completing Covariance Matrix with Cross Spectras
3841 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
3842 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
3843 ind += 1
3844 diagD = numpy.zeros(256)
3845
3635
3846 try:
3636 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
3847 Dinv=numpy.linalg.inv(D)
3637 R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0])
3848 L=numpy.linalg.cholesky(Dinv)
3638 R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1])
3849 except:
3639 R_T_d_sum = R_T_d_pos + R_T_d_neg
3850 Dinv = D*numpy.nan
3640 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
3851 L= D*numpy.nan
3641 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
3852 LT=L.T
3853
3642
3854 dp = numpy.dot(LT,d)
3643 R_T_final = R_T_spc_1# + R_T_spc_3
3855 #Initial values
3856 data_spc = dataOut.data_spc[coord,:,h]
3857 w = data_spc/data_spc
3858 if filec != None:
3859 w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
3860 if (h>6)and(error1[3]<25):
3861 p0 = dataOut.data_param[i,:,h-1].copy()
3862 else:
3863 p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, self.constants))# sin el i(data_spc, constants, i)
3864 p0[3] = fd0
3865
3644
3866 if filec != None:
3645 return R_T_final
3867 p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)
3868
3869 try:
3870 #Least Squares
3871 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,self.constants),full_output=True)
3872 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
3873 #Chi square error
3874 error0 = numpy.sum(infodict['fvec']**2)/(2*self.N)
3875 #Error with Jacobian
3876 error1 = self.library.errorFunction(minp,self.constants,LT)
3877
3646
3878 except:
3647 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
3879 minp = p0*numpy.nan
3880 error0 = numpy.nan
3881 error1 = p0*numpy.nan
3882 else :
3883 data_spc = dataOut.data_spc[coord,:,h]
3884 p0 = numpy.array(self.library.initialValuesFunction(data_spc, self.constants))
3885 minp = p0*numpy.nan
3886 error0 = numpy.nan
3887 error1 = p0*numpy.nan
3888 if dataOut.data_param is None:
3889 dataOut.data_param = numpy.zeros((self.nGroups, p0.size, self.nHeights ))*numpy.nan
3890 dataOut.data_error = numpy.zeros((self.nGroups, p0.size + 1, self.nHeights ))*numpy.nan
3891 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
3892 dataOut.data_param[i,:,h] = minp
3893
3648
3894 for ht in range(self.nHeights-1) :
3649 from scipy.stats import norm
3895 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
3650 mean,std=norm.fit(spc)
3896 dataOut.data_paramC[4*i,ht,1] = smooth
3897 signalpn0 = (clean_coh_spectra[i*2 ,1:(self.nProf-0),ht])/smooth #coh_spectra
3898 signalpn1 = (clean_coh_spectra[i*2+1,1:(self.nProf-0),ht])/smooth
3899 val0 = (signalpn0 > 0).nonzero()
3900 val0 = val0[0]
3901 if len(val0) == 0 : val0_npoints = self.nProf
3902 else : val0_npoints = len(val0)
3903
3904 val1 = (signalpn1 > 0).nonzero()
3905 val1 = val1[0]
3906 if len(val1) == 0 : val1_npoints = self.nProf
3907 else : val1_npoints = len(val1)
3908
3651
3909 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
3652 # estimate starting values from the data
3910 dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1
3653 print("A: ", A)
3911
3654 a = A-D
3912 signal0 = (signalpn0-n0)
3655 b = B
3913 vali = (signal0 < 0).nonzero()
3656 c = C#numpy.std(spc) #C
3914 vali = vali[0]
3657 d = D
3915 if len(vali) > 0 : signal0[vali] = 0
3658 #'''
3916 signal1 = (signalpn1-n1)
3659 #ippSeconds = 250*20*1.e-6/3
3917 vali = (signal1 < 0).nonzero()
3918 vali = vali[0]
3919 if len(vali) > 0 : signal1[vali] = 0
3920 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3921 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3922 doppler = self.absc[1:]
3923 if snr0 >= snrth and snr1 >= snrth and smooth :
3924 signalpn0_n0 = signalpn0
3925 signalpn0_n0[val0] = signalpn0[val0] - n0
3926 mom0 = self.moments(doppler,signalpn0-n0,self.nProf)
3927 signalpn1_n1 = signalpn1
3928 signalpn1_n1[val1] = signalpn1[val1] - n1
3929 mom1 = self.moments(doppler,signalpn1_n1,self.nProf)
3930 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
3931 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
3932 dataOut.data_spc = jspectra
3933 dataOut.spc_noise = my_noises*self.nProf*self.M
3934 if numpy.any(proc): dataOut.spc_noise = my_noises*self.nProf*self.M
3935 if getSNR:
3936 listChannels = self.groupArray.reshape((self.groupArray.size))
3937 listChannels.sort()
3938 # TEST
3939 noise_C = numpy.zeros(self.nChannels)
3940 noise_C = dataOut.getNoise()
3941 #print("noise_C",noise_C)
3942 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:],noise_C/(600.0*1.15))# PRUEBA *nProf*M
3943 #dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise_C[listChannels])# PRUEBA *nProf*M
3944 dataOut.flagNoData = False
3945 return dataOut
3946
3947 def __residFunction(self, p, dp, LT, constants):
3948
3660
3949 fm = self.library.modelFunction(p, constants)
3661 #x_t = ippSeconds * (numpy.arange(nFFTPoints) - nFFTPoints / 2.)
3950 fmp=numpy.dot(LT,fm)
3662
3951 return dp-fmp
3663 #x_t = numpy.linspace(x_t[0],x_t[-1],3200)
3664 #print("x_t: ", x_t)
3665 #print("nFFTPoints: ", nFFTPoints)
3666 x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints))
3667 #print("x_vel: ", x_vel)
3668 #x_freq = numpy.fft.fftfreq(1600,d=ippSeconds)
3669 #x_freq = numpy.fft.fftshift(x_freq)
3670 #'''
3671 # define a least squares function to optimize
3672 import matplotlib.pyplot as plt
3673 aui = R_T_spc_fun(x_vel,a,b,c,d,nFFTPoints)
3674 print("aux_max: ", numpy.nanmax(aui))
3675 #print(dataOut.heightList[hei])
3676 plt.figure()
3677 plt.plot(x,spc,marker='*',linestyle='--')
3678 plt.plot(x,gaussian(x, a, b, c, d),color='b',marker='^',linestyle='')
3679 plt.plot(x,aui,color='k')
3680 #plt.title(dataOut.heightList[hei])
3681 plt.show()
3682
3683 def minfunc(params):
3684 #print("y.shape: ", numpy.shape(y))
3685 return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2)
3686
3687 # fit
3688
3689 popt_full = fmin(minfunc,[a,b,c,d], disp=False)
3690 #print("nIter", popt_full[2])
3691 popt = popt_full#[0]
3692
3693 fun = gaussian(x, popt[0], popt[1], popt[2], popt[3])
3694 print("pop1[0]: ", popt[0])
3695 #return R_T_spc_fun(x_t,popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]), popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
3696 return fun, popt[0], popt[1], popt[2], popt[3]
3697
3698 def windowing_single(self,spc,x,A,B,C,D,nFFTPoints):
3699 '''
3700 Written by R. Flores
3701 '''
3702 from scipy.optimize import curve_fit,fmin
3703
3704 def gaussian(x, a, b, c, d):
3705 val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d
3706 return val
3707
3708 def R_gaussian(x, a, b, c):
3709 N = int(numpy.shape(x)[0])
3710
3711 val = (a*numpy.exp((-(1/2)*x*(x*c**2 + 2*1.j*b)))/numpy.sqrt(1/c**2))
3712
3713 return val
3714
3715 def T(x,N):
3716 T = 1-abs(x)/N
3717 return T
3718
3719 def R_T_spc_fun(x, a, id_dop, c, d, nFFTPoints):
3720
3721 N = int(numpy.shape(x)[0])
3722 b = 0
3723 x_max = x[-1]
3724
3725 x_pos = x[nFFTPoints:]
3726 x_neg = x[:nFFTPoints]
3727 R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0])
3728 R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1])
3729
3730 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
3731 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
3732 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
3733 max_val_1 = numpy.max(R_T_spc_1)
3734 R_T_spc_1 = R_T_spc_1*a/max_val_1
3735 #raise NotImplementedError
3736 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
3737 R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0])
3738 R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1])
3739 R_T_d_sum = R_T_d_pos + R_T_d_neg
3740 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
3741 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
3742
3743 R_T_final = R_T_spc_1 + R_T_spc_3
3744
3745 id_dop = int(id_dop)
3746
3747 R_T_final = numpy.roll(R_T_final,-id_dop)
3748
3749 return R_T_final
3750
3751 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
3752
3753 from scipy.stats import norm
3754 mean,std=norm.fit(spc)
3755
3756 # estimate starting values from the data
3757 a = A-D
3758 b = B
3759 c = C#numpy.std(spc) #C
3760 d = D
3761
3762 id_dop = numpy.argmax(spc)
3763 id_dop = int(spc.shape[0]/2 - id_dop)
3764
3765 x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints))
3766
3767 # define a least squares function to optimize
3768
3769 def minfunc(params):
3770 #print("y.shape: ", numpy.shape(y))
3771 return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2)
3772
3773 # fit
3774 popt_full = fmin(minfunc,[a,id_dop,c,d], disp=False)
3775 popt = popt_full#[0]
3776
3777 fun = gaussian(x, a, 0, popt[2], popt[3])
3778 fun = numpy.roll(fun,-int(popt[1]))
3779
3780 return fun, popt[0], popt[1], popt[2], popt[3]
3781
3782 def windowing_single_direct(self,spc_mod,x,A,B,C,D,nFFTPoints,timeInterval):
3783 '''
3784 Written by R. Flores
3785 '''
3786 from scipy.optimize import curve_fit,fmin
3787
3788 def gaussian(x, a, b, c, d):
3789 val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d
3790 return val
3791
3792 def R_gaussian(x, a, b, c, d):
3793 N = int(numpy.shape(x)[0])
3794 val = (a*numpy.exp(-2*c**2*x**2 + 2*x*1.j*b))*(numpy.sqrt(2*numpy.pi)*c)/((numpy.pi)) + d*signal.unit_impulse(N)*numpy.shape(x)[0]/2
3795
3796 return 2*val/numpy.shape(val)[0]
3797
3798 def T(x,N):
3799 T = 1-abs(x)/N
3800 return T
3801
3802 def R_T_spc_fun(x, a, b, c, d, nFFTPoints, timeInterval): #"x" should be time
3803
3804 #timeInterval = 2
3805 x_double = numpy.linspace(0,timeInterval,nFFTPoints)
3806 x_double_m = numpy.flip(x_double)
3807 x_double_aux = numpy.linspace(0,x_double[-2],nFFTPoints)
3808 x_double_t = numpy.concatenate((x_double_m,x_double_aux))
3809 x_double_t /= max(x_double_t)
3810
3811
3812 R_T_sum_1 = R_gaussian(x, a, b, c, d)
3813
3814 R_T_sum_1_flip = numpy.copy(numpy.flip(R_T_sum_1))
3815 R_T_sum_1_flip[-1] = R_T_sum_1_flip[0]
3816 R_T_sum_1_flip = numpy.roll(R_T_sum_1_flip,1)
3817
3818 R_T_sum_1_flip.imag *= -1
3819
3820 R_T_sum_1_total = numpy.concatenate((R_T_sum_1,R_T_sum_1_flip))
3821 R_T_sum_1_total *= x_double_t #times trian_fun
3822
3823 R_T_sum_1_total = R_T_sum_1_total[:nFFTPoints] + R_T_sum_1_total[nFFTPoints:]
3824
3825 R_T_spc_1 = numpy.fft.fft(R_T_sum_1_total).real
3826 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
3827
3828 freq = numpy.fft.fftfreq(nFFTPoints, d=timeInterval/nFFTPoints)
3829
3830 freq = numpy.fft.fftshift(freq)
3831
3832 freq *= 6/2 #lambda/2
3833
3834 return R_T_spc_1
3835
3836 y = spc_mod
3837
3838 #from scipy.stats import norm
3839
3840 # estimate starting values from the data
3841
3842 a = A-D
3843 b = B
3844 c = C
3845 d = D
3846
3847 # define a least squares function to optimize
3848 import matplotlib.pyplot as plt
3849 #ippSeconds = 2
3850 t_range = numpy.linspace(0,timeInterval,nFFTPoints)
3851 #aui = R_T_spc_fun(t_range,a,b,c,d,nFFTPoints,timeInterval)
3852
3853 def minfunc(params):
3854 return sum((y-R_T_spc_fun(t_range,params[0],params[1],params[2],params[3],nFFTPoints,timeInterval))**2/1)#y**2)
3855
3856 # fit
3857 popt_full = fmin(minfunc,[a,b,c,d], disp=False)
3858 popt = popt_full
3859
3860 fun = R_T_spc_fun(t_range,popt[0],popt[1],popt[2],popt[3],nFFTPoints,timeInterval)
3861
3862 return fun, popt[0], popt[1], popt[2], popt[3]
3863 # **********************************************************************************************
3864 index = 0
3865 fint = 0
3866 buffer = 0
3867 buffer2 = 0
3868 buffer3 = 0
3869
3870 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None, filec=None,coh_th=None, hei_th=None,taver=None,proc=None,nhei=None,nprofs=None,ipp=None,channelList=None,Gaussian_Windowed=0):
3871
3872 if not numpy.any(proc):
3873 nChannels = dataOut.nChannels
3874 nHeights= dataOut.heightList.size
3875 nProf = dataOut.nProfiles
3876 if numpy.any(taver): taver=int(taver)
3877 else : taver = 5
3878 tini=time.localtime(dataOut.utctime)
3879 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3880 self.index = 0
3881 jspc = self.buffer
3882 jcspc = self.buffer2
3883 jnoise = self.buffer3
3884 self.buffer = dataOut.data_spc
3885 self.buffer2 = dataOut.data_cspc
3886 self.buffer3 = dataOut.noise
3887 self.fint = 1
3888 if numpy.any(jspc) :
3889 jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights))
3890 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights))
3891 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels))
3892 else:
3893 dataOut.flagNoData = True
3894 return dataOut
3895 else:
3896 if (tini.tm_min % taver) == 0 : self.fint = 1
3897 else : self.fint = 0
3898 self.index += 1
3899 if numpy.any(self.buffer):
3900 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
3901 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
3902 self.buffer3 = numpy.concatenate((self.buffer3,dataOut.noise), axis=0)
3903 else:
3904 self.buffer = dataOut.data_spc
3905 self.buffer2 = dataOut.data_cspc
3906 self.buffer3 = dataOut.noise
3907 dataOut.flagNoData = True
3908 return dataOut
3909 if path != None:
3910 sys.path.append(path)
3911 self.library = importlib.import_module(file)
3912 if filec != None:
3913 self.weightf = importlib.import_module(filec)
3914 #self.weightf = importlib.import_module('weightfit')
3915
3916 #To be inserted as a parameter
3917 groupArray = numpy.array(groupList)
3918 #groupArray = numpy.array([[0,1],[2,3]])
3919 dataOut.groupList = groupArray
3920
3921 nGroups = groupArray.shape[0]
3922 nChannels = dataOut.nChannels
3923 nHeights = dataOut.heightList.size
3924
3925 #Parameters Array
3926 dataOut.data_param = None
3927 dataOut.data_paramC = None
3928 dataOut.clean_num_aver = None
3929 dataOut.coh_num_aver = None
3930 dataOut.tmp_spectra_i = None
3931 dataOut.tmp_cspectra_i = None
3932 dataOut.tmp_spectra_c = None
3933 dataOut.tmp_cspectra_c = None
3934 dataOut.sat_spectra = None
3935 dataOut.sat_cspectra = None
3936 dataOut.index = None
3937
3938 #Set constants
3939 constants = self.library.setConstants(dataOut)
3940 dataOut.constants = constants
3941 M = dataOut.normFactor
3942 N = dataOut.nFFTPoints
3943
3944 ippSeconds = dataOut.ippSeconds
3945 K = dataOut.nIncohInt
3946 pairsArray = numpy.array(dataOut.pairsList)
3947
3948 snrth= 15
3949 spectra = dataOut.data_spc
3950 cspectra = dataOut.data_cspc
3951 nProf = dataOut.nProfiles
3952 heights = dataOut.heightList
3953 nHei = len(heights)
3954 channels = dataOut.channelList
3955 nChan = len(channels)
3956 nIncohInt = dataOut.nIncohInt
3957 crosspairs = dataOut.groupList
3958 noise = dataOut.noise
3959 jnoise = jnoise/N
3960 noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
3961 power = numpy.sum(spectra, axis=1)
3962 nPairs = len(crosspairs)
3963 absc = dataOut.abscissaList[:-1]
3964 print('para escribir h5 ',dataOut.paramInterval)
3965 if not self.isConfig:
3966 self.isConfig = True
3967
3968 index = tini.tm_hour*12+tini.tm_min/taver
3969 dataOut.index= index
3970 jspc = jspc/N/N
3971 jcspc = jcspc/N/N
3972 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3973 jspectra = tmp_spectra*len(jspc[:,0,0,0])
3974 jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
3975 my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth,coh_th, hei_th)
3976
3977 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
3978 dataOut.data_spc = incoh_spectra
3979 dataOut.data_cspc = incoh_cspectra
3980 dataOut.sat_spectra = sat_spectra
3981 dataOut.sat_cspectra = sat_cspectra
3982 # dataOut.data_spc = tmp_spectra
3983 # dataOut.data_cspc = tmp_cspectra
3984
3985 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
3986 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
3987 # clean_num_aver = (numpy.zeros([nChan, nHei])+1)*len(jspc[:,0,0,0])
3988 # coh_num_aver = numpy.zeros([nChan, nHei])*0*len(jspc[:,0,0,0])
3989 dataOut.clean_num_aver = clean_num_aver
3990 dataOut.coh_num_aver = coh_num_aver
3991 dataOut.tmp_spectra_i = incoh_spectra
3992 dataOut.tmp_cspectra_i = incoh_cspectra
3993 dataOut.tmp_spectra_c = clean_coh_spectra
3994 dataOut.tmp_cspectra_c = clean_coh_cspectra
3995 #List of possible combinations
3996 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
3997 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3998 if Gaussian_Windowed == 1:
3999 #dataOut.data_spc = jspectra
4000 '''
4001 Written by R. Flores
4002 '''
4003 print("normFactor: ", dataOut.normFactor)
4004 data_spc_aux = numpy.copy(dataOut.data_spc)#[:,0,:]
4005 data_spc_aux[:,0,:] = (data_spc_aux[:,1,:]+data_spc_aux[:,-1,:])/2
4006 #'''
4007 from scipy.signal import medfilt
4008 import matplotlib.pyplot as plt
4009 dataOut.moments = numpy.ones((dataOut.nChannels,4,dataOut.nHeights))*numpy.NAN
4010 dataOut.VelRange = dataOut.getVelRange(0)
4011 for nChannel in range(dataOut.nChannels):
4012 for hei in range(dataOut.heightList.shape[0]):
4013 #print("ipp: ", dataOut.ippSeconds)
4014 #spc = numpy.copy(dataOut.data_spc[nChannel,:,hei])
4015 spc = data_spc_aux[nChannel,:,hei]
4016 if spc.all() == 0.:
4017 print("CONTINUE")
4018 continue
4019 #print(VelRange)
4020 #print(dataOut.getFreqRange(64))
4021 #print("Hei: ", dataOut.heightList[hei])
4022
4023 spc_mod = numpy.copy(spc)
4024 spcm = medfilt(spc_mod,11)
4025 spc_max = numpy.max(spcm)
4026 dop1_x0 = dataOut.VelRange[numpy.argmax(spcm)]
4027 #D = numpy.min(spcm)
4028 D_in = (numpy.mean(spcm[:15])+numpy.mean(spcm[-15:]))/2.
4029 #print("spc_max: ", spc_max)
4030 #print("dataOut.ippSeconds: ", dataOut.ippSeconds, dataOut.timeInterval)
4031 ##fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints)
4032 #fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints)
4033 fun, A, B, C, D = self.windowing_single_direct(spc_mod,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0/5),D_in,dataOut.nFFTPoints,dataOut.timeInterval)
4034
4035 dataOut.moments[nChannel,0,hei] = A
4036 dataOut.moments[nChannel,1,hei] = B
4037 dataOut.moments[nChannel,2,hei] = C
4038 dataOut.moments[nChannel,3,hei] = D
4039 '''
4040 if nChannel == 0:
4041 print(dataOut.heightList[hei])
4042 plt.figure()
4043 plt.plot(dataOut.VelRange,spc,marker='*',linestyle='--')
4044 plt.plot(dataOut.VelRange,fun)
4045 plt.title(dataOut.heightList[hei])
4046 plt.show()
4047 '''
4048 #plt.show()
4049 #'''
4050 dataOut.data_spc = jspectra
4051 print("SUCCESS")
4052 return dataOut
4053
4054 elif Gaussian_Windowed == 2: #Only to clean spc
4055 dataOut.VelRange = dataOut.getVelRange(0)
4056 return dataOut
4057
4058 if getSNR:
4059 listChannels = groupArray.reshape((groupArray.size))
4060 listChannels.sort()
4061 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
4062 else:
4063 if numpy.any(taver): taver=int(taver)
4064 else : taver = 5
4065 tini=time.localtime(dataOut.utctime)
4066 index = tini.tm_hour*12+tini.tm_min/taver
4067 clean_num_aver = dataOut.clean_num_aver
4068 coh_num_aver = dataOut.coh_num_aver
4069 dataOut.data_spc = dataOut.tmp_spectra_i
4070 dataOut.data_cspc = dataOut.tmp_cspectra_i
4071 clean_coh_spectra = dataOut.tmp_spectra_c
4072 clean_coh_cspectra = dataOut.tmp_cspectra_c
4073 jspectra = dataOut.data_spc+clean_coh_spectra
4074 nHeights = len(dataOut.heightList) # nhei
4075 nProf = int(dataOut.nProfiles)
4076 dataOut.nProfiles = nProf
4077 dataOut.data_param = None
4078 dataOut.data_paramC = None
4079 dataOut.code = numpy.array([[-1.,-1.,1.],[1.,1.,-1.]])
4080 #dataOut.paramInterval = 2.0
4081 #M=600
4082 #N=200
4083 dataOut.flagDecodeData=True
4084 M = int(dataOut.normFactor)
4085 N = int(dataOut.nFFTPoints)
4086 dataOut.nFFTPoints = N
4087 dataOut.nIncohInt= int(dataOut.nIncohInt)
4088 dataOut.nProfiles = int(dataOut.nProfiles)
4089 dataOut.nCohInt = int(dataOut.nCohInt)
4090 print('sale',dataOut.nProfiles,dataOut.nHeights)
4091 #dataOut.nFFTPoints=nprofs
4092 #dataOut.normFactor = nprofs
4093 dataOut.channelList = channelList
4094 nChan = len(channelList)
4095 #dataOut.ippFactor=1
4096 #ipp = ipp/150*1.e-3
4097 vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
4098 #dataOut.ippSeconds=ipp
4099 absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf
4100 print('sale 2',dataOut.ippSeconds,M,N)
4101 print('Empieza procesamiento offline')
4102 if path != None:
4103 sys.path.append(path)
4104 self.library = importlib.import_module(file)
4105 constants = self.library.setConstants(dataOut)
4106 constants['M'] = M
4107 dataOut.constants = constants
4108 if filec != None:
4109 self.weightf = importlib.import_module(filec)
4110
4111 groupArray = numpy.array(groupList)
4112 dataOut.groupList = groupArray
4113 nGroups = groupArray.shape[0]
4114 #List of possible combinations
4115 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
4116 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
4117 if dataOut.data_paramC is None:
4118 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
4119 dataOut.data_snr1_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan
4120 # dataOut.smooth_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan
4121
4122 for i in range(nGroups):
4123 coord = groupArray[i,:]
4124 #Input data array
4125 data = dataOut.data_spc[coord,:,:]/(M*N)
4126 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
4127
4128 #Cross Spectra data array for Covariance Matrixes
4129 ind = 0
4130 for pairs in listComb:
4131 pairsSel = numpy.array([coord[x],coord[y]])
4132 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
4133 ind += 1
4134 dataCross = dataOut.data_cspc[indCross,:,:]/(M*N)
4135 dataCross = dataCross**2
4136 nhei = nHeights
4137 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
4138 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
4139 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
4140 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
4141 n0 = n0i
4142 n1= n1i
4143 my_noises[2*i+0] = n0
4144 my_noises[2*i+1] = n1
4145 snrth = -13 #-13.0 # -4 -16 -25
4146 snrth = 10**(snrth/10.0)
4147 jvelr = numpy.zeros(nHeights, dtype = 'float')
4148 #snr0 = numpy.zeros(nHeights, dtype = 'float')
4149 #snr1 = numpy.zeros(nHeights, dtype = 'float')
4150 hvalid = [0]
4151
4152 coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:])
4153
4154 for h in range(nHeights):
4155 smooth = clean_num_aver[i+1,h]
4156 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
4157 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
4158 signal0 = signalpn0-n0
4159 signal1 = signalpn1-n1
4160 snr0 = numpy.sum(signal0/n0)/(nProf-1)
4161 snr1 = numpy.sum(signal1/n1)/(nProf-1)
4162 #jmax0 = MAX(signal0,maxp0)
4163 #jmax1 = MAX(signal1,maxp1)
4164 gamma = coh2[:,h]
4165
4166 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
4167
4168 if len(indxs) >0:
4169 if numpy.nanmean(gamma) > 0.07:
4170 maxp0 = numpy.argmax(signal0*gamma)
4171 maxp1 = numpy.argmax(signal1*gamma)
4172 #print('usa gamma',numpy.nanmean(gamma))
4173 else:
4174 maxp0 = numpy.argmax(signal0)
4175 maxp1 = numpy.argmax(signal1)
4176 jvelr[h] = (absc[maxp0]+absc[maxp1])/2.
4177 else: jvelr[h] = absc[0]
4178 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
4179 #print(maxp0,absc[maxp0],snr0,jvelr[h])
4180
4181 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
4182 else: fd0 = numpy.nan
4183 print(fd0)
4184 for h in range(nHeights):
4185 d = data[:,h]
4186 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
4187 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
4188 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
4189 signal0 = signalpn0-n0
4190 signal1 = signalpn1-n1
4191 snr0 = numpy.sum(signal0/n0)/(nProf-1)
4192 snr1 = numpy.sum(signal1/n1)/(nProf-1)
4193
4194 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
4195 #Covariance Matrix
4196 D = numpy.diag(d**2)
4197 ind = 0
4198 for pairs in listComb:
4199 #Coordinates in Covariance Matrix
4200 x = pairs[0]
4201 y = pairs[1]
4202 #Channel Index
4203 S12 = dataCross[ind,:,h]
4204 D12 = numpy.diag(S12)
4205 #Completing Covariance Matrix with Cross Spectras
4206 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
4207 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
4208 ind += 1
4209 diagD = numpy.zeros(256)
4210
4211 try:
4212 Dinv=numpy.linalg.inv(D)
4213 L=numpy.linalg.cholesky(Dinv)
4214 except:
4215 Dinv = D*numpy.nan
4216 L= D*numpy.nan
4217 LT=L.T
4218
4219 dp = numpy.dot(LT,d)
4220 #Initial values
4221 data_spc = dataOut.data_spc[coord,:,h]
4222 w = data_spc/data_spc
4223 if filec != None:
4224 w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
4225 if (h>6) and (error1[3]<25):
4226 p0 = dataOut.data_param[i,:,h-1].copy()
4227 else:
4228 p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i)
4229 p0[3] = fd0
4230 if filec != None:
4231 p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)
4232
4233 try:
4234 #Least Squares
4235 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
4236 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
4237 #Chi square error
4238 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
4239 #Error with Jacobian
4240 error1 = self.library.errorFunction(minp,constants,LT)
4241
4242 except:
4243 minp = p0*numpy.nan
4244 error0 = numpy.nan
4245 error1 = p0*numpy.nan
4246 else :
4247 data_spc = dataOut.data_spc[coord,:,h]
4248 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
4249 minp = p0*numpy.nan
4250 error0 = numpy.nan
4251 error1 = p0*numpy.nan
4252 if dataOut.data_param is None:
4253 dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
4254 dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
4255
4256 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
4257 dataOut.data_param[i,:,h] = minp
4258 dataOut.data_snr1_i[i*2,h] = numpy.sum(signalpn0/(nProf-1))/n0
4259 dataOut.data_snr1_i[i*2+1,h] = numpy.sum(signalpn1/(nProf-1))/n1
4260
4261 for ht in range(nHeights-1) :
4262 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
4263 dataOut.data_paramC[4*i,ht,1] = smooth
4264 signalpn0 = (clean_coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra
4265 signalpn1 = (clean_coh_spectra[i*2+1,1:(nProf-0),ht])/smooth
4266
4267 val0 = (signalpn0 > 0).nonzero()
4268 val0 = val0[0]
4269
4270 if len(val0) == 0 : val0_npoints = nProf
4271 else : val0_npoints = len(val0)
4272
4273 val1 = (signalpn1 > 0).nonzero()
4274 val1 = val1[0]
4275 if len(val1) == 0 : val1_npoints = nProf
4276 else : val1_npoints = len(val1)
4277
4278 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
4279 dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1
4280
4281 signal0 = (signalpn0-n0)
4282 vali = (signal0 < 0).nonzero()
4283 vali = vali[0]
4284 if len(vali) > 0 : signal0[vali] = 0
4285 signal1 = (signalpn1-n1)
4286 vali = (signal1 < 0).nonzero()
4287 vali = vali[0]
4288 if len(vali) > 0 : signal1[vali] = 0
4289 snr0 = numpy.sum(signal0/n0)/(nProf-1)
4290 snr1 = numpy.sum(signal1/n1)/(nProf-1)
4291 doppler = absc[1:]
4292 if snr0 >= snrth and snr1 >= snrth and smooth :
4293 signalpn0_n0 = signalpn0
4294 signalpn0_n0[val0] = signalpn0[val0] - n0
4295 mom0 = self.moments(doppler,signalpn0-n0,nProf)
4296
4297 signalpn1_n1 = signalpn1
4298 signalpn1_n1[val1] = signalpn1[val1] - n1
4299 mom1 = self.moments(doppler,signalpn1_n1,nProf)
4300 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
4301 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
4302
4303 dataOut.data_spc = jspectra
4304 dataOut.spc_noise = my_noises*nProf*M
4305
4306 if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M
4307 if 0:
4308 listChannels = groupArray.reshape((groupArray.size))
4309 listChannels.sort()
4310 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels])
4311 #print(dataOut.data_snr1_i)
4312 # Adding coherent echoes from possible satellites.
4313 #sat_spectra = numpy.zeros((nChan,nProf,nHei), dtype=float)
4314 #sat_spectra = sat_spectra[*,*,anal_header.channels]
4315 isat_spectra = numpy.zeros([2,int(nChan/2),nProf,nhei], dtype=float)
4316
4317 sat_fits = numpy.zeros([4,nhei], dtype=float)
4318 noises = my_noises/nProf
4319 #nchan2 = int(nChan/2)
4320 for beam in range(int(nChan/2)-0) :
4321 n0 = noises[2*beam]
4322 n1 = noises[2*beam+1]
4323 isat_spectra[0:2,beam,:,:] = dataOut.sat_spectra[2*beam +0:2*beam+2 ,:,:]
4324
4325 for ht in range(nhei-1) :
4326 signalpn0 = isat_spectra[0,beam,:,ht]
4327 signalpn0 = numpy.reshape(signalpn0,nProf)
4328 signalpn1 = isat_spectra[1,beam,:,ht]
4329 signalpn1 = numpy.reshape(signalpn1,nProf)
4330
4331 cval0 = len((signalpn0 > 0).nonzero()[0])
4332 if cval0 == 0 : val0_npoints = nProf
4333 else: val0_npoints = cval0
4334
4335 cval1 = len((signalpn1 > 0).nonzero()[0])
4336 if cval1 == 0 : val1_npoints = nProf
4337 else: val1_npoints = cval1
4338
4339 sat_fits[0+2*beam,ht] = numpy.sum(signalpn0/(val0_npoints*nProf))/n0
4340 sat_fits[1+2*beam,ht] = numpy.sum(signalpn1/(val1_npoints*nProf))/n1
4341
4342 dataOut.sat_fits = sat_fits
4343 return dataOut
4344
4345 def __residFunction(self, p, dp, LT, constants):
4346
4347 fm = self.library.modelFunction(p, constants)
4348 fmp=numpy.dot(LT,fm)
4349 return dp-fmp
4350
4351 def __getSNR(self, z, noise):
4352
4353 avg = numpy.average(z, axis=1)
4354 SNR = (avg.T-noise)/noise
4355 SNR = SNR.T
4356 return SNR
4357
4358 def __chisq(self, p, chindex, hindex):
4359 #similar to Resid but calculates CHI**2
4360 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
4361 dp=numpy.dot(LT,d)
4362 fmp=numpy.dot(LT,fm)
4363 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
4364 return chisq
4365
4366 class WindProfiler_V0(Operation):
4367
4368 __isConfig = False
4369
4370 __initime = None
4371 __lastdatatime = None
4372 __integrationtime = None
4373
4374 __buffer = None
4375
4376 __dataReady = False
4377
4378 __firstdata = None
4379
4380 n = None
4381
4382 def __init__(self):
4383 Operation.__init__(self)
4384
4385 def __calculateCosDir(self, elev, azim):
4386 zen = (90 - elev)*numpy.pi/180
4387 azim = azim*numpy.pi/180
4388 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
4389 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
4390
4391 signX = numpy.sign(numpy.cos(azim))
4392 signY = numpy.sign(numpy.sin(azim))
4393
4394 cosDirX = numpy.copysign(cosDirX, signX)
4395 cosDirY = numpy.copysign(cosDirY, signY)
4396 return cosDirX, cosDirY
4397
4398 def __calculateAngles(self, theta_x, theta_y, azimuth):
4399
4400 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
4401 zenith_arr = numpy.arccos(dir_cosw)
4402 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
4403
4404 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
4405 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
4406
4407 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
4408
4409 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
4410
4411 if horOnly:
4412 A = numpy.c_[dir_cosu,dir_cosv]
4413 else:
4414 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
4415 A = numpy.asmatrix(A)
4416 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
4417
4418 return A1
4419
4420 def __correctValues(self, heiRang, phi, velRadial, SNR):
4421 listPhi = phi.tolist()
4422 maxid = listPhi.index(max(listPhi))
4423 minid = listPhi.index(min(listPhi))
4424
4425 rango = list(range(len(phi)))
4426 # rango = numpy.delete(rango,maxid)
4427
4428 heiRang1 = heiRang*math.cos(phi[maxid])
4429 heiRangAux = heiRang*math.cos(phi[minid])
4430 indOut = (heiRang1 < heiRangAux[0]).nonzero()
4431 heiRang1 = numpy.delete(heiRang1,indOut)
4432
4433 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
4434 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
4435
4436 for i in rango:
4437 x = heiRang*math.cos(phi[i])
4438 y1 = velRadial[i,:]
4439 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
4440
4441 x1 = heiRang1
4442 y11 = f1(x1)
4443
4444 y2 = SNR[i,:]
4445 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
4446 y21 = f2(x1)
4447
4448 velRadial1[i,:] = y11
4449 SNR1[i,:] = y21
4450
4451 return heiRang1, velRadial1, SNR1
4452
4453 def __calculateVelUVW(self, A, velRadial):
4454
4455 #Operacion Matricial
4456 # velUVW = numpy.zeros((velRadial.shape[1],3))
4457 # for ind in range(velRadial.shape[1]):
4458 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
4459 # velUVW = velUVW.transpose()
4460 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
4461 velUVW[:,:] = numpy.dot(A,velRadial)
4462
4463
4464 return velUVW
4465
4466 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
4467
4468 def techniqueDBS(self, kwargs):
4469 """
4470 Function that implements Doppler Beam Swinging (DBS) technique.
4471
4472 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
4473 Direction correction (if necessary), Ranges and SNR
4474
4475 Output: Winds estimation (Zonal, Meridional and Vertical)
4476
4477 Parameters affected: Winds, height range, SNR
4478 """
4479 velRadial0 = kwargs['velRadial']
4480 heiRang = kwargs['heightList']
4481 SNR0 = kwargs['SNR']
4482
4483 if 'dirCosx' in kwargs and 'dirCosy' in kwargs:
4484 theta_x = numpy.array(kwargs['dirCosx'])
4485 theta_y = numpy.array(kwargs['dirCosy'])
4486 else:
4487 elev = numpy.array(kwargs['elevation'])
4488 azim = numpy.array(kwargs['azimuth'])
4489 theta_x, theta_y = self.__calculateCosDir(elev, azim)
4490 azimuth = kwargs['correctAzimuth']
4491 if 'horizontalOnly' in kwargs:
4492 horizontalOnly = kwargs['horizontalOnly']
4493 else: horizontalOnly = False
4494 if 'correctFactor' in kwargs:
4495 correctFactor = kwargs['correctFactor']
4496 else: correctFactor = 1
4497 if 'channelList' in kwargs:
4498 channelList = kwargs['channelList']
4499 if len(channelList) == 2:
4500 horizontalOnly = True
4501 arrayChannel = numpy.array(channelList)
4502 param = param[arrayChannel,:,:]
4503 theta_x = theta_x[arrayChannel]
4504 theta_y = theta_y[arrayChannel]
4505
4506 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
4507 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
4508 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
4509
4510 #Calculo de Componentes de la velocidad con DBS
4511 winds = self.__calculateVelUVW(A,velRadial1)
4512
4513 return winds, heiRang1, SNR1
4514
4515 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
4516
4517 nPairs = len(pairs_ccf)
4518 posx = numpy.asarray(posx)
4519 posy = numpy.asarray(posy)
4520
4521 #Rotacion Inversa para alinear con el azimuth
4522 if azimuth!= None:
4523 azimuth = azimuth*math.pi/180
4524 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
4525 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
4526 else:
4527 posx1 = posx
4528 posy1 = posy
4529
4530 #Calculo de Distancias
4531 distx = numpy.zeros(nPairs)
4532 disty = numpy.zeros(nPairs)
4533 dist = numpy.zeros(nPairs)
4534 ang = numpy.zeros(nPairs)
4535
4536 for i in range(nPairs):
4537 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
4538 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
4539 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
4540 ang[i] = numpy.arctan2(disty[i],distx[i])
4541
4542 return distx, disty, dist, ang
4543 #Calculo de Matrices
4544 # nPairs = len(pairs)
4545 # ang1 = numpy.zeros((nPairs, 2, 1))
4546 # dist1 = numpy.zeros((nPairs, 2, 1))
4547 #
4548 # for j in range(nPairs):
4549 # dist1[j,0,0] = dist[pairs[j][0]]
4550 # dist1[j,1,0] = dist[pairs[j][1]]
4551 # ang1[j,0,0] = ang[pairs[j][0]]
4552 # ang1[j,1,0] = ang[pairs[j][1]]
4553 #
4554 # return distx,disty, dist1,ang1
4555
4556
4557 def __calculateVelVer(self, phase, lagTRange, _lambda):
4558
4559 Ts = lagTRange[1] - lagTRange[0]
4560 velW = -_lambda*phase/(4*math.pi*Ts)
4561
4562 return velW
4563
4564 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
4565 nPairs = tau1.shape[0]
4566 nHeights = tau1.shape[1]
4567 vel = numpy.zeros((nPairs,3,nHeights))
4568 dist1 = numpy.reshape(dist, (dist.size,1))
4569
4570 angCos = numpy.cos(ang)
4571 angSin = numpy.sin(ang)
4572
4573 vel0 = dist1*tau1/(2*tau2**2)
4574 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
4575 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
4576
4577 ind = numpy.where(numpy.isinf(vel))
4578 vel[ind] = numpy.nan
4579
4580 return vel
4581
4582 # def __getPairsAutoCorr(self, pairsList, nChannels):
4583 #
4584 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
4585 #
4586 # for l in range(len(pairsList)):
4587 # firstChannel = pairsList[l][0]
4588 # secondChannel = pairsList[l][1]
4589 #
4590 # #Obteniendo pares de Autocorrelacion
4591 # if firstChannel == secondChannel:
4592 # pairsAutoCorr[firstChannel] = int(l)
4593 #
4594 # pairsAutoCorr = pairsAutoCorr.astype(int)
4595 #
4596 # pairsCrossCorr = range(len(pairsList))
4597 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
4598 #
4599 # return pairsAutoCorr, pairsCrossCorr
4600
4601 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
4602 def techniqueSA(self, kwargs):
4603
4604 """
4605 Function that implements Spaced Antenna (SA) technique.
4606
4607 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
4608 Direction correction (if necessary), Ranges and SNR
4609
4610 Output: Winds estimation (Zonal, Meridional and Vertical)
4611
4612 Parameters affected: Winds
4613 """
4614 position_x = kwargs['positionX']
4615 position_y = kwargs['positionY']
4616 azimuth = kwargs['azimuth']
4617
4618 if 'correctFactor' in kwargs:
4619 correctFactor = kwargs['correctFactor']
4620 else:
4621 correctFactor = 1
4622
4623 groupList = kwargs['groupList']
4624 pairs_ccf = groupList[1]
4625 tau = kwargs['tau']
4626 _lambda = kwargs['_lambda']
4627
4628 #Cross Correlation pairs obtained
4629 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
4630 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
4631 # pairsSelArray = numpy.array(pairsSelected)
4632 # pairs = []
4633 #
4634 # #Wind estimation pairs obtained
4635 # for i in range(pairsSelArray.shape[0]/2):
4636 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
4637 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
4638 # pairs.append((ind1,ind2))
4639
4640 indtau = tau.shape[0]/2
4641 tau1 = tau[:indtau,:]
4642 tau2 = tau[indtau:-1,:]
4643 # tau1 = tau1[pairs,:]
4644 # tau2 = tau2[pairs,:]
4645 phase1 = tau[-1,:]
4646
4647 #---------------------------------------------------------------------
4648 #Metodo Directo
4649 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
4650 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
4651 winds = stats.nanmean(winds, axis=0)
4652 #---------------------------------------------------------------------
4653 #Metodo General
4654 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
4655 # #Calculo Coeficientes de Funcion de Correlacion
4656 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
4657 # #Calculo de Velocidades
4658 # winds = self.calculateVelUV(F,G,A,B,H)
4659
4660 #---------------------------------------------------------------------
4661 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
4662 winds = correctFactor*winds
4663 return winds
4664
4665 def __checkTime(self, currentTime, paramInterval, outputInterval):
4666
4667 dataTime = currentTime + paramInterval
4668 deltaTime = dataTime - self.__initime
4669
4670 if deltaTime >= outputInterval or deltaTime < 0:
4671 self.__dataReady = True
4672 return
4673
4674 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
4675 '''
4676 Function that implements winds estimation technique with detected meteors.
4677
4678 Input: Detected meteors, Minimum meteor quantity to wind estimation
4679
4680 Output: Winds estimation (Zonal and Meridional)
4681
4682 Parameters affected: Winds
4683 '''
4684 #Settings
4685 nInt = (heightMax - heightMin)/2
4686 nInt = int(nInt)
4687 winds = numpy.zeros((2,nInt))*numpy.nan
4688
4689 #Filter errors
4690 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
4691 finalMeteor = arrayMeteor[error,:]
4692
4693 #Meteor Histogram
4694 finalHeights = finalMeteor[:,2]
4695 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
4696 nMeteorsPerI = hist[0]
4697 heightPerI = hist[1]
4698
4699 #Sort of meteors
4700 indSort = finalHeights.argsort()
4701 finalMeteor2 = finalMeteor[indSort,:]
4702
4703 # Calculating winds
4704 ind1 = 0
4705 ind2 = 0
4706
4707 for i in range(nInt):
4708 nMet = nMeteorsPerI[i]
4709 ind1 = ind2
4710 ind2 = ind1 + nMet
4711
4712 meteorAux = finalMeteor2[ind1:ind2,:]
4713
4714 if meteorAux.shape[0] >= meteorThresh:
4715 vel = meteorAux[:, 6]
4716 zen = meteorAux[:, 4]*numpy.pi/180
4717 azim = meteorAux[:, 3]*numpy.pi/180
4718
4719 n = numpy.cos(zen)
4720 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
4721 # l = m*numpy.tan(azim)
4722 l = numpy.sin(zen)*numpy.sin(azim)
4723 m = numpy.sin(zen)*numpy.cos(azim)
4724
4725 A = numpy.vstack((l, m)).transpose()
4726 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
4727 windsAux = numpy.dot(A1, vel)
4728
4729 winds[0,i] = windsAux[0]
4730 winds[1,i] = windsAux[1]
4731
4732 return winds, heightPerI[:-1]
4733
4734 def techniqueNSM_SA(self, **kwargs):
4735 metArray = kwargs['metArray']
4736 heightList = kwargs['heightList']
4737 timeList = kwargs['timeList']
4738
4739 rx_location = kwargs['rx_location']
4740 groupList = kwargs['groupList']
4741 azimuth = kwargs['azimuth']
4742 dfactor = kwargs['dfactor']
4743 k = kwargs['k']
4744
4745 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
4746 d = dist*dfactor
4747 #Phase calculation
4748 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
4749
4750 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
4751
4752 velEst = numpy.zeros((heightList.size,2))*numpy.nan
4753 azimuth1 = azimuth1*numpy.pi/180
4754
4755 for i in range(heightList.size):
4756 h = heightList[i]
4757 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
4758 metHeight = metArray1[indH,:]
4759 if metHeight.shape[0] >= 2:
4760 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
4761 iazim = metHeight[:,1].astype(int)
4762 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
4763 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
4764 A = numpy.asmatrix(A)
4765 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
4766 velHor = numpy.dot(A1,velAux)
4767
4768 velEst[i,:] = numpy.squeeze(velHor)
4769 return velEst
4770
4771 def __getPhaseSlope(self, metArray, heightList, timeList):
4772 meteorList = []
4773 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
4774 #Putting back together the meteor matrix
4775 utctime = metArray[:,0]
4776 uniqueTime = numpy.unique(utctime)
4777
4778 phaseDerThresh = 0.5
4779 ippSeconds = timeList[1] - timeList[0]
4780 sec = numpy.where(timeList>1)[0][0]
4781 nPairs = metArray.shape[1] - 6
4782 nHeights = len(heightList)
4783
4784 for t in uniqueTime:
4785 metArray1 = metArray[utctime==t,:]
4786 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
4787 tmet = metArray1[:,1].astype(int)
4788 hmet = metArray1[:,2].astype(int)
4789
4790 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
4791 metPhase[:,:] = numpy.nan
4792 metPhase[:,hmet,tmet] = metArray1[:,6:].T
4793
4794 #Delete short trails
4795 metBool = ~numpy.isnan(metPhase[0,:,:])
4796 heightVect = numpy.sum(metBool, axis = 1)
4797 metBool[heightVect<sec,:] = False
4798 metPhase[:,heightVect<sec,:] = numpy.nan
4799
4800 #Derivative
4801 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
4802 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
4803 metPhase[phDerAux] = numpy.nan
4804
4805 #--------------------------METEOR DETECTION -----------------------------------------
4806 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
4807
4808 for p in numpy.arange(nPairs):
4809 phase = metPhase[p,:,:]
4810 phDer = metDer[p,:,:]
4811
4812 for h in indMet:
4813 height = heightList[h]
4814 phase1 = phase[h,:] #82
4815 phDer1 = phDer[h,:]
4816
4817 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
4818
4819 indValid = numpy.where(~numpy.isnan(phase1))[0]
4820 initMet = indValid[0]
4821 endMet = 0
4822
4823 for i in range(len(indValid)-1):
4824
4825 #Time difference
4826 inow = indValid[i]
4827 inext = indValid[i+1]
4828 idiff = inext - inow
4829 #Phase difference
4830 phDiff = numpy.abs(phase1[inext] - phase1[inow])
4831
4832 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
4833 sizeTrail = inow - initMet + 1
4834 if sizeTrail>3*sec: #Too short meteors
4835 x = numpy.arange(initMet,inow+1)*ippSeconds
4836 y = phase1[initMet:inow+1]
4837 ynnan = ~numpy.isnan(y)
4838 x = x[ynnan]
4839 y = y[ynnan]
4840 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
4841 ylin = x*slope + intercept
4842 rsq = r_value**2
4843 if rsq > 0.5:
4844 vel = slope#*height*1000/(k*d)
4845 estAux = numpy.array([utctime,p,height, vel, rsq])
4846 meteorList.append(estAux)
4847 initMet = inext
4848 metArray2 = numpy.array(meteorList)
4849
4850 return metArray2
4851
4852 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
4853
4854 azimuth1 = numpy.zeros(len(pairslist))
4855 dist = numpy.zeros(len(pairslist))
4856
4857 for i in range(len(rx_location)):
4858 ch0 = pairslist[i][0]
4859 ch1 = pairslist[i][1]
4860
4861 diffX = rx_location[ch0][0] - rx_location[ch1][0]
4862 diffY = rx_location[ch0][1] - rx_location[ch1][1]
4863 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
4864 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
4865
4866 azimuth1 -= azimuth0
4867 return azimuth1, dist
4868
4869 def techniqueNSM_DBS(self, **kwargs):
4870 metArray = kwargs['metArray']
4871 heightList = kwargs['heightList']
4872 timeList = kwargs['timeList']
4873 azimuth = kwargs['azimuth']
4874 theta_x = numpy.array(kwargs['theta_x'])
4875 theta_y = numpy.array(kwargs['theta_y'])
4876
4877 utctime = metArray[:,0]
4878 cmet = metArray[:,1].astype(int)
4879 hmet = metArray[:,3].astype(int)
4880 SNRmet = metArray[:,4]
4881 vmet = metArray[:,5]
4882 spcmet = metArray[:,6]
4883
4884 nChan = numpy.max(cmet) + 1
4885 nHeights = len(heightList)
4886
4887 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
4888 hmet = heightList[hmet]
4889 h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights
4890
4891 velEst = numpy.zeros((heightList.size,2))*numpy.nan
4892
4893 for i in range(nHeights - 1):
4894 hmin = heightList[i]
4895 hmax = heightList[i + 1]
4896
4897 thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10)
4898 indthisH = numpy.where(thisH)
4899
4900 if numpy.size(indthisH) > 3:
4901
4902 vel_aux = vmet[thisH]
4903 chan_aux = cmet[thisH]
4904 cosu_aux = dir_cosu[chan_aux]
4905 cosv_aux = dir_cosv[chan_aux]
4906 cosw_aux = dir_cosw[chan_aux]
4907
4908 nch = numpy.size(numpy.unique(chan_aux))
4909 if nch > 1:
4910 A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True)
4911 velEst[i,:] = numpy.dot(A,vel_aux)
4912
4913 return velEst
4914
4915 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
4916
4917 param = dataOut.data_param
4918 #if dataOut.abscissaList != None:
4919 if numpy.any(dataOut.abscissaList):
4920 absc = dataOut.abscissaList[:-1]
4921 # noise = dataOut.noise
4922 heightList = dataOut.heightList
4923 SNR = dataOut.data_snr
4924
4925 if technique == 'DBS':
4926
4927 kwargs['velRadial'] = param[:,1,:] #Radial velocity
4928 kwargs['heightList'] = heightList
4929 kwargs['SNR'] = SNR
4930
4931 dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function
4932 dataOut.utctimeInit = dataOut.utctime
4933 dataOut.outputInterval = dataOut.paramInterval
4934
4935 elif technique == 'SA':
4936
4937 #Parameters
4938 # position_x = kwargs['positionX']
4939 # position_y = kwargs['positionY']
4940 # azimuth = kwargs['azimuth']
4941 #
4942 # if kwargs.has_key('crosspairsList'):
4943 # pairs = kwargs['crosspairsList']
4944 # else:
4945 # pairs = None
4946 #
4947 # if kwargs.has_key('correctFactor'):
4948 # correctFactor = kwargs['correctFactor']
4949 # else:
4950 # correctFactor = 1
4951
4952 # tau = dataOut.data_param
4953 # _lambda = dataOut.C/dataOut.frequency
4954 # pairsList = dataOut.groupList
4955 # nChannels = dataOut.nChannels
4956
4957 kwargs['groupList'] = dataOut.groupList
4958 kwargs['tau'] = dataOut.data_param
4959 kwargs['_lambda'] = dataOut.C/dataOut.frequency
4960 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
4961 dataOut.data_output = self.techniqueSA(kwargs)
4962 dataOut.utctimeInit = dataOut.utctime
4963 dataOut.outputInterval = dataOut.timeInterval
4964
4965 elif technique == 'Meteors':
4966 dataOut.flagNoData = True
4967 self.__dataReady = False
4968
4969 if 'nHours' in kwargs:
4970 nHours = kwargs['nHours']
4971 else:
4972 nHours = 1
3952
4973
3953 def __getSNR(self, z, noise):
4974 if 'meteorsPerBin' in kwargs:
4975 meteorThresh = kwargs['meteorsPerBin']
4976 else:
4977 meteorThresh = 6
3954
4978
3955 avg = numpy.average(z, axis=1)
4979 if 'hmin' in kwargs:
3956 SNR = (avg.T-noise)/noise
4980 hmin = kwargs['hmin']
3957 SNR = SNR.T
4981 else: hmin = 70
3958 return SNR
4982 if 'hmax' in kwargs:
4983 hmax = kwargs['hmax']
4984 else: hmax = 110
3959
4985
3960 def __chisq(self, p, chindex, hindex):
4986 dataOut.outputInterval = nHours*3600
3961 #similar to Resid but calculates CHI**2
4987
3962 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
4988 if self.__isConfig == False:
3963 dp=numpy.dot(LT,d)
4989 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
3964 fmp=numpy.dot(LT,fm)
4990 #Get Initial LTC time
3965 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
4991 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
3966 return chisq
4992 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
4993
4994 self.__isConfig = True
4995
4996 if self.__buffer is None:
4997 self.__buffer = dataOut.data_param
4998 self.__firstdata = copy.copy(dataOut)
4999
5000 else:
5001 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
5002
5003 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
5004
5005 if self.__dataReady:
5006 dataOut.utctimeInit = self.__initime
5007
5008 self.__initime += dataOut.outputInterval #to erase time offset
5009
5010 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
5011 dataOut.flagNoData = False
5012 self.__buffer = None
5013
5014 elif technique == 'Meteors1':
5015 dataOut.flagNoData = True
5016 self.__dataReady = False
5017
5018 if 'nMins' in kwargs:
5019 nMins = kwargs['nMins']
5020 else: nMins = 20
5021 if 'rx_location' in kwargs:
5022 rx_location = kwargs['rx_location']
5023 else: rx_location = [(0,1),(1,1),(1,0)]
5024 if 'azimuth' in kwargs:
5025 azimuth = kwargs['azimuth']
5026 else: azimuth = 51.06
5027 if 'dfactor' in kwargs:
5028 dfactor = kwargs['dfactor']
5029 if 'mode' in kwargs:
5030 mode = kwargs['mode']
5031 if 'theta_x' in kwargs:
5032 theta_x = kwargs['theta_x']
5033 if 'theta_y' in kwargs:
5034 theta_y = kwargs['theta_y']
5035 else: mode = 'SA'
5036
5037 #Borrar luego esto
5038 if dataOut.groupList is None:
5039 dataOut.groupList = [(0,1),(0,2),(1,2)]
5040 groupList = dataOut.groupList
5041 C = 3e8
5042 freq = 50e6
5043 lamb = C/freq
5044 k = 2*numpy.pi/lamb
5045
5046 timeList = dataOut.abscissaList
5047 heightList = dataOut.heightList
5048
5049 if self.__isConfig == False:
5050 dataOut.outputInterval = nMins*60
5051 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
5052 #Get Initial LTC time
5053 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
5054 minuteAux = initime.minute
5055 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
5056 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
5057
5058 self.__isConfig = True
5059
5060 if self.__buffer is None:
5061 self.__buffer = dataOut.data_param
5062 self.__firstdata = copy.copy(dataOut)
5063
5064 else:
5065 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
5066
5067 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
5068
5069 if self.__dataReady:
5070 dataOut.utctimeInit = self.__initime
5071 self.__initime += dataOut.outputInterval #to erase time offset
5072
5073 metArray = self.__buffer
5074 if mode == 'SA':
5075 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
5076 elif mode == 'DBS':
5077 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y)
5078 dataOut.data_output = dataOut.data_output.T
5079 dataOut.flagNoData = False
5080 self.__buffer = None
5081
5082 return
3967
5083
3968 class WindProfiler(Operation):
5084 class WindProfiler(Operation):
3969
5085
@@ -4137,7 +5253,6 class WindProfiler(Operation):
4137 return distx, disty, dist, ang
5253 return distx, disty, dist, ang
4138 #Calculo de Matrices
5254 #Calculo de Matrices
4139
5255
4140
4141 def __calculateVelVer(self, phase, lagTRange, _lambda):
5256 def __calculateVelVer(self, phase, lagTRange, _lambda):
4142
5257
4143 Ts = lagTRange[1] - lagTRange[0]
5258 Ts = lagTRange[1] - lagTRange[0]
@@ -4645,15 +5760,36 class EWDriftsEstimation(Operation):
4645
5760
4646 def run(self, dataOut, zenith, zenithCorrection,fileDrifts):
5761 def run(self, dataOut, zenith, zenithCorrection,fileDrifts):
4647
5762
4648 dataOut.lat=-11.95
5763 dataOut.lat = -11.95
4649 dataOut.lon=-76.87
5764 dataOut.lon = -76.87
5765 dataOut.spcst = 0.00666
5766 dataOut.pl = 0.0003
5767 dataOut.cbadn = 3
5768 dataOut.inttms = 300
5769 dataOut.azw = -115.687
5770 dataOut.elw = 86.1095
5771 dataOut.aze = 130.052
5772 dataOut.ele = 87.6558
5773 dataOut.jro14 = numpy.log10(dataOut.spc_noise[0]/dataOut.normFactor)
5774 dataOut.jro15 = numpy.log10(dataOut.spc_noise[1]/dataOut.normFactor)
5775 dataOut.jro16 = numpy.log10(dataOut.spc_noise[2]/dataOut.normFactor)
5776 dataOut.nwlos = numpy.log10(dataOut.spc_noise[3]/dataOut.normFactor)
5777
4650 heiRang = dataOut.heightList
5778 heiRang = dataOut.heightList
4651 velRadial = dataOut.data_param[:,3,:]
5779 velRadial = dataOut.data_param[:,3,:]
4652 velRadialm = dataOut.data_param[:,2:4,:]*-1
5780 velRadialm = dataOut.data_param[:,2:4,:]*-1
5781
4653 rbufc=dataOut.data_paramC[:,:,0]
5782 rbufc=dataOut.data_paramC[:,:,0]
4654 ebufc=dataOut.data_paramC[:,:,1]
5783 ebufc=dataOut.data_paramC[:,:,1]
4655 SNR = dataOut.data_snr
5784 #SNR = dataOut.data_snr
5785 SNR = dataOut.data_snr1_i
5786 rbufi = dataOut.data_snr1_i
4656 velRerr = dataOut.data_error[:,4,:]
5787 velRerr = dataOut.data_error[:,4,:]
5788 range1 = dataOut.heightList
5789 nhei = len(range1)
5790
5791 sat_fits = dataOut.sat_fits
5792
4657 channels = dataOut.channelList
5793 channels = dataOut.channelList
4658 nChan = len(channels)
5794 nChan = len(channels)
4659 my_nbeams = nChan/2
5795 my_nbeams = nChan/2
@@ -4662,18 +5798,49 class EWDriftsEstimation(Operation):
4662 else :
5798 else :
4663 moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))
5799 moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))
4664 dataOut.moments=moments
5800 dataOut.moments=moments
5801 #Incoherent
5802 smooth_w = dataOut.clean_num_aver[0,:]
5803 chisq_w = dataOut.data_error[0,0,:]
5804 p_w0 = rbufi[0,:]
5805 p_w1 = rbufi[1,:]
5806
4665 # Coherent
5807 # Coherent
4666 smooth_wC = ebufc[0,:]
5808 smooth_wC = ebufc[0,:]
4667 p_w0C = rbufc[0,:]
5809 p_w0C = rbufc[0,:]
4668 p_w1C = rbufc[1,:]
5810 p_w1C = rbufc[1,:]
4669 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
5811 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
4670 t_wC = rbufc[3,:]
5812 t_wC = rbufc[3,:]
5813 val = (numpy.isfinite(p_w0)==False).nonzero()
5814 p_w0[val]=0
5815 val = (numpy.isfinite(p_w1)==False).nonzero()
5816 p_w1[val]=0
5817 val = (numpy.isfinite(p_w0C)==False).nonzero()
5818 p_w0C[val]=0
5819 val = (numpy.isfinite(p_w1C)==False).nonzero()
5820 p_w1C[val]=0
5821 val = (numpy.isfinite(smooth_w)==False).nonzero()
5822 smooth_w[val]=0
5823 val = (numpy.isfinite(smooth_wC)==False).nonzero()
5824 smooth_wC[val]=0
5825
5826 #p_w0 = (p_w0*smooth_w+p_w0C*smooth_wC)/(smooth_w+smooth_wC)
5827 #p_w1 = (p_w1*smooth_w+p_w1C*smooth_wC)/(smooth_w+smooth_wC)
5828
5829 if len(sat_fits) >0 :
5830 p_w0C = p_w0C + sat_fits[0,:]
5831 p_w1C = p_w1C + sat_fits[1,:]
5832
4671 if my_nbeams == 1:
5833 if my_nbeams == 1:
4672 w = velRadial[0,:]
5834 w = velRadial[0,:]
4673 winds = velRadial.copy()
5835 winds = velRadial.copy()
4674 w_err = velRerr[0,:]
5836 w_err = velRerr[0,:]
4675 snr1 = 10*numpy.log10(SNR[0])
5837 u = w*numpy.nan
5838 u_err = w_err*numpy.nan
5839 p_e0 = p_w0*numpy.nan
5840 p_e1 = p_w1*numpy.nan
5841 #snr1 = 10*numpy.log10(SNR[0])
4676 if my_nbeams == 2:
5842 if my_nbeams == 2:
5843
4677 zenith = numpy.array(zenith)
5844 zenith = numpy.array(zenith)
4678 zenith -= zenithCorrection
5845 zenith -= zenithCorrection
4679 zenith *= numpy.pi/180
5846 zenith *= numpy.pi/180
@@ -4691,59 +5858,101 class EWDriftsEstimation(Operation):
4691 w_e = velRadial1[1,:]
5858 w_e = velRadial1[1,:]
4692 w_w_err = velRerr[0,:]
5859 w_w_err = velRerr[0,:]
4693 w_e_err = velRerr[1,:]
5860 w_e_err = velRerr[1,:]
4694
5861 smooth_e = dataOut.clean_num_aver[2,:]
4695 val = (numpy.isfinite(w_w)==False).nonzero()
5862 chisq_e = dataOut.data_error[1,0,:]
4696 val = val[0]
5863 p_e0 = rbufi[2,:]
4697 bad = val
5864 p_e1 = rbufi[3,:]
4698 if len(bad) > 0 :
5865
4699 w_w[bad] = w_wC[bad]
5866 tini=time.localtime(dataOut.utctime)
4700 w_w_err[bad]= numpy.nan
5867
5868 if tini[3] >= 6 and tini[3] < 18 :
5869 w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
5870 w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
5871 else:
5872 w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
5873 w_wtmp = numpy.where(range1 > 200,w_w,w_wtmp)
5874 w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
5875 w_w_errtmp = numpy.where(range1 > 200,w_w_err,w_w_errtmp)
5876 w_w = w_wtmp
5877 w_w_err = w_w_errtmp
5878
5879 #if my_nbeams == 2:
4701 smooth_eC=ebufc[4,:]
5880 smooth_eC=ebufc[4,:]
4702 p_e0C = rbufc[4,:]
5881 p_e0C = rbufc[4,:]
4703 p_e1C = rbufc[5,:]
5882 p_e1C = rbufc[5,:]
4704 w_eC = rbufc[6,:]*-1
5883 w_eC = rbufc[6,:]*-1
4705 t_eC = rbufc[7,:]
5884 t_eC = rbufc[7,:]
4706 val = (numpy.isfinite(w_e)==False).nonzero()
5885 val = (numpy.isfinite(p_e0)==False).nonzero()
4707 val = val[0]
5886 p_e0[val]=0
4708 bad = val
5887 val = (numpy.isfinite(p_e1)==False).nonzero()
4709 if len(bad) > 0 :
5888 p_e1[val]=0
4710 w_e[bad] = w_eC[bad]
5889 val = (numpy.isfinite(p_e0C)==False).nonzero()
4711 w_e_err[bad]= numpy.nan
5890 p_e0C[val]=0
4712
5891 val = (numpy.isfinite(p_e1C)==False).nonzero()
4713 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
5892 p_e1C[val]=0
4714 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
5893 val = (numpy.isfinite(smooth_e)==False).nonzero()
5894 smooth_e[val]=0
5895 val = (numpy.isfinite(smooth_eC)==False).nonzero()
5896 smooth_eC[val]=0
5897 #p_e0 = (p_e0*smooth_e+p_e0C*smooth_eC)/(smooth_e+smooth_eC)
5898 #p_e1 = (p_e1*smooth_e+p_e1C*smooth_eC)/(smooth_e+smooth_eC)
5899
5900 if len(sat_fits) >0 :
5901 p_e0C = p_e0C + sat_fits[2,:]
5902 p_e1C = p_e1C + sat_fits[3,:]
5903
5904 if tini[3] >= 6 and tini[3] < 18 :
5905 w_etmp = numpy.where(numpy.isfinite(w_eC)==True,w_eC,w_e)
5906 w_e_errtmp = numpy.where(numpy.isfinite(w_eC)==True,numpy.nan,w_e_err)
5907 else:
5908 w_etmp = numpy.where(numpy.isfinite(w_eC)==True,w_eC,w_e)
5909 w_etmp = numpy.where(range1 > 200,w_e,w_etmp)
5910 w_e_errtmp = numpy.where(numpy.isfinite(w_eC)==True,numpy.nan,w_e_err)
5911 w_e_errtmp = numpy.where(range1 > 200,w_e_err,w_e_errtmp)
5912 w_e = w_etmp
5913 w_e_err = w_e_errtmp
5914
5915 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
5916 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
4715
5917
4716 w_err = numpy.sqrt((w_w_err*numpy.sin(bet))**2.+(w_e_err*numpy.sin(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
5918 w_err = numpy.sqrt((w_w_err*numpy.sin(bet))**2.+(w_e_err*numpy.sin(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
4717 u_err = numpy.sqrt((w_w_err*numpy.cos(bet))**2.+(w_e_err*numpy.cos(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
5919 u_err = numpy.sqrt((w_w_err*numpy.cos(bet))**2.+(w_e_err*numpy.cos(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
4718
4719 winds = numpy.vstack((w,u))
4720
5920
5921 winds = numpy.vstack((w,u))
4721 dataOut.heightList = heiRang1
5922 dataOut.heightList = heiRang1
4722 snr1 = 10*numpy.log10(SNR1[0])
5923 #snr1 = 10*numpy.log10(SNR1[0])
4723 dataOut.data_output = winds
5924 dataOut.data_output = winds
4724 #snr1 = 10*numpy.log10(SNR1[0])# estaba comentado
5925 range1 = dataOut.heightList
4725 dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0]))
5926 nhei = len(range1)
4726 #print("data_snr1",dataOut.data_snr1)
5927 #print('alt ',range1*numpy.sin(86.1*numpy.pi/180))
5928 #print(numpy.min([dataOut.eldir7,dataOut.eldir8]))
5929 galt = range1*numpy.sin(numpy.min([dataOut.elw,dataOut.ele])*numpy.pi/180.)
5930 dataOut.params = numpy.vstack((range1,galt,w,w_err,u,u_err,w_w,w_w_err,w_e,w_e_err,numpy.log10(p_w0),numpy.log10(p_w0C),numpy.log10(p_w1),numpy.log10(p_w1C),numpy.log10(p_e0),numpy.log10(p_e0C),numpy.log10(p_e1),numpy.log10(p_e1C),chisq_w,chisq_e))
5931 #snr1 = 10*numpy.log10(SNR1[0])
5932 #print(min(snr1), max(snr1))
5933 snr1 = numpy.vstack((p_w0,p_w1,p_e0,p_e1))
5934 snr1db = 10*numpy.log10(snr1[0])
5935
5936 #dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0]))
5937 dataOut.data_snr1 = numpy.reshape(snr1db,(1,snr1db.shape[0]))
4727 dataOut.utctimeInit = dataOut.utctime
5938 dataOut.utctimeInit = dataOut.utctime
4728 dataOut.outputInterval = dataOut.timeInterval
5939 dataOut.outputInterval = dataOut.timeInterval
4729
5940
4730 hei_aver0 = 218
5941 hei_aver0 = 218
4731 jrange = 450 #900 para HA drifts
5942 jrange = 450 #900 para HA drifts
4732 deltah = 15.0 #dataOut.spacing(0) 25 HAD
5943 deltah = 15.0 #dataOut.spacing(0) 25 HAD
4733 h0 = 0.0 #dataOut.first_height(0)
5944 h0 = 0.0 #dataOut.first_height(0)
4734 heights = dataOut.heightList
5945
4735 nhei = len(heights)
4736
4737 range1 = numpy.arange(nhei) * deltah + h0
5946 range1 = numpy.arange(nhei) * deltah + h0
4738 jhei = (range1 >= hei_aver0).nonzero()
5947 jhei = (range1 >= hei_aver0).nonzero()
4739 if len(jhei[0]) > 0 :
5948 if len(jhei[0]) > 0 :
4740 h0_index = jhei[0][0] # Initial height for getting averages 218km
5949 h0_index = jhei[0][0] # Initial height for getting averages 218km
4741
5950
4742 mynhei = 7
5951 mynhei = 7
4743 nhei_avg = int(jrange/deltah)
5952 nhei_avg = int(jrange/deltah)
4744 h_avgs = int(nhei_avg/mynhei)
5953 h_avgs = int(nhei_avg/mynhei)
4745 nhei_avg = h_avgs*(mynhei-1)+mynhei
5954 nhei_avg = h_avgs*(mynhei-1)+mynhei
4746
5955
4747 navgs = numpy.zeros(mynhei,dtype='float')
5956 navgs = numpy.zeros(mynhei,dtype='float')
4748 delta_h = numpy.zeros(mynhei,dtype='float')
5957 delta_h = numpy.zeros(mynhei,dtype='float')
4749 range_aver = numpy.zeros(mynhei,dtype='float')
5958 range_aver = numpy.zeros(mynhei,dtype='float')
@@ -4751,11 +5960,11 class EWDriftsEstimation(Operation):
4751 range_aver[ih] = numpy.sum(range1[h0_index+h_avgs*ih:h0_index+h_avgs*(ih+1)-0])/h_avgs
5960 range_aver[ih] = numpy.sum(range1[h0_index+h_avgs*ih:h0_index+h_avgs*(ih+1)-0])/h_avgs
4752 navgs[ih] = h_avgs
5961 navgs[ih] = h_avgs
4753 delta_h[ih] = deltah*h_avgs
5962 delta_h[ih] = deltah*h_avgs
4754
5963
4755 range_aver[mynhei-1] = numpy.sum(range1[h0_index:h0_index+6*h_avgs-0])/(6*h_avgs)
5964 range_aver[mynhei-1] = numpy.sum(range1[h0_index:h0_index+6*h_avgs-0])/(6*h_avgs)
4756 navgs[mynhei-1] = 6*h_avgs
5965 navgs[mynhei-1] = 6*h_avgs
4757 delta_h[mynhei-1] = deltah*6*h_avgs
5966 delta_h[mynhei-1] = deltah*6*h_avgs
4758
5967
4759 wA = w[h0_index:h0_index+nhei_avg-0]
5968 wA = w[h0_index:h0_index+nhei_avg-0]
4760 wA_err = w_err[h0_index:h0_index+nhei_avg-0]
5969 wA_err = w_err[h0_index:h0_index+nhei_avg-0]
4761 for i in range(5) :
5970 for i in range(5) :
@@ -4765,15 +5974,14 class EWDriftsEstimation(Operation):
4765 sigma = numpy.sqrt(1./numpy.nansum(1./errs**2.))
5974 sigma = numpy.sqrt(1./numpy.nansum(1./errs**2.))
4766 wA[6*h_avgs+i] = avg
5975 wA[6*h_avgs+i] = avg
4767 wA_err[6*h_avgs+i] = sigma
5976 wA_err[6*h_avgs+i] = sigma
4768
5977
4769
4770 vals = wA[0:6*h_avgs-0]
5978 vals = wA[0:6*h_avgs-0]
4771 errs=wA_err[0:6*h_avgs-0]
5979 errs=wA_err[0:6*h_avgs-0]
4772 avg = numpy.nansum(vals/errs**2.)/numpy.nansum(1./errs**2)
5980 avg = numpy.nansum(vals/errs**2.)/numpy.nansum(1./errs**2)
4773 sigma = numpy.sqrt(1./numpy.nansum(1./errs**2.))
5981 sigma = numpy.sqrt(1./numpy.nansum(1./errs**2.))
4774 wA[nhei_avg-1] = avg
5982 wA[nhei_avg-1] = avg
4775 wA_err[nhei_avg-1] = sigma
5983 wA_err[nhei_avg-1] = sigma
4776
5984
4777 wA = wA[6*h_avgs:nhei_avg-0]
5985 wA = wA[6*h_avgs:nhei_avg-0]
4778 wA_err=wA_err[6*h_avgs:nhei_avg-0]
5986 wA_err=wA_err[6*h_avgs:nhei_avg-0]
4779 if my_nbeams == 2 :
5987 if my_nbeams == 2 :
@@ -4796,18 +6004,81 class EWDriftsEstimation(Operation):
4796 uA_err[nhei_avg-1] = sigma
6004 uA_err[nhei_avg-1] = sigma
4797 uA = uA[6*h_avgs:nhei_avg-0]
6005 uA = uA[6*h_avgs:nhei_avg-0]
4798 uA_err = uA_err[6*h_avgs:nhei_avg-0]
6006 uA_err = uA_err[6*h_avgs:nhei_avg-0]
4799 dataOut.drifts_avg = numpy.vstack((wA,uA))
6007 dataOut.drifts_avg = numpy.vstack((wA,uA))
6008
4800 if my_nbeams == 1: dataOut.drifts_avg = wA
6009 if my_nbeams == 1: dataOut.drifts_avg = wA
6010 #deltahavg= wA*0.0+deltah
6011 dataOut.range = range1
6012 galtavg = range_aver*numpy.sin(numpy.min([dataOut.elw,dataOut.ele])*numpy.pi/180.)
6013 dataOut.params_avg = numpy.vstack((wA,wA_err,uA,uA_err,range_aver,galtavg,delta_h))
6014
6015 #print('comparando dim de avg ',wA.shape,deltahavg.shape,range_aver.shape)
4801 tini=time.localtime(dataOut.utctime)
6016 tini=time.localtime(dataOut.utctime)
4802 datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
6017 datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
4803 nfile = fileDrifts+'/jro'+datefile+'drifts.txt'
6018 nfile = fileDrifts+'/jro'+datefile+'drifts_sch3.txt'
6019
4804 f1 = open(nfile,'a')
6020 f1 = open(nfile,'a')
4805 datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
6021 datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
4806 driftavgstr=str(dataOut.drifts_avg)
6022 driftavgstr=str(dataOut.drifts_avg)
4807 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
6023 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
4808 numpy.savetxt(f1,dataOut.drifts_avg,fmt='%10.2f')
6024 numpy.savetxt(f1,numpy.reshape(range_aver,(1,len(range_aver))) ,fmt='%10.2f')
6025 numpy.savetxt(f1,dataOut.drifts_avg[:,:],fmt='%10.2f')
6026 f1.close()
6027
6028 swfile = fileDrifts+'/jro'+datefile+'drifts_sw.txt'
6029 f1 = open(swfile,'a')
6030 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
6031 numpy.savetxt(f1,numpy.reshape(heiRang,(1,len(heiRang))),fmt='%10.2f')
6032 numpy.savetxt(f1,dataOut.data_param[:,0,:],fmt='%10.2f')
4809 f1.close()
6033 f1.close()
6034 dataOut.heightListtmp = dataOut.heightList
6035 '''
6036 #Envio data de drifts a mysql
6037 fechad = str(tini[0]).zfill(4)+'-'+str(tini[1]).zfill(2)+'-'+str(tini[2]).zfill(2)+' '+str(tini[3]).zfill(2)+':'+str(tini[4]).zfill(2)+':'+str(0).zfill(2)
6038 mydb = mysql.connector.connect(
6039 host="10.10.110.213",
6040 user="user_clima",
6041 password="5D.bh(B2)Y_wRNz9",
6042 database="clima_espacial"
6043 )
6044
6045 mycursor = mydb.cursor()
6046 #mycursor.execute("CREATE TABLE drifts_vertical (id INT AUTO_INCREMENT PRIMARY KEY, fecha DATETIME(6), Vertical FLOAT(10,2))")
6047
6048 sql = "INSERT INTO drifts_vertical (datetime, value) VALUES (%s, %s)"
6049 if numpy.isfinite(dataOut.drifts_avg[0,6]): vdql = dataOut.drifts_avg[0,6]
6050 else : vdql = 999
6051 val = (fechad, vdql)
6052 mycursor.execute(sql, val)
6053 mydb.commit()
6054 sql = "INSERT INTO drifts_zonal (datetime, value) VALUES (%s, %s)"
6055 if numpy.isfinite(dataOut.drifts_avg[1,6]): zdql = dataOut.drifts_avg[1,6]
6056 else : zdql = 999
6057 val = (fechad, zdql)
6058 mycursor.execute(sql, val)
6059 mydb.commit()
6060
6061 print(mycursor.rowcount, "record inserted.")
6062 '''
6063 return dataOut
4810
6064
6065 class setHeightDrifts(Operation):
6066
6067 def __init__(self):
6068 Operation.__init__(self)
6069 def run(self, dataOut):
6070 #print('h inicial ',dataOut.heightList,dataOut.heightListtmp)
6071 dataOut.heightList = dataOut.heightListtmp
6072 #print('regresa H ',dataOut.heightList)
6073 return dataOut
6074 class setHeightDriftsavg(Operation):
6075
6076 def __init__(self):
6077 Operation.__init__(self)
6078 def run(self, dataOut):
6079 #print('h inicial ',dataOut.heightList)
6080 dataOut.heightList = dataOut.params_avg[4]
6081 #print('cambia H ',dataOut.params_avg[4],dataOut.heightList)
4811 return dataOut
6082 return dataOut
4812
6083
4813 #--------------- Non Specular Meteor ----------------
6084 #--------------- Non Specular Meteor ----------------
@@ -5918,6 +7189,7 class SMOperations():
5918
7189
5919 ph0_aux = ph0 + ph1
7190 ph0_aux = ph0 + ph1
5920 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
7191 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
7192
5921 #First Estimation
7193 #First Estimation
5922 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
7194 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
5923
7195
@@ -6052,9 +7324,12 class SMOperations():
6052
7324
6053 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
7325 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
6054
7326
6055 return pairslist, distances
7327 return pairslist, distances
6056
7328
6057 class IGRFModel(Operation):
7329 class IGRFModel(Operation):
7330 '''
7331 Written by R. Flores
7332 '''
6058 """Operation to calculate Geomagnetic parameters.
7333 """Operation to calculate Geomagnetic parameters.
6059
7334
6060 Parameters:
7335 Parameters:
@@ -6077,7 +7352,7 class IGRFModel(Operation):
6077 def run(self,dataOut):
7352 def run(self,dataOut):
6078
7353
6079 try:
7354 try:
6080 from schainpy.model.proc import mkfact_short_2020
7355 from schainpy.model.proc import mkfact_short_2020_2
6081 except:
7356 except:
6082 log.warning('You should install "mkfact_short_2020" module to process IGRF Model')
7357 log.warning('You should install "mkfact_short_2020" module to process IGRF Model')
6083
7358
@@ -6091,8 +7366,9 class IGRFModel(Operation):
6091 dataOut.ut=dataOut.bd_time.tm_hour+dataOut.bd_time.tm_min/60.0+dataOut.bd_time.tm_sec/3600.0
7366 dataOut.ut=dataOut.bd_time.tm_hour+dataOut.bd_time.tm_min/60.0+dataOut.bd_time.tm_sec/3600.0
6092
7367
6093 self.aux=0
7368 self.aux=0
6094
7369 dh = dataOut.heightList[1]-dataOut.heightList[0]
6095 dataOut.h=numpy.arange(0.0,15.0*dataOut.MAXNRANGENDT,15.0,dtype='float32')
7370 #dataOut.h=numpy.arange(0.0,15.0*dataOut.MAXNRANGENDT,15.0,dtype='float32')
7371 dataOut.h=numpy.arange(0.0,dh*dataOut.MAXNRANGENDT,dh,dtype='float32')
6096 dataOut.bfm=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
7372 dataOut.bfm=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
6097 dataOut.bfm=numpy.array(dataOut.bfm,order='F')
7373 dataOut.bfm=numpy.array(dataOut.bfm,order='F')
6098 dataOut.thb=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
7374 dataOut.thb=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
@@ -6100,7 +7376,7 class IGRFModel(Operation):
6100 dataOut.bki=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
7376 dataOut.bki=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
6101 dataOut.bki=numpy.array(dataOut.bki,order='F')
7377 dataOut.bki=numpy.array(dataOut.bki,order='F')
6102
7378
6103 mkfact_short_2020.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
7379 mkfact_short_2020_2.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
6104
7380
6105 return dataOut
7381 return dataOut
6106
7382
@@ -6113,9 +7389,7 class MergeProc(ProcessingUnit):
6113
7389
6114 self.dataOut = getattr(self, self.inputs[0])
7390 self.dataOut = getattr(self, self.inputs[0])
6115 data_inputs = [getattr(self, attr) for attr in self.inputs]
7391 data_inputs = [getattr(self, attr) for attr in self.inputs]
6116 #print(self.inputs)
7392
6117 #print(numpy.shape([getattr(data, attr_data) for data in data_inputs][1]))
6118 #exit(1)
6119 if mode==0:
7393 if mode==0:
6120 data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs])
7394 data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs])
6121 setattr(self.dataOut, attr_data, data)
7395 setattr(self.dataOut, attr_data, data)
@@ -6181,39 +7455,6 class MergeProc(ProcessingUnit):
6181 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
7455 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
6182 setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
7456 setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
6183 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
7457 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
6184 #exit(1)
6185 #print(self.dataOut.data_spc_LP.shape)
6186 #print("Exit")
6187 #exit(1)
6188 #setattr(self.dataOut, 'dataLag_cspc', [getattr(data, attr_data_2) for data in data_inputs][0])
6189 #setattr(self.dataOut, 'dataLag_cspc_LP', [getattr(data, attr_data_2) for data in data_inputs][1])
6190 #setattr(self.dataOut, 'nIncohInt', [getattr(data, attr_data_3) for data in data_inputs][0])
6191 #setattr(self.dataOut, 'nIncohInt_LP', [getattr(data, attr_data_3) for data in data_inputs][1])
6192 '''
6193 print(self.dataOut.dataLag_spc_LP.shape)
6194 print(self.dataOut.dataLag_cspc_LP.shape)
6195 exit(1)
6196 '''
6197 '''
6198 print(self.dataOut.dataLag_spc_LP[0,:,100])
6199 print(self.dataOut.dataLag_spc_LP[1,:,100])
6200 exit(1)
6201 '''
6202 #self.dataOut.dataLag_spc_LP = numpy.transpose(self.dataOut.dataLag_spc_LP[0],(2,0,1))
6203 #self.dataOut.dataLag_cspc_LP = numpy.transpose(self.dataOut.dataLag_cspc_LP,(3,1,2,0))
6204 '''
6205 print("Merge")
6206 print(numpy.shape(self.dataOut.dataLag_spc))
6207 print(numpy.shape(self.dataOut.dataLag_spc_LP))
6208 print(numpy.shape(self.dataOut.dataLag_cspc))
6209 print(numpy.shape(self.dataOut.dataLag_cspc_LP))
6210 exit(1)
6211 '''
6212 #print(numpy.sum(self.dataOut.dataLag_spc_LP[2,:,164])/128)
6213 #print(numpy.sum(self.dataOut.dataLag_cspc_LP[0,:,30,1])/128)
6214 #exit(1)
6215 #print(self.dataOut.NDP)
6216 #print(self.dataOut.nNoiseProfiles)
6217
7458
6218 #self.dataOut.nIncohInt_LP = 128
7459 #self.dataOut.nIncohInt_LP = 128
6219 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
7460 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
@@ -6224,6 +7465,8 class MergeProc(ProcessingUnit):
6224 #print("sahpi",self.dataOut.nIncohInt_LP)
7465 #print("sahpi",self.dataOut.nIncohInt_LP)
6225 #exit(1)
7466 #exit(1)
6226 self.dataOut.NLAG = 16
7467 self.dataOut.NLAG = 16
7468 self.dataOut.NLAG = self.dataOut.data_acf.shape[1]
7469
6227 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
7470 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
6228
7471
6229 #print(numpy.shape(self.dataOut.data_spc))
7472 #print(numpy.shape(self.dataOut.data_spc))
@@ -6234,6 +7477,97 class MergeProc(ProcessingUnit):
6234 setattr(self.dataOut, attr_data, data)
7477 setattr(self.dataOut, attr_data, data)
6235 data = numpy.concatenate([getattr(data, attr_data_2) for data in data_inputs])
7478 data = numpy.concatenate([getattr(data, attr_data_2) for data in data_inputs])
6236 setattr(self.dataOut, attr_data_2, data)
7479 setattr(self.dataOut, attr_data_2, data)
6237 #data = numpy.concatenate([getattr(data, attr_data_3) for data in data_inputs])
7480
6238 #setattr(self.dataOut, attr_data_3, data)
7481 if mode==6: #Hybrid Spectra-Voltage
6239 #print(self.dataOut.moments.shape,self.dataOut.data_snr.shape,self.dataOut.heightList.shape) No newline at end of file
7482 #data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs],axis=1)
7483 #setattr(self.dataOut, attr_data, data)
7484 setattr(self.dataOut, 'dataLag_spc', getattr(data_inputs[1], attr_data)) #DP
7485 setattr(self.dataOut, 'dataLag_cspc', getattr(data_inputs[1], attr_data_2)) #DP
7486 setattr(self.dataOut, 'output_LP_integrated', getattr(data_inputs[0], attr_data_3)) #LP
7487 #setattr(self.dataOut, 'dataLag_cspc_LP', getattr(data_inputs[1], attr_data_4)) #LP
7488 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
7489 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
7490 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
7491 #print(self.dataOut.NSCAN)
7492 self.dataOut.nIncohInt = int(self.dataOut.NAVG * self.dataOut.nint)
7493 #print(self.dataOut.dataLag_spc.shape)
7494 self.dataOut.nProfiles = self.dataOut.nProfiles_DP = self.dataOut.dataLag_spc.shape[1]
7495 '''
7496 #self.dataOut.nIncohInt_LP = 128
7497 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
7498 self.dataOut.nProfiles_LP = 16#28#self.dataOut.nIncohInt_LP
7499 self.dataOut.nProfiles_LP = self.dataOut.data_acf.shape[1]#28#self.dataOut.nIncohInt_LP
7500 self.dataOut.NSCAN = 128
7501 self.dataOut.nIncohInt_LP = self.dataOut.nIncohInt*self.dataOut.NSCAN
7502 #print("sahpi",self.dataOut.nIncohInt_LP)
7503 #exit(1)
7504 self.dataOut.NLAG = 16
7505 self.dataOut.NLAG = self.dataOut.data_acf.shape[1]
7506 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
7507 '''
7508 #print(numpy.shape(self.dataOut.data_spc))
7509 #print("*************************GOOD*************************")
7510 #exit(1)
7511
7512 if mode==11: #MST ISR
7513 #data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs],axis=1)
7514 #setattr(self.dataOut, attr_data, data)
7515 #setattr(self.dataOut, 'ph2', [getattr(data, attr_data) for data in data_inputs][1])
7516 #setattr(self.dataOut, 'dphi', [getattr(data, attr_data_2) for data in data_inputs][1])
7517 #setattr(self.dataOut, 'sdp2', [getattr(data, attr_data_3) for data in data_inputs][1])
7518
7519 setattr(self.dataOut, 'ph2', getattr(data_inputs[1], attr_data)) #DP
7520 setattr(self.dataOut, 'dphi', getattr(data_inputs[1], attr_data_2)) #DP
7521 setattr(self.dataOut, 'sdp2', getattr(data_inputs[1], attr_data_3)) #DP
7522
7523 print("MST Density", numpy.shape(self.dataOut.ph2))
7524 print("cf MST: ", self.dataOut.cf)
7525 #exit(1)
7526 #print("MST Density", self.dataOut.ph2[116:283])
7527 print("MST Density", self.dataOut.ph2[80:120])
7528 print("MST dPhi", self.dataOut.dphi[80:120])
7529 self.dataOut.ph2 *= self.dataOut.cf#0.0008136899
7530 #print("MST Density", self.dataOut.ph2[116:283])
7531 self.dataOut.sdp2 *= 0#self.dataOut.cf#0.0008136899
7532 #print("MST Density", self.dataOut.ph2[116:283])
7533 print("MST Density", self.dataOut.ph2[80:120])
7534 self.dataOut.NSHTS = int(numpy.shape(self.dataOut.ph2)[0])
7535 dH = self.dataOut.heightList[1]-self.dataOut.heightList[0]
7536 dH /= self.dataOut.windowOfFilter
7537 self.dataOut.heightList = numpy.arange(0,self.dataOut.NSHTS)*dH + dH
7538 #print("heightList: ", self.dataOut.heightList)
7539 self.dataOut.NDP = self.dataOut.NSHTS
7540 #exit(1)
7541 #print(self.dataOut.heightList)
7542
7543 class MST_Den_Conv(Operation):
7544 '''
7545 Written by R. Flores
7546 '''
7547 """Operation to calculate Geomagnetic parameters.
7548
7549 Parameters:
7550 -----------
7551 None
7552
7553 Example
7554 --------
7555
7556 op = proc_unit.addOperation(name='MST_Den_Conv', optype='other')
7557
7558 """
7559
7560 def __init__(self, **kwargs):
7561
7562 Operation.__init__(self, **kwargs)
7563
7564 def run(self,dataOut):
7565
7566 dataOut.PowDen = numpy.zeros((1,dataOut.NDP))
7567 dataOut.PowDen[0] = numpy.copy(dataOut.ph2[:dataOut.NDP])
7568
7569 dataOut.FarDen = numpy.zeros((1,dataOut.NDP))
7570 dataOut.FarDen[0] = numpy.copy(dataOut.dphi[:dataOut.NDP])
7571 print("pow den shape", numpy.shape(dataOut.PowDen))
7572 print("far den shape", numpy.shape(dataOut.FarDen))
7573 return dataOut
General Comments 0
You need to be logged in to leave comments. Login now