##// END OF EJS Templates
schain
RSS
Update Noise para calulo de Z
avaldezp -
r1474:d5178fc6fb6a
parent child
Show More
Side by Side Unified
Context file:
r1474:d5178fc6fb6a Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,377
1 # SOPHY PROC script
2 import os, sys, json, argparse
3 import datetime
4 import time
5
6 PATH = '/DATA_RM/DATA'
7 # PATH = '/Users/jespinoza/workspace/data/'
8 PATH = '/home/soporte/Documents/HUANCAYO'
9 PARAM = {
10 'P': {'name': 'dataPP_POWER', 'zmin': -40, 'zmax': -5, 'colormap': 'jet', 'label': 'Power', 'wrname': 'Pow','cb_label': 'dB', 'ch':0},
11 'V': {'name': 'dataPP_DOP', 'zmin': -20, 'zmax': 20, 'colormap': 'seismic', 'label': 'Velocity', 'wrname': 'Vel', 'cb_label': 'm/s', 'ch':0},
12 'RH': {'name': 'RhoHV_R', 'zmin': 0, 'zmax': 1, 'colormap': 'jet', 'label': 'Coef.Correlacion', 'wrname':'R', 'cb_label': '*', 'ch':0},
13 'FD': {'name': 'PhiD_P', 'zmin': -180,'zmax': 180,'colormap': 'RdBu_r', 'label': 'Fase Diferencial', 'wrname':'P' , 'cb_label': 'ΒΊ', 'ch':0},
14 'ZD': {'name': 'Zdb_D', 'zmin': -20, 'zmax': 60, 'colormap': 'viridis','label': 'Reflect.Diferencial','wrname':'D' , 'cb_label': 'dBz','ch':0},
15 'Z': {'name': 'Zdb', 'zmin': -20, 'zmax': 70, 'colormap': 'gist_ncar','label': 'Reflectividad', 'wrname':'Z', 'cb_label': 'dBz','ch':1},
16 'W': {'name': 'Sigmav_W', 'zmin': 0, 'zmax':5, 'colormap': 'viridis','label': 'AnchoEspectral', 'wrname':'S', 'cb_label': 'hz', 'ch':1}
17 }
18
19 #
20 def max_index(r, sample_rate, ipp):
21
22 return int(sample_rate*ipp*1e6 * r / 60) + int(sample_rate*ipp*1e6 * 1.2 / 60)
23
24 def main(args):
25
26 experiment = args.experiment
27 fp = open(os.path.join(PATH, experiment, 'experiment.conf'))
28 conf = json.loads(fp.read())
29
30 ipp_km = conf['usrp_tx']['ipp']
31 ipp = ipp_km * 2 /300000
32 sample_rate = conf['usrp_rx']['sample_rate']
33 axis = ['0' if x=='elevation' else '1' for x in conf['pedestal']['axis']] # AZIMUTH 1 ELEVACION 0
34 speed_axis = conf['pedestal']['speed']
35 steeps = conf['pedestal']['table']
36 time_offset = args.time_offset
37 parameters = args.parameters
38 start_date = experiment.split('@')[1].split('T')[0].replace('-', '/')
39 end_date = start_date
40 start_time = experiment.split('@')[1].split('T')[1].replace('-', ':')
41 end_time = '23:59:59'
42 N = int(1/(speed_axis[0]*ipp)) # 1 GRADO DE RESOLUCION
43 path = os.path.join(PATH, experiment, 'rawdata')
44 path_ped = os.path.join(PATH, experiment, 'position')
45 path_plots = os.path.join(PATH, experiment, 'plots_ch0')
46 path_save = os.path.join(PATH, experiment, 'param')
47 RMIX = 10
48
49 from schainpy.controller import Project
50
51 project = Project()
52 project.setup(id='1', name='Sophy', description='sophy proc')
53
54 reader = project.addReadUnit(datatype='DigitalRFReader',
55 path=path,
56 startDate=start_date,
57 endDate=end_date,
58 startTime=start_time,
59 endTime=end_time,
60 delay=30,
61 channelList='0',
62 online=args.online,
63 walk=1,
64 ippKm = ipp_km,
65 getByBlock = 1,
66 nProfileBlocks = N,
67 )
68
69 if not conf['usrp_tx']['enable_2']: # One Pulse
70 voltage = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
71
72 if conf['usrp_tx']['code_type_1']:
73 code = [c.split() for c in conf['usrp']['code_1']]
74 op = voltage.addOperation(name='Decoder', optype='other')
75 op.addParameter(name='code', value=code)
76 op.addParameter(name='nCode', value=len(code), format='int')
77 op.addParameter(name='nBaud', value=len(code[0]), format='int')
78
79 op = voltage.addOperation(name='setH0')
80 op.addParameter(name='h0', value='-1.2')
81
82 if args.range >= 0:
83 op = voltage.addOperation(name='selectHeights')
84 op.addParameter(name='minIndex', value='0', format='int')
85 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
86
87 code=[[1]]
88 opObj11 = voltage.addOperation(name='Decoder', optype='other')
89 opObj11.addParameter(name='code', value=code)
90 opObj11.addParameter(name='nCode', value='1', format='int')
91 opObj11.addParameter(name='nBaud', value='1', format='int')
92
93 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
94 op.addParameter(name='n', value=2*len(code), format='int')
95
96 #op = voltage.addOperation(name='PulsePair_vRF', optype='other')
97 #op.addParameter(name='n', value=int(N), format='int')
98
99 if args.range >= 0:
100 op = voltage.addOperation(name='selectHeights')
101 op.addParameter(name='minIndex', value='0', format='int')
102 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
103
104
105 op = voltage.addOperation(name='PulsePair_vRF', optype='other')
106 op.addParameter(name='n', value=125, format='int')
107
108
109 proc = project.addProcUnit(datatype='ParametersProc', inputId=voltage.getId())
110 #procUnitConfObjB.addParameter(name='runNextUnit', value=True)
111
112 opObj10 = proc.addOperation(name="WeatherRadar")
113 opObj10.addParameter(name='variableList',value='Reflectividad,VelocidadRadial,AnchoEspectral')
114
115 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
116
117 op = proc.addOperation(name='PedestalInformation')
118 op.addParameter(name='path', value=path_ped, format='str')
119 op.addParameter(name='interval', value='0.04')
120 op.addParameter(name='time_offset', value=time_offset)
121 op.addParameter(name='az_offset', value=-26.2)
122
123 for param in parameters:
124 op = proc.addOperation(name='Block360_vRF4')
125 #op.addParameter(name='axis', value=','.join(axis))
126 op.addParameter(name='attr_data', value=PARAM[param]['name'])
127 op.addParameter(name='runNextOp', value=True)
128
129 op= proc.addOperation(name='WeatherParamsPlot')
130 if args.save: op.addParameter(name='save', value=path_plots, format='str')
131 op.addParameter(name='save_period', value=-1)
132 op.addParameter(name='show', value=args.show)
133 op.addParameter(name='channels', value='0,')
134 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
135 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
136 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
137 op.addParameter(name='labels', value=[PARAM[param]['label']])
138 op.addParameter(name='save_code', value=param)
139 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
140 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
141
142 desc = {
143 'Data': {
144 PARAM[param]['name']: PARAM[param]['label'],
145 'utctime': 'time'
146 },
147 'Metadata': {
148 'heightList': 'range',
149 'data_azi': 'azimuth',
150 'data_ele': 'elevation',
151 }
152 }
153
154 if args.save:
155 opObj10 = proc.addOperation(name='HDFWriter')
156 opObj10.addParameter(name='path',value=path_save+'-{}'.format(param), format='str')
157 opObj10.addParameter(name='Reset',value=True)
158 opObj10.addParameter(name='setType',value='weather')
159 opObj10.addParameter(name='description',value='desc')
160 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
161 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
162 opObj10.addParameter(name='dataList',value='{},utctime'.format(PARAM[param]['name']))
163
164 else: #Two pulses
165
166 voltage1 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
167
168 print("repetions",conf['usrp_tx']['repetitions_1'])
169
170 op = voltage1.addOperation(name='ProfileSelector')
171 op.addParameter(name='profileRangeList', value='0,{}'.format(conf['usrp_tx']['repetitions_1']-1))
172
173
174 #op3 = voltage1.addOperation(name='ProfileSelector', optype='other')
175 #op3.addParameter(name='profileRangeList', value='1,123')
176
177 '''
178 if conf['usrp_tx']['code_type_1'] != 'None':
179 code = [c.split() for c in conf['usrp_tx']['code_1']]
180 op = voltage1.addOperation(name='Decoder', optype='other')
181 op.addParameter(name='code', value=code)
182 op.addParameter(name='nCode', value=len(code), format='int')
183 op.addParameter(name='nBaud', value=len(code[0]), format='int')
184 '''
185
186 code=[[1]]
187
188 opObj11 = voltage1.addOperation(name='Decoder', optype='other')
189 opObj11.addParameter(name='code', value=code)
190 opObj11.addParameter(name='nCode', value='1', format='int')
191 opObj11.addParameter(name='nBaud', value='1', format='int')
192
193 op = voltage1.addOperation(name='setH0')
194 op.addParameter(name='h0', value='-1.2')
195
196 if args.range >= 0:
197 op = voltage1.addOperation(name='selectHeights')
198 op.addParameter(name='minIndex', value='0', format='int')
199 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
200
201 op = voltage1.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
202 op.addParameter(name='n', value=2, format='int')
203
204 op = voltage1.addOperation(name='PulsePair_vRF', optype='other')
205 #op.addParameter(name='n', value=int(N), format='int')
206 op.addParameter(name='n', value=61, format='int')
207 #op.addParameter(name='removeDC',value=True)
208
209 '''
210 if args.range >= 0:
211 print("corto",max_index(RMIX, sample_rate, ipp))
212 op = voltage1.addOperation(name='selectHeights')
213 op.addParameter(name='minIndex', value='0', format='int')
214 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
215 '''
216 proc1 = project.addProcUnit(datatype='ParametersProc', inputId=voltage1.getId())
217 proc1.addParameter(name='runNextUnit', value=True)
218
219 opObj10 = proc1.addOperation(name="WeatherRadar")
220 opObj10.addParameter(name='variableList',value='Reflectividad,VelocidadRadial,AnchoEspectral')
221 opObj10.addParameter(name='tauW',value=0.4*1e-6)
222 opObj10.addParameter(name='Pt',value=0.2)
223
224 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
225
226 op = proc1.addOperation(name='PedestalInformation')
227 op.addParameter(name='path', value=path_ped, format='str')
228 op.addParameter(name='interval', value='0.04')
229 op.addParameter(name='time_offset', value=time_offset)
230 op.addParameter(name='az_offset', value=-26.2)
231
232 for param in parameters:
233 op = proc1.addOperation(name='Block360_vRF4')
234 op.addParameter(name='attr_data', value=PARAM[param]['name'])
235 op.addParameter(name='runNextOp', value=True)
236
237 voltage2 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
238
239 op = voltage2.addOperation(name='ProfileSelector')
240 op.addParameter(name='profileRangeList', value='{},{}'.format(conf['usrp_tx']['repetitions_1'], conf['usrp_tx']['repetitions_1']+conf['usrp_tx']['repetitions_2']-1))
241
242
243 if conf['usrp_tx']['code_type_2']:
244 print(conf['usrp_tx']['code_2'])
245 codes = [ c.strip() for c in conf['usrp_tx']['code_2'].split(',')]
246 code = []
247 for c in codes:
248 code.append([int(x) for x in c])
249 print(code)
250 print(code[0])
251 op = voltage2.addOperation(name='Decoder', optype='other')
252 op.addParameter(name='code', value=code)
253 op.addParameter(name='nCode', value=len(code), format='int')
254 op.addParameter(name='nBaud', value=len(code[0]), format='int')
255 import numpy
256 pwcode = numpy.sum(numpy.array(code[0])**2)
257 print("pwcode",pwcode)
258
259 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
260 op.addParameter(name='n', value=len(code), format='int')
261 ncode = len(code)
262 else:
263 ncode = 1
264
265 op = voltage2.addOperation(name='setH0')
266 op.addParameter(name='h0', value='-1.2')
267
268 if args.range >= 0:
269 if args.range==0:
270 args.range= ipp_km
271 op = voltage2.addOperation(name='selectHeights')
272 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
273 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
274
275 #op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
276 #op.addParameter(name='n', value=int(N)/ncode, format='int')
277 op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
278 op.addParameter(name='n', value=64, format='int')
279 #op.addParameter(name='removeDC',value=True)
280
281 '''
282
283 if args.range >= 0:
284 if args.range==0:
285 args.range= ipp_km
286 op = voltage2.addOperation(name='selectHeights')
287 print("largo",max_index(RMIX, sample_rate, ipp))
288 print("largo2",max_index(args.range, sample_rate, ipp))
289
290 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
291 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
292 '''
293
294 proc2 = project.addProcUnit(datatype='ParametersProc', inputId=voltage2.getId())
295
296 opObj10 = proc2.addOperation(name="WeatherRadar")
297 opObj10.addParameter(name='variableList',value='Reflectividad,AnchoEspectral')
298 opObj10.addParameter(name='tauW',value=3.2*1e-6)
299 opObj10.addParameter(name='Pt',value=1.6)
300
301
302 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
303
304 op = proc2.addOperation(name='PedestalInformation')
305 op.addParameter(name='path', value=path_ped, format='str')
306 op.addParameter(name='interval', value='0.04')
307 op.addParameter(name='time_offset', value=time_offset)
308 op.addParameter(name='az_offset', value=-26.2)
309
310 for param in parameters:
311 op = proc2.addOperation(name='Block360_vRF4')
312 #op.addParameter(name='axis', value=','.join(axis))
313 op.addParameter(name='attr_data', value=PARAM[param]['name'])
314 op.addParameter(name='runNextOp', value=True)
315
316 merge = project.addProcUnit(datatype='MergeProc', inputId=[proc1.getId(), proc2.getId()])
317 merge.addParameter(name='attr_data', value=PARAM[param]['name'])
318 merge.addParameter(name='mode', value='7') #RM
319
320 op= merge.addOperation(name='WeatherParamsPlot')
321 if args.save: op.addParameter(name='save', value=path_plots, format='str')
322 op.addParameter(name='save_period', value=-1)
323 op.addParameter(name='show', value=args.show)
324 op.addParameter(name='channels', value='0,')
325 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
326 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
327 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
328 op.addParameter(name='labels', value=[PARAM[param]['label']])
329 op.addParameter(name='save_code', value=param)
330 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
331 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
332
333 desc = {
334 'Data': {
335 PARAM[param]['name']: PARAM[param]['label'],
336 'utctime': 'time'
337 },
338 'Metadata': {
339 'heightList': 'range',
340 'data_azi': 'azimuth',
341 'data_ele': 'elevation',
342 }
343 }
344
345 if args.save:
346 opObj10 = merge.addOperation(name='HDFWriter')
347 opObj10.addParameter(name='path',value=path_save, format='str')
348 opObj10.addParameter(name='Reset',value=True)
349 opObj10.addParameter(name='setType',value='weather')
350 opObj10.addParameter(name='description',value='desc')
351 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
352 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
353 opObj10.addParameter(name='dataList',value='{},utctime'.format(PARAM[param]['name']))
354
355 project.start()
356
357 if __name__ == '__main__':
358
359 parser = argparse.ArgumentParser(description='Script to process SOPHy data.')
360 parser.add_argument('experiment',
361 help='Experiment name')
362 parser.add_argument('--parameters', nargs='*', default=['P'],
363 help='Variables to process: P, Z, V')
364 parser.add_argument('--time_offset', default=0,
365 help='Fix time offset')
366 parser.add_argument('--range', default=0, type=float,
367 help='Max range to plot')
368 parser.add_argument('--save', action='store_true',
369 help='Create output files')
370 parser.add_argument('--show', action='store_true',
371 help='Show matplotlib plot.')
372 parser.add_argument('--online', action='store_true',
373 help='Set online mode.')
374
375 args = parser.parse_args()
376
377 main(args)
Show file before | Show file after | Hide comments
@@ -0,0 +1,376
1 # SOPHY PROC script
2 import os, sys, json, argparse
3 import datetime
4 import time
5
6 PATH = '/DATA_RM/DATA'
7 # PATH = '/Users/jespinoza/workspace/data/'
8 PATH = '/home/soporte/Documents/HUANCAYO'
9 PARAM = {
10 'P': {'name': 'dataPP_POWER', 'zmin': -50, 'zmax': -15, 'colormap': 'jet', 'label': 'Power', 'wrname': 'Pow','cb_label': 'dB', 'ch':0},
11 'V': {'name': 'dataPP_DOP', 'zmin': -20, 'zmax': 20, 'colormap': 'seismic', 'label': 'Velocity', 'wrname': 'Vel', 'cb_label': 'm/s', 'ch':0},
12 'RH': {'name': 'RhoHV_R', 'zmin': 0, 'zmax': 1, 'colormap': 'jet', 'label': 'Coef.Correlacion', 'wrname':'R', 'cb_label': '*', 'ch':0},
13 'FD': {'name': 'PhiD_P', 'zmin': -180,'zmax': 180,'colormap': 'RdBu_r', 'label': 'Fase Diferencial', 'wrname':'P' , 'cb_label': 'ΒΊ', 'ch':0},
14 'ZD': {'name': 'Zdb_D', 'zmin': -20, 'zmax': 60, 'colormap': 'viridis','label': 'Reflect.Diferencial','wrname':'D' , 'cb_label': 'dBz','ch':0},
15 'Z': {'name': 'Zdb', 'zmin': -20, 'zmax': 70, 'colormap': 'gist_ncar','label': 'Reflectividad', 'wrname':'Z', 'cb_label': 'dBz','ch':1},
16 'W': {'name': 'Sigmav_W', 'zmin': 0, 'zmax':5, 'colormap': 'viridis','label': 'AnchoEspectral', 'wrname':'S', 'cb_label': 'hz', 'ch':1}
17 }
18
19 def max_index(r, sample_rate, ipp):
20
21 return int(sample_rate*ipp*1e6 * r / 60) + int(sample_rate*ipp*1e6 * 1.2 / 60)
22
23 def main(args):
24
25 experiment = args.experiment
26 fp = open(os.path.join(PATH, experiment, 'experiment.conf'))
27 conf = json.loads(fp.read())
28
29 ipp_km = conf['usrp_tx']['ipp']
30 ipp = ipp_km * 2 /300000
31 sample_rate = conf['usrp_rx']['sample_rate']
32 axis = ['0' if x=='elevation' else '1' for x in conf['pedestal']['axis']] # AZIMUTH 1 ELEVACION 0
33 speed_axis = conf['pedestal']['speed']
34 steeps = conf['pedestal']['table']
35 time_offset = args.time_offset
36 parameters = args.parameters
37 start_date = experiment.split('@')[1].split('T')[0].replace('-', '/')
38 end_date = start_date
39 start_time = experiment.split('@')[1].split('T')[1].replace('-', ':')
40 end_time = '23:59:59'
41 N = int(1/(speed_axis[0]*ipp)) # 1 GRADO DE RESOLUCION
42 path = os.path.join(PATH, experiment, 'rawdata')
43 path_ped = os.path.join(PATH, experiment, 'position')
44 path_plots = os.path.join(PATH, experiment, 'plots_ch1')
45 path_save = os.path.join(PATH, experiment, 'param')
46 RMIX = 2
47
48 from schainpy.controller import Project
49
50 project = Project()
51 project.setup(id='1', name='Sophy', description='sophy proc')
52
53 reader = project.addReadUnit(datatype='DigitalRFReader',
54 path=path,
55 startDate=start_date,
56 endDate=end_date,
57 startTime=start_time,
58 endTime=end_time,
59 delay=30,
60 channelList='0',
61 online=args.online,
62 walk=1,
63 ippKm = ipp_km,
64 getByBlock = 1,
65 nProfileBlocks = N,
66 )
67
68 if not conf['usrp_tx']['enable_2']: # One Pulse
69 voltage = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
70
71 if conf['usrp_tx']['code_type_1']:
72 code = [c.split() for c in conf['usrp']['code_1']]
73 op = voltage.addOperation(name='Decoder', optype='other')
74 op.addParameter(name='code', value=code)
75 op.addParameter(name='nCode', value=len(code), format='int')
76 op.addParameter(name='nBaud', value=len(code[0]), format='int')
77
78 op = voltage.addOperation(name='setH0')
79 op.addParameter(name='h0', value='-1.2')
80
81 if args.range >= 0:
82 op = voltage.addOperation(name='selectHeights')
83 op.addParameter(name='minIndex', value='0', format='int')
84 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
85
86 code=[[1]]
87 opObj11 = voltage.addOperation(name='Decoder', optype='other')
88 opObj11.addParameter(name='code', value=code)
89 opObj11.addParameter(name='nCode', value='1', format='int')
90 opObj11.addParameter(name='nBaud', value='1', format='int')
91
92 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
93 op.addParameter(name='n', value=2*len(code), format='int')
94
95 #op = voltage.addOperation(name='PulsePair_vRF', optype='other')
96 #op.addParameter(name='n', value=int(N), format='int')
97
98 if args.range >= 0:
99 op = voltage.addOperation(name='selectHeights')
100 op.addParameter(name='minIndex', value='0', format='int')
101 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
102
103
104 op = voltage.addOperation(name='PulsePair_vRF', optype='other')
105 op.addParameter(name='n', value=125, format='int')
106
107
108 proc = project.addProcUnit(datatype='ParametersProc', inputId=voltage.getId())
109 #procUnitConfObjB.addParameter(name='runNextUnit', value=True)
110
111 opObj10 = proc.addOperation(name="WeatherRadar")
112 opObj10.addParameter(name='variableList',value='Reflectividad,VelocidadRadial,AnchoEspectral')
113
114 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
115
116 op = proc.addOperation(name='PedestalInformation')
117 op.addParameter(name='path', value=path_ped, format='str')
118 op.addParameter(name='interval', value='0.04')
119 op.addParameter(name='time_offset', value=time_offset)
120 op.addParameter(name='az_offset', value=-26.2)
121
122 for param in parameters:
123 op = proc.addOperation(name='Block360_vRF4')
124 #op.addParameter(name='axis', value=','.join(axis))
125 op.addParameter(name='attr_data', value=PARAM[param]['name'])
126 op.addParameter(name='runNextOp', value=True)
127
128 op= proc.addOperation(name='WeatherParamsPlot')
129 if args.save: op.addParameter(name='save', value=path_plots, format='str')
130 op.addParameter(name='save_period', value=-1)
131 op.addParameter(name='show', value=args.show)
132 op.addParameter(name='channels', value='1,')
133 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
134 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
135 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
136 op.addParameter(name='labels', value=[PARAM[param]['label']])
137 op.addParameter(name='save_code', value=param)
138 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
139 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
140
141 desc = {
142 'Data': {
143 PARAM[param]['name']: PARAM[param]['label'],
144 'utctime': 'time'
145 },
146 'Metadata': {
147 'heightList': 'range',
148 'data_azi': 'azimuth',
149 'data_ele': 'elevation',
150 }
151 }
152
153 if args.save:
154 opObj10 = proc.addOperation(name='HDFWriter')
155 opObj10.addParameter(name='path',value=path_save+'-{}'.format(param), format='str')
156 opObj10.addParameter(name='Reset',value=True)
157 opObj10.addParameter(name='setType',value='weather')
158 opObj10.addParameter(name='description',value='desc')
159 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
160 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
161 opObj10.addParameter(name='dataList',value='{},utctime'.format(PARAM[param]['name']))
162
163 else: #Two pulses
164
165 voltage1 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
166
167 print("repetions",conf['usrp_tx']['repetitions_1'])
168
169 op = voltage1.addOperation(name='ProfileSelector')
170 op.addParameter(name='profileRangeList', value='0,{}'.format(conf['usrp_tx']['repetitions_1']-1))
171
172
173 #op3 = voltage1.addOperation(name='ProfileSelector', optype='other')
174 #op3.addParameter(name='profileRangeList', value='1,123')
175
176 '''
177 if conf['usrp_tx']['code_type_1'] != 'None':
178 code = [c.split() for c in conf['usrp_tx']['code_1']]
179 op = voltage1.addOperation(name='Decoder', optype='other')
180 op.addParameter(name='code', value=code)
181 op.addParameter(name='nCode', value=len(code), format='int')
182 op.addParameter(name='nBaud', value=len(code[0]), format='int')
183 '''
184
185 code=[[1]]
186
187 opObj11 = voltage1.addOperation(name='Decoder', optype='other')
188 opObj11.addParameter(name='code', value=code)
189 opObj11.addParameter(name='nCode', value='1', format='int')
190 opObj11.addParameter(name='nBaud', value='1', format='int')
191
192 op = voltage1.addOperation(name='setH0')
193 op.addParameter(name='h0', value='-1.2')
194
195 if args.range >= 0:
196 op = voltage1.addOperation(name='selectHeights')
197 op.addParameter(name='minIndex', value='0', format='int')
198 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
199
200 op = voltage1.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
201 op.addParameter(name='n', value=2, format='int')
202
203 op = voltage1.addOperation(name='PulsePair_vRF', optype='other')
204 #op.addParameter(name='n', value=int(N), format='int')
205 op.addParameter(name='n', value=61, format='int')
206 #op.addParameter(name='removeDC',value=True)
207
208 '''
209 if args.range >= 0:
210 print("corto",max_index(RMIX, sample_rate, ipp))
211 op = voltage1.addOperation(name='selectHeights')
212 op.addParameter(name='minIndex', value='0', format='int')
213 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
214 '''
215 proc1 = project.addProcUnit(datatype='ParametersProc', inputId=voltage1.getId())
216 proc1.addParameter(name='runNextUnit', value=True)
217
218 opObj10 = proc1.addOperation(name="WeatherRadar")
219 opObj10.addParameter(name='variableList',value='Reflectividad,VelocidadRadial,AnchoEspectral')
220 opObj10.addParameter(name='tauW',value=0.4*1e-6)
221 opObj10.addParameter(name='Pt',value=0.2)
222
223 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
224
225 op = proc1.addOperation(name='PedestalInformation')
226 op.addParameter(name='path', value=path_ped, format='str')
227 op.addParameter(name='interval', value='0.04')
228 op.addParameter(name='time_offset', value=time_offset)
229 op.addParameter(name='az_offset', value=-26.2)
230
231 for param in parameters:
232 op = proc1.addOperation(name='Block360_vRF4')
233 op.addParameter(name='attr_data', value=PARAM[param]['name'])
234 op.addParameter(name='runNextOp', value=True)
235
236 voltage2 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
237
238 op = voltage2.addOperation(name='ProfileSelector')
239 op.addParameter(name='profileRangeList', value='{},{}'.format(conf['usrp_tx']['repetitions_1'], conf['usrp_tx']['repetitions_1']+conf['usrp_tx']['repetitions_2']-1))
240
241
242 if conf['usrp_tx']['code_type_2']:
243 print(conf['usrp_tx']['code_2'])
244 codes = [ c.strip() for c in conf['usrp_tx']['code_2'].split(',')]
245 code = []
246 for c in codes:
247 code.append([int(x) for x in c])
248 print(code)
249 print(code[0])
250 op = voltage2.addOperation(name='Decoder', optype='other')
251 op.addParameter(name='code', value=code)
252 op.addParameter(name='nCode', value=len(code), format='int')
253 op.addParameter(name='nBaud', value=len(code[0]), format='int')
254 import numpy
255 pwcode = numpy.sum(numpy.array(code[0])**2)
256 print("pwcode",pwcode)
257
258 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
259 op.addParameter(name='n', value=len(code), format='int')
260 ncode = len(code)
261 else:
262 ncode = 1
263
264 op = voltage2.addOperation(name='setH0')
265 op.addParameter(name='h0', value='-1.2')
266
267 if args.range >= 0:
268 if args.range==0:
269 args.range= ipp_km
270 op = voltage2.addOperation(name='selectHeights')
271 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
272 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
273
274 #op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
275 #op.addParameter(name='n', value=int(N)/ncode, format='int')
276 op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
277 op.addParameter(name='n', value=64, format='int')
278 #op.addParameter(name='removeDC',value=True)
279
280 '''
281
282 if args.range >= 0:
283 if args.range==0:
284 args.range= ipp_km
285 op = voltage2.addOperation(name='selectHeights')
286 print("largo",max_index(RMIX, sample_rate, ipp))
287 print("largo2",max_index(args.range, sample_rate, ipp))
288
289 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
290 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
291 '''
292
293 proc2 = project.addProcUnit(datatype='ParametersProc', inputId=voltage2.getId())
294
295 opObj10 = proc2.addOperation(name="WeatherRadar")
296 opObj10.addParameter(name='variableList',value='Reflectividad,AnchoEspectral')
297 opObj10.addParameter(name='tauW',value=3.2*1e-6)
298 opObj10.addParameter(name='Pt',value=1.6)
299
300
301 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
302
303 op = proc2.addOperation(name='PedestalInformation')
304 op.addParameter(name='path', value=path_ped, format='str')
305 op.addParameter(name='interval', value='0.04')
306 op.addParameter(name='time_offset', value=time_offset)
307 op.addParameter(name='az_offset', value=-26.2)
308
309 for param in parameters:
310 op = proc2.addOperation(name='Block360_vRF4')
311 #op.addParameter(name='axis', value=','.join(axis))
312 op.addParameter(name='attr_data', value=PARAM[param]['name'])
313 op.addParameter(name='runNextOp', value=True)
314
315 merge = project.addProcUnit(datatype='MergeProc', inputId=[proc1.getId(), proc2.getId()])
316 merge.addParameter(name='attr_data', value=PARAM[param]['name'])
317 merge.addParameter(name='mode', value='7') #RM
318
319 op= merge.addOperation(name='WeatherParamsPlot')
320 if args.save: op.addParameter(name='save', value=path_plots, format='str')
321 op.addParameter(name='save_period', value=-1)
322 op.addParameter(name='show', value=args.show)
323 op.addParameter(name='channels', value='1,')
324 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
325 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
326 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
327 op.addParameter(name='labels', value=[PARAM[param]['label']])
328 op.addParameter(name='save_code', value=param)
329 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
330 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
331
332 desc = {
333 'Data': {
334 PARAM[param]['name']: PARAM[param]['label'],
335 'utctime': 'time'
336 },
337 'Metadata': {
338 'heightList': 'range',
339 'data_azi': 'azimuth',
340 'data_ele': 'elevation',
341 }
342 }
343
344 if args.save:
345 opObj10 = merge.addOperation(name='HDFWriter')
346 opObj10.addParameter(name='path',value=path_save, format='str')
347 opObj10.addParameter(name='Reset',value=True)
348 opObj10.addParameter(name='setType',value='weather')
349 opObj10.addParameter(name='description',value='desc')
350 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
351 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
352 opObj10.addParameter(name='dataList',value='{},utctime'.format(PARAM[param]['name']))
353
354 project.start()
355
356 if __name__ == '__main__':
357
358 parser = argparse.ArgumentParser(description='Script to process SOPHy data.')
359 parser.add_argument('experiment',
360 help='Experiment name')
361 parser.add_argument('--parameters', nargs='*', default=['P'],
362 help='Variables to process: P, Z, V')
363 parser.add_argument('--time_offset', default=0,
364 help='Fix time offset')
365 parser.add_argument('--range', default=0, type=float,
366 help='Max range to plot')
367 parser.add_argument('--save', action='store_true',
368 help='Create output files')
369 parser.add_argument('--show', action='store_true',
370 help='Show matplotlib plot.')
371 parser.add_argument('--online', action='store_true',
372 help='Set online mode.')
373
374 args = parser.parse_args()
375
376 main(args)
Show file before | Show file after | Hide comments
@@ -3920,6 +3920,10 class WeatherRadar(Operation):
3920 Input:
3920 Input:
3921 Output:
3921 Output:
3922 Parameters affected:
3922 Parameters affected:
3923
3924 Conversion Watt
3925 Referencia
3926 https://www.tek.com/en/blog/calculating-rf-power-iq-samples
3923 '''
3927 '''
3924 isConfig = False
3928 isConfig = False
3925 variableList = None
3929 variableList = None
@@ -3950,6 +3954,7 class WeatherRadar(Operation):
3950 Numerator = ((4*numpy.pi)**3 * aL**2 * 16 *numpy.log(2)*(10**18))
3954 Numerator = ((4*numpy.pi)**3 * aL**2 * 16 *numpy.log(2)*(10**18))
3951 Denominator = (Pt *(10**(Gt/10.0))*(10**(Gr/10.0))*(10**(Glna/10.0))* lambda_**2 * SPEED_OF_LIGHT * tauW * numpy.pi*thetaT*thetaR)
3955 Denominator = (Pt *(10**(Gt/10.0))*(10**(Gr/10.0))*(10**(Glna/10.0))* lambda_**2 * SPEED_OF_LIGHT * tauW * numpy.pi*thetaT*thetaR)
3952 self.RadarConstant = Numerator/Denominator
3956 self.RadarConstant = Numerator/Denominator
3957 self.variableList = variableList
3953 if self.variableList== None:
3958 if self.variableList== None:
3954 self.variableList= ['Reflectividad','ReflectividadDiferencial','CoeficienteCorrelacion','FaseDiferencial','VelocidadRadial','AnchoEspectral']
3959 self.variableList= ['Reflectividad','ReflectividadDiferencial','CoeficienteCorrelacion','FaseDiferencial','VelocidadRadial','AnchoEspectral']
3955
3960
@@ -3961,7 +3966,7 class WeatherRadar(Operation):
3961 data_param = numpy.zeros((nCh,4,nHeis))
3966 data_param = numpy.zeros((nCh,4,nHeis))
3962 if type == "Voltage":
3967 if type == "Voltage":
3963 factor = 1
3968 factor = 1
3964 data_param[:,0,:] = dataOut.dataPP_POW/(factor)
3969 data_param[:,0,:] = dataOut.dataPP_POWER/(factor)
3965 data_param[:,1,:] = dataOut.dataPP_DOP
3970 data_param[:,1,:] = dataOut.dataPP_DOP
3966 data_param[:,2,:] = dataOut.dataPP_WIDTH
3971 data_param[:,2,:] = dataOut.dataPP_WIDTH
3967 data_param[:,3,:] = dataOut.dataPP_SNR
3972 data_param[:,3,:] = dataOut.dataPP_SNR
@@ -4010,18 +4015,36 class WeatherRadar(Operation):
4010 '''-----------------------------Potencia de Radar -Signal S-----------------------------'''
4015 '''-----------------------------Potencia de Radar -Signal S-----------------------------'''
4011
4016
4012 Pr = self.setMoments(dataOut,0)
4017 Pr = self.setMoments(dataOut,0)
4013
4018 '''---------------------------- Calculo de Noise y threshold para Reflectividad---------'''
4019 noise = numpy.zeros(self.nCh)
4020 for i in range(self.nCh):
4021 noise[i] = hildebrand_sekhon(Pr[i,:], 1)
4022 window = numpy.where(Pr[i,:]<1.3*noise[i])
4023 Pr[i,window]= 1e-10
4024 Pr = Pr/1000.0 # Conversion Watt
4014 '''-----------2 Reflectividad del Radar y Factor de Reflectividad------'''
4025 '''-----------2 Reflectividad del Radar y Factor de Reflectividad------'''
4015 self.n_radar = numpy.zeros((self.nCh,self.nHeis))
4026 self.n_radar = numpy.zeros((self.nCh,self.nHeis))
4016 self.Z_radar = numpy.zeros((self.nCh,self.nHeis))
4027 self.Z_radar = numpy.zeros((self.nCh,self.nHeis))
4028
4017 for R in range(self.nHeis):
4029 for R in range(self.nHeis):
4018 self.n_radar[:,R] = self.RadarConstant*Pr[:,R]* (self.Range[:,R])**2*(10**-10.246)
4030 self.n_radar[:,R] = self.RadarConstant*Pr[:,R]* (self.Range[:,R]*(10**3))**2
4019
4031
4020 self.Z_radar[:,R] = self.n_radar[:,R]* self.lambda_**4/( numpy.pi**5 * self.Km**2)
4032 self.Z_radar[:,R] = self.n_radar[:,R]* self.lambda_**4/( numpy.pi**5 * self.Km**2)
4021
4033
4022 '''----------- Factor de Reflectividad Equivalente lamda_ < 10 cm , lamda_= 3.2cm-------'''
4034 '''----------- Factor de Reflectividad Equivalente lamda_ < 10 cm , lamda_= 3.2cm-------'''
4023 Zeh = self.Z_radar
4035 Zeh = self.Z_radar
4024 dBZeh = 10*numpy.log10(Zeh)
4036 #print("---------------------------------------------------------------------")
4037 #print("RangedBz",10*numpy.log10((self.Range[0,-10:]*(10**3))**2))
4038 #print("CTE",10*numpy.log10(self.RadarConstant))
4039 #print("Pr first10",10*numpy.log10(Pr[0,:20]))
4040 #print("Pr last10",10*numpy.log10(Pr[0,-20:]))
4041 #print("LCTE",10*numpy.log10(self.lambda_**4/( numpy.pi**5 * self.Km**2)))
4042 if self.Pt<0.3:
4043 factor=-20.0
4044 else:
4045 factor=0
4046
4047 dBZeh = 10*numpy.log10(Zeh) + factor
4025 if type=='N':
4048 if type=='N':
4026 return dBZeh
4049 return dBZeh
4027 elif type=='D':
4050 elif type=='D':
General Comments 0
You need to be logged in to leave comments. Login now