##// END OF EJS Templates
schain
RSS
LAST UPDATE INTEGRACION
avaldez -
r1388:ec0a56a3ee09
parent child
Show More
Side by Side Unified
Context file:
r1388:ec0a56a3ee09 Loading diff...:

The requested changes are too big and content was truncated. Show full diff

Show file before | Show file after | Hide comments
@@ -0,0 +1,213
1 # Ing. AVP
2 # 06/10/2021
3 # ARCHIVO DE LECTURA
4 import os, sys
5 import datetime
6 import time
7 from schainpy.controller import Project
8 #### NOTA###########################################
9 # INPUT :
10 # VELOCIDAD PARAMETRO : V = 2Β°/seg
11 # MODO PULSE PAIR O MOMENTOS: 0 : Pulse Pair ,1 : Momentos
12 ######################################################
13 ##### PROCESAMIENTO ##################################
14 ##### OJO TENER EN CUENTA EL n= para el Pulse Pair ##
15 ##### O EL n= nFFTPoints ###
16 ######################################################
17 ######## BUSCAMOS EL numero de IPP equivalente 1Β°#####
18 ######## Sea V la velocidad del Pedestal en Β°/seg#####
19 ######## 1Β° sera Recorrido en un tiempo de 1/V ######
20 ######## IPP del Radar 400 useg --> 60 Km ############
21 ######## n = 1/(V(Β°/seg)*IPP(Km)) , NUMERO DE IPP ##
22 ######## n = 1/(V*IPP) #############################
23 ######## VELOCIDAD DEL PEDESTAL ######################
24 print("SETUP- RADAR METEOROLOGICO")
25 V = 10
26 mode = 1
27 #path = '/DATA_RM/23/6v'
28 path = '/DATA_RM/TEST_INTEGRACION_2M'
29 path_ped='/DATA_RM/TEST_PEDESTAL/P20211012-082745'
30 figpath_pp = "/home/soporte/Pictures/TEST_PP"
31 figpath_mom = "/home/soporte/Pictures/TEST_MOM"
32 plot = 0
33 integration = 1
34 save = 0
35 if save == 1:
36 if mode==0:
37 path_save = '/DATA_RM/TEST_HDF5_PP_23/6v'
38 path_save = '/DATA_RM/TEST_HDF5_PP'
39 path_save = '/DATA_RM/TEST_HDF5_PP_100'
40 else:
41 path_save = '/DATA_RM/TEST_HDF5_SPEC_23_V2/6v'
42
43 print("* PATH data ADQ :", path)
44 print("* Velocidad Pedestal :",V,"Β°/seg")
45 ############################ NRO Perfiles PROCESAMIENTO ###################
46 V=V
47 IPP=400*1e-6
48 n= int(1/(V*IPP))
49 print("* n - NRO Perfiles Proc:", n )
50 ################################## MODE ###################################
51 print("* Modo de Operacion :",mode)
52 if mode ==0:
53 print("* Met. Seleccionado : Pulse Pair")
54 else:
55 print("* Met. Momentos : Momentos")
56
57 ################################## MODE ###################################
58 print("* Grabado de datos :",save)
59 if save ==1:
60 if mode==0:
61 ope= "Pulse Pair"
62 else:
63 ope= "Momentos"
64 print("* Path-Save Data -", ope , path_save)
65
66 print("* Integracion de datos :",integration)
67
68 time.sleep(15)
69 #remotefolder = "/home/wmaster/graficos"
70 #######################################################################
71 ################# RANGO DE PLOTEO######################################
72 dBmin = '1'
73 dBmax = '85'
74 xmin = '15'
75 xmax = '15.25'
76 ymin = '0'
77 ymax = '600'
78 #######################################################################
79 ########################FECHA##########################################
80 str = datetime.date.today()
81 today = str.strftime("%Y/%m/%d")
82 str2 = str - datetime.timedelta(days=1)
83 yesterday = str2.strftime("%Y/%m/%d")
84 #######################################################################
85 ########################SIGNAL CHAIN ##################################
86 #######################################################################
87 desc = "USRP_test"
88 filename = "USRP_processing.xml"
89 controllerObj = Project()
90 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
91 #######################################################################
92 ######################## UNIDAD DE LECTURA#############################
93 #######################################################################
94 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
95 path=path,
96 startDate="2021/01/01",#today,
97 endDate="2021/12/30",#today,
98 startTime='00:00:00',
99 endTime='23:59:59',
100 delay=0,
101 #set=0,
102 online=0,
103 walk=1,
104 ippKm = 60)
105
106 opObj11 = readUnitConfObj.addOperation(name='printInfo')
107
108 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
109
110 if mode ==0:
111 ####################### METODO PULSE PAIR ######################################################################
112 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
113 opObj11.addParameter(name='n', value=int(n), format='int')#10 VOY A USAR 250 DADO QUE LA VELOCIDAD ES 10 GRADOS
114 #opObj11.addParameter(name='removeDC', value=1, format='int')
115 ####################### METODO Parametros ######################################################################
116 procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
117 if plot==1:
118 opObj11 = procUnitConfObjB.addOperation(name='GenericRTIPlot',optype='external')
119 opObj11.addParameter(name='attr_data', value='dataPP_POW')
120 opObj11.addParameter(name='colormap', value='jet')
121 opObj11.addParameter(name='xmin', value=xmin)
122 opObj11.addParameter(name='xmax', value=xmax)
123 opObj11.addParameter(name='zmin', value=dBmin)
124 opObj11.addParameter(name='zmax', value=dBmax)
125 opObj11.addParameter(name='save', value=figpath_pp)
126 opObj11.addParameter(name='showprofile', value=0)
127 opObj11.addParameter(name='save_period', value=50)
128
129 ####################### METODO ESCRITURA #######################################################################
130 if save==1:
131 opObj10 = procUnitConfObjB.addOperation(name='HDFWriter')
132 opObj10.addParameter(name='path',value=path_save)
133 #opObj10.addParameter(name='mode',value=0)
134 opObj10.addParameter(name='blocksPerFile',value='100',format='int')
135 opObj10.addParameter(name='metadataList',value='utctimeInit,timeZone,paramInterval,profileIndex,channelList,heightList,flagDataAsBlock',format='list')
136 opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,utctime',format='list')#,format='list'
137 if integration==1:
138 V=10
139 blocksPerfile=360
140 print("* Velocidad del Pedestal:",V)
141 tmp_blocksPerfile = 100
142 f_a_p= int(tmp_blocksPerfile/V)
143
144 opObj11 = procUnitConfObjB.addOperation(name='PedestalInformation')
145 opObj11.addParameter(name='path_ped', value=path_ped)
146 #opObj11.addParameter(name='path_adq', value=path_adq)
147 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
148 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
149 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
150 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
151 opObj11.addParameter(name='online', value='0', format='int')
152
153 opObj11 = procUnitConfObjB.addOperation(name='Block360')
154 opObj11.addParameter(name='n', value='10', format='int')
155 opObj11.addParameter(name='mode', value=mode, format='int')
156
157 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
158
159 opObj11= procUnitConfObjB.addOperation(name='WeatherPlot',optype='other')
160
161
162 else:
163 ####################### METODO SPECTROS ######################################################################
164 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
165 procUnitConfObjB.addParameter(name='nFFTPoints', value=n, format='int')
166 procUnitConfObjB.addParameter(name='nProfiles' , value=n, format='int')
167
168 procUnitConfObjC = controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjB.getId())
169 procUnitConfObjC.addOperation(name='SpectralMoments')
170 if plot==1:
171 dBmin = '1'
172 dBmax = '65'
173 opObj11 = procUnitConfObjC.addOperation(name='PowerPlot',optype='external')
174 opObj11.addParameter(name='xmin', value=xmin)
175 opObj11.addParameter(name='xmax', value=xmax)
176 opObj11.addParameter(name='zmin', value=dBmin)
177 opObj11.addParameter(name='zmax', value=dBmax)
178 opObj11.addParameter(name='save', value=figpath_mom)
179 opObj11.addParameter(name='showprofile', value=0)
180 opObj11.addParameter(name='save_period', value=100)
181
182 if save==1:
183 opObj10 = procUnitConfObjC.addOperation(name='HDFWriter')
184 opObj10.addParameter(name='path',value=path_save)
185 #opObj10.addParameter(name='mode',value=0)
186 opObj10.addParameter(name='blocksPerFile',value='360',format='int')
187 #opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
188 opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
189 opObj10.addParameter(name='dataList',value='data_pow,data_dop,utctime',format='list')#,format='list'
190
191 if integration==1:
192 V=10
193 blocksPerfile=360
194 print("* Velocidad del Pedestal:",V)
195 tmp_blocksPerfile = 100
196 f_a_p= int(tmp_blocksPerfile/V)
197
198 opObj11 = procUnitConfObjC.addOperation(name='PedestalInformation')
199 opObj11.addParameter(name='path_ped', value=path_ped)
200 #opObj11.addParameter(name='path_adq', value=path_adq)
201 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
202 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
203 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
204 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
205 opObj11.addParameter(name='online', value='0', format='int')
206
207 opObj11 = procUnitConfObjC.addOperation(name='Block360')
208 opObj11.addParameter(name='n', value='10', format='int')
209 opObj11.addParameter(name='mode', value=mode, format='int')
210
211 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
212 opObj11= procUnitConfObjC.addOperation(name='WeatherPlot',optype='other')
213 controllerObj.start()
Show file before | Show file after | Hide comments
@@ -0,0 +1,112
1 # Ing-AlexanderValdez
2 # Monitoreo de Pedestal
3
4 ############## IMPORTA LIBRERIAS ###################
5 import os,numpy,h5py
6 import sys,time
7 import matplotlib.pyplot as plt
8 ####################################################
9 path_ped = '/DATA_RM/TEST_PEDESTAL/P20211012-082745'
10 # Metodo para verificar numero
11 def isNumber(str):
12 try:
13 float(str)
14 return True
15 except:
16 return False
17 # Metodo para extraer el arreglo
18 def getDatavaluefromDirFilename(path,file,value):
19 dir_file= path+"/"+file
20 fp = h5py.File(dir_file,'r')
21 array = fp['Data'].get(value)[()]
22 fp.close()
23 return array
24
25 # LISTA COMPLETA DE ARCHIVOS HDF5 Pedestal
26 LIST= sorted(os.listdir(path_ped))
27 m=len(LIST)
28 print("TOTAL DE ARCHIVOS DE PEDESTAL:",m)
29 # Contadores temporales
30 k= 0
31 l= 0
32 t= 0
33 # Marca de tiempo temporal
34 time_ = numpy.zeros([m])
35 # creacion de
36 for i in range(m):
37 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_pos")
38 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="ele_pos")
39 tmp_azi_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_vel")
40 tmp_ele_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_vel")# nuevo :D
41
42 time_[i] = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="utc")
43
44 k=k +tmp_azi_pos.shape[0]
45 l=l +tmp_ele_pos.shape[0]
46 t=t +tmp_azi_vel.shape[0]
47
48 print("TOTAL DE MUESTRAS, ARCHIVOS X100:",k)
49 time.sleep(5)
50 ######CREACION DE ARREGLOS CANTIDAD DE VALORES POR MUESTRA#################
51 azi_pos = numpy.zeros([k])
52 ele_pos = numpy.zeros([l])
53 time_azi_pos= numpy.zeros([k])
54 # Contadores temporales
55 p=0
56 r=0
57 z=0
58 # VARIABLES TMP para almacenar azimuth, elevacion y tiempo
59
60 #for filename in sorted(os.listdir(path_ped)):
61 # CONDICION POR LEER EN TIEMPO REAL NO OFFLINE
62
63 for filename in LIST:
64 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="azi_pos")
65 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="ele_pos")
66 # CONDICION POR LEER EN TIEMPO REAL NO OFFLINE
67
68 if z==(m-1):
69 tmp_azi_time=numpy.arange(time_[z],time_[z]+1,1/(tmp_azi_pos.shape[0]))
70 else:
71 tmp_azi_time=numpy.arange(time_[z],time_[z+1],(time_[z+1]-time_[z])/(tmp_azi_pos.shape[0]))
72
73 print(filename,time_[z])
74 print(z,tmp_azi_pos.shape[0])
75
76 i=0
77 for i in range(tmp_azi_pos.shape[0]):
78 index=p+i
79 azi_pos[index]=tmp_azi_pos[i]
80 time_azi_pos[index]=tmp_azi_time[i]
81 p=p+tmp_azi_pos.shape[0]
82 i=0
83 for i in range(tmp_ele_pos.shape[0]):
84 index=r+i
85 ele_pos[index]=tmp_ele_pos[i]
86 r=r+tmp_ele_pos.shape[0]
87
88
89 z+=1
90
91
92 ######## GRAFIQUEMOS Y VEAMOS LOS DATOS DEL Pedestal
93 fig, ax = plt.subplots(figsize=(16,8))
94 print(time_azi_pos.shape)
95 print(azi_pos.shape)
96 t=numpy.arange(time_azi_pos.shape[0])*0.01/(60.0)
97 plt.plot(t,azi_pos,label='AZIMUTH_POS',color='blue')
98
99 # AQUI ESTOY ADICIONANDO LA POSICION EN elevaciont=numpy.arange(len(ele_pos))*0.01/60.0
100 t=numpy.arange(len(ele_pos))*0.01/60.0
101 plt.plot(t,ele_pos,label='ELEVATION_POS',color='red')#*10
102
103 #ax.set_xlim(0, 9)
104 ax.set_ylim(-5, 400)
105 plt.ylabel("Azimuth Position")
106 plt.xlabel("Muestra")
107 plt.title('Azimuth Position vs Muestra ', fontsize=20)
108 axes = plt.gca()
109 axes.yaxis.grid()
110 plt.xticks(fontsize=16)
111 plt.yticks(fontsize=16)
112 plt.show()
Show file before | Show file after | Hide comments
@@ -0,0 +1,90
1 import os,sys,json
2 import datetime
3 import time
4 from schainpy.controller import Project
5 '''
6 NOTA:
7 Este script de prueba.
8 - Unidad del lectura 'HDFReader'.
9 - Unidad de procesamiento ParametersProc
10 - Operacion SpectralMomentsPlot
11
12 '''
13
14 #######################################################################
15 ################# RANGO DE PLOTEO######################################
16 #######################################################################
17 dBmin = '1'
18 dBmax = '65'
19 xmin = '0'
20 xmax ='24'
21 #tmmin = 16.2
22 #tmmax = 16.25
23 tmmin =15
24 tmmax =15.5
25 ymin = '0'
26 ymax = '600'
27 #######################################################################
28 #######################################################################
29 #######################################################################
30 #path = '/DATA_RM/TEST_HDF5_SPEC'
31 path = '/DATA_RM/TEST_HDF5_SPEC_23/6v/'
32 figpath = '/home/soporte/Downloads/23/6v'
33 path="/home/soporte/Downloads/params-20211015T174046Z-001/params"
34 desc = "Simulator Test"
35 desc_data = {
36 'Data': {
37 'data_spc': ['Data/data_spc/channel00','Data/data_spc/channel01'\
38 ,'Data/data_spc/channel02','Data/data_spc/channel03'\
39 ,'Data/data_spc/channel04','Data/data_spc/channel05'\
40 ,'Data/data_spc/channel06','Data/data_spc/channel07'\
41 ,'Data/data_spc/channel08','Data/data_spc/channel09'],
42 'utctime':'Data/utctime'
43 },
44 'Metadata': {
45 'type' :'Metadata/type',
46 'channelList' :'Metadata/channelList',
47 'heightList' :'Metadata/heightList',
48 'ippSeconds' :'Metadata/ippSeconds',
49 'nProfiles' :'Metadata/nProfiles',
50 'codeList' :'Metadata/codeList',
51 'timeZone' :'Metadata/timeZone',
52 'azimuthList' :'Metadata/azimuthList',
53 'elevationList' :'Metadata/elevationList',
54 'nCohInt' :'Metadata/nCohInt',
55 'nIncohInt' :'Metadata/nIncohInt',
56 'nFFTPoints' :'Metadata/nFFTPoints'
57
58 }
59 }
60
61 controllerObj = Project()
62
63 controllerObj.setup(id='10',name='Test Simulator',description=desc)
64
65 readUnitConfObj = controllerObj.addReadUnit(datatype='HDFReader',
66 path=path,
67 startDate="2021/01/01", #"2020/01/01",#today,
68 endDate= "2021/12/01", #"2020/12/30",#today,
69 startTime='00:00:00',
70 endTime='23:59:59',
71 delay=0,
72 #set=0,
73 online=0,
74 walk=0,
75 description= json.dumps(desc_data))#1
76
77 procUnitConfObjA = controllerObj.addProcUnit(datatype='ParametersProc',inputId=readUnitConfObj.getId())
78 procUnitConfObjA.addOperation(name='SpectralMoments')
79
80 '''
81 opObj11 = readUnitConfObj.addOperation(name='SpectraPlot',optype='external')
82 opObj11.addParameter(name='xmin', value=tmmin)
83 opObj11.addParameter(name='xmax', value=tmmax)
84 opObj11.addParameter(name='zmin', value=dBmin)
85 opObj11.addParameter(name='zmax', value=dBmax)
86 opObj11.addParameter(name='save', value=figpath)
87 opObj11.addParameter(name='showprofile', value=0)
88 opObj11.addParameter(name='save_period', value=10)
89 '''
90 controllerObj.start()
Show file before | Show file after | Hide comments
@@ -1,518 +1,519
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
7 7 from schainpy.utils import log
8 8 # libreria wradlib
9 9 import wradlib as wrl
10 10
11 11 EARTH_RADIUS = 6.3710e3
12 12
13 13
14 14 def ll2xy(lat1, lon1, lat2, lon2):
15 15
16 16 p = 0.017453292519943295
17 17 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
18 18 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
19 19 r = 12742 * numpy.arcsin(numpy.sqrt(a))
20 20 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
21 21 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
22 22 theta = -theta + numpy.pi/2
23 23 return r*numpy.cos(theta), r*numpy.sin(theta)
24 24
25 25
26 26 def km2deg(km):
27 27 '''
28 28 Convert distance in km to degrees
29 29 '''
30 30
31 31 return numpy.rad2deg(km/EARTH_RADIUS)
32 32
33 33
34 34
35 35 class SpectralMomentsPlot(SpectraPlot):
36 36 '''
37 37 Plot for Spectral Moments
38 38 '''
39 39 CODE = 'spc_moments'
40 40 # colormap = 'jet'
41 41 # plot_type = 'pcolor'
42 42
43 43 class DobleGaussianPlot(SpectraPlot):
44 44 '''
45 45 Plot for Double Gaussian Plot
46 46 '''
47 47 CODE = 'gaussian_fit'
48 48 # colormap = 'jet'
49 49 # plot_type = 'pcolor'
50 50
51 51 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
52 52 '''
53 53 Plot SpectraCut with Double Gaussian Fit
54 54 '''
55 55 CODE = 'cut_gaussian_fit'
56 56
57 57 class SnrPlot(RTIPlot):
58 58 '''
59 59 Plot for SNR Data
60 60 '''
61 61
62 62 CODE = 'snr'
63 63 colormap = 'jet'
64 64
65 65 def update(self, dataOut):
66 66
67 67 data = {
68 68 'snr': 10*numpy.log10(dataOut.data_snr)
69 69 }
70 70
71 71 return data, {}
72 72
73 73 class DopplerPlot(RTIPlot):
74 74 '''
75 75 Plot for DOPPLER Data (1st moment)
76 76 '''
77 77
78 78 CODE = 'dop'
79 79 colormap = 'jet'
80 80
81 81 def update(self, dataOut):
82 82
83 83 data = {
84 84 'dop': 10*numpy.log10(dataOut.data_dop)
85 85 }
86 86
87 87 return data, {}
88 88
89 89 class PowerPlot(RTIPlot):
90 90 '''
91 91 Plot for Power Data (0 moment)
92 92 '''
93 93
94 94 CODE = 'pow'
95 95 colormap = 'jet'
96 96
97 97 def update(self, dataOut):
98 98
99 99 data = {
100 100 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
101 101 }
102 102
103 103 return data, {}
104 104
105 105 class SpectralWidthPlot(RTIPlot):
106 106 '''
107 107 Plot for Spectral Width Data (2nd moment)
108 108 '''
109 109
110 110 CODE = 'width'
111 111 colormap = 'jet'
112 112
113 113 def update(self, dataOut):
114 114
115 115 data = {
116 116 'width': dataOut.data_width
117 117 }
118 118
119 119 return data, {}
120 120
121 121 class SkyMapPlot(Plot):
122 122 '''
123 123 Plot for meteors detection data
124 124 '''
125 125
126 126 CODE = 'param'
127 127
128 128 def setup(self):
129 129
130 130 self.ncols = 1
131 131 self.nrows = 1
132 132 self.width = 7.2
133 133 self.height = 7.2
134 134 self.nplots = 1
135 135 self.xlabel = 'Zonal Zenith Angle (deg)'
136 136 self.ylabel = 'Meridional Zenith Angle (deg)'
137 137 self.polar = True
138 138 self.ymin = -180
139 139 self.ymax = 180
140 140 self.colorbar = False
141 141
142 142 def plot(self):
143 143
144 144 arrayParameters = numpy.concatenate(self.data['param'])
145 145 error = arrayParameters[:, -1]
146 146 indValid = numpy.where(error == 0)[0]
147 147 finalMeteor = arrayParameters[indValid, :]
148 148 finalAzimuth = finalMeteor[:, 3]
149 149 finalZenith = finalMeteor[:, 4]
150 150
151 151 x = finalAzimuth * numpy.pi / 180
152 152 y = finalZenith
153 153
154 154 ax = self.axes[0]
155 155
156 156 if ax.firsttime:
157 157 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
158 158 else:
159 159 ax.plot.set_data(x, y)
160 160
161 161 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
162 162 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
163 163 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
164 164 dt2,
165 165 len(x))
166 166 self.titles[0] = title
167 167
168 168
169 169 class GenericRTIPlot(Plot):
170 170 '''
171 171 Plot for data_xxxx object
172 172 '''
173 173
174 174 CODE = 'param'
175 175 colormap = 'viridis'
176 176 plot_type = 'pcolorbuffer'
177 177
178 178 def setup(self):
179 179 self.xaxis = 'time'
180 180 self.ncols = 1
181 181 self.nrows = self.data.shape('param')[0]
182 182 self.nplots = self.nrows
183 183 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
184 184
185 185 if not self.xlabel:
186 186 self.xlabel = 'Time'
187 187
188 188 self.ylabel = 'Range [km]'
189 189 if not self.titles:
190 190 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
191 191
192 192 def update(self, dataOut):
193 193
194 194 data = {
195 195 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
196 196 }
197 197
198 198 meta = {}
199 199
200 200 return data, meta
201 201
202 202 def plot(self):
203 203 # self.data.normalize_heights()
204 204 self.x = self.data.times
205 205 self.y = self.data.yrange
206 206 self.z = self.data['param']
207 207
208 208 self.z = 10*numpy.log10(self.z)
209 209
210 210 self.z = numpy.ma.masked_invalid(self.z)
211 211
212 212 if self.decimation is None:
213 213 x, y, z = self.fill_gaps(self.x, self.y, self.z)
214 214 else:
215 215 x, y, z = self.fill_gaps(*self.decimate())
216 216
217 217 for n, ax in enumerate(self.axes):
218 218
219 219 self.zmax = self.zmax if self.zmax is not None else numpy.max(
220 220 self.z[n])
221 221 self.zmin = self.zmin if self.zmin is not None else numpy.min(
222 222 self.z[n])
223 223
224 224 if ax.firsttime:
225 225 if self.zlimits is not None:
226 226 self.zmin, self.zmax = self.zlimits[n]
227 227
228 228 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
229 229 vmin=self.zmin,
230 230 vmax=self.zmax,
231 231 cmap=self.cmaps[n]
232 232 )
233 233 else:
234 234 if self.zlimits is not None:
235 235 self.zmin, self.zmax = self.zlimits[n]
236 236 ax.collections.remove(ax.collections[0])
237 237 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
238 238 vmin=self.zmin,
239 239 vmax=self.zmax,
240 240 cmap=self.cmaps[n]
241 241 )
242 242
243 243
244 244 class PolarMapPlot(Plot):
245 245 '''
246 246 Plot for weather radar
247 247 '''
248 248
249 249 CODE = 'param'
250 250 colormap = 'seismic'
251 251
252 252 def setup(self):
253 253 self.ncols = 1
254 254 self.nrows = 1
255 255 self.width = 9
256 256 self.height = 8
257 257 self.mode = self.data.meta['mode']
258 258 if self.channels is not None:
259 259 self.nplots = len(self.channels)
260 260 self.nrows = len(self.channels)
261 261 else:
262 262 self.nplots = self.data.shape(self.CODE)[0]
263 263 self.nrows = self.nplots
264 264 self.channels = list(range(self.nplots))
265 265 if self.mode == 'E':
266 266 self.xlabel = 'Longitude'
267 267 self.ylabel = 'Latitude'
268 268 else:
269 269 self.xlabel = 'Range (km)'
270 270 self.ylabel = 'Height (km)'
271 271 self.bgcolor = 'white'
272 272 self.cb_labels = self.data.meta['units']
273 273 self.lat = self.data.meta['latitude']
274 274 self.lon = self.data.meta['longitude']
275 275 self.xmin, self.xmax = float(
276 276 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
277 277 self.ymin, self.ymax = float(
278 278 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
279 279 # self.polar = True
280 280
281 281 def plot(self):
282 282
283 283 for n, ax in enumerate(self.axes):
284 284 data = self.data['param'][self.channels[n]]
285 285
286 286 zeniths = numpy.linspace(
287 287 0, self.data.meta['max_range'], data.shape[1])
288 288 if self.mode == 'E':
289 289 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
290 290 r, theta = numpy.meshgrid(zeniths, azimuths)
291 291 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
292 292 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
293 293 x = km2deg(x) + self.lon
294 294 y = km2deg(y) + self.lat
295 295 else:
296 296 azimuths = numpy.radians(self.data.yrange)
297 297 r, theta = numpy.meshgrid(zeniths, azimuths)
298 298 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
299 299 self.y = zeniths
300 300
301 301 if ax.firsttime:
302 302 if self.zlimits is not None:
303 303 self.zmin, self.zmax = self.zlimits[n]
304 304 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
305 305 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
306 306 vmin=self.zmin,
307 307 vmax=self.zmax,
308 308 cmap=self.cmaps[n])
309 309 else:
310 310 if self.zlimits is not None:
311 311 self.zmin, self.zmax = self.zlimits[n]
312 312 ax.collections.remove(ax.collections[0])
313 313 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
314 314 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
315 315 vmin=self.zmin,
316 316 vmax=self.zmax,
317 317 cmap=self.cmaps[n])
318 318
319 319 if self.mode == 'A':
320 320 continue
321 321
322 322 # plot district names
323 323 f = open('/data/workspace/schain_scripts/distrito.csv')
324 324 for line in f:
325 325 label, lon, lat = [s.strip() for s in line.split(',') if s]
326 326 lat = float(lat)
327 327 lon = float(lon)
328 328 # ax.plot(lon, lat, '.b', ms=2)
329 329 ax.text(lon, lat, label.decode('utf8'), ha='center',
330 330 va='bottom', size='8', color='black')
331 331
332 332 # plot limites
333 333 limites = []
334 334 tmp = []
335 335 for line in open('/data/workspace/schain_scripts/lima.csv'):
336 336 if '#' in line:
337 337 if tmp:
338 338 limites.append(tmp)
339 339 tmp = []
340 340 continue
341 341 values = line.strip().split(',')
342 342 tmp.append((float(values[0]), float(values[1])))
343 343 for points in limites:
344 344 ax.add_patch(
345 345 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
346 346
347 347 # plot Cuencas
348 348 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
349 349 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
350 350 values = [line.strip().split(',') for line in f]
351 351 points = [(float(s[0]), float(s[1])) for s in values]
352 352 ax.add_patch(Polygon(points, ec='b', fc='none'))
353 353
354 354 # plot grid
355 355 for r in (15, 30, 45, 60):
356 356 ax.add_artist(plt.Circle((self.lon, self.lat),
357 357 km2deg(r), color='0.6', fill=False, lw=0.2))
358 358 ax.text(
359 359 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
360 360 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
361 361 '{}km'.format(r),
362 362 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
363 363
364 364 if self.mode == 'E':
365 365 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
366 366 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
367 367 else:
368 368 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
369 369 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
370 370
371 371 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
372 372 self.titles = ['{} {}'.format(
373 373 self.data.parameters[x], title) for x in self.channels]
374 374
375 375 class WeatherPlot(Plot):
376 376 CODE = 'weather'
377 377 plot_name = 'weather'
378 378 plot_type = 'ppistyle'
379 379 buffering = False
380 380
381 381 def setup(self):
382 382 self.ncols = 1
383 383 self.nrows = 1
384 384 self.nplots= 1
385 385 self.ylabel= 'Range [Km]'
386 386 self.titles= ['Weather']
387 387 self.colorbar=False
388 388 self.width =8
389 389 self.height =8
390 390 self.ini =0
391 391 self.len_azi =0
392 392 self.buffer_ini = None
393 393 self.buffer_azi = None
394 394 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
395 395 self.flag =0
396 396 self.indicador= 0
397 397
398 398 def update(self, dataOut):
399 399
400 400 data = {}
401 401 meta = {}
402 data['weather'] = 10*numpy.log10(dataOut.data_360[0]/(250**2))
403 print(data['weather'])
402 print("aprox",dataOut.data_360[0])
403 data['weather'] = 10*numpy.log10(dataOut.data_360[0]/(250.0))
404 #print(data['weather'])
404 405 data['azi'] = dataOut.data_azi
405 406 print("UPDATE",data['azi'])
406 407 return data, meta
407 408
408 409 def const_ploteo(self,data_weather,data_azi,step,res):
409 410 #print("data_weather",data_weather)
410 411 print("data_azi",data_azi)
411 412 print("step",step)
412 413 if self.ini==0:
413 414 #------- AZIMUTH
414 415 n = (360/res)-len(data_azi)
415 416 start = data_azi[-1] + res
416 417 end = data_azi[0] - res
417 418 if start>end:
418 419 end = end + 360
419 420 azi_vacia = numpy.linspace(start,end,int(n))
420 421 azi_vacia = numpy.where(azi_vacia>360,azi_vacia-360,azi_vacia)
421 422 data_azi = numpy.hstack((data_azi,azi_vacia))
422 423 # RADAR
423 424 val_mean = numpy.mean(data_weather[:,0])
424 425 data_weather_cmp = numpy.ones([(360-data_weather.shape[0]),data_weather.shape[1]])*val_mean
425 426 data_weather = numpy.vstack((data_weather,data_weather_cmp))
426 427 else:
427 428 # azimuth
428 429 flag=0
429 430 start_azi = self.res_azi[0]
430 431 start = data_azi[0]
431 432 end = data_azi[-1]
432 433 print("start",start)
433 434 print("end",end)
434 435 if start< start_azi:
435 436 start = start +360
436 437 if end <start_azi:
437 438 end = end +360
438 439
439 440 print("start",start)
440 441 print("end",end)
441 442 #### AQUI SERA LA MAGIA
442 443 pos_ini = int((start-start_azi)/res)
443 444 len_azi = len(data_azi)
444 445 if (360-pos_ini)<len_azi:
445 446 if pos_ini+1==360:
446 447 pos_ini=0
447 448 else:
448 449 flag=1
449 450 dif= 360-pos_ini
450 451 comp= len_azi-dif
451 452
452 453 print(pos_ini)
453 454 print(len_azi)
454 455 print("shape",self.res_azi.shape)
455 456 if flag==0:
456 457 # AZIMUTH
457 458 self.res_azi[pos_ini:pos_ini+len_azi] = data_azi
458 459 # RADAR
459 460 self.res_weather[pos_ini:pos_ini+len_azi,:] = data_weather
460 461 else:
461 462 # AZIMUTH
462 463 self.res_azi[pos_ini:pos_ini+dif] = data_azi[0:dif]
463 464 self.res_azi[0:comp] = data_azi[dif:]
464 465 # RADAR
465 466 self.res_weather[pos_ini:pos_ini+dif,:] = data_weather[0:dif,:]
466 467 self.res_weather[0:comp,:] = data_weather[dif:,:]
467 468 flag=0
468 469 data_azi = self.res_azi
469 470 data_weather = self.res_weather
470 471
471 472 return data_weather,data_azi
472 473
473 474 def plot(self):
474 475 print("--------------------------------------",self.ini,"-----------------------------------")
475 476 #numpy.set_printoptions(suppress=True)
476 477 #print(self.data.times)
477 478 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
478 479 data = self.data[-1]
479 480 # ALTURA altura_tmp_h
480 481 altura_h = (data['weather'].shape[1])/10.0
481 482 stoprange = float(altura_h*1.5)#stoprange = float(33*1.5) por ahora 400
482 483 rangestep = float(0.15)
483 484 r = numpy.arange(0, stoprange, rangestep)
484 485 self.y = 2*r
485 486 # RADAR
486 487 #data_weather = data['weather']
487 488 # PEDESTAL
488 489 #data_azi = data['azi']
489 490 res = 1
490 491 # STEP
491 492 step = (360/(res*data['weather'].shape[0]))
492 493 #print("shape wr_data", wr_data.shape)
493 494 #print("shape wr_azi",wr_azi.shape)
494 495 #print("step",step)
495 496 print("Time---->",self.data.times[-1],thisDatetime)
496 497 #print("alturas", len(self.y))
497 498 self.res_weather, self.res_azi = self.const_ploteo(data_weather=data['weather'],data_azi=data['azi'],step=step,res=res)
498 499 #numpy.set_printoptions(suppress=True)
499 500 #print("resultado",self.res_azi)
500 501 ##########################################################
501 502 ################# PLOTEO ###################
502 503 ##########################################################
503 504
504 505 for i,ax in enumerate(self.axes):
505 506 if ax.firsttime:
506 507 plt.clf()
507 508 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=1, vmax=60)
508 509 else:
509 510 plt.clf()
510 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=1, vmax=60)
511 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=0, vmax=60)
511 512 caax = cgax.parasites[0]
512 513 paax = cgax.parasites[1]
513 514 cbar = plt.gcf().colorbar(pm, pad=0.075)
514 515 caax.set_xlabel('x_range [km]')
515 516 caax.set_ylabel('y_range [km]')
516 517 plt.text(1.0, 1.05, 'azimuth '+str(thisDatetime)+"step"+str(self.ini), transform=caax.transAxes, va='bottom',ha='right')
517 518
518 519 self.ini= self.ini+1
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,75 +1,81
1 1 import os,sys
2 2 import datetime
3 3 import time
4 4 from schainpy.controller import Project
5 5 #path='/DATA_RM/TEST_HDF5/d2021200'
6 6 #path='/DATA_RM/TEST_HDF5/d2021200'
7 7 #path='/DATA_RM/TEST_HDF5/d2021214'
8 8 #path='/DATA_RM/TEST_HDF5/d2021229'
9 9
10 10 #path='/DATA_RM/TEST_HDF5/d2021231'
11 11 #path='/DATA_RM/TEST_HDF5/ADQ_OFFLINE/d2021231'
12 12 path='/DATA_RM/TEST_HDF5/d2021231'
13 13 #path='/DATA_RM/TEST_14_HDF5/d2021257'
14 14 ## TEST ULTIMA PRUEBA 22 DE SEPTIEMBRE
15 15 path = '/DATA_RM/TEST_HDF5_PP_22/d2021265'
16 #path = '/DATA_RM/TEST_HDF5_PP_100/d2021285'
17 path = '/DATA_RM/TEST_HDF5_PP/d2021285'
18
19
16 20 path_adq=path
17 21 #path_ped='/DATA_RM/TEST_PEDESTAL/P2021200'
18 22 #path_ped='/DATA_RM/TEST_PEDESTAL/P2021214'
19 23 #path_ped='/DATA_RM/TEST_PEDESTAL/P2021230'
20 24 #path_ped='/DATA_RM/TEST_PEDESTAL/P20210819'
21 25 #path_ped='/DATA_RM/TEST_PEDESTAL/P20210819-154315'
22 26 #path_ped='/DATA_RM/TEST_PEDESTAL/P20210914-162434'
23 27 #path_ped='/DATA_RM/TEST_PEDESTAL/PEDESTAL_OFFLINE/P20210819-161524'
24 28 #pruebas con perdida de datos
25 29 #path_ped='/DATA_RM/TEST_PEDESTAL/PEDESTAL_OFFLINE/P20210819-161524_TEST'
26 30 ## TEST ULTIMA PRUEBA 22 DE SEPTIEMBRE
27 path_ped='/DATA_RM/TEST_PEDESTAL/P20210922-122731'
31 path_ped='/DATA_RM/TEST_PEDESTAL/P20211012-082745'
28 32
29 33
30 34 figpath = '/home/soporte/Pictures'
31 35 desc = "Simulator Test"
32 36
33 37 controllerObj = Project()
34 38 controllerObj.setup(id='10',name='Test Simulator',description=desc)
35 39 readUnitConfObj = controllerObj.addReadUnit(datatype='HDFReader',
36 40 path=path,
37 41 startDate="2021/01/01", #"2020/01/01",#today,
38 42 endDate= "2021/12/01", #"2020/12/30",#today,
39 43 startTime='00:00:00',
40 44 endTime='23:59:59',
41 45 t_Interval_p=0.01,
42 46 n_Muestras_p=100,
43 47 delay=30,
44 48 #set=0,
45 49 online=0,
46 50 walk=0,
47 51 nTries=6)#1
48 52
49 53 procUnitConfObjA = controllerObj.addProcUnit(datatype='ParametersProc',inputId=readUnitConfObj.getId())
50 V=2
51 blocksPerfile=100
54 V=10
55 blocksPerfile=360
52 56 print("Velocidad del Pedestal",V)
53 f_a_p= int(blocksPerfile/V)
57 tmp_blocksPerfile=100
58 f_a_p= int(tmp_blocksPerfile/V)
54 59
55 60 opObj11 = procUnitConfObjA.addOperation(name='PedestalInformation')
56 61 opObj11.addParameter(name='path_ped', value=path_ped)
57 opObj11.addParameter(name='path_adq', value=path_adq)
62 #opObj11.addParameter(name='path_adq', value=path_adq)
58 63 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
64 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
59 65 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
60 opObj11.addParameter(name='blocksPerfile', value='100', format='int')
61 66 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
62 67 opObj11.addParameter(name='online', value='0', format='int')# habilitar el enable aqui tambien
63 68
64 69
65 70 opObj11 = procUnitConfObjA.addOperation(name='Block360')
66 71 opObj11.addParameter(name='n', value='10', format='int')
72 opObj11.addParameter(name='mode', value=0, format='int')
67 73 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
68 74
69 75 opObj11= procUnitConfObjA.addOperation(name='WeatherPlot',optype='other')
70 76 #opObj11.addParameter(name='save', value=figpath)
71 77 #opObj11.addParameter(name='save_period', value=1)
72 78
73 79 controllerObj.start()
74 80 #online 1 utc_adq 1617490240.48
75 81 #online 0 utc_adq 1617489815.4804
General Comments 0
You need to be logged in to leave comments. Login now