The requested changes are too big and content was truncated. Show full diff
| @@ -0,0 +1,386 | |||
|
|
1 | #!python | |
|
|
2 | ''' | |
|
|
3 | ''' | |
|
|
4 | ||
|
|
5 | import os, sys, json, argparse | |
|
|
6 | import datetime | |
|
|
7 | import time | |
|
|
8 | ||
|
|
9 | from schainpy.controller import Project | |
|
|
10 | ||
|
|
11 | PATH = "/home/pc-igp-179/Documentos/SOPHy" | |
|
|
12 | ||
|
|
13 | ||
|
|
14 | PARAM = { | |
|
|
15 | 'S': {'zmin': -80, 'zmax':-45, 'colormap': 'jet' , 'label': 'Power', 'wrname': 'power','cb_label': 'dBm', 'ch':0}, | |
|
|
16 | 'SNR':{'zmin': -10, 'zmax': 15, 'colormap': 'jet' , 'label': 'SNR', 'wrname': 'snr','cb_label': 'dB', 'ch':0}, | |
|
|
17 | 'V': {'zmin': -12, 'zmax': 12, 'colormap': 'sophy_v', 'label': 'Velocity', 'wrname': 'velocity', 'cb_label': 'm/s', 'ch':0}, | |
|
|
18 | 'R': {'zmin': 0.5, 'zmax': 1 , 'colormap': 'sophy_r', 'label': 'RhoHV', 'wrname':'rhoHV', 'cb_label': '', 'ch':0}, | |
|
|
19 | 'P': {'zmin': -180,'zmax': 180,'colormap': 'sophy_p', 'label': 'PhiDP', 'wrname':'phiDP' , 'cb_label': 'degrees', 'ch':0}, | |
|
|
20 | 'D': {'zmin': -9 , 'zmax': 12, 'colormap': 'sophy_d', 'label': 'ZDR','wrname':'differential_reflectivity' , 'cb_label': 'dB','ch':0}, | |
|
|
21 | 'Z': {'zmin': -20, 'zmax': 80, 'colormap': 'sophy_z', 'label': 'Reflectivity ', 'wrname':'reflectivity', 'cb_label': 'dBz','ch':0}, | |
|
|
22 | 'W': {'zmin': 0 , 'zmax': 12, 'colormap': 'sophy_w', 'label': 'Spectral Width', 'wrname':'spectral_width', 'cb_label': 'm/s', 'ch':0} | |
|
|
23 | } | |
|
|
24 | ||
|
|
25 | META = ['heightList', 'data_azi', 'data_ele', 'mode_op', 'latitude', 'longitude', 'altitude', 'heading', 'radar_name', | |
|
|
26 | 'institution', 'contact', 'h0', 'range_unit', 'prf', 'prf_unit', 'variable', 'variable_unit', 'n_pulses', | |
|
|
27 | 'pulse1_range', 'pulse1_width', 'pulse2_width', 'pulse1_repetitions', 'pulse2_repetitions', 'pulse_width_unit', | |
|
|
28 | 'snr_threshold', 'data_noise'] | |
|
|
29 | ||
|
|
30 | ||
|
|
31 | def max_index(r, sample_rate, ipp, h0,ipp_km): | |
|
|
32 | ||
|
|
33 | return int(sample_rate*ipp*1e6 * r / ipp_km) + int(sample_rate*ipp*1e6 * -h0 / ipp_km) | |
|
|
34 | ||
|
|
35 | ||
|
|
36 | def main(args): | |
|
|
37 | experiment = args.experiment | |
|
|
38 | ||
|
|
39 | fp = open(os.path.join(PATH, experiment, 'experiment.json')) | |
|
|
40 | conf = json.loads(fp.read()) | |
|
|
41 | ||
|
|
42 | ipp_km = conf['usrp_tx']['ipp'] | |
|
|
43 | bottom = conf['pedestal']['bottom'] | |
|
|
44 | ipp = ipp_km * 2 /300000 | |
|
|
45 | sample_rate = conf['usrp_rx']['sample_rate'] | |
|
|
46 | speed_axis = conf['pedestal']['speed'] | |
|
|
47 | ||
|
|
48 | if args.angles: | |
|
|
49 | angles = args.angles | |
|
|
50 | else: | |
|
|
51 | angles = conf['pedestal']['table'] | |
|
|
52 | ||
|
|
53 | ||
|
|
54 | ||
|
|
55 | start_date = conf['name'].split('@')[1].split('T')[0].replace('-', '/') | |
|
|
56 | end_date = start_date | |
|
|
57 | if args.start_time: | |
|
|
58 | start_time = args.start_time | |
|
|
59 | else: | |
|
|
60 | start_time = conf['name'].split('@')[1].split('T')[1].replace('-', ':') | |
|
|
61 | ||
|
|
62 | if args.end_time: | |
|
|
63 | end_time = args.end_time | |
|
|
64 | else: | |
|
|
65 | end_time = '23:59:59' | |
|
|
66 | ||
|
|
67 | ||
|
|
68 | if args.label: | |
|
|
69 | label = '-{}'.format(args.label) | |
|
|
70 | else: | |
|
|
71 | label = '' | |
|
|
72 | ||
|
|
73 | path_plots = os.path.join(PATH, experiment, 'plots{}'.format(label)) | |
|
|
74 | path_save = os.path.join(PATH, experiment, 'param{}'.format(label)) | |
|
|
75 | ||
|
|
76 | ||
|
|
77 | time_offset = args.time_offset | |
|
|
78 | parameters = args.parameters | |
|
|
79 | ||
|
|
80 | N = int(1.0/(abs(speed_axis[0])*ipp)) | |
|
|
81 | RMIX = 6.0 | |
|
|
82 | H0 = -1.33 | |
|
|
83 | MASK = args.mask | |
|
|
84 | ||
|
|
85 | path = os.path.join(PATH, experiment, 'rawdata') | |
|
|
86 | path_ped = os.path.join(PATH, experiment, 'position') | |
|
|
87 | ||
|
|
88 | proyecto = Project() | |
|
|
89 | proyecto.setup(id = '2', name='Test_2025', description="PRUEBA") | |
|
|
90 | ||
|
|
91 | lectura = proyecto.addReadUnit(datatype='DigitalRFReader', | |
|
|
92 | path=path, | |
|
|
93 | startDate= "2025/01/06",#today, | |
|
|
94 | endDate= "2025/01/06",#today, | |
|
|
95 | startTime= start_time, | |
|
|
96 | endTime= end_time, | |
|
|
97 | delay=30, | |
|
|
98 | #set=0, | |
|
|
99 | online=0, | |
|
|
100 | walk=1, | |
|
|
101 | ippKm = ipp_km, | |
|
|
102 | getByBlock = 1, | |
|
|
103 | nProfileBlocks = N | |
|
|
104 | ) | |
|
|
105 | ||
|
|
106 | n_pulses = 1 | |
|
|
107 | pulse_1_width = conf['usrp_tx']['pulse_1'] | |
|
|
108 | pulse_1_repetitions = conf['usrp_tx']['repetitions_1'] | |
|
|
109 | pulse_2_width = conf['usrp_tx']['pulse_2'] | |
|
|
110 | pulse_2_repetitions = conf['usrp_tx']['repetitions_2'] | |
|
|
111 | ||
|
|
112 | if '1' in args.pulses: | |
|
|
113 | voltage1 = proyecto.addProcUnit( | |
|
|
114 | datatype='VoltageProc', | |
|
|
115 | inputId=lectura.getId() | |
|
|
116 | ) | |
|
|
117 | ||
|
|
118 | op = voltage1.addOperation( | |
|
|
119 | name='ProfileSelector' | |
|
|
120 | ) | |
|
|
121 | op.addParameter( | |
|
|
122 | name='profileRangeList', | |
|
|
123 | value='0,{}'.format(conf['usrp_tx']['repetitions_1']-1) | |
|
|
124 | ) | |
|
|
125 | ||
|
|
126 | if conf['usrp_tx']['code_type_1'] != 'None': | |
|
|
127 | codes = [ c.strip() for c in conf['usrp_tx']['code_1'].split(',')] | |
|
|
128 | code = [] | |
|
|
129 | for c in codes: | |
|
|
130 | code.append([int(x) for x in c]) | |
|
|
131 | op = voltage1.addOperation(name='Decoder', optype='other') | |
|
|
132 | op.addParameter(name='code', value=code) | |
|
|
133 | op.addParameter(name='nCode', value=len(code), format='int') | |
|
|
134 | op.addParameter(name='nBaud', value=len(code[0]), format='int') | |
|
|
135 | ncode = len(code) | |
|
|
136 | else: | |
|
|
137 | ncode = 1 | |
|
|
138 | code = ['0'] | |
|
|
139 | ||
|
|
140 | op = voltage1.addOperation(name='CohInt', optype='other') | |
|
|
141 | op.addParameter(name='n', value=ncode, format='int') | |
|
|
142 | ||
|
|
143 | op = voltage1.addOperation(name='setH0') | |
|
|
144 | op.addParameter(name='h0', value=H0, format='float') | |
|
|
145 | ||
|
|
146 | if args.range > 0: | |
|
|
147 | op = voltage1.addOperation(name='selectHeights') | |
|
|
148 | op.addParameter(name='minIndex', value=max_index(0, sample_rate, ipp, H0,ipp_km), format='int') | |
|
|
149 | op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp, H0,ipp_km), format='int') | |
|
|
150 | ||
|
|
151 | ||
|
|
152 | op = voltage1.addOperation(name='PulsePair_vRF', optype='other') | |
|
|
153 | op.addParameter(name='n', value=int(conf['usrp_tx']['repetitions_1'])/ncode, format='int') | |
|
|
154 | ||
|
|
155 | if args.rmDC: | |
|
|
156 | op.addParameter(name='removeDC', value=1, format='int') | |
|
|
157 | ||
|
|
158 | proc1 = proyecto.addProcUnit(datatype='ParametersProc', inputId=voltage1.getId()) | |
|
|
159 | proc1.addParameter(name='runNextUnit', value=True) | |
|
|
160 | ||
|
|
161 | opObj10 = proc1.addOperation(name="WeatherRadar") | |
|
|
162 | opObj10.addParameter(name='CR_Flag',value=True) | |
|
|
163 | print(1, len(code[0])) | |
|
|
164 | opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0])) | |
|
|
165 | opObj10.addParameter(name='Pt',value=200) | |
|
|
166 | opObj10.addParameter(name='min_index',value=max_index(0, sample_rate, ipp, H0,ipp_km)) | |
|
|
167 | ||
|
|
168 | ||
|
|
169 | op = proc1.addOperation(name='PedestalInformation') | |
|
|
170 | op.addParameter(name='path', value=path_ped, format='str') | |
|
|
171 | op.addParameter(name='interval', value='0.04') | |
|
|
172 | op.addParameter(name='time_offset', value=time_offset) | |
|
|
173 | op.addParameter(name='mode', value=args.mode) | |
|
|
174 | ||
|
|
175 | op = proc1.addOperation(name='Block360') | |
|
|
176 | op.addParameter(name='attr_data', value='data_param') | |
|
|
177 | op.addParameter(name='runNextOp', value=True) | |
|
|
178 | op.addParameter(name='angles', value=angles) | |
|
|
179 | op.addParameter(name='heading', value=conf['heading']) | |
|
|
180 | ||
|
|
181 | ||
|
|
182 | if '2' in args.pulses: | |
|
|
183 | voltage2 = proyecto.addProcUnit( | |
|
|
184 | datatype='VoltageProc', | |
|
|
185 | inputId=lectura.getId() | |
|
|
186 | ) | |
|
|
187 | ||
|
|
188 | op = voltage2.addOperation( | |
|
|
189 | name='ProfileSelector' | |
|
|
190 | ) | |
|
|
191 | op.addParameter( | |
|
|
192 | name='profileRangeList', | |
|
|
193 | value='{},{}'.format(conf['usrp_tx']['repetitions_1'], conf['usrp_tx']['repetitions_1']+conf['usrp_tx']['repetitions_2']-1) | |
|
|
194 | ) | |
|
|
195 | ||
|
|
196 | if conf['usrp_tx']['code_type_2']: | |
|
|
197 | codes = [ c.strip() for c in conf['usrp_tx']['code_2'].split(',')] | |
|
|
198 | code = [] | |
|
|
199 | for c in codes: | |
|
|
200 | code.append([int(x) for x in c]) | |
|
|
201 | op = voltage2.addOperation(name='Decoder', optype='other') | |
|
|
202 | op.addParameter(name='code', value=code) | |
|
|
203 | op.addParameter(name='nCode', value=len(code), format='int') | |
|
|
204 | op.addParameter(name='nBaud', value=len(code[0]), format='int') | |
|
|
205 | ||
|
|
206 | op = voltage2.addOperation(name='CohInt', optype='other') | |
|
|
207 | op.addParameter(name='n', value=len(code), format='int') | |
|
|
208 | ||
|
|
209 | ncode = len(code) | |
|
|
210 | else: | |
|
|
211 | ncode = 1 | |
|
|
212 | ||
|
|
213 | op = voltage2.addOperation(name='setH0') | |
|
|
214 | op.addParameter(name='h0', value=H0, format='float') | |
|
|
215 | ||
|
|
216 | if args.range > 0: | |
|
|
217 | op = voltage2.addOperation(name='selectHeights') | |
|
|
218 | op.addParameter(name='minIndex', value=max_index(0, sample_rate, ipp, H0,ipp_km), format='int') | |
|
|
219 | op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp, H0,ipp_km), format='int') | |
|
|
220 | ||
|
|
221 | ||
|
|
222 | op = voltage2.addOperation(name='PulsePair_vRF', optype='other') | |
|
|
223 | op.addParameter(name='n', value=int(conf['usrp_tx']['repetitions_1'])/ncode, format='int') | |
|
|
224 | ||
|
|
225 | proc2 = proyecto.addProcUnit(datatype='ParametersProc', inputId=voltage2.getId()) | |
|
|
226 | proc2.addParameter(name='runNextUnit', value=True) | |
|
|
227 | ||
|
|
228 | opObj10 = proc2.addOperation(name="WeatherRadar") | |
|
|
229 | opObj10.addParameter(name='CR_Flag',value=True,format='bool') | |
|
|
230 | print(2, len(code[0])) | |
|
|
231 | opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0])) | |
|
|
232 | opObj10.addParameter(name='Pt',value=200) | |
|
|
233 | opObj10.addParameter(name='min_index',value=max_index(RMIX, sample_rate, ipp, H0,ipp_km)) | |
|
|
234 | ||
|
|
235 | op = proc2.addOperation(name='PedestalInformation') | |
|
|
236 | op.addParameter(name='path', value=path_ped, format='str') | |
|
|
237 | op.addParameter(name='interval', value='0.04') | |
|
|
238 | op.addParameter(name='time_offset', value=time_offset) | |
|
|
239 | op.addParameter(name='mode', value=args.mode) | |
|
|
240 | op.addParameter(name='heading', value=conf['heading']) | |
|
|
241 | ||
|
|
242 | op = proc2.addOperation(name='Block360') | |
|
|
243 | op.addParameter(name='attr_data', value='data_param') | |
|
|
244 | op.addParameter(name='runNextOp', value=True) | |
|
|
245 | op.addParameter(name='angles', value=angles) | |
|
|
246 | op.addParameter(name='heading', value=conf['heading']) | |
|
|
247 | ||
|
|
248 | if '1' in args.pulses and '2' in args.pulses: | |
|
|
249 | merge = proyecto.addProcUnit(datatype='MergeProc', inputId=[proc1.getId(), proc2.getId()]) | |
|
|
250 | merge.addParameter(name='attr_data', value='data_param') | |
|
|
251 | merge.addParameter(name='mode', value='7') | |
|
|
252 | merge.addParameter(name='index', value=max_index(RMIX, sample_rate, ipp, H0,ipp_km)) | |
|
|
253 | ||
|
|
254 | elif '1' in args.pulses: | |
|
|
255 | merge = proc1 | |
|
|
256 | elif '2' in args.pulses: | |
|
|
257 | merge = proc2 | |
|
|
258 | ||
|
|
259 | for param in parameters: | |
|
|
260 | ||
|
|
261 | if args.plot: | |
|
|
262 | op= merge.addOperation(name='WeatherParamsPlot') | |
|
|
263 | if args.save: | |
|
|
264 | op.addParameter(name='save', value=path_plots, format='str') | |
|
|
265 | op.addParameter(name='save_period', value=-1) | |
|
|
266 | op.addParameter(name='show', value=args.show) | |
|
|
267 | op.addParameter(name='channels', value='0,') | |
|
|
268 | op.addParameter(name='zmin', value=PARAM[param]['zmin'], format='int') | |
|
|
269 | op.addParameter(name='zmax', value=PARAM[param]['zmax'], format='int') | |
|
|
270 | op.addParameter(name='yrange', value=20, format='int') | |
|
|
271 | op.addParameter(name='xrange', value=args.range, format='int') | |
|
|
272 | op.addParameter(name='attr_data', value=param, format='str') | |
|
|
273 | op.addParameter(name='labels', value=[[PARAM[param]['label']], [PARAM[param]['label']]]) | |
|
|
274 | op.addParameter(name='save_code', value=param) | |
|
|
275 | op.addParameter(name='cb_label', value=PARAM[param]['cb_label']) | |
|
|
276 | op.addParameter(name='colormap', value=PARAM[param]['colormap']) | |
|
|
277 | op.addParameter(name='bgcolor', value='black') | |
|
|
278 | op.addParameter(name='localtime', value=False) | |
|
|
279 | op.addParameter(name='shapes', value='./shapes') | |
|
|
280 | op.addParameter(name='latitude', value=conf['latitude'], format='float') | |
|
|
281 | op.addParameter(name='longitude', value=conf['longitude'], format='float') | |
|
|
282 | op.addParameter(name='map', value=True) | |
|
|
283 | ||
|
|
284 | if MASK: op.addParameter(name='mask', value=MASK, format='float') | |
|
|
285 | if args.server: | |
|
|
286 | op.addParameter(name='server', value='190.187.237.239:4444') | |
|
|
287 | op.addParameter(name='exp_code', value='400') | |
|
|
288 | ||
|
|
289 | desc = { | |
|
|
290 | 'Data': { | |
|
|
291 | 'data_param': {PARAM[param]['wrname']: ['H', 'V']}, | |
|
|
292 | 'utctime': 'time' | |
|
|
293 | }, | |
|
|
294 | 'Metadata': { | |
|
|
295 | 'heightList': 'range', | |
|
|
296 | 'data_azi': 'azimuth', | |
|
|
297 | 'data_ele': 'elevation', | |
|
|
298 | 'mode_op': 'scan_type', | |
|
|
299 | 'h0': 'range_correction', | |
|
|
300 | 'dataPP_NOISE': 'noise', | |
|
|
301 | } | |
|
|
302 | } | |
|
|
303 | ||
|
|
304 | if args.save: | |
|
|
305 | writer = merge.addOperation(name='HDFWriter') | |
|
|
306 | writer.addParameter(name='path', value=path_save, format='str') | |
|
|
307 | writer.addParameter(name='Reset', value=True) | |
|
|
308 | writer.addParameter(name='setType', value='weather') | |
|
|
309 | writer.addParameter(name='setChannel', value='0') #new parameter choose ch 0 H or ch 1 V | |
|
|
310 | writer.addParameter(name='description', value=json.dumps(desc)) | |
|
|
311 | writer.addParameter(name='blocksPerFile', value='1',format='int') | |
|
|
312 | writer.addParameter(name='metadataList', value=','.join(META)) | |
|
|
313 | writer.addParameter(name='dataList', value='data_param,utctime') | |
|
|
314 | writer.addParameter(name='weather_var', value=param) | |
|
|
315 | writer.addParameter(name='mask', value=MASK, format='float') | |
|
|
316 | writer.addParameter(name='localtime', value=False) | |
|
|
317 | # meta | |
|
|
318 | writer.addParameter(name='latitude', value=conf['latitude']) | |
|
|
319 | writer.addParameter(name='longitude', value=conf['longitude']) | |
|
|
320 | writer.addParameter(name='altitude', value=conf['altitude']) | |
|
|
321 | writer.addParameter(name='heading', value=conf['heading']) | |
|
|
322 | writer.addParameter(name='radar_name', value='SOPHy') | |
|
|
323 | writer.addParameter(name='institution', value='IGP') | |
|
|
324 | writer.addParameter(name='contact', value='dscipion@igp.gob.pe') | |
|
|
325 | writer.addParameter(name='created_by', value='Signal Chain (https://pypi.org/project/schainpy/)') | |
|
|
326 | writer.addParameter(name='range_unit', value='km') | |
|
|
327 | writer.addParameter(name='prf', value=1/ipp) | |
|
|
328 | writer.addParameter(name='prf_unit', value='hertz') | |
|
|
329 | writer.addParameter(name='variable', value=PARAM[param]['label']) | |
|
|
330 | writer.addParameter(name='variable_unit', value=PARAM[param]['cb_label']) | |
|
|
331 | writer.addParameter(name='n_pulses', value=n_pulses) | |
|
|
332 | writer.addParameter(name='pulse1_range', value=RMIX) | |
|
|
333 | writer.addParameter(name='pulse1_width', value=pulse_1_width) | |
|
|
334 | writer.addParameter(name='pulse2_width', value=pulse_2_width) | |
|
|
335 | writer.addParameter(name='pulse1_repetitions', value=pulse_1_repetitions) | |
|
|
336 | writer.addParameter(name='pulse2_repetitions', value=pulse_2_repetitions) | |
|
|
337 | writer.addParameter(name='pulse_width_unit', value='microseconds') | |
|
|
338 | writer.addParameter(name='snr_threshold', value=MASK) | |
|
|
339 | writer.addParameter(name='cr_hv', value=[67.41,67.17]) #new parameter | |
|
|
340 | ||
|
|
341 | ||
|
|
342 | return proyecto | |
|
|
343 | ||
|
|
344 | ||
|
|
345 | ||
|
|
346 | ||
|
|
347 | if __name__ == '__main__': | |
|
|
348 | ||
|
|
349 | parser = argparse.ArgumentParser(description='Script to process SOPHy data.') | |
|
|
350 | parser.add_argument('experiment', | |
|
|
351 | help='Experiment name') | |
|
|
352 | parser.add_argument('--parameters', nargs='*', default=['S'], | |
|
|
353 | help='Variables to process: P, Z, V') | |
|
|
354 | parser.add_argument('--pulses', nargs='*', default=['1', '2'], | |
|
|
355 | help='Variables to process: 1, 2') | |
|
|
356 | parser.add_argument('--range', default=60, type=float, | |
|
|
357 | help='Max range to plot') | |
|
|
358 | parser.add_argument('--server', action='store_true', | |
|
|
359 | help='Send to realtime') | |
|
|
360 | parser.add_argument('--start_time', default='', | |
|
|
361 | help='Set start time.') | |
|
|
362 | parser.add_argument('--end_time', default='', | |
|
|
363 | help='Set end time.') | |
|
|
364 | parser.add_argument('--rmDC', action='store_true', | |
|
|
365 | help='Apply remove DC.') | |
|
|
366 | parser.add_argument('--time_offset', default=0, | |
|
|
367 | help='Fix time offset') | |
|
|
368 | parser.add_argument('--mode', default=None, | |
|
|
369 | help='Type of scan') | |
|
|
370 | parser.add_argument('--angles', nargs='*', default=[], type=int, | |
|
|
371 | help='Angles to process') | |
|
|
372 | parser.add_argument('--plot', action='store_true', | |
|
|
373 | help='Create plot files') | |
|
|
374 | parser.add_argument('--save', action='store_true', | |
|
|
375 | help='Create output files') | |
|
|
376 | parser.add_argument('--show', action='store_true', | |
|
|
377 | help='Show matplotlib plot.') | |
|
|
378 | parser.add_argument('--mask', default=0.36, type=float, | |
|
|
379 | help='Filter mask over SNR') | |
|
|
380 | parser.add_argument('--label', default='', | |
|
|
381 | help='Label for plot & param folder') | |
|
|
382 | ||
|
|
383 | args = parser.parse_args() | |
|
|
384 | ||
|
|
385 | proyecto= main(args) | |
|
|
386 | proyecto.start() No newline at end of file | |
| @@ -0,0 +1,381 | |||
|
|
1 | import os, sys, json, argparse | |
|
|
2 | import multiprocessing | |
|
|
3 | import datetime | |
|
|
4 | import time | |
|
|
5 | ||
|
|
6 | PATH = "/home/pc-igp-179/Documentos/SOPHy" | |
|
|
7 | ||
|
|
8 | PARAM = { | |
|
|
9 | 'S': {'zmin': -80, 'zmax':-45, 'colormap': 'jet' , 'label': 'Power', 'wrname': 'power','cb_label': 'dBm', 'ch':0}, | |
|
|
10 | 'SNR':{'zmin': -10, 'zmax': 15, 'colormap': 'jet' , 'label': 'SNR', 'wrname': 'snr','cb_label': 'dB', 'ch':0}, | |
|
|
11 | 'V': {'zmin': -12, 'zmax': 12, 'colormap': 'sophy_v', 'label': 'Velocity', 'wrname': 'velocity', 'cb_label': 'm/s', 'ch':0}, | |
|
|
12 | 'R': {'zmin': 0.5, 'zmax': 1 , 'colormap': 'sophy_r', 'label': 'RhoHV', 'wrname':'rhoHV', 'cb_label': '', 'ch':0}, | |
|
|
13 | 'P': {'zmin': -180,'zmax': 180,'colormap': 'sophy_p', 'label': 'PhiDP', 'wrname':'phiDP' , 'cb_label': 'degrees', 'ch':0}, | |
|
|
14 | 'D': {'zmin': -9 , 'zmax': 12, 'colormap': 'sophy_d', 'label': 'ZDR','wrname':'differential_reflectivity' , 'cb_label': 'dB','ch':0}, | |
|
|
15 | 'Z': {'zmin': -20, 'zmax': 80, 'colormap': 'sophy_z', 'label': 'Reflectivity ', 'wrname':'reflectivity', 'cb_label': 'dBz','ch':0}, | |
|
|
16 | 'W': {'zmin': 0 , 'zmax': 12, 'colormap': 'sophy_w', 'label': 'Spectral Width', 'wrname':'spectral_width', 'cb_label': 'm/s', 'ch':0} | |
|
|
17 | } | |
|
|
18 | ||
|
|
19 | META = ['heightList', 'data_azi', 'data_ele', 'mode_op', 'latitude', 'longitude', 'altitude', 'heading', 'radar_name', | |
|
|
20 | 'institution', 'contact', 'h0', 'range_unit', 'prf', 'prf_unit', 'variable', 'variable_unit', 'n_pulses', | |
|
|
21 | 'pulse1_range', 'pulse1_width', 'pulse2_width', 'pulse1_repetitions', 'pulse2_repetitions', 'pulse_width_unit', | |
|
|
22 | 'snr_threshold', 'data_noise'] | |
|
|
23 | ||
|
|
24 | ||
|
|
25 | def max_index(r, sample_rate, ipp, h0,ipp_km): | |
|
|
26 | ||
|
|
27 | return int(sample_rate*ipp*1e6 * r / ipp_km) + int(sample_rate*ipp*1e6 * -h0 / ipp_km) | |
|
|
28 | ||
|
|
29 | def main(args): | |
|
|
30 | ||
|
|
31 | experiment = args.experiment | |
|
|
32 | fp = open(os.path.join(PATH, experiment, 'experiment.json')) | |
|
|
33 | conf = json.loads(fp.read()) | |
|
|
34 | ||
|
|
35 | ipp_km = conf['usrp_tx']['ipp'] | |
|
|
36 | bottom = conf['pedestal']['bottom'] | |
|
|
37 | ipp = ipp_km * 2 /300000 | |
|
|
38 | sample_rate = conf['usrp_rx']['sample_rate'] | |
|
|
39 | speed_axis = conf['pedestal']['speed'] | |
|
|
40 | if args.angles: | |
|
|
41 | angles = args.angles | |
|
|
42 | else: | |
|
|
43 | angles = conf['pedestal']['table'] | |
|
|
44 | time_offset = args.time_offset | |
|
|
45 | parameters = args.parameters | |
|
|
46 | start_date = conf['name'].split('@')[1].split('T')[0].replace('-', '/') | |
|
|
47 | end_date = start_date | |
|
|
48 | if args.start_time: | |
|
|
49 | start_time = args.start_time | |
|
|
50 | else: | |
|
|
51 | start_time = conf['name'].split('@')[1].split('T')[1].replace('-', ':') | |
|
|
52 | ||
|
|
53 | if args.end_time: | |
|
|
54 | end_time = args.end_time | |
|
|
55 | else: | |
|
|
56 | end_time = '23:59:59' | |
|
|
57 | ||
|
|
58 | N = int(1.0/(abs(speed_axis[0])*ipp)) | |
|
|
59 | ||
|
|
60 | path = os.path.join(PATH, experiment, 'rawdata') | |
|
|
61 | ||
|
|
62 | path_ped = os.path.join(PATH, experiment, 'position') | |
|
|
63 | if args.label: | |
|
|
64 | label = '-{}'.format(args.label) | |
|
|
65 | else: | |
|
|
66 | label = '' | |
|
|
67 | path_plots = os.path.join(PATH, experiment, 'plots{}'.format(label)) | |
|
|
68 | path_save = os.path.join(PATH, experiment, 'param{}'.format(label)) | |
|
|
69 | RMIX = 6.0 | |
|
|
70 | H0 = -1.33 | |
|
|
71 | MASK = args.mask | |
|
|
72 | ||
|
|
73 | from schainpy.controller import Project | |
|
|
74 | ||
|
|
75 | project = Project() | |
|
|
76 | project.setup(id='1', name='Sophy', description='sophy proc') | |
|
|
77 | ||
|
|
78 | reader = project.addReadUnit(datatype='DigitalRFReader', | |
|
|
79 | path=path, | |
|
|
80 | startDate=start_date, | |
|
|
81 | endDate=end_date, | |
|
|
82 | startTime=start_time, | |
|
|
83 | endTime=end_time, | |
|
|
84 | delay=0, | |
|
|
85 | online=args.online, | |
|
|
86 | walk=0, | |
|
|
87 | ippKm = ipp_km, | |
|
|
88 | getByBlock = 1, | |
|
|
89 | nProfileBlocks = N, | |
|
|
90 | ) | |
|
|
91 | ||
|
|
92 | ||
|
|
93 | n_pulses = 1 | |
|
|
94 | pulse_1_width = conf['usrp_tx']['pulse_1'] | |
|
|
95 | pulse_1_repetitions = conf['usrp_tx']['repetitions_1'] | |
|
|
96 | pulse_2_width = conf['usrp_tx']['pulse_2'] | |
|
|
97 | pulse_2_repetitions = conf['usrp_tx']['repetitions_2'] | |
|
|
98 | ||
|
|
99 | if '1' in args.pulses: | |
|
|
100 | voltage1 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId()) | |
|
|
101 | ||
|
|
102 | op = voltage1.addOperation(name='ProfileSelector') | |
|
|
103 | op.addParameter(name='profileRangeList', value='0,{}'.format(conf['usrp_tx']['repetitions_1']-1)) | |
|
|
104 | ||
|
|
105 | if conf['usrp_tx']['code_type_1'] != 'None': | |
|
|
106 | codes = [ c.strip() for c in conf['usrp_tx']['code_1'].split(',')] | |
|
|
107 | code = [] | |
|
|
108 | for c in codes: | |
|
|
109 | code.append([int(x) for x in c]) | |
|
|
110 | op = voltage1.addOperation(name='Decoder', optype='other') | |
|
|
111 | op.addParameter(name='code', value=code) | |
|
|
112 | op.addParameter(name='nCode', value=len(code), format='int') | |
|
|
113 | op.addParameter(name='nBaud', value=len(code[0]), format='int') | |
|
|
114 | ncode = len(code) | |
|
|
115 | else: | |
|
|
116 | ncode = 1 | |
|
|
117 | code = ['0'] | |
|
|
118 | ||
|
|
119 | ||
|
|
120 | op = voltage1.addOperation(name='CohInt', optype='other') | |
|
|
121 | op.addParameter(name='n', value=ncode, format='int') | |
|
|
122 | ||
|
|
123 | op = voltage1.addOperation(name='setH0') | |
|
|
124 | op.addParameter(name='h0', value=H0, format='float') | |
|
|
125 | ||
|
|
126 | if args.range > 0: | |
|
|
127 | op = voltage1.addOperation(name='selectHeights') | |
|
|
128 | op.addParameter(name='minIndex', value=max_index(0, sample_rate, ipp, H0,ipp_km), format='int') | |
|
|
129 | op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp, H0,ipp_km), format='int') | |
|
|
130 | ||
|
|
131 | ||
|
|
132 | ||
|
|
133 | #op = voltage1.addOperation(name='PulsePair_vRF', optype='other') | |
|
|
134 | #op.addParameter(name='n', value=int(conf['usrp_tx']['repetitions_1'])/ncode, format='int') | |
|
|
135 | ||
|
|
136 | procUnitConfObjA = project.addProcUnit(datatype='SpectraProc', inputId=voltage1.getId()) | |
|
|
137 | procUnitConfObjA.addParameter(name='nFFTPoints', value=int(conf['usrp_tx']['repetitions_1'])/ncode, format='int') | |
|
|
138 | procUnitConfObjA.addParameter(name='nProfiles', value=int(conf['usrp_tx']['repetitions_1'])/ncode, format='int') | |
|
|
139 | ||
|
|
140 | ||
|
|
141 | #opObj11 = procUnitConfObjA.addOperation(name='setRadarFrequency') | |
|
|
142 | #opObj11.addParameter(name='frequency', value='9.345e9', format='float') | |
|
|
143 | #procUnitConfObjA.addOperation(name='removeDC') | |
|
|
144 | #if args.rmDC: | |
|
|
145 | # op.addParameter(name='removeDC', value=1, format='int') | |
|
|
146 | ||
|
|
147 | proc1 = project.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjA.getId()) | |
|
|
148 | proc1.addOperation(name='SpectralMoments') | |
|
|
149 | ||
|
|
150 | proc1.addParameter(name='runNextUnit', value=True) | |
|
|
151 | ||
|
|
152 | opObj10 = proc1.addOperation(name="WeatherRadar") | |
|
|
153 | opObj10.addParameter(name='CR_Flag',value=True) | |
|
|
154 | #print(1, len(code[0])) | |
|
|
155 | opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0])) | |
|
|
156 | opObj10.addParameter(name='Pt',value=200) | |
|
|
157 | opObj10.addParameter(name='min_index',value=max_index(0, sample_rate, ipp, H0,ipp_km)) | |
|
|
158 | ||
|
|
159 | ||
|
|
160 | op = proc1.addOperation(name='PedestalInformation') | |
|
|
161 | op.addParameter(name='path', value=path_ped, format='str') | |
|
|
162 | op.addParameter(name='interval', value='0.04') | |
|
|
163 | op.addParameter(name='time_offset', value=time_offset) | |
|
|
164 | op.addParameter(name='mode', value=args.mode) | |
|
|
165 | ||
|
|
166 | op = proc1.addOperation(name='Block360') | |
|
|
167 | op.addParameter(name='attr_data', value='data_param') | |
|
|
168 | op.addParameter(name='runNextOp', value=True) | |
|
|
169 | op.addParameter(name='angles', value=angles) | |
|
|
170 | op.addParameter(name='heading', value=conf['heading']) | |
|
|
171 | ||
|
|
172 | ||
|
|
173 | if '2' in args.pulses: | |
|
|
174 | voltage2 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId()) | |
|
|
175 | ||
|
|
176 | op = voltage2.addOperation(name='ProfileSelector') | |
|
|
177 | op.addParameter(name='profileRangeList', value='{},{}'.format(conf['usrp_tx']['repetitions_1'], conf['usrp_tx']['repetitions_1']+conf['usrp_tx']['repetitions_2']-1)) | |
|
|
178 | ||
|
|
179 | if conf['usrp_tx']['code_type_2']: | |
|
|
180 | codes = [ c.strip() for c in conf['usrp_tx']['code_2'].split(',')] | |
|
|
181 | code = [] | |
|
|
182 | for c in codes: | |
|
|
183 | code.append([int(x) for x in c]) | |
|
|
184 | op = voltage2.addOperation(name='Decoder', optype='other') | |
|
|
185 | op.addParameter(name='code', value=code) | |
|
|
186 | op.addParameter(name='nCode', value=len(code), format='int') | |
|
|
187 | op.addParameter(name='nBaud', value=len(code[0]), format='int') | |
|
|
188 | ||
|
|
189 | op = voltage2.addOperation(name='CohInt', optype='other') | |
|
|
190 | op.addParameter(name='n', value=len(code), format='int') | |
|
|
191 | ||
|
|
192 | ncode = len(code) | |
|
|
193 | else: | |
|
|
194 | ncode = 1 | |
|
|
195 | ||
|
|
196 | op = voltage2.addOperation(name='setH0') | |
|
|
197 | op.addParameter(name='h0', value=H0, format='float') | |
|
|
198 | ||
|
|
199 | if args.range > 0: | |
|
|
200 | op = voltage2.addOperation(name='selectHeights') | |
|
|
201 | op.addParameter(name='minIndex', value=max_index(0, sample_rate, ipp, H0,ipp_km), format='int') | |
|
|
202 | op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp, H0,ipp_km), format='int') | |
|
|
203 | ||
|
|
204 | ||
|
|
205 | ||
|
|
206 | ||
|
|
207 | procUnitConfObjB = project.addProcUnit(datatype='SpectraProc', inputId=voltage2.getId()) | |
|
|
208 | procUnitConfObjB.addParameter(name='nFFTPoints', value=int(conf['usrp_tx']['repetitions_2'])/ncode, format='int') | |
|
|
209 | procUnitConfObjB.addParameter(name='nProfiles', value=int(conf['usrp_tx']['repetitions_2'])/ncode, format='int') | |
|
|
210 | ||
|
|
211 | ||
|
|
212 | ||
|
|
213 | #opObj11 = procUnitConfObjB.addOperation(name='setRadarFrequency') | |
|
|
214 | #opObj11.addParameter(name='frequency', value='9.345e9', format='float') | |
|
|
215 | #procUnitConfObjB.addOperation(name='removeDC') | |
|
|
216 | ||
|
|
217 | #if args.rmDC: | |
|
|
218 | # op.addParameter(name='removeDC', value=1, format='int') | |
|
|
219 | ||
|
|
220 | proc2 = project.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjB.getId()) | |
|
|
221 | proc2.addOperation(name='SpectralMoments') | |
|
|
222 | proc2.addParameter(name='runNextUnit', value=True) | |
|
|
223 | ||
|
|
224 | opObj10 = proc2.addOperation(name="WeatherRadar") | |
|
|
225 | opObj10.addParameter(name='CR_Flag',value=True,format='bool') | |
|
|
226 | opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0])) | |
|
|
227 | opObj10.addParameter(name='Pt',value=200) | |
|
|
228 | opObj10.addParameter(name='min_index',value=max_index(RMIX, sample_rate, ipp, H0,ipp_km)) | |
|
|
229 | ||
|
|
230 | op = proc2.addOperation(name='PedestalInformation') | |
|
|
231 | op.addParameter(name='path', value=path_ped, format='str') | |
|
|
232 | op.addParameter(name='interval', value='0.04') | |
|
|
233 | op.addParameter(name='time_offset', value=time_offset) | |
|
|
234 | op.addParameter(name='mode', value=args.mode) | |
|
|
235 | op.addParameter(name='heading', value=conf['heading']) | |
|
|
236 | ||
|
|
237 | op = proc2.addOperation(name='Block360') | |
|
|
238 | op.addParameter(name='attr_data', value='data_param') | |
|
|
239 | op.addParameter(name='runNextOp', value=True) | |
|
|
240 | op.addParameter(name='angles', value=angles) | |
|
|
241 | op.addParameter(name='heading', value=conf['heading']) | |
|
|
242 | ||
|
|
243 | ||
|
|
244 | if '1' in args.pulses and '2' in args.pulses: | |
|
|
245 | merge = project.addProcUnit(datatype='MergeProc', inputId=[proc1.getId(), proc2.getId()]) | |
|
|
246 | merge.addParameter(name='attr_data', value='data_param') | |
|
|
247 | merge.addParameter(name='mode', value='7') | |
|
|
248 | merge.addParameter(name='index', value=max_index(RMIX, sample_rate, ipp, H0,ipp_km)) | |
|
|
249 | ||
|
|
250 | elif '1' in args.pulses: | |
|
|
251 | merge = proc1 | |
|
|
252 | elif '2' in args.pulses: | |
|
|
253 | merge = proc2 | |
|
|
254 | ||
|
|
255 | ||
|
|
256 | for param in parameters: | |
|
|
257 | ||
|
|
258 | if args.plot: | |
|
|
259 | op= merge.addOperation(name='WeatherParamsPlot') | |
|
|
260 | if args.save: | |
|
|
261 | op.addParameter(name='save', value=path_plots, format='str') | |
|
|
262 | op.addParameter(name='save_period', value=-1) | |
|
|
263 | op.addParameter(name='show', value=args.show) | |
|
|
264 | op.addParameter(name='channels', value='0,') | |
|
|
265 | op.addParameter(name='zmin', value=PARAM[param]['zmin'], format='int') | |
|
|
266 | op.addParameter(name='zmax', value=PARAM[param]['zmax'], format='int') | |
|
|
267 | op.addParameter(name='yrange', value=20, format='int') | |
|
|
268 | op.addParameter(name='xrange', value=args.range, format='int') | |
|
|
269 | op.addParameter(name='attr_data', value=param, format='str') | |
|
|
270 | op.addParameter(name='labels', value=[[PARAM[param]['label']], [PARAM[param]['label']]]) | |
|
|
271 | op.addParameter(name='save_code', value=param) | |
|
|
272 | op.addParameter(name='cb_label', value=PARAM[param]['cb_label']) | |
|
|
273 | op.addParameter(name='colormap', value=PARAM[param]['colormap']) | |
|
|
274 | op.addParameter(name='bgcolor', value='black') | |
|
|
275 | op.addParameter(name='localtime', value=False) | |
|
|
276 | op.addParameter(name='shapes', value='./shapes') | |
|
|
277 | op.addParameter(name='latitude', value=conf['latitude'], format='float') | |
|
|
278 | op.addParameter(name='longitude', value=conf['longitude'], format='float') | |
|
|
279 | op.addParameter(name='map', value=True) | |
|
|
280 | ||
|
|
281 | if MASK: op.addParameter(name='mask', value=MASK, format='float') | |
|
|
282 | #if args.server: | |
|
|
283 | # op.addParameter(name='server', value='190.187.237.239:4444') | |
|
|
284 | # op.addParameter(name='exp_code', value='400') | |
|
|
285 | ||
|
|
286 | desc = { | |
|
|
287 | 'Data': { | |
|
|
288 | 'data_param': {PARAM[param]['wrname']: ['H', 'V']}, | |
|
|
289 | 'utctime': 'time' | |
|
|
290 | }, | |
|
|
291 | 'Metadata': { | |
|
|
292 | 'heightList': 'range', | |
|
|
293 | 'data_azi': 'azimuth', | |
|
|
294 | 'data_ele': 'elevation', | |
|
|
295 | 'mode_op': 'scan_type', | |
|
|
296 | 'h0': 'range_correction', | |
|
|
297 | 'dataPP_NOISE': 'noise', | |
|
|
298 | } | |
|
|
299 | } | |
|
|
300 | ||
|
|
301 | if args.save: | |
|
|
302 | writer = merge.addOperation(name='HDFWriter') | |
|
|
303 | writer.addParameter(name='path', value=path_save, format='str') | |
|
|
304 | writer.addParameter(name='Reset', value=True) | |
|
|
305 | writer.addParameter(name='setType', value='weather') | |
|
|
306 | writer.addParameter(name='setChannel', value='0') #new parameter choose ch 0 H or ch 1 V | |
|
|
307 | writer.addParameter(name='description', value=json.dumps(desc)) | |
|
|
308 | writer.addParameter(name='blocksPerFile', value='1',format='int') | |
|
|
309 | writer.addParameter(name='metadataList', value=','.join(META)) | |
|
|
310 | writer.addParameter(name='dataList', value='data_param,utctime') | |
|
|
311 | writer.addParameter(name='weather_var', value=param) | |
|
|
312 | writer.addParameter(name='mask', value=MASK, format='float') | |
|
|
313 | writer.addParameter(name='localtime', value=False) | |
|
|
314 | # meta | |
|
|
315 | writer.addParameter(name='latitude', value=conf['latitude']) | |
|
|
316 | writer.addParameter(name='longitude', value=conf['longitude']) | |
|
|
317 | writer.addParameter(name='altitude', value=conf['altitude']) | |
|
|
318 | writer.addParameter(name='heading', value=conf['heading']) | |
|
|
319 | writer.addParameter(name='radar_name', value='SOPHy') | |
|
|
320 | writer.addParameter(name='institution', value='IGP') | |
|
|
321 | writer.addParameter(name='contact', value='dscipion@igp.gob.pe') | |
|
|
322 | writer.addParameter(name='created_by', value='Signal Chain (https://pypi.org/project/schainpy/)') | |
|
|
323 | writer.addParameter(name='range_unit', value='km') | |
|
|
324 | writer.addParameter(name='prf', value=1/ipp) | |
|
|
325 | writer.addParameter(name='prf_unit', value='hertz') | |
|
|
326 | writer.addParameter(name='variable', value=PARAM[param]['label']) | |
|
|
327 | writer.addParameter(name='variable_unit', value=PARAM[param]['cb_label']) | |
|
|
328 | writer.addParameter(name='n_pulses', value=n_pulses) | |
|
|
329 | writer.addParameter(name='pulse1_range', value=RMIX) | |
|
|
330 | writer.addParameter(name='pulse1_width', value=pulse_1_width) | |
|
|
331 | writer.addParameter(name='pulse2_width', value=pulse_2_width) | |
|
|
332 | writer.addParameter(name='pulse1_repetitions', value=pulse_1_repetitions) | |
|
|
333 | writer.addParameter(name='pulse2_repetitions', value=pulse_2_repetitions) | |
|
|
334 | writer.addParameter(name='pulse_width_unit', value='microseconds') | |
|
|
335 | writer.addParameter(name='snr_threshold', value=MASK) | |
|
|
336 | writer.addParameter(name='cr_hv', value=[67.41,67.17]) #new parameter | |
|
|
337 | ||
|
|
338 | return project | |
|
|
339 | ||
|
|
340 | if __name__ == '__main__': | |
|
|
341 | ||
|
|
342 | parser = argparse.ArgumentParser(description='Script to process SOPHy data.') | |
|
|
343 | parser.add_argument('experiment', | |
|
|
344 | help='Experiment name') | |
|
|
345 | parser.add_argument('--parameters', nargs='*', default=['S'], | |
|
|
346 | help='Variables to process: P, Z, V') | |
|
|
347 | parser.add_argument('--pulses', nargs='*', default=['1', '2'], | |
|
|
348 | help='Variables to process: 1, 2') | |
|
|
349 | parser.add_argument('--angles', nargs='*', default=[], type=int, | |
|
|
350 | help='Angles to process') | |
|
|
351 | parser.add_argument('--time_offset', default=0, | |
|
|
352 | help='Fix time offset') | |
|
|
353 | parser.add_argument('--range', default=60, type=float, | |
|
|
354 | help='Max range to plot') | |
|
|
355 | ||
|
|
356 | parser.add_argument('--save', action='store_true', | |
|
|
357 | help='Create output files') | |
|
|
358 | parser.add_argument('--plot', action='store_true', | |
|
|
359 | help='Create plot files') | |
|
|
360 | parser.add_argument('--show', action='store_true', | |
|
|
361 | help='Show matplotlib plot.') | |
|
|
362 | parser.add_argument('--online', action='store_true', | |
|
|
363 | help='Set online mode.') | |
|
|
364 | parser.add_argument('--server', action='store_true', | |
|
|
365 | help='Send to realtime') | |
|
|
366 | parser.add_argument('--start_time', default='', | |
|
|
367 | help='Set start time.') | |
|
|
368 | parser.add_argument('--end_time', default='', | |
|
|
369 | help='Set end time.') | |
|
|
370 | parser.add_argument('--label', default='', | |
|
|
371 | help='Label for plot & param folder') | |
|
|
372 | parser.add_argument('--mode', default=None, | |
|
|
373 | help='Type of scan') | |
|
|
374 | #parser.add_argument('--rmDC', action='store_true', | |
|
|
375 | # help='Apply remove DC.') | |
|
|
376 | parser.add_argument('--mask', default=0, type=float, | |
|
|
377 | help='Set SNR threshold.') | |
|
|
378 | args = parser.parse_args() | |
|
|
379 | ||
|
|
380 | project = main(args) | |
|
|
381 | project.start() No newline at end of file | |
| @@ -1,769 +1,784 | |||
|
|
1 | 1 | import os |
|
|
2 | 2 | import datetime |
|
|
3 | 3 | import warnings |
|
|
4 | 4 | import numpy |
|
|
5 | 5 | from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter |
|
|
6 | 6 | from matplotlib.patches import Circle |
|
|
7 | 7 | from cartopy.feature import ShapelyFeature |
|
|
8 | 8 | import cartopy.io.shapereader as shpreader |
|
|
9 | 9 | |
|
|
10 | 10 | from schainpy.model.graphics.jroplot_base import Plot, plt, ccrs |
|
|
11 | 11 | from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot |
|
|
12 | 12 | from schainpy.utils import log |
|
|
13 | 13 | from schainpy.model.graphics.plotting_codes import cb_tables |
|
|
14 | 14 | |
|
|
15 | 15 | |
|
|
16 | 16 | EARTH_RADIUS = 6.3710e3 |
|
|
17 | 17 | |
|
|
18 | 18 | |
|
|
19 | 19 | def antenna_to_cartesian(ranges, azimuths, elevations): |
|
|
20 | 20 | """ |
|
|
21 | 21 | Return Cartesian coordinates from antenna coordinates. |
|
|
22 | 22 | |
|
|
23 | 23 | Parameters |
|
|
24 | 24 | ---------- |
|
|
25 | 25 | ranges : array |
|
|
26 | 26 | Distances to the center of the radar gates (bins) in kilometers. |
|
|
27 | 27 | azimuths : array |
|
|
28 | 28 | Azimuth angle of the radar in degrees. |
|
|
29 | 29 | elevations : array |
|
|
30 | 30 | Elevation angle of the radar in degrees. |
|
|
31 | 31 | |
|
|
32 | 32 | Returns |
|
|
33 | 33 | ------- |
|
|
34 | 34 | x, y, z : array |
|
|
35 | 35 | Cartesian coordinates in meters from the radar. |
|
|
36 | 36 | |
|
|
37 | 37 | Notes |
|
|
38 | 38 | ----- |
|
|
39 | 39 | The calculation for Cartesian coordinate is adapted from equations |
|
|
40 | 40 | 2.28(b) and 2.28(c) of Doviak and Zrnic [1]_ assuming a |
|
|
41 | 41 | standard atmosphere (4/3 Earth's radius model). |
|
|
42 | 42 | |
|
|
43 | 43 | .. math:: |
|
|
44 | 44 | |
|
|
45 | 45 | z = \\sqrt{r^2+R^2+2*r*R*sin(\\theta_e)} - R |
|
|
46 | 46 | |
|
|
47 | 47 | s = R * arcsin(\\frac{r*cos(\\theta_e)}{R+z}) |
|
|
48 | 48 | |
|
|
49 | 49 | x = s * sin(\\theta_a) |
|
|
50 | 50 | |
|
|
51 | 51 | y = s * cos(\\theta_a) |
|
|
52 | 52 | |
|
|
53 | 53 | Where r is the distance from the radar to the center of the gate, |
|
|
54 | 54 | :math:`\\theta_a` is the azimuth angle, :math:`\\theta_e` is the |
|
|
55 | 55 | elevation angle, s is the arc length, and R is the effective radius |
|
|
56 | 56 | of the earth, taken to be 4/3 the mean radius of earth (6371 km). |
|
|
57 | 57 | |
|
|
58 | 58 | References |
|
|
59 | 59 | ---------- |
|
|
60 | 60 | .. [1] Doviak and Zrnic, Doppler Radar and Weather Observations, Second |
|
|
61 | 61 | Edition, 1993, p. 21. |
|
|
62 | 62 | |
|
|
63 | 63 | """ |
|
|
64 | 64 | theta_e = numpy.deg2rad(elevations) # elevation angle in radians. |
|
|
65 | 65 | theta_a = numpy.deg2rad(azimuths) # azimuth angle in radians. |
|
|
66 | 66 | R = 6371.0 * 1000.0 * 4.0 / 3.0 # effective radius of earth in meters. |
|
|
67 | 67 | r = ranges * 1000.0 # distances to gates in meters. |
|
|
68 | 68 | |
|
|
69 | 69 | z = (r ** 2 + R ** 2 + 2.0 * r * R * numpy.sin(theta_e)) ** 0.5 - R |
|
|
70 | 70 | s = R * numpy.arcsin(r * numpy.cos(theta_e) / (R + z)) # arc length in m. |
|
|
71 | 71 | x = s * numpy.sin(theta_a) |
|
|
72 | 72 | y = s * numpy.cos(theta_a) |
|
|
73 | 73 | return x, y, z |
|
|
74 | 74 | |
|
|
75 | 75 | def cartesian_to_geographic_aeqd(x, y, lon_0, lat_0, R=EARTH_RADIUS): |
|
|
76 | 76 | """ |
|
|
77 | 77 | Azimuthal equidistant Cartesian to geographic coordinate transform. |
|
|
78 | 78 | |
|
|
79 | 79 | Transform a set of Cartesian/Cartographic coordinates (x, y) to |
|
|
80 | 80 | geographic coordinate system (lat, lon) using a azimuthal equidistant |
|
|
81 | 81 | map projection [1]_. |
|
|
82 | 82 | |
|
|
83 | 83 | .. math:: |
|
|
84 | 84 | |
|
|
85 | 85 | lat = \\arcsin(\\cos(c) * \\sin(lat_0) + |
|
|
86 | 86 | (y * \\sin(c) * \\cos(lat_0) / \\rho)) |
|
|
87 | 87 | |
|
|
88 | 88 | lon = lon_0 + \\arctan2( |
|
|
89 | 89 | x * \\sin(c), |
|
|
90 | 90 | \\rho * \\cos(lat_0) * \\cos(c) - y * \\sin(lat_0) * \\sin(c)) |
|
|
91 | 91 | |
|
|
92 | 92 | \\rho = \\sqrt(x^2 + y^2) |
|
|
93 | 93 | |
|
|
94 | 94 | c = \\rho / R |
|
|
95 | 95 | |
|
|
96 | 96 | Where x, y are the Cartesian position from the center of projection; |
|
|
97 | 97 | lat, lon the corresponding latitude and longitude; lat_0, lon_0 are the |
|
|
98 | 98 | latitude and longitude of the center of the projection; R is the radius of |
|
|
99 | 99 | the earth (defaults to ~6371 km). lon is adjusted to be between -180 and |
|
|
100 | 100 | 180. |
|
|
101 | 101 | |
|
|
102 | 102 | Parameters |
|
|
103 | 103 | ---------- |
|
|
104 | 104 | x, y : array-like |
|
|
105 | 105 | Cartesian coordinates in the same units as R, typically meters. |
|
|
106 | 106 | lon_0, lat_0 : float |
|
|
107 | 107 | Longitude and latitude, in degrees, of the center of the projection. |
|
|
108 | 108 | R : float, optional |
|
|
109 | 109 | Earth radius in the same units as x and y. The default value is in |
|
|
110 | 110 | units of meters. |
|
|
111 | 111 | |
|
|
112 | 112 | Returns |
|
|
113 | 113 | ------- |
|
|
114 | 114 | lon, lat : array |
|
|
115 | 115 | Longitude and latitude of Cartesian coordinates in degrees. |
|
|
116 | 116 | |
|
|
117 | 117 | References |
|
|
118 | 118 | ---------- |
|
|
119 | 119 | .. [1] Snyder, J. P. Map Projections--A Working Manual. U. S. Geological |
|
|
120 | 120 | Survey Professional Paper 1395, 1987, pp. 191-202. |
|
|
121 | 121 | |
|
|
122 | 122 | """ |
|
|
123 | 123 | x = numpy.atleast_1d(numpy.asarray(x)) |
|
|
124 | 124 | y = numpy.atleast_1d(numpy.asarray(y)) |
|
|
125 | 125 | |
|
|
126 | 126 | lat_0_rad = numpy.deg2rad(lat_0) |
|
|
127 | 127 | lon_0_rad = numpy.deg2rad(lon_0) |
|
|
128 | 128 | |
|
|
129 | 129 | rho = numpy.sqrt(x*x + y*y) |
|
|
130 | 130 | c = rho / R |
|
|
131 | 131 | |
|
|
132 | 132 | with warnings.catch_warnings(): |
|
|
133 | 133 | # division by zero may occur here but is properly addressed below so |
|
|
134 | 134 | # the warnings can be ignored |
|
|
135 | 135 | warnings.simplefilter("ignore", RuntimeWarning) |
|
|
136 | 136 | lat_rad = numpy.arcsin(numpy.cos(c) * numpy.sin(lat_0_rad) + |
|
|
137 | 137 | y * numpy.sin(c) * numpy.cos(lat_0_rad) / rho) |
|
|
138 | 138 | lat_deg = numpy.rad2deg(lat_rad) |
|
|
139 | 139 | # fix cases where the distance from the center of the projection is zero |
|
|
140 | 140 | lat_deg[rho == 0] = lat_0 |
|
|
141 | 141 | |
|
|
142 | 142 | x1 = x * numpy.sin(c) |
|
|
143 | 143 | x2 = rho*numpy.cos(lat_0_rad)*numpy.cos(c) - y*numpy.sin(lat_0_rad)*numpy.sin(c) |
|
|
144 | 144 | lon_rad = lon_0_rad + numpy.arctan2(x1, x2) |
|
|
145 | 145 | lon_deg = numpy.rad2deg(lon_rad) |
|
|
146 | 146 | # Longitudes should be from -180 to 180 degrees |
|
|
147 | 147 | lon_deg[lon_deg > 180] -= 360. |
|
|
148 | 148 | lon_deg[lon_deg < -180] += 360. |
|
|
149 | 149 | |
|
|
150 | 150 | return lon_deg, lat_deg |
|
|
151 | 151 | |
|
|
152 | 152 | def antenna_to_geographic(ranges, azimuths, elevations, site): |
|
|
153 | 153 | |
|
|
154 | 154 | x, y, z = antenna_to_cartesian(numpy.array(ranges), numpy.array(azimuths), numpy.array(elevations)) |
|
|
155 | 155 | lon, lat = cartesian_to_geographic_aeqd(x, y, site[0], site[1], R=6370997.) |
|
|
156 | 156 | |
|
|
157 | 157 | return lon, lat |
|
|
158 | 158 | |
|
|
159 | 159 | def ll2xy(lat1, lon1, lat2, lon2): |
|
|
160 | 160 | |
|
|
161 | 161 | p = 0.017453292519943295 |
|
|
162 | 162 | a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \ |
|
|
163 | 163 | numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2 |
|
|
164 | 164 | r = 12742 * numpy.arcsin(numpy.sqrt(a)) |
|
|
165 | 165 | theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p) |
|
|
166 | 166 | * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p)) |
|
|
167 | 167 | theta = -theta + numpy.pi/2 |
|
|
168 | 168 | return r*numpy.cos(theta), r*numpy.sin(theta) |
|
|
169 | 169 | |
|
|
170 | 170 | |
|
|
171 | 171 | def km2deg(km): |
|
|
172 | 172 | ''' |
|
|
173 | 173 | Convert distance in km to degrees |
|
|
174 | 174 | ''' |
|
|
175 | 175 | |
|
|
176 | 176 | return numpy.rad2deg(km/EARTH_RADIUS) |
|
|
177 | 177 | |
|
|
178 | 178 | |
|
|
179 | 179 | |
|
|
180 | 180 | class SpectralMomentsPlot(SpectraPlot): |
|
|
181 | 181 | ''' |
|
|
182 | 182 | Plot for Spectral Moments |
|
|
183 | 183 | ''' |
|
|
184 | 184 | CODE = 'spc_moments' |
|
|
185 | 185 | # colormap = 'jet' |
|
|
186 | 186 | # plot_type = 'pcolor' |
|
|
187 | 187 | |
|
|
188 | 188 | class DobleGaussianPlot(SpectraPlot): |
|
|
189 | 189 | ''' |
|
|
190 | 190 | Plot for Double Gaussian Plot |
|
|
191 | 191 | ''' |
|
|
192 | 192 | CODE = 'gaussian_fit' |
|
|
193 | 193 | # colormap = 'jet' |
|
|
194 | 194 | # plot_type = 'pcolor' |
|
|
195 | 195 | |
|
|
196 | 196 | class DoubleGaussianSpectraCutPlot(SpectraCutPlot): |
|
|
197 | 197 | ''' |
|
|
198 | 198 | Plot SpectraCut with Double Gaussian Fit |
|
|
199 | 199 | ''' |
|
|
200 | 200 | CODE = 'cut_gaussian_fit' |
|
|
201 | 201 | |
|
|
202 | 202 | class SnrPlot(RTIPlot): |
|
|
203 | 203 | ''' |
|
|
204 | 204 | Plot for SNR Data |
|
|
205 | 205 | ''' |
|
|
206 | 206 | |
|
|
207 | 207 | CODE = 'snr' |
|
|
208 | 208 | colormap = 'jet' |
|
|
209 | 209 | |
|
|
210 | 210 | def update(self, dataOut): |
|
|
211 | 211 | |
|
|
212 | 212 | data = { |
|
|
213 | 213 | 'snr': 10*numpy.log10(dataOut.data_snr) |
|
|
214 | 214 | } |
|
|
215 | 215 | |
|
|
216 | 216 | return data, {} |
|
|
217 | 217 | |
|
|
218 | 218 | class DopplerPlot(RTIPlot): |
|
|
219 | 219 | ''' |
|
|
220 | 220 | Plot for DOPPLER Data (1st moment) |
|
|
221 | 221 | ''' |
|
|
222 | 222 | |
|
|
223 | 223 | CODE = 'dop' |
|
|
224 | 224 | colormap = 'jet' |
|
|
225 | 225 | |
|
|
226 | 226 | def update(self, dataOut): |
|
|
227 | 227 | |
|
|
228 | 228 | data = { |
|
|
229 | 229 | 'dop': 10*numpy.log10(dataOut.data_dop) |
|
|
230 | 230 | } |
|
|
231 | 231 | |
|
|
232 | 232 | return data, {} |
|
|
233 | 233 | |
|
|
234 | 234 | class PowerPlot(RTIPlot): |
|
|
235 | 235 | ''' |
|
|
236 | 236 | Plot for Power Data (0 moment) |
|
|
237 | 237 | ''' |
|
|
238 | 238 | |
|
|
239 | 239 | CODE = 'pow' |
|
|
240 | 240 | colormap = 'jet' |
|
|
241 | 241 | |
|
|
242 | 242 | def update(self, dataOut): |
|
|
243 | 243 | data = { |
|
|
244 | 244 | 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor) |
|
|
245 | 245 | } |
|
|
246 | 246 | return data, {} |
|
|
247 | 247 | |
|
|
248 | 248 | class SpectralWidthPlot(RTIPlot): |
|
|
249 | 249 | ''' |
|
|
250 | 250 | Plot for Spectral Width Data (2nd moment) |
|
|
251 | 251 | ''' |
|
|
252 | 252 | |
|
|
253 | 253 | CODE = 'width' |
|
|
254 | 254 | colormap = 'jet' |
|
|
255 | 255 | |
|
|
256 | 256 | def update(self, dataOut): |
|
|
257 | 257 | |
|
|
258 | 258 | data = { |
|
|
259 | 259 | 'width': dataOut.data_width |
|
|
260 | 260 | } |
|
|
261 | 261 | |
|
|
262 | 262 | return data, {} |
|
|
263 | 263 | |
|
|
264 | 264 | class SkyMapPlot(Plot): |
|
|
265 | 265 | ''' |
|
|
266 | 266 | Plot for meteors detection data |
|
|
267 | 267 | ''' |
|
|
268 | 268 | |
|
|
269 | 269 | CODE = 'param' |
|
|
270 | 270 | |
|
|
271 | 271 | def setup(self): |
|
|
272 | 272 | |
|
|
273 | 273 | self.ncols = 1 |
|
|
274 | 274 | self.nrows = 1 |
|
|
275 | 275 | self.width = 7.2 |
|
|
276 | 276 | self.height = 7.2 |
|
|
277 | 277 | self.nplots = 1 |
|
|
278 | 278 | self.xlabel = 'Zonal Zenith Angle (deg)' |
|
|
279 | 279 | self.ylabel = 'Meridional Zenith Angle (deg)' |
|
|
280 | 280 | self.polar = True |
|
|
281 | 281 | self.ymin = -180 |
|
|
282 | 282 | self.ymax = 180 |
|
|
283 | 283 | self.colorbar = False |
|
|
284 | 284 | |
|
|
285 | 285 | def plot(self): |
|
|
286 | 286 | |
|
|
287 | 287 | arrayParameters = numpy.concatenate(self.data['param']) |
|
|
288 | 288 | error = arrayParameters[:, -1] |
|
|
289 | 289 | indValid = numpy.where(error == 0)[0] |
|
|
290 | 290 | finalMeteor = arrayParameters[indValid, :] |
|
|
291 | 291 | finalAzimuth = finalMeteor[:, 3] |
|
|
292 | 292 | finalZenith = finalMeteor[:, 4] |
|
|
293 | 293 | |
|
|
294 | 294 | x = finalAzimuth * numpy.pi / 180 |
|
|
295 | 295 | y = finalZenith |
|
|
296 | 296 | |
|
|
297 | 297 | ax = self.axes[0] |
|
|
298 | 298 | |
|
|
299 | 299 | if ax.firsttime: |
|
|
300 | 300 | ax.plot = ax.plot(x, y, 'bo', markersize=5)[0] |
|
|
301 | 301 | else: |
|
|
302 | 302 | ax.plot.set_data(x, y) |
|
|
303 | 303 | |
|
|
304 | 304 | dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S') |
|
|
305 | 305 | dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S') |
|
|
306 | 306 | title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1, |
|
|
307 | 307 | dt2, |
|
|
308 | 308 | len(x)) |
|
|
309 | 309 | self.titles[0] = title |
|
|
310 | 310 | |
|
|
311 | 311 | |
|
|
312 | 312 | class GenericRTIPlot(Plot): |
|
|
313 | 313 | ''' |
|
|
314 | 314 | Plot for data_xxxx object |
|
|
315 | 315 | ''' |
|
|
316 | 316 | |
|
|
317 | 317 | CODE = 'param' |
|
|
318 | 318 | colormap = 'viridis' |
|
|
319 | 319 | plot_type = 'pcolorbuffer' |
|
|
320 | 320 | |
|
|
321 | 321 | def setup(self): |
|
|
322 | 322 | self.xaxis = 'time' |
|
|
323 | 323 | self.ncols = 1 |
|
|
324 | 324 | self.nrows = self.data.shape('param')[0] |
|
|
325 | 325 | self.nplots = self.nrows |
|
|
326 | 326 | self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95}) |
|
|
327 | 327 | |
|
|
328 | 328 | if not self.xlabel: |
|
|
329 | 329 | self.xlabel = 'Time' |
|
|
330 | 330 | |
|
|
331 | 331 | self.ylabel = 'Range [km]' |
|
|
332 | 332 | if not self.titles: |
|
|
333 | 333 | self.titles = ['Param {}'.format(x) for x in range(self.nrows)] |
|
|
334 | 334 | |
|
|
335 | 335 | def update(self, dataOut): |
|
|
336 | 336 | |
|
|
337 | 337 | data = { |
|
|
338 | 338 | 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0) |
|
|
339 | 339 | } |
|
|
340 | 340 | |
|
|
341 | 341 | meta = {} |
|
|
342 | 342 | |
|
|
343 | 343 | return data, meta |
|
|
344 | 344 | |
|
|
345 | 345 | def plot(self): |
|
|
346 | 346 | # self.data.normalize_heights() |
|
|
347 | 347 | self.x = self.data.times |
|
|
348 | 348 | self.y = self.data.yrange |
|
|
349 | 349 | self.z = self.data['param'] |
|
|
350 | 350 | self.z = 10*numpy.log10(self.z) |
|
|
351 | 351 | self.z = numpy.ma.masked_invalid(self.z) |
|
|
352 | 352 | |
|
|
353 | 353 | if self.decimation is None: |
|
|
354 | 354 | x, y, z = self.fill_gaps(self.x, self.y, self.z) |
|
|
355 | 355 | else: |
|
|
356 | 356 | x, y, z = self.fill_gaps(*self.decimate()) |
|
|
357 | 357 | |
|
|
358 | 358 | for n, ax in enumerate(self.axes): |
|
|
359 | 359 | |
|
|
360 | 360 | self.zmax = self.zmax if self.zmax is not None else numpy.max( |
|
|
361 | 361 | self.z[n]) |
|
|
362 | 362 | self.zmin = self.zmin if self.zmin is not None else numpy.min( |
|
|
363 | 363 | self.z[n]) |
|
|
364 | 364 | |
|
|
365 | 365 | if ax.firsttime: |
|
|
366 | 366 | if self.zlimits is not None: |
|
|
367 | 367 | self.zmin, self.zmax = self.zlimits[n] |
|
|
368 | 368 | |
|
|
369 | 369 | ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n], |
|
|
370 | 370 | vmin=self.zmin, |
|
|
371 | 371 | vmax=self.zmax, |
|
|
372 | 372 | cmap=self.cmaps[n] |
|
|
373 | 373 | ) |
|
|
374 | 374 | else: |
|
|
375 | 375 | if self.zlimits is not None: |
|
|
376 | 376 | self.zmin, self.zmax = self.zlimits[n] |
|
|
377 | 377 | ax.collections.remove(ax.collections[0]) |
|
|
378 | 378 | ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n], |
|
|
379 | 379 | vmin=self.zmin, |
|
|
380 | 380 | vmax=self.zmax, |
|
|
381 | 381 | cmap=self.cmaps[n] |
|
|
382 | 382 | ) |
|
|
383 | 383 | |
|
|
384 | 384 | |
|
|
385 | 385 | class PolarMapPlot(Plot): |
|
|
386 | 386 | ''' |
|
|
387 | 387 | Plot for weather radar |
|
|
388 | 388 | ''' |
|
|
389 | 389 | |
|
|
390 | 390 | CODE = 'param' |
|
|
391 | 391 | colormap = 'seismic' |
|
|
392 | 392 | |
|
|
393 | 393 | def setup(self): |
|
|
394 | 394 | self.ncols = 1 |
|
|
395 | 395 | self.nrows = 1 |
|
|
396 | 396 | self.width = 9 |
|
|
397 | 397 | self.height = 8 |
|
|
398 | 398 | self.mode = self.data.meta['mode'] |
|
|
399 | 399 | if self.channels is not None: |
|
|
400 | 400 | self.nplots = len(self.channels) |
|
|
401 | 401 | self.nrows = len(self.channels) |
|
|
402 | 402 | else: |
|
|
403 | 403 | self.nplots = self.data.shape(self.CODE)[0] |
|
|
404 | 404 | self.nrows = self.nplots |
|
|
405 | 405 | self.channels = list(range(self.nplots)) |
|
|
406 | 406 | if self.mode == 'E': |
|
|
407 | 407 | self.xlabel = 'Longitude' |
|
|
408 | 408 | self.ylabel = 'Latitude' |
|
|
409 | 409 | else: |
|
|
410 | 410 | self.xlabel = 'Range (km)' |
|
|
411 | 411 | self.ylabel = 'Height (km)' |
|
|
412 | 412 | self.bgcolor = 'white' |
|
|
413 | 413 | self.cb_labels = self.data.meta['units'] |
|
|
414 | 414 | self.lat = self.data.meta['latitude'] |
|
|
415 | 415 | self.lon = self.data.meta['longitude'] |
|
|
416 | 416 | self.xmin, self.xmax = float( |
|
|
417 | 417 | km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon) |
|
|
418 | 418 | self.ymin, self.ymax = float( |
|
|
419 | 419 | km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat) |
|
|
420 | 420 | # self.polar = True |
|
|
421 | 421 | |
|
|
422 | 422 | def plot(self): |
|
|
423 | 423 | |
|
|
424 | 424 | for n, ax in enumerate(self.axes): |
|
|
425 | 425 | data = self.data['param'][self.channels[n]] |
|
|
426 | 426 | |
|
|
427 | 427 | zeniths = numpy.linspace( |
|
|
428 | 428 | 0, self.data.meta['max_range'], data.shape[1]) |
|
|
429 | 429 | if self.mode == 'E': |
|
|
430 | 430 | azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2 |
|
|
431 | 431 | r, theta = numpy.meshgrid(zeniths, azimuths) |
|
|
432 | 432 | x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin( |
|
|
433 | 433 | theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])) |
|
|
434 | 434 | x = km2deg(x) + self.lon |
|
|
435 | 435 | y = km2deg(y) + self.lat |
|
|
436 | 436 | else: |
|
|
437 | 437 | azimuths = numpy.radians(self.data.yrange) |
|
|
438 | 438 | r, theta = numpy.meshgrid(zeniths, azimuths) |
|
|
439 | 439 | x, y = r*numpy.cos(theta), r*numpy.sin(theta) |
|
|
440 | 440 | self.y = zeniths |
|
|
441 | 441 | |
|
|
442 | 442 | if ax.firsttime: |
|
|
443 | 443 | if self.zlimits is not None: |
|
|
444 | 444 | self.zmin, self.zmax = self.zlimits[n] |
|
|
445 | 445 | ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)), |
|
|
446 | 446 | x, y, numpy.ma.array(data, mask=numpy.isnan(data)), |
|
|
447 | 447 | vmin=self.zmin, |
|
|
448 | 448 | vmax=self.zmax, |
|
|
449 | 449 | cmap=self.cmaps[n]) |
|
|
450 | 450 | else: |
|
|
451 | 451 | if self.zlimits is not None: |
|
|
452 | 452 | self.zmin, self.zmax = self.zlimits[n] |
|
|
453 | 453 | ax.collections.remove(ax.collections[0]) |
|
|
454 | 454 | ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)), |
|
|
455 | 455 | x, y, numpy.ma.array(data, mask=numpy.isnan(data)), |
|
|
456 | 456 | vmin=self.zmin, |
|
|
457 | 457 | vmax=self.zmax, |
|
|
458 | 458 | cmap=self.cmaps[n]) |
|
|
459 | 459 | |
|
|
460 | 460 | if self.mode == 'A': |
|
|
461 | 461 | continue |
|
|
462 | 462 | |
|
|
463 | 463 | # plot district names |
|
|
464 | 464 | f = open('/data/workspace/schain_scripts/distrito.csv') |
|
|
465 | 465 | for line in f: |
|
|
466 | 466 | label, lon, lat = [s.strip() for s in line.split(',') if s] |
|
|
467 | 467 | lat = float(lat) |
|
|
468 | 468 | lon = float(lon) |
|
|
469 | 469 | # ax.plot(lon, lat, '.b', ms=2) |
|
|
470 | 470 | ax.text(lon, lat, label.decode('utf8'), ha='center', |
|
|
471 | 471 | va='bottom', size='8', color='black') |
|
|
472 | 472 | |
|
|
473 | 473 | # plot limites |
|
|
474 | 474 | limites = [] |
|
|
475 | 475 | tmp = [] |
|
|
476 | 476 | for line in open('/data/workspace/schain_scripts/lima.csv'): |
|
|
477 | 477 | if '#' in line: |
|
|
478 | 478 | if tmp: |
|
|
479 | 479 | limites.append(tmp) |
|
|
480 | 480 | tmp = [] |
|
|
481 | 481 | continue |
|
|
482 | 482 | values = line.strip().split(',') |
|
|
483 | 483 | tmp.append((float(values[0]), float(values[1]))) |
|
|
484 | 484 | for points in limites: |
|
|
485 | 485 | ax.add_patch( |
|
|
486 | 486 | Polygon(points, ec='k', fc='none', ls='--', lw=0.5)) |
|
|
487 | 487 | |
|
|
488 | 488 | # plot Cuencas |
|
|
489 | 489 | for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'): |
|
|
490 | 490 | f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca)) |
|
|
491 | 491 | values = [line.strip().split(',') for line in f] |
|
|
492 | 492 | points = [(float(s[0]), float(s[1])) for s in values] |
|
|
493 | 493 | ax.add_patch(Polygon(points, ec='b', fc='none')) |
|
|
494 | 494 | |
|
|
495 | 495 | # plot grid |
|
|
496 | 496 | for r in (15, 30, 45, 60): |
|
|
497 | 497 | ax.add_artist(plt.Circle((self.lon, self.lat), |
|
|
498 | 498 | km2deg(r), color='0.6', fill=False, lw=0.2)) |
|
|
499 | 499 | ax.text( |
|
|
500 | 500 | self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180), |
|
|
501 | 501 | self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180), |
|
|
502 | 502 | '{}km'.format(r), |
|
|
503 | 503 | ha='center', va='bottom', size='8', color='0.6', weight='heavy') |
|
|
504 | 504 | |
|
|
505 | 505 | if self.mode == 'E': |
|
|
506 | 506 | title = 'El={}$^\circ$'.format(self.data.meta['elevation']) |
|
|
507 | 507 | label = 'E{:02d}'.format(int(self.data.meta['elevation'])) |
|
|
508 | 508 | else: |
|
|
509 | 509 | title = 'Az={}$^\circ$'.format(self.data.meta['azimuth']) |
|
|
510 | 510 | label = 'A{:02d}'.format(int(self.data.meta['azimuth'])) |
|
|
511 | 511 | |
|
|
512 | 512 | self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels] |
|
|
513 | 513 | self.titles = ['{} {}'.format( |
|
|
514 | 514 | self.data.parameters[x], title) for x in self.channels] |
|
|
515 | 515 | |
|
|
516 | ||
|
|
517 | ||
|
|
516 | 518 | class WeatherParamsPlot(Plot): |
|
|
517 | ||
|
|
519 | ||
|
|
518 | 520 | plot_type = 'scattermap' |
|
|
519 | 521 | buffering = False |
|
|
520 | 522 | |
|
|
521 | 523 | def setup(self): |
|
|
522 | 524 | |
|
|
523 | 525 | self.ncols = 1 |
|
|
524 | 526 | self.nrows = 1 |
|
|
525 | 527 | self.nplots= 1 |
|
|
526 | 528 | |
|
|
527 | 529 | if self.channels is not None: |
|
|
528 | 530 | self.nplots = len(self.channels) |
|
|
529 | 531 | self.ncols = len(self.channels) |
|
|
530 | 532 | else: |
|
|
531 | 533 | self.nplots = self.data.shape(self.CODE)[0] |
|
|
532 | 534 | self.ncols = self.nplots |
|
|
533 | 535 | self.channels = list(range(self.nplots)) |
|
|
534 | 536 | |
|
|
535 | 537 | self.colorbar=True |
|
|
536 | 538 | if len(self.channels)>1: |
|
|
537 | 539 | self.width = 12 |
|
|
538 | 540 | else: |
|
|
539 | 541 | self.width =8 |
|
|
540 | 542 | self.height =7 |
|
|
541 | 543 | self.ini =0 |
|
|
542 | 544 | self.len_azi =0 |
|
|
543 | 545 | self.buffer_ini = None |
|
|
544 | 546 | self.buffer_ele = None |
|
|
545 | 547 | self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.1}) |
|
|
546 | 548 | self.flag =0 |
|
|
547 | 549 | self.indicador= 0 |
|
|
548 | 550 | self.last_data_ele = None |
|
|
549 | 551 | self.val_mean = None |
|
|
550 | 552 | |
|
|
551 | 553 | def update(self, dataOut): |
|
|
552 | 554 | |
|
|
553 | 555 | vars = { |
|
|
554 | 556 | 'S' : 0, |
|
|
555 | 557 | 'V' : 1, |
|
|
556 | 558 | 'W' : 2, |
|
|
557 | 559 | 'SNR' : 3, |
|
|
558 | 560 | 'Z' : 4, |
|
|
559 | 561 | 'D' : 5, |
|
|
560 | 562 | 'P' : 6, |
|
|
561 | 563 | 'R' : 7, |
|
|
562 | 564 | } |
|
|
563 | 565 | |
|
|
564 | 566 | data = {} |
|
|
565 | 567 | meta = {} |
|
|
566 | 568 | |
|
|
567 | if hasattr(dataOut, 'nFFTPoints'): | |
|
|
568 | factor = dataOut.normFactor | |
|
|
569 | else: | |
|
|
570 | factor = 1 | |
|
|
569 | ##if hasattr(dataOut, 'nFFTPoints'): | |
|
|
570 | ## factor = dataOut.normFactor*10.0 # CONSIDERACION ENTRE PULSE PAIR Y FFT | |
|
|
571 | ##else: | |
|
|
572 | ## factor = 1 | |
|
|
571 | 573 | |
|
|
572 | 574 | if hasattr(dataOut, 'dparam'): |
|
|
573 | 575 | tmp = getattr(dataOut, 'data_param') |
|
|
574 | 576 | else: |
|
|
575 | 577 | #print("-------------------self.attr_data[0]",self.attr_data[0]) |
|
|
576 | 578 | if 'S' in self.attr_data[0]: |
|
|
577 | 579 | if self.attr_data[0]=='S': |
|
|
578 |
tmp = 10*numpy.log10(10.0*getattr(dataOut, 'data_param')[:,0,:] |
|
|
|
580 | tmp = 10*numpy.log10(10.0*getattr(dataOut, 'data_param')[:,0,:]) ## /(factor)) ya no considerar factor se aplica factor jroproc_parametrs | |
|
|
579 | 581 | if self.attr_data[0]=='SNR': |
|
|
580 | 582 | tmp = 10*numpy.log10(getattr(dataOut, 'data_param')[:,3,:]) |
|
|
581 | 583 | else: |
|
|
582 | 584 | tmp = getattr(dataOut, 'data_param')[:,vars[self.attr_data[0]],:] |
|
|
583 | 585 | |
|
|
584 | 586 | if self.mask: |
|
|
585 | 587 | mask = dataOut.data_param[:,3,:] < self.mask |
|
|
586 | 588 | tmp[mask] = numpy.nan |
|
|
587 | 589 | mask = numpy.nansum((tmp, numpy.roll(tmp, 1),numpy.roll(tmp, -1)), axis=0) == tmp |
|
|
588 | 590 | tmp[mask] = numpy.nan |
|
|
591 | ||
|
|
592 | #################################################################### | |
|
|
593 | #SE GUARDAN LOS DATOS DE LOS PARAMETROS YA SEA PP O SPECTRA EN UN ARCHIVO .npy | |
|
|
594 | ##elapsed_time = time.time() - self.start_time | |
|
|
595 | ##filename = f'{dataOut.inputUnit}_{self.attr_data[0]}_{elapsed_time:.0f}.npy' # Nombre único con timestamp | |
|
|
596 | ||
|
|
597 | # Guardar el array en el nuevo archivo | |
|
|
598 | ##with open(filename, 'wb') as f: | |
|
|
599 | ## numpy.save(f, tmp) | |
|
|
600 | ||
|
|
601 | ##print("Se creó el archivo:", filename) | |
|
|
602 | ||
|
|
589 | 603 | |
|
|
604 | ##################################################################### | |
|
|
590 | 605 | r = dataOut.heightList |
|
|
591 | 606 | delta_height = r[1]-r[0] |
|
|
592 | 607 | valid = numpy.where(r>=0)[0] |
|
|
593 | 608 | data['r'] = numpy.arange(len(valid))*delta_height |
|
|
594 | 609 | |
|
|
595 | 610 | data['data'] = [0, 0] |
|
|
596 | 611 | |
|
|
597 | 612 | try: |
|
|
598 | 613 | data['data'][0] = tmp[0][:,valid] |
|
|
599 | 614 | data['data'][1] = tmp[1][:,valid] |
|
|
600 | 615 | except: |
|
|
601 | 616 | data['data'][0] = tmp[0][:,valid] |
|
|
602 | 617 | data['data'][1] = tmp[0][:,valid] |
|
|
603 | 618 | |
|
|
604 | 619 | if dataOut.mode_op == 'PPI': |
|
|
605 | 620 | self.CODE = 'PPI' |
|
|
606 | 621 | self.title = self.CODE |
|
|
607 | 622 | elif dataOut.mode_op == 'RHI': |
|
|
608 | 623 | self.CODE = 'RHI' |
|
|
609 | 624 | self.title = self.CODE |
|
|
610 | 625 | |
|
|
611 | 626 | data['azi'] = dataOut.data_azi |
|
|
612 | 627 | data['ele'] = dataOut.data_ele |
|
|
613 | 628 | |
|
|
614 | 629 | if isinstance(dataOut.mode_op, bytes): |
|
|
615 | 630 | try: |
|
|
616 | 631 | dataOut.mode_op = dataOut.mode_op.decode() |
|
|
617 | 632 | except: |
|
|
618 | 633 | dataOut.mode_op = str(dataOut.mode_op, 'utf-8') |
|
|
619 | 634 | data['mode_op'] = dataOut.mode_op |
|
|
620 | 635 | self.mode = dataOut.mode_op |
|
|
621 | 636 | |
|
|
622 | 637 | return data, meta |
|
|
623 | 638 | |
|
|
624 | 639 | def plot(self): |
|
|
625 | 640 | data = self.data[-1] |
|
|
626 | 641 | z = data['data'] |
|
|
627 | 642 | r = data['r'] |
|
|
628 | 643 | self.titles = [] |
|
|
629 | 644 | |
|
|
630 | 645 | self.zmax = self.zmax if self.zmax else numpy.nanmax(z) |
|
|
631 | 646 | self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z) |
|
|
632 | 647 | |
|
|
633 | 648 | if isinstance(data['mode_op'], bytes): |
|
|
634 | 649 | data['mode_op'] = data['mode_op'].decode() |
|
|
635 | 650 | |
|
|
636 | 651 | if data['mode_op'] == 'RHI': |
|
|
637 | 652 | r, theta = numpy.meshgrid(r, numpy.radians(data['ele'])) |
|
|
638 | 653 | len_aux = int(data['azi'].shape[0]/4) |
|
|
639 | 654 | mean = numpy.mean(data['azi'][len_aux:-len_aux]) |
|
|
640 | 655 | x, y = r*numpy.cos(theta), r*numpy.sin(theta) |
|
|
641 | 656 | if self.yrange: |
|
|
642 | 657 | self.ylabel= 'Height [km]' |
|
|
643 | 658 | self.xlabel= 'Distance from radar [km]' |
|
|
644 | 659 | self.ymax = self.yrange |
|
|
645 | 660 | self.ymin = 0 |
|
|
646 | 661 | self.xmax = self.xrange if self.xrange else numpy.nanmax(r) |
|
|
647 | 662 | self.xmin = -self.xrange if self.xrange else -numpy.nanmax(r) |
|
|
648 | 663 | self.setrhilimits = False |
|
|
649 | 664 | else: |
|
|
650 | 665 | self.ymin = 0 |
|
|
651 | 666 | self.ymax = numpy.nanmax(r) |
|
|
652 | 667 | self.xmin = -numpy.nanmax(r) |
|
|
653 | 668 | self.xmax = numpy.nanmax(r) |
|
|
654 | 669 | |
|
|
655 | 670 | elif data['mode_op'] == 'PPI': |
|
|
656 | 671 | r, theta = numpy.meshgrid(r, -numpy.radians(data['azi'])+numpy.pi/2) |
|
|
657 | 672 | len_aux = int(data['ele'].shape[0]/4) |
|
|
658 | 673 | mean = numpy.mean(data['ele'][len_aux:-len_aux]) |
|
|
659 | 674 | x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(mean)), r*numpy.sin( |
|
|
660 | 675 | theta)*numpy.cos(numpy.radians(mean)) |
|
|
661 | 676 | x = km2deg(x) + self.longitude |
|
|
662 | 677 | y = km2deg(y) + self.latitude |
|
|
663 | 678 | if self.xrange: |
|
|
664 | 679 | self.ylabel= 'Latitude' |
|
|
665 | 680 | self.xlabel= 'Longitude' |
|
|
666 | 681 | |
|
|
667 | 682 | self.xmin = km2deg(-self.xrange) + self.longitude |
|
|
668 | 683 | self.xmax = km2deg(self.xrange) + self.longitude |
|
|
669 | 684 | |
|
|
670 | 685 | self.ymin = km2deg(-self.xrange) + self.latitude |
|
|
671 | 686 | self.ymax = km2deg(self.xrange) + self.latitude |
|
|
672 | 687 | else: |
|
|
673 | 688 | self.xmin = km2deg(-numpy.nanmax(r)) + self.longitude |
|
|
674 | 689 | self.xmax = km2deg(numpy.nanmax(r)) + self.longitude |
|
|
675 | 690 | |
|
|
676 | 691 | self.ymin = km2deg(-numpy.nanmax(r)) + self.latitude |
|
|
677 | 692 | self.ymax = km2deg(numpy.nanmax(r)) + self.latitude |
|
|
678 | 693 | |
|
|
679 | 694 | self.clear_figures() |
|
|
680 | 695 | |
|
|
681 | 696 | if data['mode_op'] == 'PPI': |
|
|
682 | 697 | axes = self.axes['PPI'] |
|
|
683 | 698 | else: |
|
|
684 | 699 | axes = self.axes['RHI'] |
|
|
685 | 700 | |
|
|
686 | 701 | if self.colormap in cb_tables: |
|
|
687 | 702 | norm = cb_tables[self.colormap]['norm'] |
|
|
688 | 703 | else: |
|
|
689 | 704 | norm = None |
|
|
690 | 705 | |
|
|
691 | 706 | for i, ax in enumerate(axes): |
|
|
692 | 707 | |
|
|
693 | 708 | if norm is None: |
|
|
694 | 709 | ax.plt = ax.pcolormesh(x, y, z[i], cmap=self.colormap, vmin=self.zmin, vmax=self.zmax) |
|
|
695 | 710 | else: |
|
|
696 | 711 | ax.plt = ax.pcolormesh(x, y, z[i], cmap=self.colormap, norm=norm) |
|
|
697 | 712 | |
|
|
698 | 713 | if data['mode_op'] == 'RHI': |
|
|
699 | 714 | len_aux = int(data['azi'].shape[0]/4) |
|
|
700 | 715 | mean = numpy.mean(data['azi'][len_aux:-len_aux]) |
|
|
701 | 716 | if len(self.channels) !=1: |
|
|
702 | 717 | self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels] |
|
|
703 | 718 | else: |
|
|
704 | 719 | self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])] |
|
|
705 | 720 | elif data['mode_op'] == 'PPI': |
|
|
706 | 721 | len_aux = int(data['ele'].shape[0]/4) |
|
|
707 | 722 | mean = numpy.mean(data['ele'][len_aux:-len_aux]) |
|
|
708 | 723 | if len(self.channels) !=1: |
|
|
709 | 724 | self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels] |
|
|
710 | 725 | else: |
|
|
711 | 726 | self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])] |
|
|
712 | 727 | self.mode_value = round(mean,1) |
|
|
713 | 728 | |
|
|
714 | 729 | if data['mode_op'] == 'PPI': |
|
|
715 | 730 | if self.map: |
|
|
716 | 731 | gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, |
|
|
717 | 732 | linewidth=1, color='gray', alpha=0.5, linestyle='--') |
|
|
718 | 733 | gl.xlabel_style = {'size': 8} |
|
|
719 | 734 | gl.ylabel_style = {'size': 8} |
|
|
720 | 735 | gl.xlabels_top = False |
|
|
721 | 736 | gl.ylabels_right = False |
|
|
722 | 737 | shape_d = os.path.join(self.shapes,'Distritos/PER_adm3.shp') |
|
|
723 | 738 | shape_p = os.path.join(self.shapes,'PER_ADM2/PER_ADM2.shp') |
|
|
724 | 739 | capitales = os.path.join(self.shapes,'CAPITALES/cap_distrito.shp') |
|
|
725 | 740 | vias = os.path.join(self.shapes,'Carreteras/VIAS_NACIONAL_250000.shp') |
|
|
726 | 741 | reader_d = shpreader.BasicReader(shape_d, encoding='latin1') |
|
|
727 | 742 | reader_p = shpreader.BasicReader(shape_p, encoding='latin1') |
|
|
728 | 743 | reader_c = shpreader.BasicReader(capitales, encoding='latin1') |
|
|
729 | 744 | reader_v = shpreader.BasicReader(vias, encoding='latin1') |
|
|
730 | 745 | caps = [x for x in reader_c.records() if x.attributes['DEPARTA']=='PIURA' and x.attributes['CATEGORIA']=='CIUDAD'] |
|
|
731 | 746 | districts = [x for x in reader_d.records() if x.attributes['NAME_1']=='Piura'] |
|
|
732 | 747 | provs = [x for x in reader_p.records()] |
|
|
733 | 748 | vias = [x for x in reader_v.records()] |
|
|
734 | 749 | |
|
|
735 | 750 | # Display limits and streets |
|
|
736 | 751 | shape_feature = ShapelyFeature([x.geometry for x in districts], ccrs.PlateCarree(), facecolor="none", edgecolor='grey', lw=0.5) |
|
|
737 | 752 | ax.add_feature(shape_feature) |
|
|
738 | 753 | shape_feature = ShapelyFeature([x.geometry for x in provs], ccrs.PlateCarree(), facecolor="none", edgecolor='white', lw=1) |
|
|
739 | 754 | ax.add_feature(shape_feature) |
|
|
740 | 755 | shape_feature = ShapelyFeature([x.geometry for x in vias], ccrs.PlateCarree(), facecolor="none", edgecolor='yellow', lw=1) |
|
|
741 | 756 | ax.add_feature(shape_feature) |
|
|
742 | 757 | |
|
|
743 | 758 | for cap in caps: |
|
|
744 | 759 | if cap.attributes['NOMBRE'] in ('PIURA', 'SULLANA', 'PAITA', 'SECHURA', 'TALARA'): |
|
|
745 | 760 | ax.text(cap.attributes['X'], cap.attributes['Y'], cap.attributes['NOMBRE'], size=8, color='white', weight='bold') |
|
|
746 | 761 | elif cap.attributes['NOMBRE'] in ('NEGRITOS', 'SAN LUCAS', 'QUERECOTILLO', 'TAMBO GRANDE', 'CHULUCANAS', 'CATACAOS', 'LA UNION'): |
|
|
747 | 762 | ax.text(cap.attributes['X'], cap.attributes['Y'], cap.attributes['NOMBRE'].title(), size=7, color='white') |
|
|
748 | 763 | else: |
|
|
749 | 764 | ax.grid(color='grey', alpha=0.5, linestyle='--', linewidth=1) |
|
|
750 | 765 | |
|
|
751 | 766 | if self.xrange<=10: |
|
|
752 | 767 | ranges = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] |
|
|
753 | 768 | elif self.xrange<=30: |
|
|
754 | 769 | ranges = [5, 10, 15, 20, 25, 30, 35] |
|
|
755 | 770 | elif self.xrange<=60: |
|
|
756 | 771 | ranges = [10, 20, 30, 40, 50, 60] |
|
|
757 | 772 | elif self.xrange<=100: |
|
|
758 | 773 | ranges = [15, 30, 45, 60, 75, 90] |
|
|
759 | 774 | |
|
|
760 | 775 | for R in ranges: |
|
|
761 | 776 | if R <= self.xrange: |
|
|
762 | 777 | circle = Circle((self.longitude, self.latitude), km2deg(R), facecolor='none', |
|
|
763 | 778 | edgecolor='skyblue', linewidth=1, alpha=0.5) |
|
|
764 | 779 | ax.add_patch(circle) |
|
|
765 | 780 | ax.text(km2deg(R)*numpy.cos(numpy.radians(45))+self.longitude, |
|
|
766 | 781 | km2deg(R)*numpy.sin(numpy.radians(45))+self.latitude, |
|
|
767 | 782 | '{}km'.format(R), color='skyblue', size=7) |
|
|
768 | 783 | elif data['mode_op'] == 'RHI': |
|
|
769 | 784 | ax.grid(color='grey', alpha=0.5, linestyle='--', linewidth=1) |
| @@ -1,863 +1,865 | |||
|
|
1 | 1 | ''' |
|
|
2 | 2 | Created on Jul 3, 2014 |
|
|
3 | 3 | |
|
|
4 | 4 | @author: roj-idl71 |
|
|
5 | 5 | ''' |
|
|
6 | 6 | # SUBCHANNELS EN VEZ DE CHANNELS |
|
|
7 | 7 | # BENCHMARKS -> PROBLEMAS CON ARCHIVOS GRANDES -> INCONSTANTE EN EL TIEMPO |
|
|
8 | 8 | # ACTUALIZACION DE VERSION |
|
|
9 | 9 | # HEADERS |
|
|
10 | 10 | # MODULO DE ESCRITURA |
|
|
11 | 11 | # METADATA |
|
|
12 | 12 | |
|
|
13 | 13 | import os |
|
|
14 | 14 | import time |
|
|
15 | 15 | import datetime |
|
|
16 | 16 | import numpy |
|
|
17 | 17 | import timeit |
|
|
18 | 18 | from fractions import Fraction |
|
|
19 | 19 | from time import time |
|
|
20 | 20 | from time import sleep |
|
|
21 | 21 | |
|
|
22 | 22 | import schainpy.admin |
|
|
23 | 23 | from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader |
|
|
24 | 24 | from schainpy.model.data.jrodata import Voltage |
|
|
25 | 25 | from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator |
|
|
26 | 26 | |
|
|
27 | 27 | import pickle |
|
|
28 | 28 | try: |
|
|
29 | 29 | os.environ["HDF5_USE_FILE_LOCKING"] = "FALSE" |
|
|
30 | 30 | import digital_rf |
|
|
31 | 31 | except: |
|
|
32 | 32 | pass |
|
|
33 | 33 | |
|
|
34 | 34 | |
|
|
35 | 35 | class DigitalRFReader(ProcessingUnit): |
|
|
36 | 36 | ''' |
|
|
37 | 37 | classdocs |
|
|
38 | 38 | ''' |
|
|
39 | 39 | |
|
|
40 | 40 | def __init__(self): |
|
|
41 | 41 | ''' |
|
|
42 | 42 | Constructor |
|
|
43 | 43 | ''' |
|
|
44 | 44 | |
|
|
45 | 45 | ProcessingUnit.__init__(self) |
|
|
46 | 46 | |
|
|
47 | 47 | self.dataOut = Voltage() |
|
|
48 | 48 | self.__printInfo = True |
|
|
49 | 49 | self.__flagDiscontinuousBlock = False |
|
|
50 | 50 | self.__bufferIndex = 9999999 |
|
|
51 | 51 | self.__codeType = 0 |
|
|
52 | 52 | self.__ippKm = None |
|
|
53 | 53 | self.__nCode = None |
|
|
54 | 54 | self.__nBaud = None |
|
|
55 | 55 | self.__code = None |
|
|
56 | 56 | self.dtype = None |
|
|
57 | 57 | self.oldAverage = None |
|
|
58 | 58 | self.path = None |
|
|
59 | 59 | |
|
|
60 | 60 | def close(self): |
|
|
61 | 61 | print('Average of writing to digital rf format is ', self.oldAverage * 1000) |
|
|
62 | 62 | return |
|
|
63 | 63 | |
|
|
64 | 64 | def __getCurrentSecond(self): |
|
|
65 | 65 | |
|
|
66 | 66 | return self.__thisUnixSample / self.__sample_rate |
|
|
67 | 67 | |
|
|
68 | 68 | thisSecond = property(__getCurrentSecond, "I'm the 'thisSecond' property.") |
|
|
69 | 69 | |
|
|
70 | 70 | def __setFileHeader(self): |
|
|
71 | 71 | ''' |
|
|
72 | 72 | In this method will be initialized every parameter of dataOut object (header, no data) |
|
|
73 | 73 | ''' |
|
|
74 | 74 | ippSeconds = 1.0 * self.__nSamples / self.__sample_rate |
|
|
75 | 75 | if not self.getByBlock: |
|
|
76 | 76 | nProfiles = 1.0 / ippSeconds # Number of profiles in one second |
|
|
77 | 77 | else: |
|
|
78 | 78 | nProfiles = self.nProfileBlocks # Number of profiles in one block |
|
|
79 | 79 | |
|
|
80 | 80 | try: |
|
|
81 | 81 | self.dataOut.radarControllerHeaderObj = RadarControllerHeader( |
|
|
82 | 82 | self.__radarControllerHeader) |
|
|
83 | 83 | except: |
|
|
84 | 84 | self.dataOut.radarControllerHeaderObj = RadarControllerHeader( |
|
|
85 | 85 | txA=0, |
|
|
86 | 86 | txB=0, |
|
|
87 | 87 | nWindows=1, |
|
|
88 | 88 | nHeights=self.__nSamples, |
|
|
89 | 89 | firstHeight=self.__firstHeigth, |
|
|
90 | 90 | deltaHeight=self.__deltaHeigth, |
|
|
91 | 91 | codeType=self.__codeType, |
|
|
92 | 92 | nCode=self.__nCode, nBaud=self.__nBaud, |
|
|
93 | 93 | code=self.__code) |
|
|
94 | 94 | |
|
|
95 | 95 | try: |
|
|
96 | 96 | self.dataOut.systemHeaderObj = SystemHeader(self.__systemHeader) |
|
|
97 | 97 | except: |
|
|
98 | 98 | self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples, |
|
|
99 | 99 | nProfiles=nProfiles, |
|
|
100 | 100 | nChannels=len( |
|
|
101 | 101 | self.__channelList), |
|
|
102 | 102 | adcResolution=14) |
|
|
103 | 103 | self.dataOut.type = "Voltage" |
|
|
104 | 104 | |
|
|
105 | 105 | self.dataOut.data = None |
|
|
106 | 106 | |
|
|
107 | 107 | self.dataOut.dtype = self.dtype |
|
|
108 | 108 | |
|
|
109 | 109 | # self.dataOut.nChannels = 0 |
|
|
110 | 110 | |
|
|
111 | 111 | # self.dataOut.nHeights = 0 |
|
|
112 | 112 | |
|
|
113 | 113 | self.dataOut.nProfiles = int(nProfiles) |
|
|
114 | 114 | |
|
|
115 | 115 | self.dataOut.heightList = self.__firstHeigth + \ |
|
|
116 | 116 | numpy.arange(self.__nSamples, dtype=numpy.float_) * \ |
|
|
117 | 117 | self.__deltaHeigth |
|
|
118 | 118 | |
|
|
119 | 119 | #self.dataOut.channelList = list(range(self.__num_subchannels)) |
|
|
120 | 120 | self.dataOut.channelList = list(range(len(self.__channelList))) |
|
|
121 | 121 | if not self.getByBlock: |
|
|
122 | 122 | |
|
|
123 | 123 | self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights |
|
|
124 | 124 | else: |
|
|
125 | 125 | self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights*self.nProfileBlocks |
|
|
126 | 126 | |
|
|
127 | 127 | # self.dataOut.channelIndexList = None |
|
|
128 | 128 | |
|
|
129 | 129 | self.dataOut.flagNoData = True |
|
|
130 | 130 | if not self.getByBlock: |
|
|
131 | 131 | self.dataOut.flagDataAsBlock = False |
|
|
132 | 132 | else: |
|
|
133 | 133 | self.dataOut.flagDataAsBlock = True |
|
|
134 | 134 | # Set to TRUE if the data is discontinuous |
|
|
135 | 135 | self.dataOut.flagDiscontinuousBlock = False |
|
|
136 | 136 | |
|
|
137 | 137 | self.dataOut.utctime = None |
|
|
138 | 138 | |
|
|
139 | 139 | # timezone like jroheader, difference in minutes between UTC and localtime |
|
|
140 | 140 | self.dataOut.timeZone = self.__timezone / 60 |
|
|
141 | 141 | |
|
|
142 | 142 | self.dataOut.dstFlag = 0 |
|
|
143 | 143 | |
|
|
144 | 144 | self.dataOut.errorCount = 0 |
|
|
145 | 145 | |
|
|
146 | 146 | try: |
|
|
147 | 147 | self.dataOut.nCohInt = self.fixed_metadata_dict.get( |
|
|
148 | 148 | 'nCohInt', self.nCohInt) |
|
|
149 | 149 | |
|
|
150 | 150 | # asumo que la data esta decodificada |
|
|
151 | 151 | self.dataOut.flagDecodeData = self.fixed_metadata_dict.get( |
|
|
152 | 152 | 'flagDecodeData', self.flagDecodeData) |
|
|
153 | 153 | |
|
|
154 | 154 | # asumo que la data esta sin flip |
|
|
155 | 155 | self.dataOut.flagDeflipData = self.fixed_metadata_dict['flagDeflipData'] |
|
|
156 | 156 | |
|
|
157 | 157 | self.dataOut.flagShiftFFT = self.fixed_metadata_dict['flagShiftFFT'] |
|
|
158 | 158 | |
|
|
159 | 159 | self.dataOut.useLocalTime = self.fixed_metadata_dict['useLocalTime'] |
|
|
160 | 160 | except: |
|
|
161 | 161 | pass |
|
|
162 | 162 | |
|
|
163 | 163 | self.dataOut.ippSeconds = ippSeconds |
|
|
164 | 164 | |
|
|
165 | 165 | # Time interval between profiles |
|
|
166 | 166 | # self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt |
|
|
167 | 167 | |
|
|
168 | 168 | self.dataOut.frequency = self.__frequency |
|
|
169 | 169 | |
|
|
170 | 170 | self.dataOut.realtime = self.__online |
|
|
171 | 171 | |
|
|
172 | 172 | def findDatafiles(self, path, startDate=None, endDate=None): |
|
|
173 | 173 | |
|
|
174 | 174 | if not os.path.isdir(path): |
|
|
175 | 175 | return [] |
|
|
176 | 176 | |
|
|
177 | 177 | try: |
|
|
178 | 178 | digitalReadObj = digital_rf.DigitalRFReader( |
|
|
179 | 179 | path, load_all_metadata=True) |
|
|
180 | 180 | except: |
|
|
181 | 181 | digitalReadObj = digital_rf.DigitalRFReader(path) |
|
|
182 | 182 | |
|
|
183 | 183 | channelNameList = digitalReadObj.get_channels() |
|
|
184 | 184 | |
|
|
185 | 185 | if not channelNameList: |
|
|
186 | 186 | return [] |
|
|
187 | 187 | |
|
|
188 | 188 | metadata_dict = digitalReadObj.get_rf_file_metadata(channelNameList[0]) |
|
|
189 | 189 | |
|
|
190 | 190 | sample_rate = metadata_dict['sample_rate'][0] |
|
|
191 | 191 | |
|
|
192 | 192 | this_metadata_file = digitalReadObj.get_metadata(channelNameList[0]) |
|
|
193 | 193 | |
|
|
194 | 194 | try: |
|
|
195 | 195 | timezone = this_metadata_file['timezone'].value |
|
|
196 | 196 | except: |
|
|
197 | 197 | timezone = 0 |
|
|
198 | 198 | |
|
|
199 | 199 | startUTCSecond, endUTCSecond = digitalReadObj.get_bounds( |
|
|
200 | 200 | channelNameList[0]) / sample_rate - timezone |
|
|
201 | 201 | |
|
|
202 | 202 | startDatetime = datetime.datetime.utcfromtimestamp(startUTCSecond) |
|
|
203 | 203 | endDatatime = datetime.datetime.utcfromtimestamp(endUTCSecond) |
|
|
204 | 204 | |
|
|
205 | 205 | if not startDate: |
|
|
206 | 206 | startDate = startDatetime.date() |
|
|
207 | 207 | |
|
|
208 | 208 | if not endDate: |
|
|
209 | 209 | endDate = endDatatime.date() |
|
|
210 | 210 | |
|
|
211 | 211 | dateList = [] |
|
|
212 | 212 | |
|
|
213 | 213 | thisDatetime = startDatetime |
|
|
214 | 214 | |
|
|
215 | 215 | while(thisDatetime <= endDatatime): |
|
|
216 | 216 | |
|
|
217 | 217 | thisDate = thisDatetime.date() |
|
|
218 | 218 | |
|
|
219 | 219 | if thisDate < startDate: |
|
|
220 | 220 | continue |
|
|
221 | 221 | |
|
|
222 | 222 | if thisDate > endDate: |
|
|
223 | 223 | break |
|
|
224 | 224 | |
|
|
225 | 225 | dateList.append(thisDate) |
|
|
226 | 226 | thisDatetime += datetime.timedelta(1) |
|
|
227 | 227 | |
|
|
228 | 228 | return dateList |
|
|
229 | 229 | |
|
|
230 | 230 | def setup(self, path=None, |
|
|
231 | 231 | startDate=None, |
|
|
232 | 232 | endDate=None, |
|
|
233 | 233 | startTime=datetime.time(0, 0, 0), |
|
|
234 | 234 | endTime=datetime.time(23, 59, 59), |
|
|
235 | 235 | channelList=None, |
|
|
236 | 236 | nSamples=None, |
|
|
237 | 237 | online=False, |
|
|
238 | 238 | delay=60, |
|
|
239 | 239 | buffer_size=1024, |
|
|
240 | 240 | ippKm=None, |
|
|
241 | 241 | nCohInt=1, |
|
|
242 | 242 | nCode=1, |
|
|
243 | 243 | nBaud=1, |
|
|
244 | 244 | flagDecodeData=False, |
|
|
245 | 245 | code=numpy.ones((1, 1), dtype=int), |
|
|
246 | 246 | getByBlock=0, |
|
|
247 | 247 | nProfileBlocks=1, |
|
|
248 | 248 | **kwargs): |
|
|
249 | 249 | ''' |
|
|
250 | 250 | In this method we should set all initial parameters. |
|
|
251 | 251 | |
|
|
252 | 252 | Inputs: |
|
|
253 | 253 | path |
|
|
254 | 254 | startDate |
|
|
255 | 255 | endDate |
|
|
256 | 256 | startTime |
|
|
257 | 257 | endTime |
|
|
258 | 258 | set |
|
|
259 | 259 | expLabel |
|
|
260 | 260 | ext |
|
|
261 | 261 | online |
|
|
262 | 262 | delay |
|
|
263 | 263 | ''' |
|
|
264 | 264 | self.path = path |
|
|
265 | 265 | self.nCohInt = nCohInt |
|
|
266 | 266 | self.flagDecodeData = flagDecodeData |
|
|
267 | 267 | self.i = 0 |
|
|
268 | 268 | |
|
|
269 | 269 | self.getByBlock = getByBlock |
|
|
270 | 270 | self.nProfileBlocks = nProfileBlocks |
|
|
271 | 271 | if online: |
|
|
272 | 272 | print('Waiting for RF data..') |
|
|
273 | 273 | sleep(40) |
|
|
274 | 274 | |
|
|
275 | 275 | if not os.path.isdir(path): |
|
|
276 | 276 | raise ValueError("[Reading] Directory %s does not exist" % path) |
|
|
277 | 277 | |
|
|
278 | 278 | #print("path",path) |
|
|
279 | 279 | try: |
|
|
280 | 280 | self.digitalReadObj = digital_rf.DigitalRFReader( |
|
|
281 | 281 | path, load_all_metadata=True) |
|
|
282 | 282 | except: |
|
|
283 | 283 | self.digitalReadObj = digital_rf.DigitalRFReader(path) |
|
|
284 | 284 | |
|
|
285 | 285 | channelNameList = self.digitalReadObj.get_channels() |
|
|
286 | 286 | |
|
|
287 | 287 | if not channelNameList: |
|
|
288 | 288 | raise ValueError("[Reading] Directory %s does not have any files" % path) |
|
|
289 | 289 | |
|
|
290 | 290 | if not channelList: |
|
|
291 | 291 | channelList = list(range(len(channelNameList))) |
|
|
292 | 292 | |
|
|
293 | 293 | ########## Reading metadata ###################### |
|
|
294 | 294 | |
|
|
295 | 295 | top_properties = self.digitalReadObj.get_properties( |
|
|
296 | 296 | channelNameList[channelList[0]]) |
|
|
297 | 297 | |
|
|
298 | 298 | self.__num_subchannels = top_properties['num_subchannels'] |
|
|
299 | 299 | self.__sample_rate = 1.0 * \ |
|
|
300 | 300 | top_properties['sample_rate_numerator'] / \ |
|
|
301 | 301 | top_properties['sample_rate_denominator'] |
|
|
302 | 302 | # self.__samples_per_file = top_properties['samples_per_file'][0] |
|
|
303 | 303 | self.__deltaHeigth = 1e6 * 0.15 / self.__sample_rate # why 0.15? |
|
|
304 | 304 | |
|
|
305 | 305 | this_metadata_file = self.digitalReadObj.get_digital_metadata( |
|
|
306 | 306 | channelNameList[channelList[0]]) |
|
|
307 | 307 | metadata_bounds = this_metadata_file.get_bounds() |
|
|
308 | 308 | self.fixed_metadata_dict = this_metadata_file.read( |
|
|
309 | 309 | metadata_bounds[0])[metadata_bounds[0]] # GET FIRST HEADER |
|
|
310 | 310 | |
|
|
311 | 311 | try: |
|
|
312 | 312 | self.__processingHeader = self.fixed_metadata_dict['processingHeader'] |
|
|
313 | 313 | self.__radarControllerHeader = self.fixed_metadata_dict['radarControllerHeader'] |
|
|
314 | 314 | self.__systemHeader = self.fixed_metadata_dict['systemHeader'] |
|
|
315 | 315 | self.dtype = pickle.loads(self.fixed_metadata_dict['dtype']) |
|
|
316 | 316 | except: |
|
|
317 | 317 | pass |
|
|
318 | 318 | |
|
|
319 | 319 | self.__frequency = None |
|
|
320 | 320 | |
|
|
321 | 321 | self.__frequency = self.fixed_metadata_dict.get('frequency', 1) |
|
|
322 | 322 | |
|
|
323 | self.__frequency = 9.345e9 | |
|
|
324 | ||
|
|
323 | 325 | self.__timezone = self.fixed_metadata_dict.get('timezone', 18000) |
|
|
324 | 326 | |
|
|
325 | 327 | try: |
|
|
326 | 328 | nSamples = self.fixed_metadata_dict['nSamples'] |
|
|
327 | 329 | except: |
|
|
328 | 330 | nSamples = None |
|
|
329 | 331 | |
|
|
330 | 332 | self.__firstHeigth = 0 |
|
|
331 | 333 | |
|
|
332 | 334 | try: |
|
|
333 | 335 | codeType = self.__radarControllerHeader['codeType'] |
|
|
334 | 336 | except: |
|
|
335 | 337 | codeType = 0 |
|
|
336 | 338 | |
|
|
337 | 339 | try: |
|
|
338 | 340 | if codeType: |
|
|
339 | 341 | nCode = self.__radarControllerHeader['nCode'] |
|
|
340 | 342 | nBaud = self.__radarControllerHeader['nBaud'] |
|
|
341 | 343 | code = self.__radarControllerHeader['code'] |
|
|
342 | 344 | except: |
|
|
343 | 345 | pass |
|
|
344 | 346 | |
|
|
345 | 347 | if not ippKm: |
|
|
346 | 348 | try: |
|
|
347 | 349 | # seconds to km |
|
|
348 | 350 | ippKm = self.__radarControllerHeader['ipp'] |
|
|
349 | 351 | except: |
|
|
350 | 352 | ippKm = None |
|
|
351 | 353 | #################################################### |
|
|
352 | 354 | self.__ippKm = ippKm |
|
|
353 | 355 | startUTCSecond = None |
|
|
354 | 356 | endUTCSecond = None |
|
|
355 | 357 | |
|
|
356 | 358 | if startDate: |
|
|
357 | 359 | startDatetime = datetime.datetime.combine(startDate, startTime) |
|
|
358 | 360 | startUTCSecond = ( |
|
|
359 | 361 | startDatetime - datetime.datetime(1970, 1, 1)).total_seconds()# + self.__timezone |
|
|
360 | 362 | |
|
|
361 | 363 | if endDate: |
|
|
362 | 364 | endDatetime = datetime.datetime.combine(endDate, endTime) |
|
|
363 | 365 | endUTCSecond = (endDatetime - datetime.datetime(1970, |
|
|
364 | 366 | 1, 1)).total_seconds()# + self.__timezone |
|
|
365 | 367 | start_index, end_index = self.digitalReadObj.get_bounds(channelNameList[channelList[0]]) |
|
|
366 | 368 | if start_index==None or end_index==None: |
|
|
367 | 369 | print("Check error No data, start_index: ",start_index,",end_index: ",end_index) |
|
|
368 | 370 | #return 0 |
|
|
369 | 371 | if not startUTCSecond: |
|
|
370 | 372 | startUTCSecond = start_index / self.__sample_rate |
|
|
371 | 373 | if start_index > startUTCSecond * self.__sample_rate: |
|
|
372 | 374 | startUTCSecond = start_index / self.__sample_rate |
|
|
373 | 375 | |
|
|
374 | 376 | if not endUTCSecond: |
|
|
375 | 377 | endUTCSecond = end_index / self.__sample_rate |
|
|
376 | 378 | |
|
|
377 | 379 | if end_index < endUTCSecond * self.__sample_rate: |
|
|
378 | 380 | endUTCSecond = end_index / self.__sample_rate #Check UTC and LT time |
|
|
379 | 381 | |
|
|
380 | 382 | if not nSamples: |
|
|
381 | 383 | if not ippKm: |
|
|
382 | 384 | raise ValueError("[Reading] nSamples or ippKm should be defined") |
|
|
383 | 385 | nSamples = int(ippKm / (1e6 * 0.15 / self.__sample_rate)) |
|
|
384 | 386 | |
|
|
385 | 387 | channelBoundList = [] |
|
|
386 | 388 | channelNameListFiltered = [] |
|
|
387 | 389 | |
|
|
388 | 390 | for thisIndexChannel in channelList: |
|
|
389 | 391 | thisChannelName = channelNameList[thisIndexChannel] |
|
|
390 | 392 | start_index, end_index = self.digitalReadObj.get_bounds( |
|
|
391 | 393 | thisChannelName) |
|
|
392 | 394 | channelBoundList.append((start_index, end_index)) |
|
|
393 | 395 | channelNameListFiltered.append(thisChannelName) |
|
|
394 | 396 | |
|
|
395 | 397 | self.profileIndex = 0 |
|
|
396 | 398 | self.i = 0 |
|
|
397 | 399 | self.__delay = delay |
|
|
398 | 400 | |
|
|
399 | 401 | self.__codeType = codeType |
|
|
400 | 402 | self.__nCode = nCode |
|
|
401 | 403 | self.__nBaud = nBaud |
|
|
402 | 404 | self.__code = code |
|
|
403 | 405 | |
|
|
404 | 406 | self.__datapath = path |
|
|
405 | 407 | self.__online = online |
|
|
406 | 408 | self.__channelList = channelList |
|
|
407 | 409 | self.__channelNameList = channelNameListFiltered |
|
|
408 | 410 | self.__channelBoundList = channelBoundList |
|
|
409 | 411 | self.__nSamples = nSamples |
|
|
410 | 412 | if self.getByBlock: |
|
|
411 | 413 | nSamples = nSamples*nProfileBlocks |
|
|
412 | 414 | |
|
|
413 | 415 | |
|
|
414 | 416 | self.__samples_to_read = int(nSamples) # FIJO: AHORA 40 |
|
|
415 | 417 | self.__nChannels = len(self.__channelList) |
|
|
416 | 418 | #print("------------------------------------------") |
|
|
417 | 419 | #print("self.__samples_to_read",self.__samples_to_read) |
|
|
418 | 420 | #print("self.__nSamples",self.__nSamples) |
|
|
419 | 421 | # son iguales y el buffer_index da 0 |
|
|
420 | 422 | self.__startUTCSecond = startUTCSecond |
|
|
421 | 423 | self.__endUTCSecond = endUTCSecond |
|
|
422 | 424 | |
|
|
423 | 425 | self.__timeInterval = 1.0 * self.__samples_to_read / \ |
|
|
424 | 426 | self.__sample_rate # Time interval |
|
|
425 | 427 | |
|
|
426 | 428 | if online: |
|
|
427 | 429 | # self.__thisUnixSample = int(endUTCSecond*self.__sample_rate - 4*self.__samples_to_read) |
|
|
428 | 430 | startUTCSecond = numpy.floor(endUTCSecond) |
|
|
429 | 431 | |
|
|
430 | 432 | # por que en el otro metodo lo primero q se hace es sumar samplestoread |
|
|
431 | 433 | self.__thisUnixSample = int(startUTCSecond * self.__sample_rate) - self.__samples_to_read |
|
|
432 | 434 | |
|
|
433 | 435 | #self.__data_buffer = numpy.zeros( |
|
|
434 | 436 | # (self.__num_subchannels, self.__samples_to_read), dtype=numpy.complex) |
|
|
435 | 437 | print("samplestoread",self.__samples_to_read) |
|
|
436 | 438 | self.__data_buffer = numpy.zeros((int(len(channelList)), self.__samples_to_read), dtype=numpy.complex_) |
|
|
437 | 439 | |
|
|
438 | 440 | |
|
|
439 | 441 | self.__setFileHeader() |
|
|
440 | 442 | self.isConfig = True |
|
|
441 | 443 | |
|
|
442 | 444 | print("[Reading] Digital RF Data was found from %s to %s " % ( |
|
|
443 | 445 | datetime.datetime.utcfromtimestamp( |
|
|
444 | 446 | self.__startUTCSecond - self.__timezone), |
|
|
445 | 447 | datetime.datetime.utcfromtimestamp( |
|
|
446 | 448 | self.__endUTCSecond - self.__timezone) |
|
|
447 | 449 | )) |
|
|
448 | 450 | |
|
|
449 | 451 | print("[Reading] Starting process from %s to %s" % (datetime.datetime.utcfromtimestamp(startUTCSecond - self.__timezone), |
|
|
450 | 452 | datetime.datetime.utcfromtimestamp(endUTCSecond - self.__timezone))) |
|
|
451 | 453 | self.oldAverage = None |
|
|
452 | 454 | self.count = 0 |
|
|
453 | 455 | self.executionTime = 0 |
|
|
454 | 456 | |
|
|
455 | 457 | def __reload(self): |
|
|
456 | 458 | |
|
|
457 | 459 | # print "%s not in range [%s, %s]" %( |
|
|
458 | 460 | # datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), |
|
|
459 | 461 | # datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone), |
|
|
460 | 462 | # datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone) |
|
|
461 | 463 | # ) |
|
|
462 | 464 | print("[Reading] reloading metadata ...") |
|
|
463 | 465 | |
|
|
464 | 466 | try: |
|
|
465 | 467 | self.digitalReadObj.reload(complete_update=True) |
|
|
466 | 468 | except: |
|
|
467 | 469 | self.digitalReadObj = digital_rf.DigitalRFReader(self.path) |
|
|
468 | 470 | |
|
|
469 | 471 | start_index, end_index = self.digitalReadObj.get_bounds( |
|
|
470 | 472 | self.__channelNameList[self.__channelList[0]]) |
|
|
471 | 473 | |
|
|
472 | 474 | if start_index > self.__startUTCSecond * self.__sample_rate: |
|
|
473 | 475 | self.__startUTCSecond = 1.0 * start_index / self.__sample_rate |
|
|
474 | 476 | |
|
|
475 | 477 | if end_index > self.__endUTCSecond * self.__sample_rate: |
|
|
476 | 478 | self.__endUTCSecond = 1.0 * end_index / self.__sample_rate |
|
|
477 | 479 | print() |
|
|
478 | 480 | print("[Reading] New timerange found [%s, %s] " % ( |
|
|
479 | 481 | datetime.datetime.utcfromtimestamp( |
|
|
480 | 482 | self.__startUTCSecond - self.__timezone), |
|
|
481 | 483 | datetime.datetime.utcfromtimestamp( |
|
|
482 | 484 | self.__endUTCSecond - self.__timezone) |
|
|
483 | 485 | )) |
|
|
484 | 486 | |
|
|
485 | 487 | return True |
|
|
486 | 488 | |
|
|
487 | 489 | return False |
|
|
488 | 490 | |
|
|
489 | 491 | def timeit(self, toExecute): |
|
|
490 | 492 | t0 = time.time() |
|
|
491 | 493 | toExecute() |
|
|
492 | 494 | self.executionTime = time.time() - t0 |
|
|
493 | 495 | if self.oldAverage is None: |
|
|
494 | 496 | self.oldAverage = self.executionTime |
|
|
495 | 497 | self.oldAverage = (self.executionTime + self.count * |
|
|
496 | 498 | self.oldAverage) / (self.count + 1.0) |
|
|
497 | 499 | self.count = self.count + 1.0 |
|
|
498 | 500 | return |
|
|
499 | 501 | |
|
|
500 | 502 | def __readNextBlock(self, seconds=30, volt_scale=1/20000.0): |
|
|
501 | 503 | ''' |
|
|
502 | 504 | NOTA: APLICACION RADAR METEOROLOGICO |
|
|
503 | 505 | VALORES OBTENIDOS CON LA USRP, volt_scale = 1,conexion directa al Ch Rx. |
|
|
504 | 506 | |
|
|
505 | 507 | MAXIMO |
|
|
506 | 508 | 9886 -> 0.980 Voltiospp |
|
|
507 | 509 | 4939 -> 0.480 Voltiospp |
|
|
508 | 510 | 14825 -> 1.440 Voltiospp |
|
|
509 | 511 | 18129 -> 1.940 Voltiospp |
|
|
510 | 512 | Para llevar al valor correspondiente de Voltaje, debemos dividir por 20000 |
|
|
511 | 513 | y obtenemos la Amplitud correspondiente de entrada IQ. |
|
|
512 | 514 | volt_scale = (1/20000.0) |
|
|
513 | 515 | ''' |
|
|
514 | 516 | # Set the next data |
|
|
515 | 517 | self.__flagDiscontinuousBlock = False |
|
|
516 | 518 | self.__thisUnixSample += self.__samples_to_read |
|
|
517 | 519 | |
|
|
518 | 520 | if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate: |
|
|
519 | 521 | print ("[Reading] There are no more data into selected time-range") |
|
|
520 | 522 | if self.__online: |
|
|
521 | 523 | sleep(3) |
|
|
522 | 524 | self.__reload() |
|
|
523 | 525 | else: |
|
|
524 | 526 | return False |
|
|
525 | 527 | |
|
|
526 | 528 | if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate: |
|
|
527 | 529 | return False |
|
|
528 | 530 | self.__thisUnixSample -= self.__samples_to_read |
|
|
529 | 531 | |
|
|
530 | 532 | indexChannel = 0 |
|
|
531 | 533 | |
|
|
532 | 534 | dataOk = False |
|
|
533 | 535 | |
|
|
534 | 536 | for thisChannelName in self.__channelNameList: # TODO VARIOS CHANNELS? |
|
|
535 | 537 | for indexSubchannel in range(self.__num_subchannels): |
|
|
536 | 538 | try: |
|
|
537 | 539 | t0 = time() |
|
|
538 | 540 | #print("thisUNixSample",self.__thisUnixSample) |
|
|
539 | 541 | result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample, |
|
|
540 | 542 | self.__samples_to_read, |
|
|
541 | 543 | thisChannelName, sub_channel=indexSubchannel) |
|
|
542 | 544 | #print("result--------------",result) |
|
|
543 | 545 | self.executionTime = time() - t0 |
|
|
544 | 546 | if self.oldAverage is None: |
|
|
545 | 547 | self.oldAverage = self.executionTime |
|
|
546 | 548 | self.oldAverage = ( |
|
|
547 | 549 | self.executionTime + self.count * self.oldAverage) / (self.count + 1.0) |
|
|
548 | 550 | self.count = self.count + 1.0 |
|
|
549 | 551 | |
|
|
550 | 552 | except IOError as e: |
|
|
551 | 553 | # read next profile |
|
|
552 | 554 | self.__flagDiscontinuousBlock = True |
|
|
553 | 555 | print("[Reading] %s" % datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), e) |
|
|
554 | 556 | bot = 0 |
|
|
555 | 557 | while(self.__flagDiscontinuousBlock): |
|
|
556 | 558 | bot +=1 |
|
|
557 | 559 | self.__thisUnixSample += self.__samples_to_read |
|
|
558 | 560 | try: |
|
|
559 | 561 | result = result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample,self.__samples_to_read,thisChannelName, sub_channel=indexSubchannel) |
|
|
560 | 562 | self.__flagDiscontinuousBlock=False |
|
|
561 | 563 | print("Searching.. N°: ",bot,"Success",self.__thisUnixSample) |
|
|
562 | 564 | except: |
|
|
563 | 565 | print("Searching...N°: ",bot,"Fail", self.__thisUnixSample) |
|
|
564 | 566 | if self.__flagDiscontinuousBlock==True: |
|
|
565 | 567 | break |
|
|
566 | 568 | else: |
|
|
567 | 569 | print("New data index found...",self.__thisUnixSample) |
|
|
568 | 570 | #break |
|
|
569 | 571 | |
|
|
570 | 572 | if result.shape[0] != self.__samples_to_read: |
|
|
571 | 573 | self.__flagDiscontinuousBlock = True |
|
|
572 | 574 | print("[Reading] %s: Too few samples were found, just %d/%d samples" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), |
|
|
573 | 575 | result.shape[0], |
|
|
574 | 576 | self.__samples_to_read)) |
|
|
575 | 577 | break |
|
|
576 | 578 | |
|
|
577 | 579 | self.__data_buffer[indexChannel, :] = result * volt_scale |
|
|
578 | 580 | indexChannel+=1 |
|
|
579 | 581 | |
|
|
580 | 582 | dataOk = True |
|
|
581 | 583 | |
|
|
582 | 584 | self.__utctime = self.__thisUnixSample / self.__sample_rate |
|
|
583 | 585 | |
|
|
584 | 586 | if not dataOk: |
|
|
585 | 587 | return False |
|
|
586 | 588 | |
|
|
587 | 589 | print("[Reading] %s: %d samples <> %f sec" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), |
|
|
588 | 590 | self.__samples_to_read, |
|
|
589 | 591 | self.__timeInterval)) |
|
|
590 | 592 | |
|
|
591 | 593 | self.__bufferIndex = 0 |
|
|
592 | 594 | |
|
|
593 | 595 | return True |
|
|
594 | 596 | |
|
|
595 | 597 | def __isBufferEmpty(self): |
|
|
596 | 598 | |
|
|
597 | 599 | return self.__bufferIndex > self.__samples_to_read - self.__nSamples # 40960 - 40 |
|
|
598 | 600 | |
|
|
599 | 601 | def getData(self, seconds=30, nTries=5): |
|
|
600 | 602 | ''' |
|
|
601 | 603 | This method gets the data from files and put the data into the dataOut object |
|
|
602 | 604 | |
|
|
603 | 605 | In addition, increase el the buffer counter in one. |
|
|
604 | 606 | |
|
|
605 | 607 | Return: |
|
|
606 | 608 | data : retorna un perfil de voltages (alturas * canales) copiados desde el |
|
|
607 | 609 | buffer. Si no hay mas archivos a leer retorna None. |
|
|
608 | 610 | |
|
|
609 | 611 | Affected: |
|
|
610 | 612 | self.dataOut |
|
|
611 | 613 | self.profileIndex |
|
|
612 | 614 | self.flagDiscontinuousBlock |
|
|
613 | 615 | self.flagIsNewBlock |
|
|
614 | 616 | ''' |
|
|
615 | 617 | #print("getdata") |
|
|
616 | 618 | err_counter = 0 |
|
|
617 | 619 | self.dataOut.flagNoData = True |
|
|
618 | 620 | |
|
|
619 | 621 | |
|
|
620 | 622 | if self.__isBufferEmpty(): |
|
|
621 | 623 | #print("hi") |
|
|
622 | 624 | self.__flagDiscontinuousBlock = False |
|
|
623 | 625 | |
|
|
624 | 626 | while True: |
|
|
625 | 627 | if self.__readNextBlock(): |
|
|
626 | 628 | break |
|
|
627 | 629 | if self.__thisUnixSample > self.__endUTCSecond * self.__sample_rate: |
|
|
628 | 630 | raise schainpy.admin.SchainError('Error') |
|
|
629 | 631 | return |
|
|
630 | 632 | |
|
|
631 | 633 | if self.__flagDiscontinuousBlock: |
|
|
632 | 634 | raise schainpy.admin.SchainError('discontinuous block found') |
|
|
633 | 635 | return |
|
|
634 | 636 | |
|
|
635 | 637 | if not self.__online: |
|
|
636 | 638 | raise schainpy.admin.SchainError('Online?') |
|
|
637 | 639 | return |
|
|
638 | 640 | |
|
|
639 | 641 | err_counter += 1 |
|
|
640 | 642 | if err_counter > nTries: |
|
|
641 | 643 | raise schainpy.admin.SchainError('Max retrys reach') |
|
|
642 | 644 | return |
|
|
643 | 645 | |
|
|
644 | 646 | print('[Reading] waiting %d seconds to read a new block' % seconds) |
|
|
645 | 647 | sleep(seconds) |
|
|
646 | 648 | |
|
|
647 | 649 | |
|
|
648 | 650 | if not self.getByBlock: |
|
|
649 | 651 | |
|
|
650 | 652 | #print("self.__bufferIndex",self.__bufferIndex)# este valor siempre es cero aparentemente |
|
|
651 | 653 | self.dataOut.data = self.__data_buffer[:, self.__bufferIndex:self.__bufferIndex + self.__nSamples] |
|
|
652 | 654 | self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate |
|
|
653 | 655 | self.dataOut.flagNoData = False |
|
|
654 | 656 | self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock |
|
|
655 | 657 | self.dataOut.profileIndex = self.profileIndex |
|
|
656 | 658 | |
|
|
657 | 659 | self.__bufferIndex += self.__nSamples |
|
|
658 | 660 | self.profileIndex += 1 |
|
|
659 | 661 | |
|
|
660 | 662 | if self.profileIndex == self.dataOut.nProfiles: |
|
|
661 | 663 | self.profileIndex = 0 |
|
|
662 | 664 | else: |
|
|
663 | 665 | # ojo debo anadir el readNextBLock y el __isBufferEmpty( |
|
|
664 | 666 | self.dataOut.flagNoData = False |
|
|
665 | 667 | buffer = self.__data_buffer[:,self.__bufferIndex:self.__bufferIndex + self.__samples_to_read] |
|
|
666 | 668 | buffer = buffer.reshape((self.__nChannels, self.nProfileBlocks, int(self.__samples_to_read/self.nProfileBlocks))) |
|
|
667 | 669 | self.dataOut.nProfileBlocks = self.nProfileBlocks |
|
|
668 | 670 | self.dataOut.data = buffer |
|
|
669 | 671 | self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate |
|
|
670 | 672 | self.profileIndex += self.__samples_to_read |
|
|
671 | 673 | self.__bufferIndex += self.__samples_to_read |
|
|
672 | 674 | self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock |
|
|
673 | 675 | return True |
|
|
674 | 676 | |
|
|
675 | 677 | |
|
|
676 | 678 | def printInfo(self): |
|
|
677 | 679 | ''' |
|
|
678 | 680 | ''' |
|
|
679 | 681 | if self.__printInfo == False: |
|
|
680 | 682 | return |
|
|
681 | 683 | |
|
|
682 | 684 | # self.systemHeaderObj.printInfo() |
|
|
683 | 685 | # self.radarControllerHeaderObj.printInfo() |
|
|
684 | 686 | |
|
|
685 | 687 | self.__printInfo = False |
|
|
686 | 688 | |
|
|
687 | 689 | def printNumberOfBlock(self): |
|
|
688 | 690 | ''' |
|
|
689 | 691 | ''' |
|
|
690 | 692 | return |
|
|
691 | 693 | # print self.profileIndex |
|
|
692 | 694 | |
|
|
693 | 695 | def run(self, **kwargs): |
|
|
694 | 696 | ''' |
|
|
695 | 697 | This method will be called many times so here you should put all your code |
|
|
696 | 698 | ''' |
|
|
697 | 699 | |
|
|
698 | 700 | if not self.isConfig: |
|
|
699 | 701 | self.setup(**kwargs) |
|
|
700 | 702 | |
|
|
701 | 703 | self.getData(seconds=self.__delay) |
|
|
702 | 704 | |
|
|
703 | 705 | return |
|
|
704 | 706 | |
|
|
705 | 707 | @MPDecorator |
|
|
706 | 708 | class DigitalRFWriter(Operation): |
|
|
707 | 709 | ''' |
|
|
708 | 710 | classdocs |
|
|
709 | 711 | ''' |
|
|
710 | 712 | |
|
|
711 | 713 | def __init__(self, **kwargs): |
|
|
712 | 714 | ''' |
|
|
713 | 715 | Constructor |
|
|
714 | 716 | ''' |
|
|
715 | 717 | Operation.__init__(self, **kwargs) |
|
|
716 | 718 | self.metadata_dict = {} |
|
|
717 | 719 | self.dataOut = None |
|
|
718 | 720 | self.dtype = None |
|
|
719 | 721 | self.oldAverage = 0 |
|
|
720 | 722 | |
|
|
721 | 723 | def setHeader(self): |
|
|
722 | 724 | |
|
|
723 | 725 | self.metadata_dict['frequency'] = self.dataOut.frequency |
|
|
724 | 726 | self.metadata_dict['timezone'] = self.dataOut.timeZone |
|
|
725 | 727 | self.metadata_dict['dtype'] = pickle.dumps(self.dataOut.dtype) |
|
|
726 | 728 | self.metadata_dict['nProfiles'] = self.dataOut.nProfiles |
|
|
727 | 729 | self.metadata_dict['heightList'] = self.dataOut.heightList |
|
|
728 | 730 | self.metadata_dict['channelList'] = self.dataOut.channelList |
|
|
729 | 731 | self.metadata_dict['flagDecodeData'] = self.dataOut.flagDecodeData |
|
|
730 | 732 | self.metadata_dict['flagDeflipData'] = self.dataOut.flagDeflipData |
|
|
731 | 733 | self.metadata_dict['flagShiftFFT'] = self.dataOut.flagShiftFFT |
|
|
732 | 734 | self.metadata_dict['useLocalTime'] = self.dataOut.useLocalTime |
|
|
733 | 735 | self.metadata_dict['nCohInt'] = self.dataOut.nCohInt |
|
|
734 | 736 | self.metadata_dict['type'] = self.dataOut.type |
|
|
735 | 737 | self.metadata_dict['flagDataAsBlock']= getattr( |
|
|
736 | 738 | self.dataOut, 'flagDataAsBlock', None) # chequear |
|
|
737 | 739 | |
|
|
738 | 740 | def setup(self, dataOut, path, frequency, fileCadence, dirCadence, metadataCadence, set=0, metadataFile='metadata', ext='.h5'): |
|
|
739 | 741 | ''' |
|
|
740 | 742 | In this method we should set all initial parameters. |
|
|
741 | 743 | Input: |
|
|
742 | 744 | dataOut: Input data will also be outputa data |
|
|
743 | 745 | ''' |
|
|
744 | 746 | self.setHeader() |
|
|
745 | 747 | self.__ippSeconds = dataOut.ippSeconds |
|
|
746 | 748 | self.__deltaH = dataOut.getDeltaH() |
|
|
747 | 749 | self.__sample_rate = 1e6 * 0.15 / self.__deltaH |
|
|
748 | 750 | self.__dtype = dataOut.dtype |
|
|
749 | 751 | if len(dataOut.dtype) == 2: |
|
|
750 | 752 | self.__dtype = dataOut.dtype[0] |
|
|
751 | 753 | self.__nSamples = dataOut.systemHeaderObj.nSamples |
|
|
752 | 754 | self.__nProfiles = dataOut.nProfiles |
|
|
753 | 755 | |
|
|
754 | 756 | if self.dataOut.type != 'Voltage': |
|
|
755 | 757 | raise 'Digital RF cannot be used with this data type' |
|
|
756 | 758 | self.arr_data = numpy.ones((1, dataOut.nFFTPoints * len( |
|
|
757 | 759 | self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)]) |
|
|
758 | 760 | else: |
|
|
759 | 761 | self.arr_data = numpy.ones((self.__nSamples, len( |
|
|
760 | 762 | self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)]) |
|
|
761 | 763 | |
|
|
762 | 764 | file_cadence_millisecs = 1000 |
|
|
763 | 765 | |
|
|
764 | 766 | sample_rate_fraction = Fraction(self.__sample_rate).limit_denominator() |
|
|
765 | 767 | sample_rate_numerator = int(sample_rate_fraction.numerator) |
|
|
766 | 768 | sample_rate_denominator = int(sample_rate_fraction.denominator) |
|
|
767 | 769 | start_global_index = dataOut.utctime * self.__sample_rate |
|
|
768 | 770 | |
|
|
769 | 771 | uuid = 'prueba' |
|
|
770 | 772 | compression_level = 0 |
|
|
771 | 773 | checksum = False |
|
|
772 | 774 | is_complex = True |
|
|
773 | 775 | num_subchannels = len(dataOut.channelList) |
|
|
774 | 776 | is_continuous = True |
|
|
775 | 777 | marching_periods = False |
|
|
776 | 778 | |
|
|
777 | 779 | self.digitalWriteObj = digital_rf.DigitalRFWriter(path, self.__dtype, dirCadence, |
|
|
778 | 780 | fileCadence, start_global_index, |
|
|
779 | 781 | sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum, |
|
|
780 | 782 | is_complex, num_subchannels, is_continuous, marching_periods) |
|
|
781 | 783 | metadata_dir = os.path.join(path, 'metadata') |
|
|
782 | 784 | os.system('mkdir %s' % (metadata_dir)) |
|
|
783 | 785 | self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, # 236, file_cadence_millisecs / 1000 |
|
|
784 | 786 | sample_rate_numerator, sample_rate_denominator, |
|
|
785 | 787 | metadataFile) |
|
|
786 | 788 | self.isConfig = True |
|
|
787 | 789 | self.currentSample = 0 |
|
|
788 | 790 | self.oldAverage = 0 |
|
|
789 | 791 | self.count = 0 |
|
|
790 | 792 | return |
|
|
791 | 793 | |
|
|
792 | 794 | def writeMetadata(self): |
|
|
793 | 795 | start_idx = self.__sample_rate * self.dataOut.utctime |
|
|
794 | 796 | |
|
|
795 | 797 | self.metadata_dict['processingHeader'] = self.dataOut.processingHeaderObj.getAsDict( |
|
|
796 | 798 | ) |
|
|
797 | 799 | self.metadata_dict['radarControllerHeader'] = self.dataOut.radarControllerHeaderObj.getAsDict( |
|
|
798 | 800 | ) |
|
|
799 | 801 | self.metadata_dict['systemHeader'] = self.dataOut.systemHeaderObj.getAsDict( |
|
|
800 | 802 | ) |
|
|
801 | 803 | self.digitalMetadataWriteObj.write(start_idx, self.metadata_dict) |
|
|
802 | 804 | return |
|
|
803 | 805 | |
|
|
804 | 806 | def timeit(self, toExecute): |
|
|
805 | 807 | t0 = time() |
|
|
806 | 808 | toExecute() |
|
|
807 | 809 | self.executionTime = time() - t0 |
|
|
808 | 810 | if self.oldAverage is None: |
|
|
809 | 811 | self.oldAverage = self.executionTime |
|
|
810 | 812 | self.oldAverage = (self.executionTime + self.count * |
|
|
811 | 813 | self.oldAverage) / (self.count + 1.0) |
|
|
812 | 814 | self.count = self.count + 1.0 |
|
|
813 | 815 | return |
|
|
814 | 816 | |
|
|
815 | 817 | def writeData(self): |
|
|
816 | 818 | if self.dataOut.type != 'Voltage': |
|
|
817 | 819 | raise 'Digital RF cannot be used with this data type' |
|
|
818 | 820 | for channel in self.dataOut.channelList: |
|
|
819 | 821 | for i in range(self.dataOut.nFFTPoints): |
|
|
820 | 822 | self.arr_data[1][channel * self.dataOut.nFFTPoints + |
|
|
821 | 823 | i]['r'] = self.dataOut.data[channel][i].real |
|
|
822 | 824 | self.arr_data[1][channel * self.dataOut.nFFTPoints + |
|
|
823 | 825 | i]['i'] = self.dataOut.data[channel][i].imag |
|
|
824 | 826 | else: |
|
|
825 | 827 | for i in range(self.dataOut.systemHeaderObj.nSamples): |
|
|
826 | 828 | for channel in self.dataOut.channelList: |
|
|
827 | 829 | self.arr_data[i][channel]['r'] = self.dataOut.data[channel][i].real |
|
|
828 | 830 | self.arr_data[i][channel]['i'] = self.dataOut.data[channel][i].imag |
|
|
829 | 831 | |
|
|
830 | 832 | def f(): return self.digitalWriteObj.rf_write(self.arr_data) |
|
|
831 | 833 | self.timeit(f) |
|
|
832 | 834 | |
|
|
833 | 835 | return |
|
|
834 | 836 | |
|
|
835 | 837 | def run(self, dataOut, frequency=49.92e6, path=None, fileCadence=1000, dirCadence=36000, metadataCadence=1, **kwargs): |
|
|
836 | 838 | ''' |
|
|
837 | 839 | This method will be called many times so here you should put all your code |
|
|
838 | 840 | Inputs: |
|
|
839 | 841 | dataOut: object with the data |
|
|
840 | 842 | ''' |
|
|
841 | 843 | # print dataOut.__dict__ |
|
|
842 | 844 | self.dataOut = dataOut |
|
|
843 | 845 | if not self.isConfig: |
|
|
844 | 846 | self.setup(dataOut, path, frequency, fileCadence, |
|
|
845 | 847 | dirCadence, metadataCadence, **kwargs) |
|
|
846 | 848 | self.writeMetadata() |
|
|
847 | 849 | |
|
|
848 | 850 | self.writeData() |
|
|
849 | 851 | |
|
|
850 | 852 | ## self.currentSample += 1 |
|
|
851 | 853 | # if self.dataOut.flagDataAsBlock or self.currentSample == 1: |
|
|
852 | 854 | # self.writeMetadata() |
|
|
853 | 855 | ## if self.currentSample == self.__nProfiles: self.currentSample = 0 |
|
|
854 | 856 | |
|
|
855 | 857 | return dataOut# en la version 2.7 no aparece este return |
|
|
856 | 858 | |
|
|
857 | 859 | def close(self): |
|
|
858 | 860 | print('[Writing] - Closing files ') |
|
|
859 | 861 | print('Average of writing to digital rf format is ', self.oldAverage * 1000) |
|
|
860 | 862 | try: |
|
|
861 | 863 | self.digitalWriteObj.close() |
|
|
862 | 864 | except: |
|
|
863 | 865 | pass |
|
|
1 | NO CONTENT: modified file | |
| The requested commit or file is too big and content was truncated. Show full diff | |||
General Comments 0
You need to be logged in to leave comments.
Login now