##// END OF EJS Templates
update fix axes, test de prueba de canales
jespinoza -
r1476:d08bdce68e00
parent child
Show More
@@ -0,0 +1,362
1 #!python
2 '''
3 '''
4
5 import os, sys
6 import datetime
7 import time
8
9 #path = os.path.dirname(os.getcwd())
10 #path = os.path.dirname(path)
11 #sys.path.insert(0, path)
12
13 from schainpy.controller import Project
14
15 desc = "USRP_test"
16 filename = "USRP_processing.xml"
17 controllerObj = Project()
18 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
19
20 ############## USED TO PLOT IQ VOLTAGE, POWER AND SPECTRA #############
21
22 #######################################################################
23 ######PATH DE LECTURA, ESCRITURA, GRAFICOS Y ENVIO WEB#################
24 #######################################################################
25 #path = '/media/data/data/vientos/57.2063km/echoes/NCO_Woodman'
26 #path = '/DATA_RM/TEST_INTEGRACION'
27 #path = '/DATA_RM/PRUEBA_USRP_RP'
28 #path = '/DATA_RM/PRUEBA_USRP_RP'
29
30 path = '/DATA_RM/TEST_2M'
31 path = '/DATA_RM/TEST_2M_UD'
32 path = '/DATA_RM/2MHZ17022022'
33 path = '/DATA_RM/10MHZTEST/'
34 path = '/DATA_RM/10MHZDRONE/'
35
36
37 path= '/home/soporte/TEST_500mVPP'
38 path= '/home/soporte/TEST_1VPP+500mVDC'
39 path = '/home/soporte/TEST_500mVPP+500mVDC'
40 path = '/home/soporte/TEST_1.5VPP'
41 path = '/home/soporte/TEST_2VPP'
42 path= '/home/soporte/TEST_1VPP'
43 path = '/home/soporte/Documents/HUANCAYO/TEST_HYO_PM@2022-05-14T11-28-19/rawdata'
44
45 #HYO_PM@2022-05-28T00-00-17
46 path = '/DATA_RM/DATA/HYO_PM@2022-05-28T00-00-17/rawdata'
47
48 #figpath = '/home/soporte/Pictures/TEST_RP_0001'
49 #figpath = '/home/soporte/Pictures/TEST_RP_6000'
50 figpath = '/home/soporte/Pictures/USRP_TEST_2M'
51 figpath = '/home/soporte/Pictures/USRP_TEST_2M_UD'
52 figpath = '/home/soporte/Pictures/10MHZDRONE'
53 figpath = '/home/soporte/Pictures/500mVPP'
54 figpath = '/home/soporte/Pictures/1VPP+500mVDC'
55 figpath = '/home/soporte/Pictures/TEST_500mVPP+500mVDC'
56 figpath = '/home/soporte/Pictures/TEST_1.5VPP'
57 figpath = '/home/soporte/Pictures/TEST_2VPP'
58 figpath = '/home/soporte/Pictures/TEST_1VPP'
59
60
61
62
63 #remotefolder = "/home/wmaster/graficos"
64 #######################################################################
65 ################# RANGO DE PLOTEO######################################
66 #######################################################################
67 dBmin = '-60'#'-20'
68 dBmax = '-5'#'-85'
69 xmin = '0'
70 xmax ='24'
71 ymin = '0'
72 ymax = '10'
73 #######################################################################
74 ########################FECHA##########################################
75 #######################################################################
76 str = datetime.date.today()
77 today = str.strftime("%Y/%m/%d")
78 str2 = str - datetime.timedelta(days=1)
79 yesterday = str2.strftime("%Y/%m/%d")
80 #######################################################################
81 ######################## UNIDAD DE LECTURA#############################
82 #######################################################################
83 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
84 path=path,
85 startDate="2022/05/28",#today,
86 endDate="2022/05/28",#today,
87 startTime='00:00:00',# inicio libre
88 #startTime='00:00:00',
89 endTime='23:59:59',
90 delay=0,
91 #set=0,
92 online=0,
93 walk=1,
94 ippKm = 60)
95
96 opObj11 = readUnitConfObj.addOperation(name='printInfo')
97 #opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
98 #######################################################################
99 ################ OPERACIONES DOMINIO DEL TIEMPO########################
100 #######################################################################
101
102 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
103
104 op3 = procUnitConfObjA.addOperation(name='ProfileSelector', optype='other')
105 op3.addParameter(name='profileRangeList', value='1,123')
106
107
108 code=[[1]]
109
110 opObj11 = procUnitConfObjA.addOperation(name='Decoder', optype='other')
111 opObj11.addParameter(name='code', value=code)
112 opObj11.addParameter(name='nCode', value='1', format='int')
113 opObj11.addParameter(name='nBaud', value='1', format='int')
114
115
116 '''
117 op3 = procUnitConfObjA.addOperation(name='ProfileSelector', optype='other')
118 op3.addParameter(name='profileRangeList', value='122,249')
119 code8=[[1,1,1,0,1,1,0,1],[1,1,1,0,0,0,1,0]]
120
121 opObj11 = procUnitConfObjA.addOperation(name='Decoder', optype='other')
122 opObj11.addParameter(name='code', value=code8)
123 opObj11.addParameter(name='nCode', value='2', format='int')
124 opObj11.addParameter(name='nBaud', value='8', format='int')
125 '''
126 op = procUnitConfObjA.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
127 op.addParameter(name='n', value=2, format='int')
128
129
130 '''
131
132 # OJO SCOPE
133 opObj10 = procUnitConfObjA.addOperation(name='ScopePlot', optype='external')
134 opObj10.addParameter(name='id', value='10', format='int')
135 opObj10.addParameter(name='xmin', value='0', format='int')
136 opObj10.addParameter(name='xmax', value='60', format='int')
137 opObj10.addParameter(name='type', value='iq')
138 #opObj10.addParameter(name='ymin', value='-0.20000', format='int')
139 #opObj10.addParameter(name='ymax', value='0.20000', format='int')
140 opObj10.addParameter(name='save', value=figpath, format='str')
141 opObj10.addParameter(name='save_period', value=1, format='int')
142 '''
143 '''
144 opObj11 = procUnitConfObjA.addOperation(name='selectHeights')
145 opObj11.addParameter(name='minIndex', value='1', format='int')
146 # opObj11.addParameter(name='maxIndex', value='10000', format='int')
147 opObj11.addParameter(name='maxIndex', value='200', format='int')
148 '''
149 #
150 # codigo64='1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1,1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0,'+\
151 # '1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1,0,0,0,1,0,0,1,0,0,0,0,1,1,1,0,1,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1'
152
153 #opObj11 = procUnitConfObjA.addOperation(name='setRadarFrequency')
154 #opObj11.addParameter(name='frequency', value='49920000')
155
156 '''
157 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
158 opObj11.addParameter(name='n', value='625', format='int')#10
159 opObj11.addParameter(name='removeDC', value=1, format='int')
160 '''
161
162 # Ploteo TEST
163 '''
164 opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')
165 opObj11 = procUnitConfObjA.addOperation(name='PulsepairSignalPlot', optype='other')
166 opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')
167 #opObj11.addParameter(name='xmax', value=8)
168 opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')
169 '''
170 # OJO SCOPE
171 #opObj10 = procUnitConfObjA.addOperation(name='ScopePlot', optype='external')
172 #opObj10.addParameter(name='id', value='10', format='int')
173 ##opObj10.addParameter(name='xmin', value='0', format='int')
174 ##opObj10.addParameter(name='xmax', value='50', format='int')
175 #opObj10.addParameter(name='type', value='iq')
176 ##opObj10.addParameter(name='ymin', value='-5000', format='int')
177 ##opObj10.addParameter(name='ymax', value='8500', format='int')
178 #opObj11.addParameter(name='save', value=figpath, format='str')
179 #opObj11.addParameter(name='save_period', value=10, format='int')
180
181 #opObj10 = procUnitConfObjA.addOperation(name='setH0')
182 #opObj10.addParameter(name='h0', value='-5000', format='float')
183
184 #opObj11 = procUnitConfObjA.addOperation(name='filterByHeights')
185 #opObj11.addParameter(name='window', value='1', format='int')
186
187 #codigo='1,1,-1,1,1,-1,1,-1,-1,1,-1,-1,-1,1,-1,-1,-1,1,-1,-1,-1,1,1,1,1,-1,-1,-1'
188 #opObj11 = procUnitConfObjSousy.addOperation(name='Decoder', optype='other')
189 #opObj11.addParameter(name='code', value=codigo, format='floatlist')
190 #opObj11.addParameter(name='nCode', value='1', format='int')
191 #opObj11.addParameter(name='nBaud', value='28', format='int')
192
193 #opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other')
194 #opObj11.addParameter(name='n', value='100', format='int')
195
196 #######################################################################
197 ########## OPERACIONES ParametersProc########################
198 #######################################################################
199 ###procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
200 '''
201
202 opObj11 = procUnitConfObjA.addOperation(name='PedestalInformation')
203 opObj11.addParameter(name='path_ped', value=path_ped)
204 opObj11.addParameter(name='path_adq', value=path_adq)
205 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
206 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
207 opObj11.addParameter(name='blocksPerfile', value='100', format='int')
208 opObj11.addParameter(name='f_a_p', value='25', format='int')
209 opObj11.addParameter(name='online', value='0', format='int')
210
211 opObj11 = procUnitConfObjA.addOperation(name='Block360')
212 opObj11.addParameter(name='n', value='40', format='int')
213
214 opObj11= procUnitConfObjA.addOperation(name='WeatherPlot',optype='other')
215 opObj11.addParameter(name='save', value=figpath)
216 opObj11.addParameter(name='save_period', value=1)
217
218 8
219 '''
220
221
222 '''
223 opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other')
224 opObj11.addParameter(name='n', value='250', format='int')
225 '''
226 #######################################################################
227 ########## OPERACIONES DOMINIO DE LA FRECUENCIA########################
228 #######################################################################
229 '''
230 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
231 procUnitConfObjB.addParameter(name='nFFTPoints', value='64', format='int')
232 procUnitConfObjB.addParameter(name='nProfiles', value='64', format='int')
233 '''
234
235 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
236 procUnitConfObjB.addParameter(name='nFFTPoints', value='61', format='int')
237 procUnitConfObjB.addParameter(name='nProfiles', value='61', format='int')
238
239 '''
240 procUnitConfObjC = controllerObj.addProcUnit(datatype='SpectraHeisProc', inputId=procUnitConfObjA.getId())
241 #procUnitConfObjB.addParameter(name='nFFTPoints', value='64', format='int')
242 #procUnitConfObjB.addParameter(name='nProfiles', value='64', format='int')
243 opObj11 = procUnitConfObjC.addOperation(name='IncohInt4SpectraHeis', optype='other')
244 #opObj11.addParameter(name='timeInterval', value='4', format='int')
245 opObj11.addParameter(name='n', value='100', format='int')
246
247 #procUnitConfObjB.addParameter(name='pairsList', value='(0,0),(1,1),(0,1)', format='pairsList')
248
249 #opObj13 = procUnitConfObjB.addOperation(name='removeDC')
250 #opObj13.addParameter(name='mode', value='2', format='int')
251
252 #opObj11 = procUnitConfObjB.addOperation(name='IncohInt', optype='other')
253 #opObj11.addParameter(name='n', value='8', format='float')
254 #######################################################################
255 ########## PLOTEO DOMINIO DE LA FRECUENCIA#############################
256 #######################################################################
257 #----
258 '''
259 '''
260 opObj11 = procUnitConfObjC.addOperation(name='SpectraHeisPlot')
261 opObj11.addParameter(name='id', value='10', format='int')
262 opObj11.addParameter(name='wintitle', value='Spectra_Alturas', format='str')
263 #opObj11.addParameter(name='xmin', value=-100000, format='float')
264 #opObj11.addParameter(name='xmax', value=100000, format='float')
265 opObj11.addParameter(name='oneFigure', value=False,format='bool')
266 #opObj11.addParameter(name='zmin', value=-10, format='int')
267 #opObj11.addParameter(name='zmax', value=40, format='int')
268 opObj11.addParameter(name='ymin', value=10, format='int')
269 opObj11.addParameter(name='ymax', value=55, format='int')
270 opObj11.addParameter(name='grid', value=True, format
271 [Reading] 2022-05-23 12:27:32.732775: 21333 samples <> 0.010667 sec
272 ='bool')
273 #opObj11.addParameter(name='showprofile', value='1', format='int')
274 opObj11.addParameter(name='save', value=figpath, format='str')
275 #opObj11.addParameter(name='save_period', value=10, format='int')
276 '''
277 '''
278 opObj11 = procUnitConfObjC.addOperation(name='RTIHeisPlot')
279 opObj11.addParameter(name='id', value='10', format='int')
280 opObj11.addParameter(name='wintitle', value='RTI_Alturas', format='str')
281 opObj11.addParameter(name='xmin', value=11.0, format='float')
282 opObj11.addParameter(name='xmax', value=18.0, format='float')
283 opObj11.addParameter(name='zmin', value=10, format='int')
284 opObj11.addParameter(name='zmax', value=30, format='int')
285 opObj11.addParameter(name='ymin', value=5, format='int')
286 opObj11.addParameter(name='ymax', value=28, format='int')
287 opObj11.addParameter(name='showprofile', value='1', format='int')
288 opObj11.addParameter(name='save', value=figpath, format='str')
289 opObj11.addParameter(name='save_period', value=10, format='int')
290 '''
291
292 #SpectraPlot
293
294 opObj11 = procUnitConfObjB.addOperation(name='SpectraPlot', optype='external')
295 opObj11.addParameter(name='id', value='1', format='int')
296 opObj11.addParameter(name='wintitle', value='Spectra', format='str')
297 #opObj11.addParameter(name='xmin', value=-0.01, format='float')
298 #opObj11.addParameter(name='xmax', value=0.01, format='float')
299 opObj11.addParameter(name='zmin', value=dBmin, format='int')
300 opObj11.addParameter(name='zmax', value=dBmax, format='int')
301 opObj11.addParameter(name='ymin', value=ymin, format='int')
302 opObj11.addParameter(name='ymax', value=ymax, format='int')
303 opObj11.addParameter(name='showprofile', value='1', format='int')
304 opObj11.addParameter(name='save', value=figpath, format='str')
305 opObj11.addParameter(name='save_period', value=10, format='int')
306
307
308 #RTIPLOT
309 '''
310 opObj11 = procUnitConfObjB.addOperation(name='RTIPlot', optype='external')
311 opObj11.addParameter(name='id', value='2', format='int')
312 opObj11.addParameter(name='wintitle', value='RTIPlot', format='str')
313 opObj11.addParameter(name='zmin', value=dBmin, format='int')
314 opObj11.addParameter(name='zmax', value=dBmax, format='int')
315 #opObj11.addParameter(name='ymin', value=ymin, format='int')
316 #opObj11.addParameter(name='ymax', value=ymax, format='int')
317 #opObj11.addParameter(name='xmin', value=15, format='int')
318 #opObj11.addParameter(name='xmax', value=16, format='int')
319 opObj11.addParameter(name='zmin', value=dBmin, format='int')
320 opObj11.addParameter(name='zmax', value=dBmax, format='int')
321
322 opObj11.addParameter(name='showprofile', value='1', format='int')
323 opObj11.addParameter(name='save', value=figpath, format='str')
324 opObj11.addParameter(name='save_period', value=10, format='int')
325 '''
326 '''
327 # opObj11 = procUnitConfObjB.addOperation(name='CrossSpectraPlot', optype='other')
328 # opObj11.addParameter(name='id', value='3', format='int')
329 # opObj11.addParameter(name='wintitle', value='CrossSpectraPlot', format='str')
330 # opObj11.addParameter(name='ymin', value=ymin, format='int')
331 # opObj11.addParameter(name='ymax', value=ymax, format='int')
332 # opObj11.addParameter(name='phase_cmap', value='jet', format='str')
333 # opObj11.addParameter(name='zmin', value=dBmin, format='int')
334 # opObj11.addParameter(name='zmax', value=dBmax, format='int')
335 # opObj11.addParameter(name='figpath', value=figures_path, format='str')
336 # opObj11.addParameter(name='save', value=0, format='bool')
337 # opObj11.addParameter(name='pairsList', value='(0,1)', format='pairsList')
338 # #
339 # opObj11 = procUnitConfObjB.addOperation(name='CoherenceMap', optype='other')
340 # opObj11.addParameter(name='id', value='4', format='int')
341 # opObj11.addParameter(name='wintitle', value='Coherence', format='str')
342 # opObj11.addParameter(name='phase_cmap', value='jet', format='str')
343 # opObj11.addParameter(name='xmin', value=xmin, format='float')
344 # opObj11.addParameter(name='xmax', value=xmax, format='float')
345 # opObj11.addParameter(name='figpath', value=figures_path, format='str')
346 # opObj11.addParameter(name='save', value=0, format='bool')
347 # opObj11.addParameter(name='pairsList', value='(0,1)', format='pairsList')
348 #
349 '''
350 '''
351 #######################################################################
352 ############### UNIDAD DE ESCRITURA ###################################
353 #######################################################################
354 #opObj11 = procUnitConfObjB.addOperation(name='SpectraWriter', optype='other')
355 #opObj11.addParameter(name='path', value=wr_path)
356 #opObj11.addParameter(name='blocksPerFile', value='50', format='int')
357 print ("Escribiendo el archivo XML")
358 print ("Leyendo el archivo XML")
359 '''
360
361
362 controllerObj.start()
@@ -1,716 +1,716
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Base class to create plot operations
6 6
7 7 """
8 8
9 9 import os
10 10 import sys
11 11 import zmq
12 12 import time
13 13 import numpy
14 14 import datetime
15 15 from collections import deque
16 16 from functools import wraps
17 17 from threading import Thread
18 18 import matplotlib
19 19
20 20 if 'BACKEND' in os.environ:
21 21 matplotlib.use(os.environ['BACKEND'])
22 22 elif 'linux' in sys.platform:
23 23 matplotlib.use("TkAgg")
24 24 elif 'darwin' in sys.platform:
25 25 matplotlib.use('MacOSX')
26 26 else:
27 27 from schainpy.utils import log
28 28 log.warning('Using default Backend="Agg"', 'INFO')
29 29 matplotlib.use('Agg')
30 30
31 31 import matplotlib.pyplot as plt
32 32 from matplotlib.patches import Polygon
33 33 from mpl_toolkits.axes_grid1 import make_axes_locatable
34 34 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
35 35
36 36 from schainpy.model.data.jrodata import PlotterData
37 37 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
38 38 from schainpy.utils import log
39 39
40 40 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
41 41 blu_values = matplotlib.pyplot.get_cmap(
42 42 'seismic_r', 20)(numpy.arange(20))[10:15]
43 43 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
44 44 'jro', numpy.vstack((blu_values, jet_values)))
45 45 matplotlib.pyplot.register_cmap(cmap=ncmap)
46 46
47 47 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
48 48 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
49 49
50 50 EARTH_RADIUS = 6.3710e3
51 51
52 52 def ll2xy(lat1, lon1, lat2, lon2):
53 53
54 54 p = 0.017453292519943295
55 55 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
56 56 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
57 57 r = 12742 * numpy.arcsin(numpy.sqrt(a))
58 58 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
59 59 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
60 60 theta = -theta + numpy.pi/2
61 61 return r*numpy.cos(theta), r*numpy.sin(theta)
62 62
63 63
64 64 def km2deg(km):
65 65 '''
66 66 Convert distance in km to degrees
67 67 '''
68 68
69 69 return numpy.rad2deg(km/EARTH_RADIUS)
70 70
71 71
72 72 def figpause(interval):
73 73 backend = plt.rcParams['backend']
74 74 if backend in matplotlib.rcsetup.interactive_bk:
75 75 figManager = matplotlib._pylab_helpers.Gcf.get_active()
76 76 if figManager is not None:
77 77 canvas = figManager.canvas
78 78 if canvas.figure.stale:
79 79 canvas.draw()
80 80 try:
81 81 canvas.start_event_loop(interval)
82 82 except:
83 83 pass
84 84 return
85 85
86 86 def popup(message):
87 87 '''
88 88 '''
89 89
90 90 fig = plt.figure(figsize=(12, 8), facecolor='r')
91 91 text = '\n'.join([s.strip() for s in message.split(':')])
92 92 fig.text(0.01, 0.5, text, ha='left', va='center',
93 93 size='20', weight='heavy', color='w')
94 94 fig.show()
95 95 figpause(1000)
96 96
97 97
98 98 class Throttle(object):
99 99 '''
100 100 Decorator that prevents a function from being called more than once every
101 101 time period.
102 102 To create a function that cannot be called more than once a minute, but
103 103 will sleep until it can be called:
104 104 @Throttle(minutes=1)
105 105 def foo():
106 106 pass
107 107
108 108 for i in range(10):
109 109 foo()
110 110 print "This function has run %s times." % i
111 111 '''
112 112
113 113 def __init__(self, seconds=0, minutes=0, hours=0):
114 114 self.throttle_period = datetime.timedelta(
115 115 seconds=seconds, minutes=minutes, hours=hours
116 116 )
117 117
118 118 self.time_of_last_call = datetime.datetime.min
119 119
120 120 def __call__(self, fn):
121 121 @wraps(fn)
122 122 def wrapper(*args, **kwargs):
123 123 coerce = kwargs.pop('coerce', None)
124 124 if coerce:
125 125 self.time_of_last_call = datetime.datetime.now()
126 126 return fn(*args, **kwargs)
127 127 else:
128 128 now = datetime.datetime.now()
129 129 time_since_last_call = now - self.time_of_last_call
130 130 time_left = self.throttle_period - time_since_last_call
131 131
132 132 if time_left > datetime.timedelta(seconds=0):
133 133 return
134 134
135 135 self.time_of_last_call = datetime.datetime.now()
136 136 return fn(*args, **kwargs)
137 137
138 138 return wrapper
139 139
140 140 def apply_throttle(value):
141 141
142 142 @Throttle(seconds=value)
143 143 def fnThrottled(fn):
144 144 fn()
145 145
146 146 return fnThrottled
147 147
148 148
149 149 @MPDecorator
150 150 class Plot(Operation):
151 151 """Base class for Schain plotting operations
152 152
153 153 This class should never be use directtly you must subclass a new operation,
154 154 children classes must be defined as follow:
155 155
156 156 ExamplePlot(Plot):
157 157
158 158 CODE = 'code'
159 159 colormap = 'jet'
160 160 plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
161 161
162 162 def setup(self):
163 163 pass
164 164
165 165 def plot(self):
166 166 pass
167 167
168 168 """
169 169
170 170 CODE = 'Figure'
171 171 colormap = 'jet'
172 172 bgcolor = 'white'
173 173 buffering = True
174 174 __missing = 1E30
175 175
176 176 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
177 177 'showprofile']
178 178
179 179 def __init__(self):
180 180
181 181 Operation.__init__(self)
182 182 self.isConfig = False
183 183 self.isPlotConfig = False
184 184 self.save_time = 0
185 185 self.sender_time = 0
186 186 self.data = None
187 187 self.firsttime = True
188 188 self.sender_queue = deque(maxlen=10)
189 189 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
190 190
191 191 def __fmtTime(self, x, pos):
192 192 '''
193 193 '''
194 194
195 195 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
196 196
197 197 def __setup(self, **kwargs):
198 198 '''
199 199 Initialize variables
200 200 '''
201 201
202 202 self.figures = []
203 203 self.axes = []
204 204 self.cb_axes = []
205 205 self.localtime = kwargs.pop('localtime', True)
206 206 self.show = kwargs.get('show', True)
207 207 self.save = kwargs.get('save', False)
208 208 self.save_period = kwargs.get('save_period', 0)
209 209 self.colormap = kwargs.get('colormap', self.colormap)
210 210 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
211 211 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
212 212 self.colormaps = kwargs.get('colormaps', None)
213 213 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
214 214 self.showprofile = kwargs.get('showprofile', False)
215 215 self.title = kwargs.get('wintitle', self.CODE.upper())
216 216 self.cb_label = kwargs.get('cb_label', None)
217 217 self.cb_labels = kwargs.get('cb_labels', None)
218 218 self.labels = kwargs.get('labels', None)
219 219 self.xaxis = kwargs.get('xaxis', 'frequency')
220 220 self.zmin = kwargs.get('zmin', None)
221 221 self.zmax = kwargs.get('zmax', None)
222 222 self.zlimits = kwargs.get('zlimits', None)
223 223 self.xmin = kwargs.get('xmin', None)
224 224 self.xmax = kwargs.get('xmax', None)
225 225 self.xrange = kwargs.get('xrange', 12)
226 226 self.xscale = kwargs.get('xscale', None)
227 227 self.ymin = kwargs.get('ymin', None)
228 228 self.ymax = kwargs.get('ymax', None)
229 229 self.yscale = kwargs.get('yscale', None)
230 230 self.xlabel = kwargs.get('xlabel', None)
231 231 self.attr_time = kwargs.get('attr_time', 'utctime')
232 232 self.attr_data = kwargs.get('attr_data', 'data_param')
233 233 self.decimation = kwargs.get('decimation', None)
234 234 self.oneFigure = kwargs.get('oneFigure', True)
235 235 self.width = kwargs.get('width', None)
236 236 self.height = kwargs.get('height', None)
237 237 self.colorbar = kwargs.get('colorbar', True)
238 238 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
239 239 self.channels = kwargs.get('channels', None)
240 240 self.titles = kwargs.get('titles', [])
241 241 self.polar = False
242 242 self.type = kwargs.get('type', 'iq')
243 243 self.grid = kwargs.get('grid', False)
244 244 self.pause = kwargs.get('pause', False)
245 245 self.save_code = kwargs.get('save_code', self.CODE)
246 246 self.throttle = kwargs.get('throttle', 0)
247 247 self.exp_code = kwargs.get('exp_code', None)
248 248 self.server = kwargs.get('server', False)
249 249 self.sender_period = kwargs.get('sender_period', 60)
250 250 self.tag = kwargs.get('tag', '')
251 251 self.height_index = kwargs.get('height_index', None)
252 252 self.__throttle_plot = apply_throttle(self.throttle)
253 253 code = self.attr_data if self.attr_data else self.CODE
254 254 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
255 255 self.ang_min = kwargs.get('ang_min', None)
256 256 self.ang_max = kwargs.get('ang_max', None)
257 257 self.mode = kwargs.get('mode', None)
258 258
259 259
260 260
261 261 if self.server:
262 262 if not self.server.startswith('tcp://'):
263 263 self.server = 'tcp://{}'.format(self.server)
264 264 log.success(
265 265 'Sending to server: {}'.format(self.server),
266 266 self.name
267 267 )
268 268
269 269 if isinstance(self.attr_data, str):
270 270 self.attr_data = [self.attr_data]
271 271
272 272 def __setup_plot(self):
273 273 '''
274 274 Common setup for all figures, here figures and axes are created
275 275 '''
276 276
277 277 self.setup()
278 278
279 279 self.time_label = 'LT' if self.localtime else 'UTC'
280 280
281 281 if self.width is None:
282 282 self.width = 8
283 283
284 284 self.figures = []
285 285 self.axes = []
286 286 self.cb_axes = []
287 287 self.pf_axes = []
288 288 self.cmaps = []
289 289
290 290 size = '15%' if self.ncols == 1 else '30%'
291 291 pad = '4%' if self.ncols == 1 else '8%'
292 292
293 293 if self.oneFigure:
294 294 if self.height is None:
295 295 self.height = 1.4 * self.nrows + 1
296 296 fig = plt.figure(figsize=(self.width, self.height),
297 297 edgecolor='k',
298 298 facecolor='w')
299 299 self.figures.append(fig)
300 300 for n in range(self.nplots):
301 301 ax = fig.add_subplot(self.nrows, self.ncols,
302 302 n + 1, polar=self.polar)
303 303 ax.tick_params(labelsize=8)
304 304 ax.firsttime = True
305 305 ax.index = 0
306 306 ax.press = None
307 307 self.axes.append(ax)
308 308 if self.showprofile:
309 309 cax = self.__add_axes(ax, size=size, pad=pad)
310 310 cax.tick_params(labelsize=8)
311 311 self.pf_axes.append(cax)
312 312 else:
313 313 if self.height is None:
314 314 self.height = 3
315 315 for n in range(self.nplots):
316 316 fig = plt.figure(figsize=(self.width, self.height),
317 317 edgecolor='k',
318 318 facecolor='w')
319 319 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
320 320 ax.tick_params(labelsize=8)
321 321 ax.firsttime = True
322 322 ax.index = 0
323 323 ax.press = None
324 324 self.figures.append(fig)
325 325 self.axes.append(ax)
326 326 if self.showprofile:
327 327 cax = self.__add_axes(ax, size=size, pad=pad)
328 328 cax.tick_params(labelsize=8)
329 329 self.pf_axes.append(cax)
330 330
331 331 for n in range(self.nrows):
332 332 if self.colormaps is not None:
333 333 cmap = plt.get_cmap(self.colormaps[n])
334 334 else:
335 335 cmap = plt.get_cmap(self.colormap)
336 336 cmap.set_bad(self.bgcolor, 1.)
337 337 self.cmaps.append(cmap)
338 338
339 339 def __add_axes(self, ax, size='30%', pad='8%'):
340 340 '''
341 341 Add new axes to the given figure
342 342 '''
343 343 divider = make_axes_locatable(ax)
344 344 nax = divider.new_horizontal(size=size, pad=pad)
345 345 ax.figure.add_axes(nax)
346 346 return nax
347 347
348 348 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
349 349 '''
350 350 Create a masked array for missing data
351 351 '''
352 352 if x_buffer.shape[0] < 2:
353 353 return x_buffer, y_buffer, z_buffer
354 354
355 355 deltas = x_buffer[1:] - x_buffer[0:-1]
356 356 x_median = numpy.median(deltas)
357 357
358 358 index = numpy.where(deltas > 5 * x_median)
359 359
360 360 if len(index[0]) != 0:
361 361 z_buffer[::, index[0], ::] = self.__missing
362 362 z_buffer = numpy.ma.masked_inside(z_buffer,
363 363 0.99 * self.__missing,
364 364 1.01 * self.__missing)
365 365
366 366 return x_buffer, y_buffer, z_buffer
367 367
368 368 def decimate(self):
369 369
370 370 # dx = int(len(self.x)/self.__MAXNUMX) + 1
371 371 dy = int(len(self.y) / self.decimation) + 1
372 372
373 373 # x = self.x[::dx]
374 374 x = self.x
375 375 y = self.y[::dy]
376 376 z = self.z[::, ::, ::dy]
377 377
378 378 return x, y, z
379 379
380 380 def format(self):
381 381 '''
382 382 Set min and max values, labels, ticks and titles
383 383 '''
384 384
385 385 for n, ax in enumerate(self.axes):
386 386 if ax.firsttime:
387 387 if self.xaxis != 'time':
388 388 xmin = self.xmin
389 389 xmax = self.xmax
390 390 else:
391 391 xmin = self.tmin
392 392 xmax = self.tmin + self.xrange*60*60
393 393 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
394 394 ax.xaxis.set_major_locator(LinearLocator(9))
395 395 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
396 396 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
397 397 ax.set_facecolor(self.bgcolor)
398 398 if self.xscale:
399 399 ax.xaxis.set_major_formatter(FuncFormatter(
400 400 lambda x, pos: '{0:g}'.format(x*self.xscale)))
401 401 if self.yscale:
402 402 ax.yaxis.set_major_formatter(FuncFormatter(
403 403 lambda x, pos: '{0:g}'.format(x*self.yscale)))
404 404 if self.xlabel is not None:
405 405 ax.set_xlabel(self.xlabel)
406 406 if self.ylabel is not None:
407 407 ax.set_ylabel(self.ylabel)
408 408 if self.showprofile:
409 409 self.pf_axes[n].set_ylim(ymin, ymax)
410 410 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
411 411 self.pf_axes[n].set_xlabel('dB')
412 412 self.pf_axes[n].grid(b=True, axis='x')
413 413 [tick.set_visible(False)
414 414 for tick in self.pf_axes[n].get_yticklabels()]
415 415 if self.colorbar:
416 416 ax.cbar = plt.colorbar(
417 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
417 ax.plt, ax=ax, fraction=0.05, pad=0.06, aspect=10)
418 418 ax.cbar.ax.tick_params(labelsize=8)
419 419 ax.cbar.ax.press = None
420 420 if self.cb_label:
421 421 ax.cbar.set_label(self.cb_label, size=8)
422 422 elif self.cb_labels:
423 423 ax.cbar.set_label(self.cb_labels[n], size=8)
424 424 else:
425 425 ax.cbar = None
426 426 ax.set_xlim(xmin, xmax)
427 427 ax.set_ylim(ymin, ymax)
428 428 ax.firsttime = False
429 429 if self.grid:
430 430 ax.grid(True)
431 431 if not self.polar:
432 432 ax.set_title('{} {} {}'.format(
433 433 self.titles[n],
434 434 self.getDateTime(self.data.max_time).strftime(
435 435 '%Y-%m-%d %H:%M:%S'),
436 436 self.time_label),
437 437 size=8)
438 438 else:
439 439 #ax.set_title('{}'.format(self.titles[n]), size=8)
440 440 ax.set_title('{} {} {}'.format(
441 441 self.titles[n],
442 442 self.getDateTime(self.data.max_time).strftime(
443 443 '%Y-%m-%d %H:%M:%S'),
444 444 self.time_label),
445 445 size=8)
446 446 ax.set_ylim(0, self.ymax)
447 447 #ax.set_yticks(numpy.arange(0, self.ymax, 20))
448 ax.yaxis.labelpad = 20
448 ax.yaxis.labelpad = 28
449 449
450 450 if self.firsttime:
451 451 for n, fig in enumerate(self.figures):
452 452 fig.subplots_adjust(**self.plots_adjust)
453 453 self.firsttime = False
454 454
455 455 def clear_figures(self):
456 456 '''
457 457 Reset axes for redraw plots
458 458 '''
459 459
460 460 for ax in self.axes+self.pf_axes+self.cb_axes:
461 461 ax.clear()
462 462 ax.firsttime = True
463 463 if hasattr(ax, 'cbar') and ax.cbar:
464 464 ax.cbar.remove()
465 465
466 466 def __plot(self):
467 467 '''
468 468 Main function to plot, format and save figures
469 469 '''
470 470
471 471 self.plot()
472 472 self.format()
473 473
474 474 for n, fig in enumerate(self.figures):
475 475 if self.nrows == 0 or self.nplots == 0:
476 476 log.warning('No data', self.name)
477 477 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
478 478 fig.canvas.manager.set_window_title(self.CODE)
479 479 continue
480 480
481 481 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
482 482 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
483 483 fig.canvas.draw()
484 484 if self.show:
485 485 fig.show()
486 486 figpause(0.01)
487 487
488 488 if self.save:
489 489 self.save_figure(n)
490 490
491 491 if self.server:
492 492 self.send_to_server()
493 493
494 494 def __update(self, dataOut, timestamp):
495 495 '''
496 496 '''
497 497
498 498 metadata = {
499 499 'yrange': dataOut.heightList,
500 500 'interval': dataOut.timeInterval,
501 501 'channels': dataOut.channelList
502 502 }
503 503
504 504 data, meta = self.update(dataOut)
505 505 metadata.update(meta)
506 506 self.data.update(data, timestamp, metadata)
507 507
508 508 def save_figure(self, n):
509 509 '''
510 510 '''
511 511 if self.oneFigure:
512 512 if (self.data.max_time - self.save_time) <= self.save_period:
513 513 return
514 514
515 515 self.save_time = self.data.max_time
516 516
517 517 fig = self.figures[n]
518
518
519 519 if self.throttle == 0:
520 520 if self.oneFigure:
521 521 figname = os.path.join(
522 522 self.save,
523 523 self.save_code,
524 524 '{}_{}.png'.format(
525 525 self.save_code,
526 526 self.getDateTime(self.data.max_time).strftime(
527 527 '%Y%m%d_%H%M%S'
528 528 ),
529 529 )
530 530 )
531 531 else:
532 532 figname = os.path.join(
533 533 self.save,
534 534 self.save_code,
535 535 '{}_ch{}_{}.png'.format(
536 536 self.save_code,n,
537 537 self.getDateTime(self.data.max_time).strftime(
538 538 '%Y%m%d_%H%M%S'
539 539 ),
540 540 )
541 541 )
542 542 log.log('Saving figure: {}'.format(figname), self.name)
543 543 if not os.path.isdir(os.path.dirname(figname)):
544 544 os.makedirs(os.path.dirname(figname))
545 545 fig.savefig(figname)
546 546
547 547 figname = os.path.join(
548 548 self.save,
549 549 '{}_{}.png'.format(
550 550 self.save_code,
551 551 self.getDateTime(self.data.min_time).strftime(
552 552 '%Y%m%d'
553 553 ),
554 554 )
555 555 )
556 556
557 557 log.log('Saving figure: {}'.format(figname), self.name)
558 558 if not os.path.isdir(os.path.dirname(figname)):
559 559 os.makedirs(os.path.dirname(figname))
560 560 fig.savefig(figname)
561 561
562 562 def send_to_server(self):
563 563 '''
564 564 '''
565 565
566 566 if self.exp_code == None:
567 567 log.warning('Missing `exp_code` skipping sending to server...')
568 568
569 569 last_time = self.data.max_time
570 570 interval = last_time - self.sender_time
571 571 if interval < self.sender_period:
572 572 return
573 573
574 574 self.sender_time = last_time
575 575
576 576 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
577 577 for attr in attrs:
578 578 value = getattr(self, attr)
579 579 if value:
580 580 if isinstance(value, (numpy.float32, numpy.float64)):
581 581 value = round(float(value), 2)
582 582 self.data.meta[attr] = value
583 583 if self.colormap == 'jet':
584 584 self.data.meta['colormap'] = 'Jet'
585 585 elif 'RdBu' in self.colormap:
586 586 self.data.meta['colormap'] = 'RdBu'
587 587 else:
588 588 self.data.meta['colormap'] = 'Viridis'
589 589 self.data.meta['interval'] = int(interval)
590 590
591 591 self.sender_queue.append(last_time)
592 592
593 593 while True:
594 594 try:
595 595 tm = self.sender_queue.popleft()
596 596 except IndexError:
597 597 break
598 598 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
599 599 self.socket.send_string(msg)
600 600 socks = dict(self.poll.poll(2000))
601 601 if socks.get(self.socket) == zmq.POLLIN:
602 602 reply = self.socket.recv_string()
603 603 if reply == 'ok':
604 604 log.log("Response from server ok", self.name)
605 605 time.sleep(0.1)
606 606 continue
607 607 else:
608 608 log.warning(
609 609 "Malformed reply from server: {}".format(reply), self.name)
610 610 else:
611 611 log.warning(
612 612 "No response from server, retrying...", self.name)
613 613 self.sender_queue.appendleft(tm)
614 614 self.socket.setsockopt(zmq.LINGER, 0)
615 615 self.socket.close()
616 616 self.poll.unregister(self.socket)
617 617 self.socket = self.context.socket(zmq.REQ)
618 618 self.socket.connect(self.server)
619 619 self.poll.register(self.socket, zmq.POLLIN)
620 620 break
621 621
622 622 def setup(self):
623 623 '''
624 624 This method should be implemented in the child class, the following
625 625 attributes should be set:
626 626
627 627 self.nrows: number of rows
628 628 self.ncols: number of cols
629 629 self.nplots: number of plots (channels or pairs)
630 630 self.ylabel: label for Y axes
631 631 self.titles: list of axes title
632 632
633 633 '''
634 634 raise NotImplementedError
635 635
636 636 def plot(self):
637 637 '''
638 638 Must be defined in the child class, the actual plotting method
639 639 '''
640 640 raise NotImplementedError
641 641
642 642 def update(self, dataOut):
643 643 '''
644 644 Must be defined in the child class, update self.data with new data
645 645 '''
646 646
647 647 data = {
648 648 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
649 649 }
650 650 meta = {}
651 651
652 652 return data, meta
653 653
654 654 def run(self, dataOut, **kwargs):
655 655 '''
656 656 Main plotting routine
657 657 '''
658 658
659 659 if self.isConfig is False:
660 660 self.__setup(**kwargs)
661 661
662 662 if self.localtime:
663 663 self.getDateTime = datetime.datetime.fromtimestamp
664 664 else:
665 665 self.getDateTime = datetime.datetime.utcfromtimestamp
666 666
667 667 self.data.setup()
668 668 self.isConfig = True
669 669 if self.server:
670 670 self.context = zmq.Context()
671 671 self.socket = self.context.socket(zmq.REQ)
672 672 self.socket.connect(self.server)
673 673 self.poll = zmq.Poller()
674 674 self.poll.register(self.socket, zmq.POLLIN)
675 675
676 676 tm = getattr(dataOut, self.attr_time)
677 677
678 678 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
679 679 self.save_time = tm
680 680 self.__plot()
681 681 self.tmin += self.xrange*60*60
682 682 self.data.setup()
683 683 self.clear_figures()
684 684
685 685 self.__update(dataOut, tm)
686 686
687 687 if self.isPlotConfig is False:
688 688 self.__setup_plot()
689 689 self.isPlotConfig = True
690 690 if self.xaxis == 'time':
691 691 dt = self.getDateTime(tm)
692 692 if self.xmin is None:
693 693 self.tmin = tm
694 694 self.xmin = dt.hour
695 695 minutes = (self.xmin-int(self.xmin)) * 60
696 696 seconds = (minutes - int(minutes)) * 60
697 697 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
698 698 datetime.datetime(1970, 1, 1)).total_seconds()
699 699 if self.localtime:
700 700 self.tmin += time.timezone
701 701
702 702 if self.xmin is not None and self.xmax is not None:
703 703 self.xrange = self.xmax - self.xmin
704 704
705 705 if self.throttle == 0:
706 706 self.__plot()
707 707 else:
708 708 self.__throttle_plot(self.__plot)#, coerce=coerce)
709 709
710 710 def close(self):
711 711
712 712 if self.data and not self.data.flagNoData:
713 713 self.save_time = 0
714 714 self.__plot()
715 715 if self.data and not self.data.flagNoData and self.pause:
716 716 figpause(10)
@@ -1,1936 +1,1936
1 1 import os
2 2 import datetime
3 3 import numpy
4 4 from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
5 5
6 6 from schainpy.model.graphics.jroplot_base import Plot, plt
7 7 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
8 8 from schainpy.utils import log
9 9 # libreria wradlib
10 10 #import wradlib as wrl
11 11
12 12 EARTH_RADIUS = 6.3710e3
13 13
14 14
15 15 def ll2xy(lat1, lon1, lat2, lon2):
16 16
17 17 p = 0.017453292519943295
18 18 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
19 19 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
20 20 r = 12742 * numpy.arcsin(numpy.sqrt(a))
21 21 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
22 22 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
23 23 theta = -theta + numpy.pi/2
24 24 return r*numpy.cos(theta), r*numpy.sin(theta)
25 25
26 26
27 27 def km2deg(km):
28 28 '''
29 29 Convert distance in km to degrees
30 30 '''
31 31
32 32 return numpy.rad2deg(km/EARTH_RADIUS)
33 33
34 34
35 35
36 36 class SpectralMomentsPlot(SpectraPlot):
37 37 '''
38 38 Plot for Spectral Moments
39 39 '''
40 40 CODE = 'spc_moments'
41 41 # colormap = 'jet'
42 42 # plot_type = 'pcolor'
43 43
44 44 class DobleGaussianPlot(SpectraPlot):
45 45 '''
46 46 Plot for Double Gaussian Plot
47 47 '''
48 48 CODE = 'gaussian_fit'
49 49 # colormap = 'jet'
50 50 # plot_type = 'pcolor'
51 51
52 52 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
53 53 '''
54 54 Plot SpectraCut with Double Gaussian Fit
55 55 '''
56 56 CODE = 'cut_gaussian_fit'
57 57
58 58 class SnrPlot(RTIPlot):
59 59 '''
60 60 Plot for SNR Data
61 61 '''
62 62
63 63 CODE = 'snr'
64 64 colormap = 'jet'
65 65
66 66 def update(self, dataOut):
67 67
68 68 data = {
69 69 'snr': 10*numpy.log10(dataOut.data_snr)
70 70 }
71 71
72 72 return data, {}
73 73
74 74 class DopplerPlot(RTIPlot):
75 75 '''
76 76 Plot for DOPPLER Data (1st moment)
77 77 '''
78 78
79 79 CODE = 'dop'
80 80 colormap = 'jet'
81 81
82 82 def update(self, dataOut):
83 83
84 84 data = {
85 85 'dop': 10*numpy.log10(dataOut.data_dop)
86 86 }
87 87
88 88 return data, {}
89 89
90 90 class PowerPlot(RTIPlot):
91 91 '''
92 92 Plot for Power Data (0 moment)
93 93 '''
94 94
95 95 CODE = 'pow'
96 96 colormap = 'jet'
97 97
98 98 def update(self, dataOut):
99 99 data = {
100 100 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
101 101 }
102 102 return data, {}
103 103
104 104 class SpectralWidthPlot(RTIPlot):
105 105 '''
106 106 Plot for Spectral Width Data (2nd moment)
107 107 '''
108 108
109 109 CODE = 'width'
110 110 colormap = 'jet'
111 111
112 112 def update(self, dataOut):
113 113
114 114 data = {
115 115 'width': dataOut.data_width
116 116 }
117 117
118 118 return data, {}
119 119
120 120 class SkyMapPlot(Plot):
121 121 '''
122 122 Plot for meteors detection data
123 123 '''
124 124
125 125 CODE = 'param'
126 126
127 127 def setup(self):
128 128
129 129 self.ncols = 1
130 130 self.nrows = 1
131 131 self.width = 7.2
132 132 self.height = 7.2
133 133 self.nplots = 1
134 134 self.xlabel = 'Zonal Zenith Angle (deg)'
135 135 self.ylabel = 'Meridional Zenith Angle (deg)'
136 136 self.polar = True
137 137 self.ymin = -180
138 138 self.ymax = 180
139 139 self.colorbar = False
140 140
141 141 def plot(self):
142 142
143 143 arrayParameters = numpy.concatenate(self.data['param'])
144 144 error = arrayParameters[:, -1]
145 145 indValid = numpy.where(error == 0)[0]
146 146 finalMeteor = arrayParameters[indValid, :]
147 147 finalAzimuth = finalMeteor[:, 3]
148 148 finalZenith = finalMeteor[:, 4]
149 149
150 150 x = finalAzimuth * numpy.pi / 180
151 151 y = finalZenith
152 152
153 153 ax = self.axes[0]
154 154
155 155 if ax.firsttime:
156 156 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
157 157 else:
158 158 ax.plot.set_data(x, y)
159 159
160 160 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
161 161 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
162 162 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
163 163 dt2,
164 164 len(x))
165 165 self.titles[0] = title
166 166
167 167
168 168 class GenericRTIPlot(Plot):
169 169 '''
170 170 Plot for data_xxxx object
171 171 '''
172 172
173 173 CODE = 'param'
174 174 colormap = 'viridis'
175 175 plot_type = 'pcolorbuffer'
176 176
177 177 def setup(self):
178 178 self.xaxis = 'time'
179 179 self.ncols = 1
180 180 self.nrows = self.data.shape('param')[0]
181 181 self.nplots = self.nrows
182 182 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
183 183
184 184 if not self.xlabel:
185 185 self.xlabel = 'Time'
186 186
187 187 self.ylabel = 'Range [km]'
188 188 if not self.titles:
189 189 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
190 190
191 191 def update(self, dataOut):
192 192
193 193 data = {
194 194 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
195 195 }
196 196
197 197 meta = {}
198 198
199 199 return data, meta
200 200
201 201 def plot(self):
202 202 # self.data.normalize_heights()
203 203 self.x = self.data.times
204 204 self.y = self.data.yrange
205 205 self.z = self.data['param']
206 206 self.z = 10*numpy.log10(self.z)
207 207 self.z = numpy.ma.masked_invalid(self.z)
208 208
209 209 if self.decimation is None:
210 210 x, y, z = self.fill_gaps(self.x, self.y, self.z)
211 211 else:
212 212 x, y, z = self.fill_gaps(*self.decimate())
213 213
214 214 for n, ax in enumerate(self.axes):
215 215
216 216 self.zmax = self.zmax if self.zmax is not None else numpy.max(
217 217 self.z[n])
218 218 self.zmin = self.zmin if self.zmin is not None else numpy.min(
219 219 self.z[n])
220 220
221 221 if ax.firsttime:
222 222 if self.zlimits is not None:
223 223 self.zmin, self.zmax = self.zlimits[n]
224 224
225 225 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
226 226 vmin=self.zmin,
227 227 vmax=self.zmax,
228 228 cmap=self.cmaps[n]
229 229 )
230 230 else:
231 231 if self.zlimits is not None:
232 232 self.zmin, self.zmax = self.zlimits[n]
233 233 ax.collections.remove(ax.collections[0])
234 234 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
235 235 vmin=self.zmin,
236 236 vmax=self.zmax,
237 237 cmap=self.cmaps[n]
238 238 )
239 239
240 240
241 241 class PolarMapPlot(Plot):
242 242 '''
243 243 Plot for weather radar
244 244 '''
245 245
246 246 CODE = 'param'
247 247 colormap = 'seismic'
248 248
249 249 def setup(self):
250 250 self.ncols = 1
251 251 self.nrows = 1
252 252 self.width = 9
253 253 self.height = 8
254 254 self.mode = self.data.meta['mode']
255 255 if self.channels is not None:
256 256 self.nplots = len(self.channels)
257 257 self.nrows = len(self.channels)
258 258 else:
259 259 self.nplots = self.data.shape(self.CODE)[0]
260 260 self.nrows = self.nplots
261 261 self.channels = list(range(self.nplots))
262 262 if self.mode == 'E':
263 263 self.xlabel = 'Longitude'
264 264 self.ylabel = 'Latitude'
265 265 else:
266 266 self.xlabel = 'Range (km)'
267 267 self.ylabel = 'Height (km)'
268 268 self.bgcolor = 'white'
269 269 self.cb_labels = self.data.meta['units']
270 270 self.lat = self.data.meta['latitude']
271 271 self.lon = self.data.meta['longitude']
272 272 self.xmin, self.xmax = float(
273 273 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
274 274 self.ymin, self.ymax = float(
275 275 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
276 276 # self.polar = True
277 277
278 278 def plot(self):
279 279
280 280 for n, ax in enumerate(self.axes):
281 281 data = self.data['param'][self.channels[n]]
282 282
283 283 zeniths = numpy.linspace(
284 284 0, self.data.meta['max_range'], data.shape[1])
285 285 if self.mode == 'E':
286 286 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
287 287 r, theta = numpy.meshgrid(zeniths, azimuths)
288 288 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
289 289 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
290 290 x = km2deg(x) + self.lon
291 291 y = km2deg(y) + self.lat
292 292 else:
293 293 azimuths = numpy.radians(self.data.yrange)
294 294 r, theta = numpy.meshgrid(zeniths, azimuths)
295 295 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
296 296 self.y = zeniths
297 297
298 298 if ax.firsttime:
299 299 if self.zlimits is not None:
300 300 self.zmin, self.zmax = self.zlimits[n]
301 301 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
302 302 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
303 303 vmin=self.zmin,
304 304 vmax=self.zmax,
305 305 cmap=self.cmaps[n])
306 306 else:
307 307 if self.zlimits is not None:
308 308 self.zmin, self.zmax = self.zlimits[n]
309 309 ax.collections.remove(ax.collections[0])
310 310 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
311 311 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
312 312 vmin=self.zmin,
313 313 vmax=self.zmax,
314 314 cmap=self.cmaps[n])
315 315
316 316 if self.mode == 'A':
317 317 continue
318 318
319 319 # plot district names
320 320 f = open('/data/workspace/schain_scripts/distrito.csv')
321 321 for line in f:
322 322 label, lon, lat = [s.strip() for s in line.split(',') if s]
323 323 lat = float(lat)
324 324 lon = float(lon)
325 325 # ax.plot(lon, lat, '.b', ms=2)
326 326 ax.text(lon, lat, label.decode('utf8'), ha='center',
327 327 va='bottom', size='8', color='black')
328 328
329 329 # plot limites
330 330 limites = []
331 331 tmp = []
332 332 for line in open('/data/workspace/schain_scripts/lima.csv'):
333 333 if '#' in line:
334 334 if tmp:
335 335 limites.append(tmp)
336 336 tmp = []
337 337 continue
338 338 values = line.strip().split(',')
339 339 tmp.append((float(values[0]), float(values[1])))
340 340 for points in limites:
341 341 ax.add_patch(
342 342 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
343 343
344 344 # plot Cuencas
345 345 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
346 346 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
347 347 values = [line.strip().split(',') for line in f]
348 348 points = [(float(s[0]), float(s[1])) for s in values]
349 349 ax.add_patch(Polygon(points, ec='b', fc='none'))
350 350
351 351 # plot grid
352 352 for r in (15, 30, 45, 60):
353 353 ax.add_artist(plt.Circle((self.lon, self.lat),
354 354 km2deg(r), color='0.6', fill=False, lw=0.2))
355 355 ax.text(
356 356 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
357 357 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
358 358 '{}km'.format(r),
359 359 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
360 360
361 361 if self.mode == 'E':
362 362 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
363 363 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
364 364 else:
365 365 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
366 366 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
367 367
368 368 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
369 369 self.titles = ['{} {}'.format(
370 370 self.data.parameters[x], title) for x in self.channels]
371 371
372 372 class WeatherPlot(Plot):
373 373 CODE = 'weather'
374 374 plot_name = 'weather'
375 375 plot_type = 'ppistyle'
376 376 buffering = False
377 377
378 378 def setup(self):
379 379 self.ncols = 1
380 380 self.nrows = 1
381 381 self.width =8
382 382 self.height =8
383 383 self.nplots= 1
384 384 self.ylabel= 'Range [Km]'
385 385 self.titles= ['Weather']
386 386 self.colorbar=False
387 387 self.ini =0
388 388 self.len_azi =0
389 389 self.buffer_ini = None
390 390 self.buffer_azi = None
391 391 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
392 392 self.flag =0
393 393 self.indicador= 0
394 394 self.last_data_azi = None
395 395 self.val_mean = None
396 396
397 397 def update(self, dataOut):
398 398
399 399 data = {}
400 400 meta = {}
401 401 if hasattr(dataOut, 'dataPP_POWER'):
402 402 factor = 1
403 403 if hasattr(dataOut, 'nFFTPoints'):
404 404 factor = dataOut.normFactor
405 405 #print("DIME EL SHAPE PORFAVOR",dataOut.data_360.shape)
406 406 data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor))
407 407 data['azi'] = dataOut.data_azi
408 408 data['ele'] = dataOut.data_ele
409 409 return data, meta
410 410
411 411 def get2List(self,angulos):
412 412 list1=[]
413 413 list2=[]
414 414 for i in reversed(range(len(angulos))):
415 415 diff_ = angulos[i]-angulos[i-1]
416 416 if diff_ >1.5:
417 417 list1.append(i-1)
418 418 list2.append(diff_)
419 419 return list(reversed(list1)),list(reversed(list2))
420 420
421 421 def fixData360(self,list_,ang_):
422 422 if list_[0]==-1:
423 423 vec = numpy.where(ang_<ang_[0])
424 424 ang_[vec] = ang_[vec]+360
425 425 return ang_
426 426 return ang_
427 427
428 428 def fixData360HL(self,angulos):
429 429 vec = numpy.where(angulos>=360)
430 430 angulos[vec]=angulos[vec]-360
431 431 return angulos
432 432
433 433 def search_pos(self,pos,list_):
434 434 for i in range(len(list_)):
435 435 if pos == list_[i]:
436 436 return True,i
437 437 i=None
438 438 return False,i
439 439
440 440 def fixDataComp(self,ang_,list1_,list2_):
441 441 size = len(ang_)
442 442 size2 = 0
443 443 for i in range(len(list2_)):
444 444 size2=size2+round(list2_[i])-1
445 445 new_size= size+size2
446 446 ang_new = numpy.zeros(new_size)
447 447 ang_new2 = numpy.zeros(new_size)
448 448
449 449 tmp = 0
450 450 c = 0
451 451 for i in range(len(ang_)):
452 452 ang_new[tmp +c] = ang_[i]
453 453 ang_new2[tmp+c] = ang_[i]
454 454 condition , value = self.search_pos(i,list1_)
455 455 if condition:
456 456 pos = tmp + c + 1
457 457 for k in range(round(list2_[value])-1):
458 458 ang_new[pos+k] = ang_new[pos+k-1]+1
459 459 ang_new2[pos+k] = numpy.nan
460 460 tmp = pos +k
461 461 c = 0
462 462 c=c+1
463 463 return ang_new,ang_new2
464 464
465 465 def globalCheckPED(self,angulos):
466 466 l1,l2 = self.get2List(angulos)
467 467 if len(l1)>0:
468 468 angulos2 = self.fixData360(list_=l1,ang_=angulos)
469 469 l1,l2 = self.get2List(angulos2)
470 470
471 471 ang1_,ang2_ = self.fixDataComp(ang_=angulos2,list1_=l1,list2_=l2)
472 472 ang1_ = self.fixData360HL(ang1_)
473 473 ang2_ = self.fixData360HL(ang2_)
474 474 else:
475 475 ang1_= angulos
476 476 ang2_= angulos
477 477 return ang1_,ang2_
478 478
479 479 def analizeDATA(self,data_azi):
480 480 list1 = []
481 481 list2 = []
482 482 dat = data_azi
483 483 for i in reversed(range(1,len(dat))):
484 484 if dat[i]>dat[i-1]:
485 485 diff = int(dat[i])-int(dat[i-1])
486 486 else:
487 487 diff = 360+int(dat[i])-int(dat[i-1])
488 488 if diff > 1:
489 489 list1.append(i-1)
490 490 list2.append(diff-1)
491 491 return list1,list2
492 492
493 493 def fixDATANEW(self,data_azi,data_weather):
494 494 list1,list2 = self.analizeDATA(data_azi)
495 495 if len(list1)== 0:
496 496 return data_azi,data_weather
497 497 else:
498 498 resize = 0
499 499 for i in range(len(list2)):
500 500 resize= resize + list2[i]
501 501 new_data_azi = numpy.resize(data_azi,resize)
502 502 new_data_weather= numpy.resize(date_weather,resize)
503 503
504 504 for i in range(len(list2)):
505 505 j=0
506 506 position=list1[i]+1
507 507 for j in range(list2[i]):
508 508 new_data_azi[position+j]=new_data_azi[position+j-1]+1
509 509 return new_data_azi
510 510
511 511 def fixDATA(self,data_azi):
512 512 data=data_azi
513 513 for i in range(len(data)):
514 514 if numpy.isnan(data[i]):
515 515 data[i]=data[i-1]+1
516 516 return data
517 517
518 518 def replaceNAN(self,data_weather,data_azi,val):
519 519 data= data_azi
520 520 data_T= data_weather
521 521 if data.shape[0]> data_T.shape[0]:
522 522 data_N = numpy.ones( [data.shape[0],data_T.shape[1]])
523 523 c = 0
524 524 for i in range(len(data)):
525 525 if numpy.isnan(data[i]):
526 526 data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
527 527 else:
528 528 data_N[i,:]=data_T[c,:]
529 529 c=c+1
530 530 return data_N
531 531 else:
532 532 for i in range(len(data)):
533 533 if numpy.isnan(data[i]):
534 534 data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
535 535 return data_T
536 536
537 537 def const_ploteo(self,data_weather,data_azi,step,res):
538 538 if self.ini==0:
539 539 #-------
540 540 n = (360/res)-len(data_azi)
541 541 #--------------------- new -------------------------
542 542 data_azi_new ,data_azi_old= self.globalCheckPED(data_azi)
543 543 #------------------------
544 544 start = data_azi_new[-1] + res
545 545 end = data_azi_new[0] - res
546 546 #------ new
547 547 self.last_data_azi = end
548 548 if start>end:
549 549 end = end + 360
550 550 azi_vacia = numpy.linspace(start,end,int(n))
551 551 azi_vacia = numpy.where(azi_vacia>360,azi_vacia-360,azi_vacia)
552 552 data_azi = numpy.hstack((data_azi_new,azi_vacia))
553 553 # RADAR
554 554 val_mean = numpy.mean(data_weather[:,-1])
555 555 self.val_mean = val_mean
556 556 data_weather_cmp = numpy.ones([(360-data_weather.shape[0]),data_weather.shape[1]])*val_mean
557 557 data_weather = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean)
558 558 data_weather = numpy.vstack((data_weather,data_weather_cmp))
559 559 else:
560 560 # azimuth
561 561 flag=0
562 562 start_azi = self.res_azi[0]
563 563 #-----------new------------
564 564 data_azi ,data_azi_old= self.globalCheckPED(data_azi)
565 565 data_weather = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean)
566 566 #--------------------------
567 567 start = data_azi[0]
568 568 end = data_azi[-1]
569 569 self.last_data_azi= end
570 570 if start< start_azi:
571 571 start = start +360
572 572 if end <start_azi:
573 573 end = end +360
574 574
575 575 pos_ini = int((start-start_azi)/res)
576 576 len_azi = len(data_azi)
577 577 if (360-pos_ini)<len_azi:
578 578 if pos_ini+1==360:
579 579 pos_ini=0
580 580 else:
581 581 flag=1
582 582 dif= 360-pos_ini
583 583 comp= len_azi-dif
584 584 #-----------------
585 585 if flag==0:
586 586 # AZIMUTH
587 587 self.res_azi[pos_ini:pos_ini+len_azi] = data_azi
588 588 # RADAR
589 589 self.res_weather[pos_ini:pos_ini+len_azi,:] = data_weather
590 590 else:
591 591 # AZIMUTH
592 592 self.res_azi[pos_ini:pos_ini+dif] = data_azi[0:dif]
593 593 self.res_azi[0:comp] = data_azi[dif:]
594 594 # RADAR
595 595 self.res_weather[pos_ini:pos_ini+dif,:] = data_weather[0:dif,:]
596 596 self.res_weather[0:comp,:] = data_weather[dif:,:]
597 597 flag=0
598 598 data_azi = self.res_azi
599 599 data_weather = self.res_weather
600 600
601 601 return data_weather,data_azi
602 602
603 603 def plot(self):
604 604 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
605 605 data = self.data[-1]
606 606 r = self.data.yrange
607 607 delta_height = r[1]-r[0]
608 608 r_mask = numpy.where(r>=0)[0]
609 609 r = numpy.arange(len(r_mask))*delta_height
610 610 self.y = 2*r
611 611 # RADAR
612 612 #data_weather = data['weather']
613 613 # PEDESTAL
614 614 #data_azi = data['azi']
615 615 res = 1
616 616 # STEP
617 617 step = (360/(res*data['weather'].shape[0]))
618 618
619 619 self.res_weather, self.res_azi = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_azi=data['azi'],step=step,res=res)
620 620 self.res_ele = numpy.mean(data['ele'])
621 621 ################# PLOTEO ###################
622 622 for i,ax in enumerate(self.axes):
623 623 self.zmin = self.zmin if self.zmin else 20
624 624 self.zmax = self.zmax if self.zmax else 80
625 625 if ax.firsttime:
626 626 plt.clf()
627 627 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=self.zmin, vmax=self.zmax)
628 628 else:
629 629 plt.clf()
630 630 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=self.zmin, vmax=self.zmax)
631 631 caax = cgax.parasites[0]
632 632 paax = cgax.parasites[1]
633 633 cbar = plt.gcf().colorbar(pm, pad=0.075)
634 634 caax.set_xlabel('x_range [km]')
635 635 caax.set_ylabel('y_range [km]')
636 636 plt.text(1.0, 1.05, 'Azimuth '+str(thisDatetime)+" Step "+str(self.ini)+ " EL: "+str(round(self.res_ele, 1)), transform=caax.transAxes, va='bottom',ha='right')
637 637
638 638 self.ini= self.ini+1
639 639
640 640
641 641 class WeatherRHIPlot(Plot):
642 642 CODE = 'weather'
643 643 plot_name = 'weather'
644 644 plot_type = 'rhistyle'
645 645 buffering = False
646 646 data_ele_tmp = None
647 647
648 648 def setup(self):
649 649 print("********************")
650 650 print("********************")
651 651 print("********************")
652 652 print("SETUP WEATHER PLOT")
653 653 self.ncols = 1
654 654 self.nrows = 1
655 655 self.nplots= 1
656 656 self.ylabel= 'Range [Km]'
657 657 self.titles= ['Weather']
658 658 if self.channels is not None:
659 659 self.nplots = len(self.channels)
660 660 self.nrows = len(self.channels)
661 661 else:
662 662 self.nplots = self.data.shape(self.CODE)[0]
663 663 self.nrows = self.nplots
664 664 self.channels = list(range(self.nplots))
665 665 print("channels",self.channels)
666 666 print("que saldra", self.data.shape(self.CODE)[0])
667 667 self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)]
668 668 print("self.titles",self.titles)
669 669 self.colorbar=False
670 670 self.width =12
671 671 self.height =8
672 672 self.ini =0
673 673 self.len_azi =0
674 674 self.buffer_ini = None
675 675 self.buffer_ele = None
676 676 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
677 677 self.flag =0
678 678 self.indicador= 0
679 679 self.last_data_ele = None
680 680 self.val_mean = None
681 681
682 682 def update(self, dataOut):
683 683
684 684 data = {}
685 685 meta = {}
686 686 if hasattr(dataOut, 'dataPP_POWER'):
687 687 factor = 1
688 688 if hasattr(dataOut, 'nFFTPoints'):
689 689 factor = dataOut.normFactor
690 690 print("dataOut",dataOut.data_360.shape)
691 691 #
692 692 data['weather'] = 10*numpy.log10(dataOut.data_360/(factor))
693 693 #
694 694 #data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor))
695 695 data['azi'] = dataOut.data_azi
696 696 data['ele'] = dataOut.data_ele
697 697 #print("UPDATE")
698 698 #print("data[weather]",data['weather'].shape)
699 699 #print("data[azi]",data['azi'])
700 700 return data, meta
701 701
702 702 def get2List(self,angulos):
703 703 list1=[]
704 704 list2=[]
705 705 for i in reversed(range(len(angulos))):
706 706 if not i==0:#el caso de i=0 evalula el primero de la lista con el ultimo y no es relevante
707 707 diff_ = angulos[i]-angulos[i-1]
708 708 if abs(diff_) >1.5:
709 709 list1.append(i-1)
710 710 list2.append(diff_)
711 711 return list(reversed(list1)),list(reversed(list2))
712 712
713 713 def fixData90(self,list_,ang_):
714 714 if list_[0]==-1:
715 715 vec = numpy.where(ang_<ang_[0])
716 716 ang_[vec] = ang_[vec]+90
717 717 return ang_
718 718 return ang_
719 719
720 720 def fixData90HL(self,angulos):
721 721 vec = numpy.where(angulos>=90)
722 722 angulos[vec]=angulos[vec]-90
723 723 return angulos
724 724
725 725
726 726 def search_pos(self,pos,list_):
727 727 for i in range(len(list_)):
728 728 if pos == list_[i]:
729 729 return True,i
730 730 i=None
731 731 return False,i
732 732
733 733 def fixDataComp(self,ang_,list1_,list2_,tipo_case):
734 734 size = len(ang_)
735 735 size2 = 0
736 736 for i in range(len(list2_)):
737 737 size2=size2+round(abs(list2_[i]))-1
738 738 new_size= size+size2
739 739 ang_new = numpy.zeros(new_size)
740 740 ang_new2 = numpy.zeros(new_size)
741 741
742 742 tmp = 0
743 743 c = 0
744 744 for i in range(len(ang_)):
745 745 ang_new[tmp +c] = ang_[i]
746 746 ang_new2[tmp+c] = ang_[i]
747 747 condition , value = self.search_pos(i,list1_)
748 748 if condition:
749 749 pos = tmp + c + 1
750 750 for k in range(round(abs(list2_[value]))-1):
751 751 if tipo_case==0 or tipo_case==3:#subida
752 752 ang_new[pos+k] = ang_new[pos+k-1]+1
753 753 ang_new2[pos+k] = numpy.nan
754 754 elif tipo_case==1 or tipo_case==2:#bajada
755 755 ang_new[pos+k] = ang_new[pos+k-1]-1
756 756 ang_new2[pos+k] = numpy.nan
757 757
758 758 tmp = pos +k
759 759 c = 0
760 760 c=c+1
761 761 return ang_new,ang_new2
762 762
763 763 def globalCheckPED(self,angulos,tipo_case):
764 764 l1,l2 = self.get2List(angulos)
765 765 ##print("l1",l1)
766 766 ##print("l2",l2)
767 767 if len(l1)>0:
768 768 #angulos2 = self.fixData90(list_=l1,ang_=angulos)
769 769 #l1,l2 = self.get2List(angulos2)
770 770 ang1_,ang2_ = self.fixDataComp(ang_=angulos,list1_=l1,list2_=l2,tipo_case=tipo_case)
771 771 #ang1_ = self.fixData90HL(ang1_)
772 772 #ang2_ = self.fixData90HL(ang2_)
773 773 else:
774 774 ang1_= angulos
775 775 ang2_= angulos
776 776 return ang1_,ang2_
777 777
778 778
779 779 def replaceNAN(self,data_weather,data_ele,val):
780 780 data= data_ele
781 781 data_T= data_weather
782 782 if data.shape[0]> data_T.shape[0]:
783 783 data_N = numpy.ones( [data.shape[0],data_T.shape[1]])
784 784 c = 0
785 785 for i in range(len(data)):
786 786 if numpy.isnan(data[i]):
787 787 data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
788 788 else:
789 789 data_N[i,:]=data_T[c,:]
790 790 c=c+1
791 791 return data_N
792 792 else:
793 793 for i in range(len(data)):
794 794 if numpy.isnan(data[i]):
795 795 data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
796 796 return data_T
797 797
798 798 def check_case(self,data_ele,ang_max,ang_min):
799 799 start = data_ele[0]
800 800 end = data_ele[-1]
801 801 number = (end-start)
802 802 len_ang=len(data_ele)
803 803 print("start",start)
804 804 print("end",end)
805 805 print("number",number)
806 806
807 807 print("len_ang",len_ang)
808 808
809 809 #exit(1)
810 810
811 811 if start<end and (round(abs(number)+1)>=len_ang or (numpy.argmin(data_ele)==0)):#caso subida
812 812 return 0
813 813 #elif start>end and (round(abs(number)+1)>=len_ang or(numpy.argmax(data_ele)==0)):#caso bajada
814 814 # return 1
815 815 elif round(abs(number)+1)>=len_ang and (start>end or(numpy.argmax(data_ele)==0)):#caso bajada
816 816 return 1
817 817 elif round(abs(number)+1)<len_ang and data_ele[-2]>data_ele[-1]:# caso BAJADA CAMBIO ANG MAX
818 818 return 2
819 819 elif round(abs(number)+1)<len_ang and data_ele[-2]<data_ele[-1] :# caso SUBIDA CAMBIO ANG MIN
820 820 return 3
821 821
822 822
823 823 def const_ploteo(self,val_ch,data_weather,data_ele,step,res,ang_max,ang_min):
824 824 ang_max= ang_max
825 825 ang_min= ang_min
826 826 data_weather=data_weather
827 827 val_ch=val_ch
828 828 ##print("*********************DATA WEATHER**************************************")
829 829 ##print(data_weather)
830 830 if self.ini==0:
831 831 '''
832 832 print("**********************************************")
833 833 print("**********************************************")
834 834 print("***************ini**************")
835 835 print("**********************************************")
836 836 print("**********************************************")
837 837 '''
838 838 #print("data_ele",data_ele)
839 839 #----------------------------------------------------------
840 840 tipo_case = self.check_case(data_ele,ang_max,ang_min)
841 841 print("check_case",tipo_case)
842 842 #exit(1)
843 843 #--------------------- new -------------------------
844 844 data_ele_new ,data_ele_old= self.globalCheckPED(data_ele,tipo_case)
845 845
846 846 #-------------------------CAMBIOS RHI---------------------------------
847 847 start= ang_min
848 848 end = ang_max
849 849 n= (ang_max-ang_min)/res
850 850 #------ new
851 851 self.start_data_ele = data_ele_new[0]
852 852 self.end_data_ele = data_ele_new[-1]
853 853 if tipo_case==0 or tipo_case==3: # SUBIDA
854 854 n1= round(self.start_data_ele)- start
855 855 n2= end - round(self.end_data_ele)
856 856 print(self.start_data_ele)
857 857 print(self.end_data_ele)
858 858 if n1>0:
859 859 ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1)
860 860 ele1_nan= numpy.ones(n1)*numpy.nan
861 861 data_ele = numpy.hstack((ele1,data_ele_new))
862 862 print("ele1_nan",ele1_nan.shape)
863 863 print("data_ele_old",data_ele_old.shape)
864 864 data_ele_old = numpy.hstack((ele1_nan,data_ele_old))
865 865 if n2>0:
866 866 ele2= numpy.linspace(self.end_data_ele+1,end,n2)
867 867 ele2_nan= numpy.ones(n2)*numpy.nan
868 868 data_ele = numpy.hstack((data_ele,ele2))
869 869 print("ele2_nan",ele2_nan.shape)
870 870 print("data_ele_old",data_ele_old.shape)
871 871 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
872 872
873 873 if tipo_case==1 or tipo_case==2: # BAJADA
874 874 data_ele_new = data_ele_new[::-1] # reversa
875 875 data_ele_old = data_ele_old[::-1]# reversa
876 876 data_weather = data_weather[::-1,:]# reversa
877 877 vec= numpy.where(data_ele_new<ang_max)
878 878 data_ele_new = data_ele_new[vec]
879 879 data_ele_old = data_ele_old[vec]
880 880 data_weather = data_weather[vec[0]]
881 881 vec2= numpy.where(0<data_ele_new)
882 882 data_ele_new = data_ele_new[vec2]
883 883 data_ele_old = data_ele_old[vec2]
884 884 data_weather = data_weather[vec2[0]]
885 885 self.start_data_ele = data_ele_new[0]
886 886 self.end_data_ele = data_ele_new[-1]
887 887
888 888 n1= round(self.start_data_ele)- start
889 889 n2= end - round(self.end_data_ele)-1
890 890 print(self.start_data_ele)
891 891 print(self.end_data_ele)
892 892 if n1>0:
893 893 ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1)
894 894 ele1_nan= numpy.ones(n1)*numpy.nan
895 895 data_ele = numpy.hstack((ele1,data_ele_new))
896 896 data_ele_old = numpy.hstack((ele1_nan,data_ele_old))
897 897 if n2>0:
898 898 ele2= numpy.linspace(self.end_data_ele+1,end,n2)
899 899 ele2_nan= numpy.ones(n2)*numpy.nan
900 900 data_ele = numpy.hstack((data_ele,ele2))
901 901 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
902 902 # RADAR
903 903 # NOTA data_ele y data_weather es la variable que retorna
904 904 val_mean = numpy.mean(data_weather[:,-1])
905 905 self.val_mean = val_mean
906 906 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
907 907 self.data_ele_tmp[val_ch]= data_ele_old
908 908 else:
909 909 #print("**********************************************")
910 910 #print("****************VARIABLE**********************")
911 911 #-------------------------CAMBIOS RHI---------------------------------
912 912 #---------------------------------------------------------------------
913 913 ##print("INPUT data_ele",data_ele)
914 914 flag=0
915 915 start_ele = self.res_ele[0]
916 916 tipo_case = self.check_case(data_ele,ang_max,ang_min)
917 917 #print("TIPO DE DATA",tipo_case)
918 918 #-----------new------------
919 919 data_ele ,data_ele_old = self.globalCheckPED(data_ele,tipo_case)
920 920 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
921 921
922 922 #-------------------------------NEW RHI ITERATIVO-------------------------
923 923
924 924 if tipo_case==0 : # SUBIDA
925 925 vec = numpy.where(data_ele<ang_max)
926 926 data_ele = data_ele[vec]
927 927 data_ele_old = data_ele_old[vec]
928 928 data_weather = data_weather[vec[0]]
929 929
930 930 vec2 = numpy.where(0<data_ele)
931 931 data_ele= data_ele[vec2]
932 932 data_ele_old= data_ele_old[vec2]
933 933 ##print(data_ele_new)
934 934 data_weather= data_weather[vec2[0]]
935 935
936 936 new_i_ele = int(round(data_ele[0]))
937 937 new_f_ele = int(round(data_ele[-1]))
938 938 #print(new_i_ele)
939 939 #print(new_f_ele)
940 940 #print(data_ele,len(data_ele))
941 941 #print(data_ele_old,len(data_ele_old))
942 942 if new_i_ele< 2:
943 943 self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan
944 944 self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean)
945 945 self.data_ele_tmp[val_ch][new_i_ele:new_i_ele+len(data_ele)]=data_ele_old
946 946 self.res_ele[new_i_ele:new_i_ele+len(data_ele)]= data_ele
947 947 self.res_weather[val_ch][new_i_ele:new_i_ele+len(data_ele),:]= data_weather
948 948 data_ele = self.res_ele
949 949 data_weather = self.res_weather[val_ch]
950 950
951 951 elif tipo_case==1 : #BAJADA
952 952 data_ele = data_ele[::-1] # reversa
953 953 data_ele_old = data_ele_old[::-1]# reversa
954 954 data_weather = data_weather[::-1,:]# reversa
955 955 vec= numpy.where(data_ele<ang_max)
956 956 data_ele = data_ele[vec]
957 957 data_ele_old = data_ele_old[vec]
958 958 data_weather = data_weather[vec[0]]
959 959 vec2= numpy.where(0<data_ele)
960 960 data_ele = data_ele[vec2]
961 961 data_ele_old = data_ele_old[vec2]
962 962 data_weather = data_weather[vec2[0]]
963 963
964 964
965 965 new_i_ele = int(round(data_ele[0]))
966 966 new_f_ele = int(round(data_ele[-1]))
967 967 #print(data_ele)
968 968 #print(ang_max)
969 969 #print(data_ele_old)
970 970 if new_i_ele <= 1:
971 971 new_i_ele = 1
972 972 if round(data_ele[-1])>=ang_max-1:
973 973 self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan
974 974 self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean)
975 975 self.data_ele_tmp[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1]=data_ele_old
976 976 self.res_ele[new_i_ele-1:new_i_ele+len(data_ele)-1]= data_ele
977 977 self.res_weather[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1,:]= data_weather
978 978 data_ele = self.res_ele
979 979 data_weather = self.res_weather[val_ch]
980 980
981 981 elif tipo_case==2: #bajada
982 982 vec = numpy.where(data_ele<ang_max)
983 983 data_ele = data_ele[vec]
984 984 data_weather= data_weather[vec[0]]
985 985
986 986 len_vec = len(vec)
987 987 data_ele_new = data_ele[::-1] # reversa
988 988 data_weather = data_weather[::-1,:]
989 989 new_i_ele = int(data_ele_new[0])
990 990 new_f_ele = int(data_ele_new[-1])
991 991
992 992 n1= new_i_ele- ang_min
993 993 n2= ang_max - new_f_ele-1
994 994 if n1>0:
995 995 ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1)
996 996 ele1_nan= numpy.ones(n1)*numpy.nan
997 997 data_ele = numpy.hstack((ele1,data_ele_new))
998 998 data_ele_old = numpy.hstack((ele1_nan,data_ele_new))
999 999 if n2>0:
1000 1000 ele2= numpy.linspace(new_f_ele+1,ang_max,n2)
1001 1001 ele2_nan= numpy.ones(n2)*numpy.nan
1002 1002 data_ele = numpy.hstack((data_ele,ele2))
1003 1003 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1004 1004
1005 1005 self.data_ele_tmp[val_ch] = data_ele_old
1006 1006 self.res_ele = data_ele
1007 1007 self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1008 1008 data_ele = self.res_ele
1009 1009 data_weather = self.res_weather[val_ch]
1010 1010
1011 1011 elif tipo_case==3:#subida
1012 1012 vec = numpy.where(0<data_ele)
1013 1013 data_ele= data_ele[vec]
1014 1014 data_ele_new = data_ele
1015 1015 data_ele_old= data_ele_old[vec]
1016 1016 data_weather= data_weather[vec[0]]
1017 1017 pos_ini = numpy.argmin(data_ele)
1018 1018 if pos_ini>0:
1019 1019 len_vec= len(data_ele)
1020 1020 vec3 = numpy.linspace(pos_ini,len_vec-1,len_vec-pos_ini).astype(int)
1021 1021 #print(vec3)
1022 1022 data_ele= data_ele[vec3]
1023 1023 data_ele_new = data_ele
1024 1024 data_ele_old= data_ele_old[vec3]
1025 1025 data_weather= data_weather[vec3]
1026 1026
1027 1027 new_i_ele = int(data_ele_new[0])
1028 1028 new_f_ele = int(data_ele_new[-1])
1029 1029 n1= new_i_ele- ang_min
1030 1030 n2= ang_max - new_f_ele-1
1031 1031 if n1>0:
1032 1032 ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1)
1033 1033 ele1_nan= numpy.ones(n1)*numpy.nan
1034 1034 data_ele = numpy.hstack((ele1,data_ele_new))
1035 1035 data_ele_old = numpy.hstack((ele1_nan,data_ele_new))
1036 1036 if n2>0:
1037 1037 ele2= numpy.linspace(new_f_ele+1,ang_max,n2)
1038 1038 ele2_nan= numpy.ones(n2)*numpy.nan
1039 1039 data_ele = numpy.hstack((data_ele,ele2))
1040 1040 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1041 1041
1042 1042 self.data_ele_tmp[val_ch] = data_ele_old
1043 1043 self.res_ele = data_ele
1044 1044 self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1045 1045 data_ele = self.res_ele
1046 1046 data_weather = self.res_weather[val_ch]
1047 1047 #print("self.data_ele_tmp",self.data_ele_tmp)
1048 1048 return data_weather,data_ele
1049 1049
1050 1050
1051 1051 def plot(self):
1052 1052 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
1053 1053 data = self.data[-1]
1054 1054 r = self.data.yrange
1055 1055 delta_height = r[1]-r[0]
1056 1056 r_mask = numpy.where(r>=0)[0]
1057 1057 ##print("delta_height",delta_height)
1058 1058 #print("r_mask",r_mask,len(r_mask))
1059 1059 r = numpy.arange(len(r_mask))*delta_height
1060 1060 self.y = 2*r
1061 1061 res = 1
1062 1062 ###print("data['weather'].shape[0]",data['weather'].shape[0])
1063 1063 ang_max = self.ang_max
1064 1064 ang_min = self.ang_min
1065 1065 var_ang =ang_max - ang_min
1066 1066 step = (int(var_ang)/(res*data['weather'].shape[0]))
1067 1067 ###print("step",step)
1068 1068 #--------------------------------------------------------
1069 1069 ##print('weather',data['weather'].shape)
1070 1070 ##print('ele',data['ele'].shape)
1071 1071
1072 1072 ###self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min)
1073 1073 ###self.res_azi = numpy.mean(data['azi'])
1074 1074 ###print("self.res_ele",self.res_ele)
1075 1075 plt.clf()
1076 1076 subplots = [121, 122]
1077 1077 cg={'angular_spacing': 20.}
1078 1078 if self.ini==0:
1079 1079 self.data_ele_tmp = numpy.ones([self.nplots,int(var_ang)])*numpy.nan
1080 1080 self.res_weather= numpy.ones([self.nplots,int(var_ang),len(r_mask)])*numpy.nan
1081 1081 print("SHAPE",self.data_ele_tmp.shape)
1082 1082
1083 1083 for i,ax in enumerate(self.axes):
1084 1084 self.res_weather[i], self.res_ele = self.const_ploteo(val_ch=i, data_weather=data['weather'][i][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min)
1085 1085 self.res_azi = numpy.mean(data['azi'])
1086 1086 if i==0:
1087 1087 print("*****************************************************************************to plot**************************",self.res_weather[i].shape)
1088 1088 self.zmin = self.zmin if self.zmin else 20
1089 1089 self.zmax = self.zmax if self.zmax else 80
1090 1090 if ax.firsttime:
1091 1091 #plt.clf()
1092 1092 cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj=cg,vmin=self.zmin, vmax=self.zmax)
1093 1093 #fig=self.figures[0]
1094 1094 else:
1095 1095 #plt.clf()
1096 1096 if i==0:
1097 1097 print(self.res_weather[i])
1098 1098 print(self.res_ele)
1099 1099 cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj=cg,vmin=self.zmin, vmax=self.zmax)
1100 1100 caax = cgax.parasites[0]
1101 1101 paax = cgax.parasites[1]
1102 1102 cbar = plt.gcf().colorbar(pm, pad=0.075)
1103 1103 caax.set_xlabel('x_range [km]')
1104 1104 caax.set_ylabel('y_range [km]')
1105 1105 plt.text(1.0, 1.05, 'Elevacion '+str(thisDatetime)+" Step "+str(self.ini)+ " Azi: "+str(round(self.res_azi,2)), transform=caax.transAxes, va='bottom',ha='right')
1106 1106 print("***************************self.ini****************************",self.ini)
1107 1107 self.ini= self.ini+1
1108 1108
1109 1109 class Weather_vRF_Plot(Plot):
1110 1110 CODE = 'PPI'
1111 1111 plot_name = 'PPI'
1112 1112 #plot_type = 'ppistyle'
1113 1113 buffering = False
1114 1114
1115 1115 def setup(self):
1116 1116
1117 1117 self.ncols = 1
1118 1118 self.nrows = 1
1119 1119 self.width =8
1120 1120 self.height =8
1121 1121 self.nplots= 1
1122 1122 self.ylabel= 'Range [Km]'
1123 1123 self.xlabel= 'Range [Km]'
1124 1124 self.titles= ['PPI']
1125 1125 self.polar = True
1126 1126 if self.channels is not None:
1127 1127 self.nplots = len(self.channels)
1128 1128 self.nrows = len(self.channels)
1129 1129 else:
1130 1130 self.nplots = self.data.shape(self.CODE)[0]
1131 1131 self.nrows = self.nplots
1132 1132 self.channels = list(range(self.nplots))
1133 1133
1134 1134 if self.CODE == 'POWER':
1135 1135 self.cb_label = r'Power (dB)'
1136 1136 elif self.CODE == 'DOPPLER':
1137 1137 self.cb_label = r'Velocity (m/s)'
1138 1138 self.colorbar=True
1139 1139 self.width = 9
1140 1140 self.height =8
1141 1141 self.ini =0
1142 1142 self.len_azi =0
1143 1143 self.buffer_ini = None
1144 1144 self.buffer_ele = None
1145 1145 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.15, 'right': 0.9, 'bottom': 0.08})
1146 1146 self.flag =0
1147 1147 self.indicador= 0
1148 1148 self.last_data_ele = None
1149 1149 self.val_mean = None
1150 1150
1151 1151 def update(self, dataOut):
1152 1152
1153 1153 data = {}
1154 1154 meta = {}
1155 1155 if hasattr(dataOut, 'dataPP_POWER'):
1156 1156 factor = 1
1157 1157 if hasattr(dataOut, 'nFFTPoints'):
1158 1158 factor = dataOut.normFactor
1159 1159
1160 1160 if 'pow' in self.attr_data[0].lower():
1161 1161 data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))
1162 1162 else:
1163 1163 data['data'] = getattr(dataOut, self.attr_data[0])/(factor)
1164 1164
1165 1165 data['azi'] = dataOut.data_azi
1166 1166 data['ele'] = dataOut.data_ele
1167 1167
1168 1168 return data, meta
1169 1169
1170 1170 def plot(self):
1171 1171 data = self.data[-1]
1172 1172 r = self.data.yrange
1173 1173 delta_height = r[1]-r[0]
1174 1174 r_mask = numpy.where(r>=0)[0]
1175 1175 self.r_mask = r_mask
1176 1176 r = numpy.arange(len(r_mask))*delta_height
1177 1177 self.y = 2*r
1178 1178
1179 1179 try:
1180 1180 z = data['data'][self.channels[0]][:,r_mask]
1181 1181
1182 1182 except:
1183 1183 z = data['data'][0][:,r_mask]
1184 1184
1185 1185 self.titles = []
1186 1186
1187 1187 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
1188 1188 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
1189 1189 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1190 1190 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1191 1191 self.ang_min = self.ang_min if self.ang_min else 0
1192 1192 self.ang_max = self.ang_max if self.ang_max else 360
1193 1193
1194 1194 r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )
1195 1195
1196 1196 for i,ax in enumerate(self.axes):
1197 1197
1198 1198 if ax.firsttime:
1199 1199 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1200 1200 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1201 1201 ax.set_theta_direction(-1)
1202 1202
1203 1203 else:
1204 1204 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1205 1205 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1206 1206 ax.set_theta_direction(-1)
1207 1207
1208 1208 ax.grid(True)
1209 1209
1210 1210 if len(self.channels) !=1:
1211 1211 self.titles = ['PPI {} at EL: {} Channel {}'.format(self.self.labels[x], str(round(numpy.mean(data['ele']),1)), x) for x in range(self.nrows)]
1212 1212 else:
1213 1213 self.titles = ['PPI {} at EL: {} Channel {}'.format(self.labels[0], str(round(numpy.mean(data['ele']),1)), self.channels[0])]
1214 1214
1215 1215 class WeatherRHI_vRF2_Plot(Plot):
1216 1216 CODE = 'weather'
1217 1217 plot_name = 'weather'
1218 1218 plot_type = 'rhistyle'
1219 1219 buffering = False
1220 1220 data_ele_tmp = None
1221 1221
1222 1222 def setup(self):
1223 1223 print("********************")
1224 1224 print("********************")
1225 1225 print("********************")
1226 1226 print("SETUP WEATHER PLOT")
1227 1227 self.ncols = 1
1228 1228 self.nrows = 1
1229 1229 self.nplots= 1
1230 1230 self.ylabel= 'Range [Km]'
1231 1231 self.titles= ['Weather']
1232 1232 if self.channels is not None:
1233 1233 self.nplots = len(self.channels)
1234 1234 self.nrows = len(self.channels)
1235 1235 else:
1236 1236 self.nplots = self.data.shape(self.CODE)[0]
1237 1237 self.nrows = self.nplots
1238 1238 self.channels = list(range(self.nplots))
1239 1239 print("channels",self.channels)
1240 1240 print("que saldra", self.data.shape(self.CODE)[0])
1241 1241 self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)]
1242 1242 print("self.titles",self.titles)
1243 1243 self.colorbar=False
1244 1244 self.width =8
1245 1245 self.height =8
1246 1246 self.ini =0
1247 1247 self.len_azi =0
1248 1248 self.buffer_ini = None
1249 1249 self.buffer_ele = None
1250 1250 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
1251 1251 self.flag =0
1252 1252 self.indicador= 0
1253 1253 self.last_data_ele = None
1254 1254 self.val_mean = None
1255 1255
1256 1256 def update(self, dataOut):
1257 1257
1258 1258 data = {}
1259 1259 meta = {}
1260 1260 if hasattr(dataOut, 'dataPP_POWER'):
1261 1261 factor = 1
1262 1262 if hasattr(dataOut, 'nFFTPoints'):
1263 1263 factor = dataOut.normFactor
1264 1264 print("dataOut",dataOut.data_360.shape)
1265 1265 #
1266 1266 data['weather'] = 10*numpy.log10(dataOut.data_360/(factor))
1267 1267 #
1268 1268 #data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor))
1269 1269 data['azi'] = dataOut.data_azi
1270 1270 data['ele'] = dataOut.data_ele
1271 1271 data['case_flag'] = dataOut.case_flag
1272 1272 #print("UPDATE")
1273 1273 #print("data[weather]",data['weather'].shape)
1274 1274 #print("data[azi]",data['azi'])
1275 1275 return data, meta
1276 1276
1277 1277 def get2List(self,angulos):
1278 1278 list1=[]
1279 1279 list2=[]
1280 1280 for i in reversed(range(len(angulos))):
1281 1281 if not i==0:#el caso de i=0 evalula el primero de la lista con el ultimo y no es relevante
1282 1282 diff_ = angulos[i]-angulos[i-1]
1283 1283 if abs(diff_) >1.5:
1284 1284 list1.append(i-1)
1285 1285 list2.append(diff_)
1286 1286 return list(reversed(list1)),list(reversed(list2))
1287 1287
1288 1288 def fixData90(self,list_,ang_):
1289 1289 if list_[0]==-1:
1290 1290 vec = numpy.where(ang_<ang_[0])
1291 1291 ang_[vec] = ang_[vec]+90
1292 1292 return ang_
1293 1293 return ang_
1294 1294
1295 1295 def fixData90HL(self,angulos):
1296 1296 vec = numpy.where(angulos>=90)
1297 1297 angulos[vec]=angulos[vec]-90
1298 1298 return angulos
1299 1299
1300 1300
1301 1301 def search_pos(self,pos,list_):
1302 1302 for i in range(len(list_)):
1303 1303 if pos == list_[i]:
1304 1304 return True,i
1305 1305 i=None
1306 1306 return False,i
1307 1307
1308 1308 def fixDataComp(self,ang_,list1_,list2_,tipo_case):
1309 1309 size = len(ang_)
1310 1310 size2 = 0
1311 1311 for i in range(len(list2_)):
1312 1312 size2=size2+round(abs(list2_[i]))-1
1313 1313 new_size= size+size2
1314 1314 ang_new = numpy.zeros(new_size)
1315 1315 ang_new2 = numpy.zeros(new_size)
1316 1316
1317 1317 tmp = 0
1318 1318 c = 0
1319 1319 for i in range(len(ang_)):
1320 1320 ang_new[tmp +c] = ang_[i]
1321 1321 ang_new2[tmp+c] = ang_[i]
1322 1322 condition , value = self.search_pos(i,list1_)
1323 1323 if condition:
1324 1324 pos = tmp + c + 1
1325 1325 for k in range(round(abs(list2_[value]))-1):
1326 1326 if tipo_case==0 or tipo_case==3:#subida
1327 1327 ang_new[pos+k] = ang_new[pos+k-1]+1
1328 1328 ang_new2[pos+k] = numpy.nan
1329 1329 elif tipo_case==1 or tipo_case==2:#bajada
1330 1330 ang_new[pos+k] = ang_new[pos+k-1]-1
1331 1331 ang_new2[pos+k] = numpy.nan
1332 1332
1333 1333 tmp = pos +k
1334 1334 c = 0
1335 1335 c=c+1
1336 1336 return ang_new,ang_new2
1337 1337
1338 1338 def globalCheckPED(self,angulos,tipo_case):
1339 1339 l1,l2 = self.get2List(angulos)
1340 1340 ##print("l1",l1)
1341 1341 ##print("l2",l2)
1342 1342 if len(l1)>0:
1343 1343 #angulos2 = self.fixData90(list_=l1,ang_=angulos)
1344 1344 #l1,l2 = self.get2List(angulos2)
1345 1345 ang1_,ang2_ = self.fixDataComp(ang_=angulos,list1_=l1,list2_=l2,tipo_case=tipo_case)
1346 1346 #ang1_ = self.fixData90HL(ang1_)
1347 1347 #ang2_ = self.fixData90HL(ang2_)
1348 1348 else:
1349 1349 ang1_= angulos
1350 1350 ang2_= angulos
1351 1351 return ang1_,ang2_
1352 1352
1353 1353
1354 1354 def replaceNAN(self,data_weather,data_ele,val):
1355 1355 data= data_ele
1356 1356 data_T= data_weather
1357 1357 if data.shape[0]> data_T.shape[0]:
1358 1358 data_N = numpy.ones( [data.shape[0],data_T.shape[1]])
1359 1359 c = 0
1360 1360 for i in range(len(data)):
1361 1361 if numpy.isnan(data[i]):
1362 1362 data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
1363 1363 else:
1364 1364 data_N[i,:]=data_T[c,:]
1365 1365 c=c+1
1366 1366 return data_N
1367 1367 else:
1368 1368 for i in range(len(data)):
1369 1369 if numpy.isnan(data[i]):
1370 1370 data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan
1371 1371 return data_T
1372 1372
1373 1373 def check_case(self,data_ele,ang_max,ang_min):
1374 1374 start = data_ele[0]
1375 1375 end = data_ele[-1]
1376 1376 number = (end-start)
1377 1377 len_ang=len(data_ele)
1378 1378 print("start",start)
1379 1379 print("end",end)
1380 1380 print("number",number)
1381 1381
1382 1382 print("len_ang",len_ang)
1383 1383
1384 1384 #exit(1)
1385 1385
1386 1386 if start<end and (round(abs(number)+1)>=len_ang or (numpy.argmin(data_ele)==0)):#caso subida
1387 1387 return 0
1388 1388 #elif start>end and (round(abs(number)+1)>=len_ang or(numpy.argmax(data_ele)==0)):#caso bajada
1389 1389 # return 1
1390 1390 elif round(abs(number)+1)>=len_ang and (start>end or(numpy.argmax(data_ele)==0)):#caso bajada
1391 1391 return 1
1392 1392 elif round(abs(number)+1)<len_ang and data_ele[-2]>data_ele[-1]:# caso BAJADA CAMBIO ANG MAX
1393 1393 return 2
1394 1394 elif round(abs(number)+1)<len_ang and data_ele[-2]<data_ele[-1] :# caso SUBIDA CAMBIO ANG MIN
1395 1395 return 3
1396 1396
1397 1397
1398 1398 def const_ploteo(self,val_ch,data_weather,data_ele,step,res,ang_max,ang_min,case_flag):
1399 1399 ang_max= ang_max
1400 1400 ang_min= ang_min
1401 1401 data_weather=data_weather
1402 1402 val_ch=val_ch
1403 1403 ##print("*********************DATA WEATHER**************************************")
1404 1404 ##print(data_weather)
1405 1405 if self.ini==0:
1406 1406 '''
1407 1407 print("**********************************************")
1408 1408 print("**********************************************")
1409 1409 print("***************ini**************")
1410 1410 print("**********************************************")
1411 1411 print("**********************************************")
1412 1412 '''
1413 1413 #print("data_ele",data_ele)
1414 1414 #----------------------------------------------------------
1415 1415 tipo_case = case_flag[-1]
1416 1416 #tipo_case = self.check_case(data_ele,ang_max,ang_min)
1417 1417 print("check_case",tipo_case)
1418 1418 #exit(1)
1419 1419 #--------------------- new -------------------------
1420 1420 data_ele_new ,data_ele_old= self.globalCheckPED(data_ele,tipo_case)
1421 1421
1422 1422 #-------------------------CAMBIOS RHI---------------------------------
1423 1423 start= ang_min
1424 1424 end = ang_max
1425 1425 n= (ang_max-ang_min)/res
1426 1426 #------ new
1427 1427 self.start_data_ele = data_ele_new[0]
1428 1428 self.end_data_ele = data_ele_new[-1]
1429 1429 if tipo_case==0 or tipo_case==3: # SUBIDA
1430 1430 n1= round(self.start_data_ele)- start
1431 1431 n2= end - round(self.end_data_ele)
1432 1432 print(self.start_data_ele)
1433 1433 print(self.end_data_ele)
1434 1434 if n1>0:
1435 1435 ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1)
1436 1436 ele1_nan= numpy.ones(n1)*numpy.nan
1437 1437 data_ele = numpy.hstack((ele1,data_ele_new))
1438 1438 print("ele1_nan",ele1_nan.shape)
1439 1439 print("data_ele_old",data_ele_old.shape)
1440 1440 data_ele_old = numpy.hstack((ele1_nan,data_ele_old))
1441 1441 if n2>0:
1442 1442 ele2= numpy.linspace(self.end_data_ele+1,end,n2)
1443 1443 ele2_nan= numpy.ones(n2)*numpy.nan
1444 1444 data_ele = numpy.hstack((data_ele,ele2))
1445 1445 print("ele2_nan",ele2_nan.shape)
1446 1446 print("data_ele_old",data_ele_old.shape)
1447 1447 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1448 1448
1449 1449 if tipo_case==1 or tipo_case==2: # BAJADA
1450 1450 data_ele_new = data_ele_new[::-1] # reversa
1451 1451 data_ele_old = data_ele_old[::-1]# reversa
1452 1452 data_weather = data_weather[::-1,:]# reversa
1453 1453 vec= numpy.where(data_ele_new<ang_max)
1454 1454 data_ele_new = data_ele_new[vec]
1455 1455 data_ele_old = data_ele_old[vec]
1456 1456 data_weather = data_weather[vec[0]]
1457 1457 vec2= numpy.where(0<data_ele_new)
1458 1458 data_ele_new = data_ele_new[vec2]
1459 1459 data_ele_old = data_ele_old[vec2]
1460 1460 data_weather = data_weather[vec2[0]]
1461 1461 self.start_data_ele = data_ele_new[0]
1462 1462 self.end_data_ele = data_ele_new[-1]
1463 1463
1464 1464 n1= round(self.start_data_ele)- start
1465 1465 n2= end - round(self.end_data_ele)-1
1466 1466 print(self.start_data_ele)
1467 1467 print(self.end_data_ele)
1468 1468 if n1>0:
1469 1469 ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1)
1470 1470 ele1_nan= numpy.ones(n1)*numpy.nan
1471 1471 data_ele = numpy.hstack((ele1,data_ele_new))
1472 1472 data_ele_old = numpy.hstack((ele1_nan,data_ele_old))
1473 1473 if n2>0:
1474 1474 ele2= numpy.linspace(self.end_data_ele+1,end,n2)
1475 1475 ele2_nan= numpy.ones(n2)*numpy.nan
1476 1476 data_ele = numpy.hstack((data_ele,ele2))
1477 1477 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1478 1478 # RADAR
1479 1479 # NOTA data_ele y data_weather es la variable que retorna
1480 1480 val_mean = numpy.mean(data_weather[:,-1])
1481 1481 self.val_mean = val_mean
1482 1482 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1483 1483 print("eleold",data_ele_old)
1484 1484 print(self.data_ele_tmp[val_ch])
1485 1485 print(data_ele_old.shape[0])
1486 1486 print(self.data_ele_tmp[val_ch].shape[0])
1487 1487 if (data_ele_old.shape[0]==91 or self.data_ele_tmp[val_ch].shape[0]==91):
1488 1488 import sys
1489 1489 print("EXIT",self.ini)
1490 1490
1491 1491 sys.exit(1)
1492 1492 self.data_ele_tmp[val_ch]= data_ele_old
1493 1493 else:
1494 1494 #print("**********************************************")
1495 1495 #print("****************VARIABLE**********************")
1496 1496 #-------------------------CAMBIOS RHI---------------------------------
1497 1497 #---------------------------------------------------------------------
1498 1498 ##print("INPUT data_ele",data_ele)
1499 1499 flag=0
1500 1500 start_ele = self.res_ele[0]
1501 1501 #tipo_case = self.check_case(data_ele,ang_max,ang_min)
1502 1502 tipo_case = case_flag[-1]
1503 1503 #print("TIPO DE DATA",tipo_case)
1504 1504 #-----------new------------
1505 1505 data_ele ,data_ele_old = self.globalCheckPED(data_ele,tipo_case)
1506 1506 data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1507 1507
1508 1508 #-------------------------------NEW RHI ITERATIVO-------------------------
1509 1509
1510 1510 if tipo_case==0 : # SUBIDA
1511 1511 vec = numpy.where(data_ele<ang_max)
1512 1512 data_ele = data_ele[vec]
1513 1513 data_ele_old = data_ele_old[vec]
1514 1514 data_weather = data_weather[vec[0]]
1515 1515
1516 1516 vec2 = numpy.where(0<data_ele)
1517 1517 data_ele= data_ele[vec2]
1518 1518 data_ele_old= data_ele_old[vec2]
1519 1519 ##print(data_ele_new)
1520 1520 data_weather= data_weather[vec2[0]]
1521 1521
1522 1522 new_i_ele = int(round(data_ele[0]))
1523 1523 new_f_ele = int(round(data_ele[-1]))
1524 1524 #print(new_i_ele)
1525 1525 #print(new_f_ele)
1526 1526 #print(data_ele,len(data_ele))
1527 1527 #print(data_ele_old,len(data_ele_old))
1528 1528 if new_i_ele< 2:
1529 1529 self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan
1530 1530 self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean)
1531 1531 self.data_ele_tmp[val_ch][new_i_ele:new_i_ele+len(data_ele)]=data_ele_old
1532 1532 self.res_ele[new_i_ele:new_i_ele+len(data_ele)]= data_ele
1533 1533 self.res_weather[val_ch][new_i_ele:new_i_ele+len(data_ele),:]= data_weather
1534 1534 data_ele = self.res_ele
1535 1535 data_weather = self.res_weather[val_ch]
1536 1536
1537 1537 elif tipo_case==1 : #BAJADA
1538 1538 data_ele = data_ele[::-1] # reversa
1539 1539 data_ele_old = data_ele_old[::-1]# reversa
1540 1540 data_weather = data_weather[::-1,:]# reversa
1541 1541 vec= numpy.where(data_ele<ang_max)
1542 1542 data_ele = data_ele[vec]
1543 1543 data_ele_old = data_ele_old[vec]
1544 1544 data_weather = data_weather[vec[0]]
1545 1545 vec2= numpy.where(0<data_ele)
1546 1546 data_ele = data_ele[vec2]
1547 1547 data_ele_old = data_ele_old[vec2]
1548 1548 data_weather = data_weather[vec2[0]]
1549 1549
1550 1550
1551 1551 new_i_ele = int(round(data_ele[0]))
1552 1552 new_f_ele = int(round(data_ele[-1]))
1553 1553 #print(data_ele)
1554 1554 #print(ang_max)
1555 1555 #print(data_ele_old)
1556 1556 if new_i_ele <= 1:
1557 1557 new_i_ele = 1
1558 1558 if round(data_ele[-1])>=ang_max-1:
1559 1559 self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan
1560 1560 self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean)
1561 1561 self.data_ele_tmp[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1]=data_ele_old
1562 1562 self.res_ele[new_i_ele-1:new_i_ele+len(data_ele)-1]= data_ele
1563 1563 self.res_weather[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1,:]= data_weather
1564 1564 data_ele = self.res_ele
1565 1565 data_weather = self.res_weather[val_ch]
1566 1566
1567 1567 elif tipo_case==2: #bajada
1568 1568 vec = numpy.where(data_ele<ang_max)
1569 1569 data_ele = data_ele[vec]
1570 1570 data_weather= data_weather[vec[0]]
1571 1571
1572 1572 len_vec = len(vec)
1573 1573 data_ele_new = data_ele[::-1] # reversa
1574 1574 data_weather = data_weather[::-1,:]
1575 1575 new_i_ele = int(data_ele_new[0])
1576 1576 new_f_ele = int(data_ele_new[-1])
1577 1577
1578 1578 n1= new_i_ele- ang_min
1579 1579 n2= ang_max - new_f_ele-1
1580 1580 if n1>0:
1581 1581 ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1)
1582 1582 ele1_nan= numpy.ones(n1)*numpy.nan
1583 1583 data_ele = numpy.hstack((ele1,data_ele_new))
1584 1584 data_ele_old = numpy.hstack((ele1_nan,data_ele_new))
1585 1585 if n2>0:
1586 1586 ele2= numpy.linspace(new_f_ele+1,ang_max,n2)
1587 1587 ele2_nan= numpy.ones(n2)*numpy.nan
1588 1588 data_ele = numpy.hstack((data_ele,ele2))
1589 1589 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1590 1590
1591 1591 self.data_ele_tmp[val_ch] = data_ele_old
1592 1592 self.res_ele = data_ele
1593 1593 self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1594 1594 data_ele = self.res_ele
1595 1595 data_weather = self.res_weather[val_ch]
1596 1596
1597 1597 elif tipo_case==3:#subida
1598 1598 vec = numpy.where(0<data_ele)
1599 1599 data_ele= data_ele[vec]
1600 1600 data_ele_new = data_ele
1601 1601 data_ele_old= data_ele_old[vec]
1602 1602 data_weather= data_weather[vec[0]]
1603 1603 pos_ini = numpy.argmin(data_ele)
1604 1604 if pos_ini>0:
1605 1605 len_vec= len(data_ele)
1606 1606 vec3 = numpy.linspace(pos_ini,len_vec-1,len_vec-pos_ini).astype(int)
1607 1607 #print(vec3)
1608 1608 data_ele= data_ele[vec3]
1609 1609 data_ele_new = data_ele
1610 1610 data_ele_old= data_ele_old[vec3]
1611 1611 data_weather= data_weather[vec3]
1612 1612
1613 1613 new_i_ele = int(data_ele_new[0])
1614 1614 new_f_ele = int(data_ele_new[-1])
1615 1615 n1= new_i_ele- ang_min
1616 1616 n2= ang_max - new_f_ele-1
1617 1617 if n1>0:
1618 1618 ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1)
1619 1619 ele1_nan= numpy.ones(n1)*numpy.nan
1620 1620 data_ele = numpy.hstack((ele1,data_ele_new))
1621 1621 data_ele_old = numpy.hstack((ele1_nan,data_ele_new))
1622 1622 if n2>0:
1623 1623 ele2= numpy.linspace(new_f_ele+1,ang_max,n2)
1624 1624 ele2_nan= numpy.ones(n2)*numpy.nan
1625 1625 data_ele = numpy.hstack((data_ele,ele2))
1626 1626 data_ele_old = numpy.hstack((data_ele_old,ele2_nan))
1627 1627
1628 1628 self.data_ele_tmp[val_ch] = data_ele_old
1629 1629 self.res_ele = data_ele
1630 1630 self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean)
1631 1631 data_ele = self.res_ele
1632 1632 data_weather = self.res_weather[val_ch]
1633 1633 #print("self.data_ele_tmp",self.data_ele_tmp)
1634 1634 return data_weather,data_ele
1635 1635
1636 1636
1637 1637 def plot(self):
1638 1638 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
1639 1639 data = self.data[-1]
1640 1640 r = self.data.yrange
1641 1641 delta_height = r[1]-r[0]
1642 1642 r_mask = numpy.where(r>=0)[0]
1643 1643 ##print("delta_height",delta_height)
1644 1644 #print("r_mask",r_mask,len(r_mask))
1645 1645 r = numpy.arange(len(r_mask))*delta_height
1646 1646 self.y = 2*r
1647 1647 res = 1
1648 1648 ###print("data['weather'].shape[0]",data['weather'].shape[0])
1649 1649 ang_max = self.ang_max
1650 1650 ang_min = self.ang_min
1651 1651 var_ang =ang_max - ang_min
1652 1652 step = (int(var_ang)/(res*data['weather'].shape[0]))
1653 1653 ###print("step",step)
1654 1654 #--------------------------------------------------------
1655 1655 ##print('weather',data['weather'].shape)
1656 1656 ##print('ele',data['ele'].shape)
1657 1657
1658 1658 ###self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min)
1659 1659 ###self.res_azi = numpy.mean(data['azi'])
1660 1660 ###print("self.res_ele",self.res_ele)
1661 1661 plt.clf()
1662 1662 subplots = [121, 122]
1663 1663 try:
1664 1664 if self.data[-2]['ele'].max()<data['ele'].max():
1665 1665 self.ini=0
1666 1666 except:
1667 1667 pass
1668 1668 if self.ini==0:
1669 1669 self.data_ele_tmp = numpy.ones([self.nplots,int(var_ang)])*numpy.nan
1670 1670 self.res_weather= numpy.ones([self.nplots,int(var_ang),len(r_mask)])*numpy.nan
1671 1671 print("SHAPE",self.data_ele_tmp.shape)
1672 1672
1673 1673 for i,ax in enumerate(self.axes):
1674 1674 self.res_weather[i], self.res_ele = self.const_ploteo(val_ch=i, data_weather=data['weather'][i][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min,case_flag=self.data['case_flag'])
1675 1675 self.res_azi = numpy.mean(data['azi'])
1676 1676
1677 1677 if ax.firsttime:
1678 1678 #plt.clf()
1679 1679 print("Frist Plot")
1680 1680 cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj='cg',vmin=20, vmax=80)
1681 1681 #fig=self.figures[0]
1682 1682 else:
1683 1683 #plt.clf()
1684 1684 print("ELSE PLOT")
1685 1685 cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj='cg',vmin=20, vmax=80)
1686 1686 caax = cgax.parasites[0]
1687 1687 paax = cgax.parasites[1]
1688 1688 cbar = plt.gcf().colorbar(pm, pad=0.075)
1689 1689 caax.set_xlabel('x_range [km]')
1690 1690 caax.set_ylabel('y_range [km]')
1691 1691 plt.text(1.0, 1.05, 'Elevacion '+str(thisDatetime)+" Step "+str(self.ini)+ " Azi: "+str(round(self.res_azi,2)), transform=caax.transAxes, va='bottom',ha='right')
1692 1692 print("***************************self.ini****************************",self.ini)
1693 1693 self.ini= self.ini+1
1694 1694
1695 1695
1696 1696
1697 1697
1698 1698
1699 1699 class WeatherRHI_vRF4_Plot(Plot):
1700 1700 CODE = 'RHI'
1701 1701 plot_name = 'RHI'
1702 1702 #plot_type = 'rhistyle'
1703 1703 buffering = False
1704 1704
1705 1705 def setup(self):
1706 1706
1707 1707 self.ncols = 1
1708 1708 self.nrows = 1
1709 1709 self.nplots= 1
1710 1710 self.ylabel= 'Range [Km]'
1711 1711 self.xlabel= 'Range [Km]'
1712 1712 self.titles= ['RHI']
1713 1713 self.polar = True
1714 1714 self.grid = True
1715 1715 if self.channels is not None:
1716 1716 self.nplots = len(self.channels)
1717 1717 self.nrows = len(self.channels)
1718 1718 else:
1719 1719 self.nplots = self.data.shape(self.CODE)[0]
1720 1720 self.nrows = self.nplots
1721 1721 self.channels = list(range(self.nplots))
1722 1722
1723 1723 if self.CODE == 'Power':
1724 1724 self.cb_label = r'Power (dB)'
1725 1725 elif self.CODE == 'Doppler':
1726 1726 self.cb_label = r'Velocity (m/s)'
1727 1727 self.colorbar=True
1728 1728 self.width =8
1729 1729 self.height =8
1730 1730 self.ini =0
1731 1731 self.len_azi =0
1732 1732 self.buffer_ini = None
1733 1733 self.buffer_ele = None
1734 1734 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
1735 1735 self.flag =0
1736 1736 self.indicador= 0
1737 1737 self.last_data_ele = None
1738 1738 self.val_mean = None
1739 1739
1740 1740 def update(self, dataOut):
1741 1741
1742 1742 data = {}
1743 1743 meta = {}
1744 1744 if hasattr(dataOut, 'dataPP_POWER'):
1745 1745 factor = 1
1746 1746 if hasattr(dataOut, 'nFFTPoints'):
1747 1747 factor = dataOut.normFactor
1748 1748
1749 1749 if 'pow' in self.attr_data[0].lower():
1750 1750 data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))
1751 1751 else:
1752 1752 data['data'] = getattr(dataOut, self.attr_data[0])/(factor)
1753 1753
1754 1754 data['azi'] = dataOut.data_azi
1755 1755 data['ele'] = dataOut.data_ele
1756 1756
1757 1757 return data, meta
1758 1758
1759 1759 def plot(self):
1760 1760 data = self.data[-1]
1761 1761 r = self.data.yrange
1762 1762 delta_height = r[1]-r[0]
1763 1763 r_mask = numpy.where(r>=0)[0]
1764 1764 self.r_mask =r_mask
1765 1765 r = numpy.arange(len(r_mask))*delta_height
1766 1766 self.y = 2*r
1767 1767
1768 1768 try:
1769 1769 z = data['data'][self.channels[0]][:,r_mask]
1770 1770 except:
1771 1771 z = data['data'][0][:,r_mask]
1772 1772
1773 1773 self.titles = []
1774 1774
1775 1775 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
1776 1776 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
1777 1777 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1778 1778 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1779 1779 self.ang_min = self.ang_min if self.ang_min else 0
1780 1780 self.ang_max = self.ang_max if self.ang_max else 90
1781 1781
1782 1782 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )
1783 1783
1784 1784 for i,ax in enumerate(self.axes):
1785 1785
1786 1786 if ax.firsttime:
1787 1787 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1788 1788 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1789 1789
1790 1790 else:
1791 1791 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1792 1792 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1793 1793 ax.grid(True)
1794 1794 if len(self.channels) !=1:
1795 1795 self.titles = ['RHI {} at AZ: {} Channel {}'.format(self.labels[x], str(round(numpy.mean(data['azi']),1)), x) for x in range(self.nrows)]
1796 1796 else:
1797 1797 self.titles = ['RHI {} at AZ: {} Channel {}'.format(self.labels[0], str(round(numpy.mean(data['azi']),1)), self.channels[0])]
1798 1798
1799 1799 class WeatherParamsPlot(Plot):
1800 1800 #CODE = 'RHI'
1801 1801 #plot_name = 'RHI'
1802 1802 #plot_type = 'rhistyle'
1803 1803 buffering = False
1804 1804
1805 1805 def setup(self):
1806 1806
1807 1807 self.ncols = 1
1808 1808 self.nrows = 1
1809 1809 self.nplots= 1
1810 1810 self.ylabel= 'Range [km]'
1811 1811 self.xlabel= 'Range [km]'
1812 1812 self.polar = True
1813 1813 self.grid = True
1814 1814 if self.channels is not None:
1815 1815 self.nplots = len(self.channels)
1816 1816 self.nrows = len(self.channels)
1817 1817 else:
1818 1818 self.nplots = self.data.shape(self.CODE)[0]
1819 1819 self.nrows = self.nplots
1820 1820 self.channels = list(range(self.nplots))
1821 1821
1822 1822 self.colorbar=True
1823 1823 self.width =8
1824 1824 self.height =8
1825 1825 self.ini =0
1826 1826 self.len_azi =0
1827 1827 self.buffer_ini = None
1828 1828 self.buffer_ele = None
1829 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.08, 'right': 0.92, 'bottom': 0.01,'top':0.99})
1829 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
1830 1830 self.flag =0
1831 1831 self.indicador= 0
1832 1832 self.last_data_ele = None
1833 1833 self.val_mean = None
1834 1834
1835 1835 def update(self, dataOut):
1836 1836
1837 1837 data = {}
1838 1838 meta = {}
1839 1839 if hasattr(dataOut, 'dataPP_POWER'):
1840 1840 factor = 1
1841 1841 if hasattr(dataOut, 'nFFTPoints'):
1842 1842 factor = dataOut.normFactor
1843 1843
1844 1844 if 'pow' in self.attr_data[0].lower():
1845 1845 data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))
1846 1846 else:
1847 1847 data['data'] = getattr(dataOut, self.attr_data[0])/(factor)
1848 1848
1849 1849 if dataOut.mode_op == 'PPI':
1850 1850 self.CODE = 'PPI'
1851 1851 self.title = self.CODE
1852 1852 elif dataOut.mode_op == 'RHI':
1853 1853 self.CODE = 'RHI'
1854 1854 self.title = self.CODE
1855 1855
1856 1856 data['azi'] = dataOut.data_azi
1857 1857 data['ele'] = dataOut.data_ele
1858 1858 data['mode_op'] = dataOut.mode_op
1859 1859
1860 1860 return data, meta
1861 1861
1862 1862 def plot(self):
1863 1863 data = self.data[-1]
1864 1864 r = self.data.yrange
1865 1865 delta_height = r[1]-r[0]
1866 1866 r_mask = numpy.where(r>=0)[0]
1867 1867 self.r_mask =r_mask
1868 1868 r = numpy.arange(len(r_mask))*delta_height
1869 1869 self.y = 2*r
1870 1870
1871 1871 try:
1872 1872 z = data['data'][self.channels[0]][:,r_mask]
1873 1873 except:
1874 1874 z = data['data'][0][:,r_mask]
1875 1875
1876 1876 self.titles = []
1877 1877
1878 1878 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
1879 1879 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
1880 1880 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1881 1881 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1882 1882 print("mode inside plot",self.data['mode_op'],data['mode_op'])
1883 1883 if data['mode_op'] == 'RHI':
1884 1884 try:
1885 1885 if self.data['mode_op'][-2] == 'PPI':
1886 1886 self.ang_min = None
1887 1887 self.ang_max = None
1888 1888 except:
1889 1889 pass
1890 1890 self.ang_min = self.ang_min if self.ang_min else 0
1891 1891 self.ang_max = self.ang_max if self.ang_max else 90
1892 1892 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )
1893 1893 elif data['mode_op'] == 'PPI':
1894 1894 try:
1895 1895 if self.data['mode_op'][-2] == 'RHI':
1896 1896 self.ang_min = None
1897 1897 self.ang_max = None
1898 1898 except:
1899 1899 pass
1900 1900 self.ang_min = self.ang_min if self.ang_min else 0
1901 1901 self.ang_max = self.ang_max if self.ang_max else 360
1902 1902 r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )
1903 1903
1904 1904 self.clear_figures()
1905 1905
1906 1906 for i,ax in enumerate(self.axes):
1907 1907
1908 1908 if ax.firsttime:
1909 1909 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1910 1910 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1911 1911 if data['mode_op'] == 'PPI':
1912 1912 ax.set_theta_direction(-1)
1913 1913 ax.set_theta_offset(numpy.pi/2)
1914 1914
1915 1915 else:
1916 1916 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
1917 1917 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
1918 1918 if data['mode_op'] == 'PPI':
1919 1919 ax.set_theta_direction(-1)
1920 1920 ax.set_theta_offset(numpy.pi/2)
1921 1921
1922 1922 ax.grid(True)
1923 1923 if data['mode_op'] == 'RHI':
1924 1924 len_aux = int(data['azi'].shape[0]/4)
1925 1925 mean = numpy.mean(data['azi'][len_aux:-len_aux])
1926 1926 if len(self.channels) !=1:
1927 1927 self.titles = ['RHI {} at AZ: {} Channel {}'.format(self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
1928 1928 else:
1929 1929 self.titles = ['RHI {} at AZ: {} Channel {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
1930 1930 elif data['mode_op'] == 'PPI':
1931 1931 len_aux = int(data['ele'].shape[0]/4)
1932 1932 mean = numpy.mean(data['ele'][len_aux:-len_aux])
1933 1933 if len(self.channels) !=1:
1934 1934 self.titles = ['PPI {} at EL: {} Channel {}'.format(self.self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
1935 1935 else:
1936 1936 self.titles = ['PPI {} at EL: {} Channel {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
@@ -1,377 +1,378
1
1 2 # SOPHY PROC script
2 3 import os, sys, json, argparse
3 4 import datetime
4 5 import time
5 6
6 7 PATH = '/DATA_RM/DATA'
7 8 # PATH = '/Users/jespinoza/workspace/data/'
8 PATH = '/home/soporte/Documents/HUANCAYO'
9 #PATH = '/home/soporte/Documents/HUANCAYO'
9 10 PARAM = {
10 'P': {'name': 'dataPP_POWER', 'zmin': -40, 'zmax': -5, 'colormap': 'jet', 'label': 'Power', 'wrname': 'Pow','cb_label': 'dB', 'ch':0},
11 'P': {'name': 'dataPP_POWER', 'zmin': -45, 'zmax': -25, 'colormap': 'jet', 'label': 'Power', 'wrname': 'Pow','cb_label': 'dB', 'ch':0},
11 12 'V': {'name': 'dataPP_DOP', 'zmin': -20, 'zmax': 20, 'colormap': 'seismic', 'label': 'Velocity', 'wrname': 'Vel', 'cb_label': 'm/s', 'ch':0},
12 13 'RH': {'name': 'RhoHV_R', 'zmin': 0, 'zmax': 1, 'colormap': 'jet', 'label': 'Coef.Correlacion', 'wrname':'R', 'cb_label': '*', 'ch':0},
13 14 'FD': {'name': 'PhiD_P', 'zmin': -180,'zmax': 180,'colormap': 'RdBu_r', 'label': 'Fase Diferencial', 'wrname':'P' , 'cb_label': 'ΒΊ', 'ch':0},
14 15 'ZD': {'name': 'Zdb_D', 'zmin': -20, 'zmax': 60, 'colormap': 'viridis','label': 'Reflect.Diferencial','wrname':'D' , 'cb_label': 'dBz','ch':0},
15 16 'Z': {'name': 'Zdb', 'zmin': -20, 'zmax': 70, 'colormap': 'gist_ncar','label': 'Reflectividad', 'wrname':'Z', 'cb_label': 'dBz','ch':1},
16 17 'W': {'name': 'Sigmav_W', 'zmin': 0, 'zmax':5, 'colormap': 'viridis','label': 'AnchoEspectral', 'wrname':'S', 'cb_label': 'hz', 'ch':1}
17 18 }
18 19
19 20 #
20 21 def max_index(r, sample_rate, ipp):
21 22
22 23 return int(sample_rate*ipp*1e6 * r / 60) + int(sample_rate*ipp*1e6 * 1.2 / 60)
23 24
24 25 def main(args):
25 26
26 27 experiment = args.experiment
27 28 fp = open(os.path.join(PATH, experiment, 'experiment.conf'))
28 29 conf = json.loads(fp.read())
29 30
30 31 ipp_km = conf['usrp_tx']['ipp']
31 32 ipp = ipp_km * 2 /300000
32 33 sample_rate = conf['usrp_rx']['sample_rate']
33 34 axis = ['0' if x=='elevation' else '1' for x in conf['pedestal']['axis']] # AZIMUTH 1 ELEVACION 0
34 35 speed_axis = conf['pedestal']['speed']
35 36 steeps = conf['pedestal']['table']
36 37 time_offset = args.time_offset
37 38 parameters = args.parameters
38 39 start_date = experiment.split('@')[1].split('T')[0].replace('-', '/')
39 40 end_date = start_date
40 41 start_time = experiment.split('@')[1].split('T')[1].replace('-', ':')
41 42 end_time = '23:59:59'
42 43 N = int(1/(speed_axis[0]*ipp)) # 1 GRADO DE RESOLUCION
43 44 path = os.path.join(PATH, experiment, 'rawdata')
44 45 path_ped = os.path.join(PATH, experiment, 'position')
45 46 path_plots = os.path.join(PATH, experiment, 'plots_ch0')
46 47 path_save = os.path.join(PATH, experiment, 'param')
47 RMIX = 10
48 RMIX = 20
48 49
49 50 from schainpy.controller import Project
50 51
51 52 project = Project()
52 53 project.setup(id='1', name='Sophy', description='sophy proc')
53 54
54 55 reader = project.addReadUnit(datatype='DigitalRFReader',
55 56 path=path,
56 57 startDate=start_date,
57 58 endDate=end_date,
58 59 startTime=start_time,
59 60 endTime=end_time,
60 61 delay=30,
61 62 channelList='0',
62 63 online=args.online,
63 64 walk=1,
64 65 ippKm = ipp_km,
65 66 getByBlock = 1,
66 67 nProfileBlocks = N,
67 68 )
68 69
69 70 if not conf['usrp_tx']['enable_2']: # One Pulse
70 71 voltage = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
71 72
72 73 if conf['usrp_tx']['code_type_1']:
73 74 code = [c.split() for c in conf['usrp']['code_1']]
74 75 op = voltage.addOperation(name='Decoder', optype='other')
75 76 op.addParameter(name='code', value=code)
76 77 op.addParameter(name='nCode', value=len(code), format='int')
77 78 op.addParameter(name='nBaud', value=len(code[0]), format='int')
78 79
79 80 op = voltage.addOperation(name='setH0')
80 81 op.addParameter(name='h0', value='-1.2')
81 82
82 83 if args.range >= 0:
83 84 op = voltage.addOperation(name='selectHeights')
84 85 op.addParameter(name='minIndex', value='0', format='int')
85 86 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
86 87
87 88 code=[[1]]
88 89 opObj11 = voltage.addOperation(name='Decoder', optype='other')
89 90 opObj11.addParameter(name='code', value=code)
90 91 opObj11.addParameter(name='nCode', value='1', format='int')
91 92 opObj11.addParameter(name='nBaud', value='1', format='int')
92 93
93 94 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
94 95 op.addParameter(name='n', value=2*len(code), format='int')
95 96
96 97 #op = voltage.addOperation(name='PulsePair_vRF', optype='other')
97 98 #op.addParameter(name='n', value=int(N), format='int')
98 99
99 100 if args.range >= 0:
100 101 op = voltage.addOperation(name='selectHeights')
101 102 op.addParameter(name='minIndex', value='0', format='int')
102 103 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
103 104
104 105
105 106 op = voltage.addOperation(name='PulsePair_vRF', optype='other')
106 107 op.addParameter(name='n', value=125, format='int')
107 108
108 109
109 110 proc = project.addProcUnit(datatype='ParametersProc', inputId=voltage.getId())
110 111 #procUnitConfObjB.addParameter(name='runNextUnit', value=True)
111 112
112 113 opObj10 = proc.addOperation(name="WeatherRadar")
113 114 opObj10.addParameter(name='variableList',value='Reflectividad,VelocidadRadial,AnchoEspectral')
114 115
115 116 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
116 117
117 118 op = proc.addOperation(name='PedestalInformation')
118 119 op.addParameter(name='path', value=path_ped, format='str')
119 120 op.addParameter(name='interval', value='0.04')
120 121 op.addParameter(name='time_offset', value=time_offset)
121 122 op.addParameter(name='az_offset', value=-26.2)
122 123
123 124 for param in parameters:
124 125 op = proc.addOperation(name='Block360_vRF4')
125 126 #op.addParameter(name='axis', value=','.join(axis))
126 127 op.addParameter(name='attr_data', value=PARAM[param]['name'])
127 128 op.addParameter(name='runNextOp', value=True)
128 129
129 130 op= proc.addOperation(name='WeatherParamsPlot')
130 131 if args.save: op.addParameter(name='save', value=path_plots, format='str')
131 132 op.addParameter(name='save_period', value=-1)
132 133 op.addParameter(name='show', value=args.show)
133 134 op.addParameter(name='channels', value='0,')
134 135 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
135 136 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
136 137 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
137 138 op.addParameter(name='labels', value=[PARAM[param]['label']])
138 139 op.addParameter(name='save_code', value=param)
139 140 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
140 141 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
141 142
142 143 desc = {
143 144 'Data': {
144 145 PARAM[param]['name']: PARAM[param]['label'],
145 146 'utctime': 'time'
146 147 },
147 148 'Metadata': {
148 149 'heightList': 'range',
149 150 'data_azi': 'azimuth',
150 151 'data_ele': 'elevation',
151 152 }
152 153 }
153 154
154 155 if args.save:
155 156 opObj10 = proc.addOperation(name='HDFWriter')
156 157 opObj10.addParameter(name='path',value=path_save+'-{}'.format(param), format='str')
157 158 opObj10.addParameter(name='Reset',value=True)
158 159 opObj10.addParameter(name='setType',value='weather')
159 160 opObj10.addParameter(name='description',value='desc')
160 161 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
161 162 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
162 163 opObj10.addParameter(name='dataList',value='{},utctime'.format(PARAM[param]['name']))
163 164
164 165 else: #Two pulses
165 166
166 167 voltage1 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
167 168
168 169 print("repetions",conf['usrp_tx']['repetitions_1'])
169 170
170 171 op = voltage1.addOperation(name='ProfileSelector')
171 172 op.addParameter(name='profileRangeList', value='0,{}'.format(conf['usrp_tx']['repetitions_1']-1))
172 173
173 174
174 175 #op3 = voltage1.addOperation(name='ProfileSelector', optype='other')
175 176 #op3.addParameter(name='profileRangeList', value='1,123')
176 177
177 178 '''
178 179 if conf['usrp_tx']['code_type_1'] != 'None':
179 180 code = [c.split() for c in conf['usrp_tx']['code_1']]
180 181 op = voltage1.addOperation(name='Decoder', optype='other')
181 182 op.addParameter(name='code', value=code)
182 183 op.addParameter(name='nCode', value=len(code), format='int')
183 184 op.addParameter(name='nBaud', value=len(code[0]), format='int')
184 185 '''
185 186
186 187 code=[[1]]
187 188
188 189 opObj11 = voltage1.addOperation(name='Decoder', optype='other')
189 190 opObj11.addParameter(name='code', value=code)
190 191 opObj11.addParameter(name='nCode', value='1', format='int')
191 192 opObj11.addParameter(name='nBaud', value='1', format='int')
192 193
193 194 op = voltage1.addOperation(name='setH0')
194 195 op.addParameter(name='h0', value='-1.2')
195 196
196 197 if args.range >= 0:
197 198 op = voltage1.addOperation(name='selectHeights')
198 199 op.addParameter(name='minIndex', value='0', format='int')
199 200 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
200 201
201 202 op = voltage1.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
202 203 op.addParameter(name='n', value=2, format='int')
203 204
204 205 op = voltage1.addOperation(name='PulsePair_vRF', optype='other')
205 206 #op.addParameter(name='n', value=int(N), format='int')
206 207 op.addParameter(name='n', value=61, format='int')
207 208 #op.addParameter(name='removeDC',value=True)
208 209
209 210 '''
210 211 if args.range >= 0:
211 212 print("corto",max_index(RMIX, sample_rate, ipp))
212 213 op = voltage1.addOperation(name='selectHeights')
213 214 op.addParameter(name='minIndex', value='0', format='int')
214 215 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
215 216 '''
216 217 proc1 = project.addProcUnit(datatype='ParametersProc', inputId=voltage1.getId())
217 218 proc1.addParameter(name='runNextUnit', value=True)
218 219
219 220 opObj10 = proc1.addOperation(name="WeatherRadar")
220 221 opObj10.addParameter(name='variableList',value='Reflectividad,VelocidadRadial,AnchoEspectral')
221 222 opObj10.addParameter(name='tauW',value=0.4*1e-6)
222 223 opObj10.addParameter(name='Pt',value=0.2)
223 224
224 225 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
225 226
226 227 op = proc1.addOperation(name='PedestalInformation')
227 228 op.addParameter(name='path', value=path_ped, format='str')
228 229 op.addParameter(name='interval', value='0.04')
229 230 op.addParameter(name='time_offset', value=time_offset)
230 231 op.addParameter(name='az_offset', value=-26.2)
231 232
232 233 for param in parameters:
233 234 op = proc1.addOperation(name='Block360_vRF4')
234 235 op.addParameter(name='attr_data', value=PARAM[param]['name'])
235 236 op.addParameter(name='runNextOp', value=True)
236 237
237 238 voltage2 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
238 239
239 240 op = voltage2.addOperation(name='ProfileSelector')
240 241 op.addParameter(name='profileRangeList', value='{},{}'.format(conf['usrp_tx']['repetitions_1'], conf['usrp_tx']['repetitions_1']+conf['usrp_tx']['repetitions_2']-1))
241 242
242 243
243 244 if conf['usrp_tx']['code_type_2']:
244 245 print(conf['usrp_tx']['code_2'])
245 246 codes = [ c.strip() for c in conf['usrp_tx']['code_2'].split(',')]
246 247 code = []
247 248 for c in codes:
248 249 code.append([int(x) for x in c])
249 250 print(code)
250 251 print(code[0])
251 252 op = voltage2.addOperation(name='Decoder', optype='other')
252 253 op.addParameter(name='code', value=code)
253 254 op.addParameter(name='nCode', value=len(code), format='int')
254 255 op.addParameter(name='nBaud', value=len(code[0]), format='int')
255 256 import numpy
256 257 pwcode = numpy.sum(numpy.array(code[0])**2)
257 258 print("pwcode",pwcode)
258 259
259 260 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
260 261 op.addParameter(name='n', value=len(code), format='int')
261 262 ncode = len(code)
262 263 else:
263 264 ncode = 1
264 265
265 266 op = voltage2.addOperation(name='setH0')
266 267 op.addParameter(name='h0', value='-1.2')
267 268
268 269 if args.range >= 0:
269 270 if args.range==0:
270 271 args.range= ipp_km
271 272 op = voltage2.addOperation(name='selectHeights')
272 273 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
273 274 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
274 275
275 276 #op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
276 277 #op.addParameter(name='n', value=int(N)/ncode, format='int')
277 278 op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
278 279 op.addParameter(name='n', value=64, format='int')
279 280 #op.addParameter(name='removeDC',value=True)
280 281
281 282 '''
282 283
283 284 if args.range >= 0:
284 285 if args.range==0:
285 286 args.range= ipp_km
286 287 op = voltage2.addOperation(name='selectHeights')
287 288 print("largo",max_index(RMIX, sample_rate, ipp))
288 289 print("largo2",max_index(args.range, sample_rate, ipp))
289 290
290 291 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
291 292 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
292 293 '''
293 294
294 295 proc2 = project.addProcUnit(datatype='ParametersProc', inputId=voltage2.getId())
295 296
296 297 opObj10 = proc2.addOperation(name="WeatherRadar")
297 298 opObj10.addParameter(name='variableList',value='Reflectividad,AnchoEspectral')
298 opObj10.addParameter(name='tauW',value=3.2*1e-6)
299 opObj10.addParameter(name='Pt',value=1.6)
299 opObj10.addParameter(name='tauW',value=6.3*1e-6)
300 opObj10.addParameter(name='Pt',value=3.2)
300 301
301 302
302 303 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
303 304
304 305 op = proc2.addOperation(name='PedestalInformation')
305 306 op.addParameter(name='path', value=path_ped, format='str')
306 307 op.addParameter(name='interval', value='0.04')
307 308 op.addParameter(name='time_offset', value=time_offset)
308 309 op.addParameter(name='az_offset', value=-26.2)
309 310
310 311 for param in parameters:
311 312 op = proc2.addOperation(name='Block360_vRF4')
312 313 #op.addParameter(name='axis', value=','.join(axis))
313 314 op.addParameter(name='attr_data', value=PARAM[param]['name'])
314 315 op.addParameter(name='runNextOp', value=True)
315 316
316 317 merge = project.addProcUnit(datatype='MergeProc', inputId=[proc1.getId(), proc2.getId()])
317 318 merge.addParameter(name='attr_data', value=PARAM[param]['name'])
318 319 merge.addParameter(name='mode', value='7') #RM
319 320
320 321 op= merge.addOperation(name='WeatherParamsPlot')
321 322 if args.save: op.addParameter(name='save', value=path_plots, format='str')
322 323 op.addParameter(name='save_period', value=-1)
323 324 op.addParameter(name='show', value=args.show)
324 325 op.addParameter(name='channels', value='0,')
325 326 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
326 327 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
327 328 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
328 329 op.addParameter(name='labels', value=[PARAM[param]['label']])
329 330 op.addParameter(name='save_code', value=param)
330 331 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
331 332 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
332 333
333 334 desc = {
334 335 'Data': {
335 336 PARAM[param]['name']: PARAM[param]['label'],
336 337 'utctime': 'time'
337 338 },
338 339 'Metadata': {
339 340 'heightList': 'range',
340 341 'data_azi': 'azimuth',
341 342 'data_ele': 'elevation',
342 343 }
343 344 }
344 345
345 346 if args.save:
346 347 opObj10 = merge.addOperation(name='HDFWriter')
347 348 opObj10.addParameter(name='path',value=path_save, format='str')
348 349 opObj10.addParameter(name='Reset',value=True)
349 350 opObj10.addParameter(name='setType',value='weather')
350 351 opObj10.addParameter(name='description',value='desc')
351 352 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
352 353 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
353 354 opObj10.addParameter(name='dataList',value='{},utctime'.format(PARAM[param]['name']))
354 355
355 356 project.start()
356 357
357 358 if __name__ == '__main__':
358 359
359 360 parser = argparse.ArgumentParser(description='Script to process SOPHy data.')
360 361 parser.add_argument('experiment',
361 362 help='Experiment name')
362 363 parser.add_argument('--parameters', nargs='*', default=['P'],
363 364 help='Variables to process: P, Z, V')
364 365 parser.add_argument('--time_offset', default=0,
365 366 help='Fix time offset')
366 367 parser.add_argument('--range', default=0, type=float,
367 368 help='Max range to plot')
368 369 parser.add_argument('--save', action='store_true',
369 370 help='Create output files')
370 371 parser.add_argument('--show', action='store_true',
371 372 help='Show matplotlib plot.')
372 373 parser.add_argument('--online', action='store_true',
373 374 help='Set online mode.')
374 375
375 376 args = parser.parse_args()
376 377
377 378 main(args)
@@ -1,376 +1,376
1 1 # SOPHY PROC script
2 2 import os, sys, json, argparse
3 3 import datetime
4 4 import time
5 5
6 6 PATH = '/DATA_RM/DATA'
7 7 # PATH = '/Users/jespinoza/workspace/data/'
8 PATH = '/home/soporte/Documents/HUANCAYO'
8 #PATH = '/home/soporte/Documents/HUANCAYO'
9 9 PARAM = {
10 10 'P': {'name': 'dataPP_POWER', 'zmin': -50, 'zmax': -15, 'colormap': 'jet', 'label': 'Power', 'wrname': 'Pow','cb_label': 'dB', 'ch':0},
11 11 'V': {'name': 'dataPP_DOP', 'zmin': -20, 'zmax': 20, 'colormap': 'seismic', 'label': 'Velocity', 'wrname': 'Vel', 'cb_label': 'm/s', 'ch':0},
12 12 'RH': {'name': 'RhoHV_R', 'zmin': 0, 'zmax': 1, 'colormap': 'jet', 'label': 'Coef.Correlacion', 'wrname':'R', 'cb_label': '*', 'ch':0},
13 13 'FD': {'name': 'PhiD_P', 'zmin': -180,'zmax': 180,'colormap': 'RdBu_r', 'label': 'Fase Diferencial', 'wrname':'P' , 'cb_label': 'ΒΊ', 'ch':0},
14 14 'ZD': {'name': 'Zdb_D', 'zmin': -20, 'zmax': 60, 'colormap': 'viridis','label': 'Reflect.Diferencial','wrname':'D' , 'cb_label': 'dBz','ch':0},
15 15 'Z': {'name': 'Zdb', 'zmin': -20, 'zmax': 70, 'colormap': 'gist_ncar','label': 'Reflectividad', 'wrname':'Z', 'cb_label': 'dBz','ch':1},
16 16 'W': {'name': 'Sigmav_W', 'zmin': 0, 'zmax':5, 'colormap': 'viridis','label': 'AnchoEspectral', 'wrname':'S', 'cb_label': 'hz', 'ch':1}
17 17 }
18 18
19 19 def max_index(r, sample_rate, ipp):
20 20
21 21 return int(sample_rate*ipp*1e6 * r / 60) + int(sample_rate*ipp*1e6 * 1.2 / 60)
22 22
23 23 def main(args):
24 24
25 25 experiment = args.experiment
26 26 fp = open(os.path.join(PATH, experiment, 'experiment.conf'))
27 27 conf = json.loads(fp.read())
28 28
29 29 ipp_km = conf['usrp_tx']['ipp']
30 30 ipp = ipp_km * 2 /300000
31 31 sample_rate = conf['usrp_rx']['sample_rate']
32 32 axis = ['0' if x=='elevation' else '1' for x in conf['pedestal']['axis']] # AZIMUTH 1 ELEVACION 0
33 33 speed_axis = conf['pedestal']['speed']
34 34 steeps = conf['pedestal']['table']
35 35 time_offset = args.time_offset
36 36 parameters = args.parameters
37 37 start_date = experiment.split('@')[1].split('T')[0].replace('-', '/')
38 38 end_date = start_date
39 39 start_time = experiment.split('@')[1].split('T')[1].replace('-', ':')
40 40 end_time = '23:59:59'
41 41 N = int(1/(speed_axis[0]*ipp)) # 1 GRADO DE RESOLUCION
42 42 path = os.path.join(PATH, experiment, 'rawdata')
43 43 path_ped = os.path.join(PATH, experiment, 'position')
44 44 path_plots = os.path.join(PATH, experiment, 'plots_ch1')
45 45 path_save = os.path.join(PATH, experiment, 'param')
46 46 RMIX = 2
47 47
48 48 from schainpy.controller import Project
49 49
50 50 project = Project()
51 51 project.setup(id='1', name='Sophy', description='sophy proc')
52 52
53 53 reader = project.addReadUnit(datatype='DigitalRFReader',
54 54 path=path,
55 55 startDate=start_date,
56 56 endDate=end_date,
57 57 startTime=start_time,
58 58 endTime=end_time,
59 59 delay=30,
60 60 channelList='0',
61 61 online=args.online,
62 62 walk=1,
63 63 ippKm = ipp_km,
64 64 getByBlock = 1,
65 65 nProfileBlocks = N,
66 66 )
67 67
68 68 if not conf['usrp_tx']['enable_2']: # One Pulse
69 69 voltage = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
70 70
71 71 if conf['usrp_tx']['code_type_1']:
72 72 code = [c.split() for c in conf['usrp']['code_1']]
73 73 op = voltage.addOperation(name='Decoder', optype='other')
74 74 op.addParameter(name='code', value=code)
75 75 op.addParameter(name='nCode', value=len(code), format='int')
76 76 op.addParameter(name='nBaud', value=len(code[0]), format='int')
77 77
78 78 op = voltage.addOperation(name='setH0')
79 79 op.addParameter(name='h0', value='-1.2')
80 80
81 81 if args.range >= 0:
82 82 op = voltage.addOperation(name='selectHeights')
83 83 op.addParameter(name='minIndex', value='0', format='int')
84 84 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
85 85
86 86 code=[[1]]
87 87 opObj11 = voltage.addOperation(name='Decoder', optype='other')
88 88 opObj11.addParameter(name='code', value=code)
89 89 opObj11.addParameter(name='nCode', value='1', format='int')
90 90 opObj11.addParameter(name='nBaud', value='1', format='int')
91 91
92 92 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
93 93 op.addParameter(name='n', value=2*len(code), format='int')
94 94
95 95 #op = voltage.addOperation(name='PulsePair_vRF', optype='other')
96 96 #op.addParameter(name='n', value=int(N), format='int')
97 97
98 98 if args.range >= 0:
99 99 op = voltage.addOperation(name='selectHeights')
100 100 op.addParameter(name='minIndex', value='0', format='int')
101 101 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
102 102
103 103
104 104 op = voltage.addOperation(name='PulsePair_vRF', optype='other')
105 105 op.addParameter(name='n', value=125, format='int')
106 106
107 107
108 108 proc = project.addProcUnit(datatype='ParametersProc', inputId=voltage.getId())
109 109 #procUnitConfObjB.addParameter(name='runNextUnit', value=True)
110 110
111 111 opObj10 = proc.addOperation(name="WeatherRadar")
112 112 opObj10.addParameter(name='variableList',value='Reflectividad,VelocidadRadial,AnchoEspectral')
113 113
114 114 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
115 115
116 116 op = proc.addOperation(name='PedestalInformation')
117 117 op.addParameter(name='path', value=path_ped, format='str')
118 118 op.addParameter(name='interval', value='0.04')
119 119 op.addParameter(name='time_offset', value=time_offset)
120 120 op.addParameter(name='az_offset', value=-26.2)
121 121
122 122 for param in parameters:
123 123 op = proc.addOperation(name='Block360_vRF4')
124 124 #op.addParameter(name='axis', value=','.join(axis))
125 125 op.addParameter(name='attr_data', value=PARAM[param]['name'])
126 126 op.addParameter(name='runNextOp', value=True)
127 127
128 128 op= proc.addOperation(name='WeatherParamsPlot')
129 129 if args.save: op.addParameter(name='save', value=path_plots, format='str')
130 130 op.addParameter(name='save_period', value=-1)
131 131 op.addParameter(name='show', value=args.show)
132 132 op.addParameter(name='channels', value='1,')
133 133 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
134 134 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
135 135 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
136 136 op.addParameter(name='labels', value=[PARAM[param]['label']])
137 137 op.addParameter(name='save_code', value=param)
138 138 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
139 139 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
140 140
141 141 desc = {
142 142 'Data': {
143 143 PARAM[param]['name']: PARAM[param]['label'],
144 144 'utctime': 'time'
145 145 },
146 146 'Metadata': {
147 147 'heightList': 'range',
148 148 'data_azi': 'azimuth',
149 149 'data_ele': 'elevation',
150 150 }
151 151 }
152 152
153 153 if args.save:
154 154 opObj10 = proc.addOperation(name='HDFWriter')
155 155 opObj10.addParameter(name='path',value=path_save+'-{}'.format(param), format='str')
156 156 opObj10.addParameter(name='Reset',value=True)
157 157 opObj10.addParameter(name='setType',value='weather')
158 158 opObj10.addParameter(name='description',value='desc')
159 159 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
160 160 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
161 161 opObj10.addParameter(name='dataList',value='{},utctime'.format(PARAM[param]['name']))
162 162
163 163 else: #Two pulses
164 164
165 165 voltage1 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
166 166
167 167 print("repetions",conf['usrp_tx']['repetitions_1'])
168 168
169 169 op = voltage1.addOperation(name='ProfileSelector')
170 170 op.addParameter(name='profileRangeList', value='0,{}'.format(conf['usrp_tx']['repetitions_1']-1))
171 171
172 172
173 173 #op3 = voltage1.addOperation(name='ProfileSelector', optype='other')
174 174 #op3.addParameter(name='profileRangeList', value='1,123')
175 175
176 176 '''
177 177 if conf['usrp_tx']['code_type_1'] != 'None':
178 178 code = [c.split() for c in conf['usrp_tx']['code_1']]
179 179 op = voltage1.addOperation(name='Decoder', optype='other')
180 180 op.addParameter(name='code', value=code)
181 181 op.addParameter(name='nCode', value=len(code), format='int')
182 182 op.addParameter(name='nBaud', value=len(code[0]), format='int')
183 183 '''
184 184
185 185 code=[[1]]
186 186
187 187 opObj11 = voltage1.addOperation(name='Decoder', optype='other')
188 188 opObj11.addParameter(name='code', value=code)
189 189 opObj11.addParameter(name='nCode', value='1', format='int')
190 190 opObj11.addParameter(name='nBaud', value='1', format='int')
191 191
192 192 op = voltage1.addOperation(name='setH0')
193 193 op.addParameter(name='h0', value='-1.2')
194 194
195 195 if args.range >= 0:
196 196 op = voltage1.addOperation(name='selectHeights')
197 197 op.addParameter(name='minIndex', value='0', format='int')
198 198 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
199 199
200 200 op = voltage1.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
201 201 op.addParameter(name='n', value=2, format='int')
202 202
203 203 op = voltage1.addOperation(name='PulsePair_vRF', optype='other')
204 204 #op.addParameter(name='n', value=int(N), format='int')
205 205 op.addParameter(name='n', value=61, format='int')
206 206 #op.addParameter(name='removeDC',value=True)
207 207
208 208 '''
209 209 if args.range >= 0:
210 210 print("corto",max_index(RMIX, sample_rate, ipp))
211 211 op = voltage1.addOperation(name='selectHeights')
212 212 op.addParameter(name='minIndex', value='0', format='int')
213 213 op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
214 214 '''
215 215 proc1 = project.addProcUnit(datatype='ParametersProc', inputId=voltage1.getId())
216 216 proc1.addParameter(name='runNextUnit', value=True)
217 217
218 218 opObj10 = proc1.addOperation(name="WeatherRadar")
219 219 opObj10.addParameter(name='variableList',value='Reflectividad,VelocidadRadial,AnchoEspectral')
220 220 opObj10.addParameter(name='tauW',value=0.4*1e-6)
221 221 opObj10.addParameter(name='Pt',value=0.2)
222 222
223 223 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
224 224
225 225 op = proc1.addOperation(name='PedestalInformation')
226 226 op.addParameter(name='path', value=path_ped, format='str')
227 227 op.addParameter(name='interval', value='0.04')
228 228 op.addParameter(name='time_offset', value=time_offset)
229 229 op.addParameter(name='az_offset', value=-26.2)
230 230
231 231 for param in parameters:
232 232 op = proc1.addOperation(name='Block360_vRF4')
233 233 op.addParameter(name='attr_data', value=PARAM[param]['name'])
234 234 op.addParameter(name='runNextOp', value=True)
235 235
236 236 voltage2 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
237 237
238 238 op = voltage2.addOperation(name='ProfileSelector')
239 239 op.addParameter(name='profileRangeList', value='{},{}'.format(conf['usrp_tx']['repetitions_1'], conf['usrp_tx']['repetitions_1']+conf['usrp_tx']['repetitions_2']-1))
240 240
241 241
242 242 if conf['usrp_tx']['code_type_2']:
243 243 print(conf['usrp_tx']['code_2'])
244 244 codes = [ c.strip() for c in conf['usrp_tx']['code_2'].split(',')]
245 245 code = []
246 246 for c in codes:
247 247 code.append([int(x) for x in c])
248 248 print(code)
249 249 print(code[0])
250 250 op = voltage2.addOperation(name='Decoder', optype='other')
251 251 op.addParameter(name='code', value=code)
252 252 op.addParameter(name='nCode', value=len(code), format='int')
253 253 op.addParameter(name='nBaud', value=len(code[0]), format='int')
254 254 import numpy
255 255 pwcode = numpy.sum(numpy.array(code[0])**2)
256 256 print("pwcode",pwcode)
257 257
258 258 op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
259 259 op.addParameter(name='n', value=len(code), format='int')
260 260 ncode = len(code)
261 261 else:
262 262 ncode = 1
263 263
264 264 op = voltage2.addOperation(name='setH0')
265 265 op.addParameter(name='h0', value='-1.2')
266 266
267 267 if args.range >= 0:
268 268 if args.range==0:
269 269 args.range= ipp_km
270 270 op = voltage2.addOperation(name='selectHeights')
271 271 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
272 272 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
273 273
274 274 #op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
275 275 #op.addParameter(name='n', value=int(N)/ncode, format='int')
276 276 op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
277 277 op.addParameter(name='n', value=64, format='int')
278 278 #op.addParameter(name='removeDC',value=True)
279 279
280 280 '''
281 281
282 282 if args.range >= 0:
283 283 if args.range==0:
284 284 args.range= ipp_km
285 285 op = voltage2.addOperation(name='selectHeights')
286 286 print("largo",max_index(RMIX, sample_rate, ipp))
287 287 print("largo2",max_index(args.range, sample_rate, ipp))
288 288
289 289 op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
290 290 op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
291 291 '''
292 292
293 293 proc2 = project.addProcUnit(datatype='ParametersProc', inputId=voltage2.getId())
294 294
295 295 opObj10 = proc2.addOperation(name="WeatherRadar")
296 296 opObj10.addParameter(name='variableList',value='Reflectividad,AnchoEspectral')
297 297 opObj10.addParameter(name='tauW',value=3.2*1e-6)
298 298 opObj10.addParameter(name='Pt',value=1.6)
299 299
300 300
301 301 # {"latitude": -12.0404828587, "longitude": -75.2147483647, "altitude": 3379.2147483647}
302 302
303 303 op = proc2.addOperation(name='PedestalInformation')
304 304 op.addParameter(name='path', value=path_ped, format='str')
305 305 op.addParameter(name='interval', value='0.04')
306 306 op.addParameter(name='time_offset', value=time_offset)
307 307 op.addParameter(name='az_offset', value=-26.2)
308 308
309 309 for param in parameters:
310 310 op = proc2.addOperation(name='Block360_vRF4')
311 311 #op.addParameter(name='axis', value=','.join(axis))
312 312 op.addParameter(name='attr_data', value=PARAM[param]['name'])
313 313 op.addParameter(name='runNextOp', value=True)
314 314
315 315 merge = project.addProcUnit(datatype='MergeProc', inputId=[proc1.getId(), proc2.getId()])
316 316 merge.addParameter(name='attr_data', value=PARAM[param]['name'])
317 317 merge.addParameter(name='mode', value='7') #RM
318 318
319 319 op= merge.addOperation(name='WeatherParamsPlot')
320 320 if args.save: op.addParameter(name='save', value=path_plots, format='str')
321 321 op.addParameter(name='save_period', value=-1)
322 322 op.addParameter(name='show', value=args.show)
323 323 op.addParameter(name='channels', value='1,')
324 324 op.addParameter(name='zmin', value=PARAM[param]['zmin'])
325 325 op.addParameter(name='zmax', value=PARAM[param]['zmax'])
326 326 op.addParameter(name='attr_data', value=PARAM[param]['name'], format='str')
327 327 op.addParameter(name='labels', value=[PARAM[param]['label']])
328 328 op.addParameter(name='save_code', value=param)
329 329 op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
330 330 op.addParameter(name='colormap', value=PARAM[param]['colormap'])
331 331
332 332 desc = {
333 333 'Data': {
334 334 PARAM[param]['name']: PARAM[param]['label'],
335 335 'utctime': 'time'
336 336 },
337 337 'Metadata': {
338 338 'heightList': 'range',
339 339 'data_azi': 'azimuth',
340 340 'data_ele': 'elevation',
341 341 }
342 342 }
343 343
344 344 if args.save:
345 345 opObj10 = merge.addOperation(name='HDFWriter')
346 346 opObj10.addParameter(name='path',value=path_save, format='str')
347 347 opObj10.addParameter(name='Reset',value=True)
348 348 opObj10.addParameter(name='setType',value='weather')
349 349 opObj10.addParameter(name='description',value='desc')
350 350 opObj10.addParameter(name='blocksPerFile',value='1',format='int')
351 351 opObj10.addParameter(name='metadataList',value='heightList,data_azi,data_ele')
352 352 opObj10.addParameter(name='dataList',value='{},utctime'.format(PARAM[param]['name']))
353 353
354 354 project.start()
355 355
356 356 if __name__ == '__main__':
357 357
358 358 parser = argparse.ArgumentParser(description='Script to process SOPHy data.')
359 359 parser.add_argument('experiment',
360 360 help='Experiment name')
361 361 parser.add_argument('--parameters', nargs='*', default=['P'],
362 362 help='Variables to process: P, Z, V')
363 363 parser.add_argument('--time_offset', default=0,
364 364 help='Fix time offset')
365 365 parser.add_argument('--range', default=0, type=float,
366 366 help='Max range to plot')
367 367 parser.add_argument('--save', action='store_true',
368 368 help='Create output files')
369 369 parser.add_argument('--show', action='store_true',
370 370 help='Show matplotlib plot.')
371 371 parser.add_argument('--online', action='store_true',
372 372 help='Set online mode.')
373 373
374 374 args = parser.parse_args()
375 375
376 376 main(args)
General Comments 0
You need to be logged in to leave comments. Login now