##// END OF EJS Templates
LAST CHANGE WR
avaldez -
r1393:ccbe99dfcead
parent child
Show More

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

@@ -0,0 +1,217
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 = '/DATA_RM/WR_20_OCT'
30
31 #path_ped='/DATA_RM/TEST_PEDESTAL/P20211012-082745'
32 path_ped='/DATA_RM/TEST_PEDESTAL/P20211020-131248'
33
34 figpath_pp = "/home/soporte/Pictures/TEST_PP"
35 figpath_mom = "/home/soporte/Pictures/TEST_MOM"
36 plot = 0
37 integration = 1
38 save = 0
39 if save == 1:
40 if mode==0:
41 path_save = '/DATA_RM/TEST_HDF5_PP_23/6v'
42 path_save = '/DATA_RM/TEST_HDF5_PP'
43 path_save = '/DATA_RM/TEST_HDF5_PP_100'
44 else:
45 path_save = '/DATA_RM/TEST_HDF5_SPEC_23_V2/6v'
46
47 print("* PATH data ADQ :", path)
48 print("* Velocidad Pedestal :",V,"Β°/seg")
49 ############################ NRO Perfiles PROCESAMIENTO ###################
50 V=V
51 IPP=400*1e-6
52 n= int(1/(V*IPP))
53 print("* n - NRO Perfiles Proc:", n )
54 ################################## MODE ###################################
55 print("* Modo de Operacion :",mode)
56 if mode ==0:
57 print("* Met. Seleccionado : Pulse Pair")
58 else:
59 print("* Met. Momentos : Momentos")
60
61 ################################## MODE ###################################
62 print("* Grabado de datos :",save)
63 if save ==1:
64 if mode==0:
65 ope= "Pulse Pair"
66 else:
67 ope= "Momentos"
68 print("* Path-Save Data -", ope , path_save)
69
70 print("* Integracion de datos :",integration)
71
72 time.sleep(15)
73 #remotefolder = "/home/wmaster/graficos"
74 #######################################################################
75 ################# RANGO DE PLOTEO######################################
76 dBmin = '1'
77 dBmax = '85'
78 xmin = '15'
79 xmax = '15.25'
80 ymin = '0'
81 ymax = '600'
82 #######################################################################
83 ########################FECHA##########################################
84 str = datetime.date.today()
85 today = str.strftime("%Y/%m/%d")
86 str2 = str - datetime.timedelta(days=1)
87 yesterday = str2.strftime("%Y/%m/%d")
88 #######################################################################
89 ########################SIGNAL CHAIN ##################################
90 #######################################################################
91 desc = "USRP_test"
92 filename = "USRP_processing.xml"
93 controllerObj = Project()
94 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
95 #######################################################################
96 ######################## UNIDAD DE LECTURA#############################
97 #######################################################################
98 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
99 path=path,
100 startDate="2021/01/01",#today,
101 endDate="2021/12/30",#today,
102 startTime='00:00:00',
103 endTime='23:59:59',
104 delay=0,
105 #set=0,
106 online=0,
107 walk=1,
108 ippKm = 60)
109
110 opObj11 = readUnitConfObj.addOperation(name='printInfo')
111
112 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
113
114 if mode ==0:
115 ####################### METODO PULSE PAIR ######################################################################
116 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
117 opObj11.addParameter(name='n', value=int(n), format='int')#10 VOY A USAR 250 DADO QUE LA VELOCIDAD ES 10 GRADOS
118 #opObj11.addParameter(name='removeDC', value=1, format='int')
119 ####################### METODO Parametros ######################################################################
120 procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
121 if plot==1:
122 opObj11 = procUnitConfObjB.addOperation(name='GenericRTIPlot',optype='external')
123 opObj11.addParameter(name='attr_data', value='dataPP_POW')
124 opObj11.addParameter(name='colormap', value='jet')
125 opObj11.addParameter(name='xmin', value=xmin)
126 opObj11.addParameter(name='xmax', value=xmax)
127 opObj11.addParameter(name='zmin', value=dBmin)
128 opObj11.addParameter(name='zmax', value=dBmax)
129 opObj11.addParameter(name='save', value=figpath_pp)
130 opObj11.addParameter(name='showprofile', value=0)
131 opObj11.addParameter(name='save_period', value=50)
132
133 ####################### METODO ESCRITURA #######################################################################
134 if save==1:
135 opObj10 = procUnitConfObjB.addOperation(name='HDFWriter')
136 opObj10.addParameter(name='path',value=path_save)
137 #opObj10.addParameter(name='mode',value=0)
138 opObj10.addParameter(name='blocksPerFile',value='100',format='int')
139 opObj10.addParameter(name='metadataList',value='utctimeInit,timeZone,paramInterval,profileIndex,channelList,heightList,flagDataAsBlock',format='list')
140 opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,utctime',format='list')#,format='list'
141 if integration==1:
142 V=10
143 blocksPerfile=360
144 print("* Velocidad del Pedestal:",V)
145 tmp_blocksPerfile = 100
146 f_a_p= int(tmp_blocksPerfile/V)
147
148 opObj11 = procUnitConfObjB.addOperation(name='PedestalInformation')
149 opObj11.addParameter(name='path_ped', value=path_ped)
150 #opObj11.addParameter(name='path_adq', value=path_adq)
151 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
152 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
153 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
154 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
155 opObj11.addParameter(name='online', value='0', format='int')
156
157 opObj11 = procUnitConfObjB.addOperation(name='Block360')
158 opObj11.addParameter(name='n', value='10', format='int')
159 opObj11.addParameter(name='mode', value=mode, format='int')
160
161 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
162
163 opObj11= procUnitConfObjB.addOperation(name='WeatherPlot',optype='other')
164
165
166 else:
167 ####################### METODO SPECTROS ######################################################################
168 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
169 procUnitConfObjB.addParameter(name='nFFTPoints', value=n, format='int')
170 procUnitConfObjB.addParameter(name='nProfiles' , value=n, format='int')
171
172 procUnitConfObjC = controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjB.getId())
173 procUnitConfObjC.addOperation(name='SpectralMoments')
174 if plot==1:
175 dBmin = '1'
176 dBmax = '65'
177 opObj11 = procUnitConfObjC.addOperation(name='PowerPlot',optype='external')
178 opObj11.addParameter(name='xmin', value=xmin)
179 opObj11.addParameter(name='xmax', value=xmax)
180 opObj11.addParameter(name='zmin', value=dBmin)
181 opObj11.addParameter(name='zmax', value=dBmax)
182 opObj11.addParameter(name='save', value=figpath_mom)
183 opObj11.addParameter(name='showprofile', value=0)
184 opObj11.addParameter(name='save_period', value=100)
185
186 if save==1:
187 opObj10 = procUnitConfObjC.addOperation(name='HDFWriter')
188 opObj10.addParameter(name='path',value=path_save)
189 #opObj10.addParameter(name='mode',value=0)
190 opObj10.addParameter(name='blocksPerFile',value='360',format='int')
191 #opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
192 opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
193 opObj10.addParameter(name='dataList',value='data_pow,data_dop,utctime',format='list')#,format='list'
194
195 if integration==1:
196 V=10
197 blocksPerfile=360
198 print("* Velocidad del Pedestal:",V)
199 tmp_blocksPerfile = 100
200 f_a_p= int(tmp_blocksPerfile/V)
201
202 opObj11 = procUnitConfObjC.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='blocksPerfile', value=blocksPerfile, format='int')
207 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
208 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
209 opObj11.addParameter(name='online', value='0', format='int')
210
211 opObj11 = procUnitConfObjC.addOperation(name='Block360')
212 opObj11.addParameter(name='n', value='30', format='int')
213 opObj11.addParameter(name='mode', value=mode, format='int')
214
215 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
216 opObj11= procUnitConfObjC.addOperation(name='WeatherPlot',optype='other')
217 controllerObj.start()
@@ -0,0 +1,13
1
2 PULSE PAIR
3 help [[1.364 1.376 1.28 ... 1.352 1.372 1.332]
4 [2.912 3.012 3.06 ... 3.056 2.6 2.936]]
5 help [[1.256 1.304 1.212 ... 1.228 1.328 1.528]
6 [2.744 2.604 2.492 ... 2.544 2.916 2.644]]
7 help [[1.176 1.248 1.48 ... 1.388 1.396 1.216]
8 [2.564 2.524 2.756 ... 2.772 2.7 2.684]]
9
10 help [[1.312 1.328 1.28 ... 1.28 1.444 1.532]
11 [3. 2.964 3.092 ... 3.104 3.016 2.984]]
12
13 MOMENTOS
@@ -0,0 +1,69
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
9 print("SETUP- RADAR METEOROLOGICO")
10 V = 10
11 #######################################################################
12 ################# RANGO DE PLOTEO######################################
13 dBmin = '1'
14 dBmax = '65'
15 xmin = '13.2'
16 xmax = '13.5'
17 ymin = '0'
18 ymax = '60'
19
20 path = '/DATA_RM/WR_20_OCT'
21 figpath_pp = "/home/soporte/Pictures/TEST_PP"
22
23
24 IPP=400*1e-6
25 n= int(1/(V*IPP))
26 #n=250
27 print("* n - NRO Perfiles Proc:", n )
28
29 desc = "USRP_test"
30 filename = "USRP_processing.xml"
31 controllerObj = Project()
32 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
33 #######################################################################
34 ######################## UNIDAD DE LECTURA#############################
35 #######################################################################
36 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
37 path=path,
38 startDate="2021/01/01",#today,
39 endDate="2021/12/30",#today,
40 startTime='00:00:00',
41 endTime='23:59:59',
42 delay=0,
43 #set=0,
44 online=0,
45 walk=1,
46 ippKm = 60)
47
48 opObj11 = readUnitConfObj.addOperation(name='printInfo')
49
50 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
51
52 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
53 opObj11.addParameter(name='n', value=int(n), format='int')#10 VOY A USAR 250 DADO QUE LA VELOCIDAD ES 10 GRADOS
54 #opObj11.addParameter(name='removeDC', value=1, format='int')
55
56 procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
57
58 opObj11 = procUnitConfObjB.addOperation(name='GenericRTIPlot',optype='external')
59 opObj11.addParameter(name='attr_data', value='dataPP_POWER')
60 opObj11.addParameter(name='colormap', value='jet')
61 opObj11.addParameter(name='xmin', value=xmin)
62 opObj11.addParameter(name='xmax', value=xmax)
63 opObj11.addParameter(name='zmin', value=dBmin)
64 opObj11.addParameter(name='zmax', value=dBmax)
65 opObj11.addParameter(name='save', value=figpath_pp)
66 opObj11.addParameter(name='showprofile', value=0)
67 #opObj11.addParameter(name='save_period', value=10)
68
69 controllerObj.start()
@@ -0,0 +1,80
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
9 print("SETUP- RADAR METEOROLOGICO")
10 V = 10
11 #######################################################################
12 ################# RANGO DE PLOTEO######################################
13 dBmin = '1'
14 dBmax = '65'
15 xmin = '13.2'
16 xmax = '13.5'
17 ymin = '0'
18 ymax = '60'
19
20 path = '/DATA_RM/WR_20_OCT'
21 figpath_spec = "/home/soporte/Pictures/TEST_MOM"
22
23
24 IPP=400*1e-6
25 n= int(1/(V*IPP))
26 print("* n - NRO Perfiles Proc:", n )
27 time.sleep(5)
28 desc = "USRP_test"
29 filename = "USRP_processing.xml"
30 controllerObj = Project()
31 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
32 #######################################################################
33 ######################## UNIDAD DE LECTURA#############################
34 #######################################################################
35 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
36 path=path,
37 startDate="2021/01/01",#today,
38 endDate="2021/12/30",#today,
39 startTime='00:00:00',
40 endTime='23:59:59',
41 delay=0,
42 #set=0,
43 online=0,
44 walk=1,
45 ippKm = 60)
46
47 opObj11 = readUnitConfObj.addOperation(name='printInfo')
48
49 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
50
51 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
52 procUnitConfObjB.addParameter(name='nFFTPoints', value=n, format='int')
53 procUnitConfObjB.addParameter(name='nProfiles' , value=n, format='int')
54 '''
55 opObj11 = procUnitConfObjB.addOperation(name='RTIPlot', optype='external')
56 #.addParameter(name='id', value='2', format='int')
57 opObj11.addParameter(name='wintitle', value='RTIPlot', format='str')
58 opObj11.addParameter(name='xmin', value=xmin)
59 opObj11.addParameter(name='xmax', value=xmax)
60 opObj11.addParameter(name='zmin', value=dBmin, format='int')
61 opObj11.addParameter(name='zmax', value=dBmax, format='int')
62 '''
63 #opObj13 = procUnitConfObjB.addOperation(name='removeDC')
64 #opObj13.addParameter(name='mode', value='2', format='int')
65
66 procUnitConfObjC = controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjB.getId())
67 procUnitConfObjC.addOperation(name='SpectralMoments')
68
69 dBmin = '1'
70 dBmax = '65'
71 opObj11 = procUnitConfObjC.addOperation(name='PowerPlot',optype='external')
72 opObj11.addParameter(name='xmin', value=xmin)
73 opObj11.addParameter(name='xmax', value=xmax)
74 opObj11.addParameter(name='zmin', value=dBmin)
75 opObj11.addParameter(name='zmax', value=dBmax)
76 opObj11.addParameter(name='save', value=figpath_spec)
77 opObj11.addParameter(name='showprofile', value=0)
78 #opObj11.addParameter(name='save_period', value=10)
79
80 controllerObj.start()
@@ -1,1069 +1,1068
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 """Definition of diferent Data objects for different types of data
6 6
7 7 Here you will find the diferent data objects for the different types
8 8 of data, this data objects must be used as dataIn or dataOut objects in
9 9 processing units and operations. Currently the supported data objects are:
10 10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 11 """
12 12
13 13 import copy
14 14 import numpy
15 15 import datetime
16 16 import json
17 17
18 18 import schainpy.admin
19 19 from schainpy.utils import log
20 20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 21 from schainpy.model.data import _noise
22 22
23 23
24 24 def getNumpyDtype(dataTypeCode):
25 25
26 26 if dataTypeCode == 0:
27 27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 28 elif dataTypeCode == 1:
29 29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 30 elif dataTypeCode == 2:
31 31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 32 elif dataTypeCode == 3:
33 33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 34 elif dataTypeCode == 4:
35 35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 36 elif dataTypeCode == 5:
37 37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 38 else:
39 39 raise ValueError('dataTypeCode was not defined')
40 40
41 41 return numpyDtype
42 42
43 43
44 44 def getDataTypeCode(numpyDtype):
45 45
46 46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 47 datatype = 0
48 48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 49 datatype = 1
50 50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 51 datatype = 2
52 52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 53 datatype = 3
54 54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 55 datatype = 4
56 56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 57 datatype = 5
58 58 else:
59 59 datatype = None
60 60
61 61 return datatype
62 62
63 63
64 64 def hildebrand_sekhon(data, navg):
65 65 """
66 66 This method is for the objective determination of the noise level in Doppler spectra. This
67 67 implementation technique is based on the fact that the standard deviation of the spectral
68 68 densities is equal to the mean spectral density for white Gaussian noise
69 69
70 70 Inputs:
71 71 Data : heights
72 72 navg : numbers of averages
73 73
74 74 Return:
75 75 mean : noise's level
76 76 """
77 77
78 78 sortdata = numpy.sort(data, axis=None)
79 79 '''
80 80 lenOfData = len(sortdata)
81 81 nums_min = lenOfData*0.2
82 82
83 83 if nums_min <= 5:
84 84
85 85 nums_min = 5
86 86
87 87 sump = 0.
88 88 sumq = 0.
89 89
90 90 j = 0
91 91 cont = 1
92 92
93 93 while((cont == 1)and(j < lenOfData)):
94 94
95 95 sump += sortdata[j]
96 96 sumq += sortdata[j]**2
97 97
98 98 if j > nums_min:
99 99 rtest = float(j)/(j-1) + 1.0/navg
100 100 if ((sumq*j) > (rtest*sump**2)):
101 101 j = j - 1
102 102 sump = sump - sortdata[j]
103 103 sumq = sumq - sortdata[j]**2
104 104 cont = 0
105 105
106 106 j += 1
107 107
108 108 lnoise = sump / j
109 109 '''
110 110 return _noise.hildebrand_sekhon(sortdata, navg)
111 111
112 112
113 113 class Beam:
114 114
115 115 def __init__(self):
116 116 self.codeList = []
117 117 self.azimuthList = []
118 118 self.zenithList = []
119 119
120 120
121 121 class GenericData(object):
122 122
123 123 flagNoData = True
124 124
125 125 def copy(self, inputObj=None):
126 126
127 127 if inputObj == None:
128 128 return copy.deepcopy(self)
129 129
130 130 for key in list(inputObj.__dict__.keys()):
131 131
132 132 attribute = inputObj.__dict__[key]
133 133
134 134 # If this attribute is a tuple or list
135 135 if type(inputObj.__dict__[key]) in (tuple, list):
136 136 self.__dict__[key] = attribute[:]
137 137 continue
138 138
139 139 # If this attribute is another object or instance
140 140 if hasattr(attribute, '__dict__'):
141 141 self.__dict__[key] = attribute.copy()
142 142 continue
143 143
144 144 self.__dict__[key] = inputObj.__dict__[key]
145 145
146 146 def deepcopy(self):
147 147
148 148 return copy.deepcopy(self)
149 149
150 150 def isEmpty(self):
151 151
152 152 return self.flagNoData
153 153
154 154 def isReady(self):
155 155
156 156 return not self.flagNoData
157 157
158 158
159 159 class JROData(GenericData):
160 160
161 161 systemHeaderObj = SystemHeader()
162 162 radarControllerHeaderObj = RadarControllerHeader()
163 163 type = None
164 164 datatype = None # dtype but in string
165 165 nProfiles = None
166 166 heightList = None
167 167 channelList = None
168 168 flagDiscontinuousBlock = False
169 169 useLocalTime = False
170 170 utctime = None
171 171 timeZone = None
172 172 dstFlag = None
173 173 errorCount = None
174 174 blocksize = None
175 175 flagDecodeData = False # asumo q la data no esta decodificada
176 176 flagDeflipData = False # asumo q la data no esta sin flip
177 177 flagShiftFFT = False
178 178 nCohInt = None
179 179 windowOfFilter = 1
180 180 C = 3e8
181 181 frequency = 49.92e6
182 182 realtime = False
183 183 beacon_heiIndexList = None
184 184 last_block = None
185 185 blocknow = None
186 186 azimuth = None
187 187 zenith = None
188 188 beam = Beam()
189 189 profileIndex = None
190 190 error = None
191 191 data = None
192 192 nmodes = None
193 193 metadata_list = ['heightList', 'timeZone', 'type']
194 194
195 195 def __str__(self):
196 196
197 197 return '{} - {}'.format(self.type, self.datatime())
198 198
199 199 def getNoise(self):
200 200
201 201 raise NotImplementedError
202 202
203 203 @property
204 204 def nChannels(self):
205 205
206 206 return len(self.channelList)
207 207
208 208 @property
209 209 def channelIndexList(self):
210 210
211 211 return list(range(self.nChannels))
212 212
213 213 @property
214 214 def nHeights(self):
215 215
216 216 return len(self.heightList)
217 217
218 218 def getDeltaH(self):
219 219
220 220 return self.heightList[1] - self.heightList[0]
221 221
222 222 @property
223 223 def ltctime(self):
224 224
225 225 if self.useLocalTime:
226 226 return self.utctime - self.timeZone * 60
227 227
228 228 return self.utctime
229 229
230 230 @property
231 231 def datatime(self):
232 232
233 233 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
234 234 return datatimeValue
235 235
236 236 def getTimeRange(self):
237 237
238 238 datatime = []
239 239
240 240 datatime.append(self.ltctime)
241 241 datatime.append(self.ltctime + self.timeInterval + 1)
242 242
243 243 datatime = numpy.array(datatime)
244 244
245 245 return datatime
246 246
247 247 def getFmaxTimeResponse(self):
248 248
249 249 period = (10**-6) * self.getDeltaH() / (0.15)
250 250
251 251 PRF = 1. / (period * self.nCohInt)
252 252
253 253 fmax = PRF
254 254
255 255 return fmax
256 256
257 257 def getFmax(self):
258 258 PRF = 1. / (self.ippSeconds * self.nCohInt)
259 259
260 260 fmax = PRF
261 261 return fmax
262 262
263 263 def getVmax(self):
264 264
265 265 _lambda = self.C / self.frequency
266 266
267 267 vmax = self.getFmax() * _lambda / 2
268 268
269 269 return vmax
270 270
271 271 @property
272 272 def ippSeconds(self):
273 273 '''
274 274 '''
275 275 return self.radarControllerHeaderObj.ippSeconds
276 276
277 277 @ippSeconds.setter
278 278 def ippSeconds(self, ippSeconds):
279 279 '''
280 280 '''
281 281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
282 282
283 283 @property
284 284 def code(self):
285 285 '''
286 286 '''
287 287 return self.radarControllerHeaderObj.code
288 288
289 289 @code.setter
290 290 def code(self, code):
291 291 '''
292 292 '''
293 293 self.radarControllerHeaderObj.code = code
294 294
295 295 @property
296 296 def nCode(self):
297 297 '''
298 298 '''
299 299 return self.radarControllerHeaderObj.nCode
300 300
301 301 @nCode.setter
302 302 def nCode(self, ncode):
303 303 '''
304 304 '''
305 305 self.radarControllerHeaderObj.nCode = ncode
306 306
307 307 @property
308 308 def nBaud(self):
309 309 '''
310 310 '''
311 311 return self.radarControllerHeaderObj.nBaud
312 312
313 313 @nBaud.setter
314 314 def nBaud(self, nbaud):
315 315 '''
316 316 '''
317 317 self.radarControllerHeaderObj.nBaud = nbaud
318 318
319 319 @property
320 320 def ipp(self):
321 321 '''
322 322 '''
323 323 return self.radarControllerHeaderObj.ipp
324 324
325 325 @ipp.setter
326 326 def ipp(self, ipp):
327 327 '''
328 328 '''
329 329 self.radarControllerHeaderObj.ipp = ipp
330 330
331 331 @property
332 332 def metadata(self):
333 333 '''
334 334 '''
335 335
336 336 return {attr: getattr(self, attr) for attr in self.metadata_list}
337 337
338 338
339 339 class Voltage(JROData):
340 340
341 341 dataPP_POW = None
342 342 dataPP_DOP = None
343 343 dataPP_WIDTH = None
344 344 dataPP_SNR = None
345 345
346 346 def __init__(self):
347 347 '''
348 348 Constructor
349 349 '''
350 350
351 351 self.useLocalTime = True
352 352 self.radarControllerHeaderObj = RadarControllerHeader()
353 353 self.systemHeaderObj = SystemHeader()
354 354 self.type = "Voltage"
355 355 self.data = None
356 356 self.nProfiles = None
357 357 self.heightList = None
358 358 self.channelList = None
359 359 self.flagNoData = True
360 360 self.flagDiscontinuousBlock = False
361 361 self.utctime = None
362 362 self.timeZone = 0
363 363 self.dstFlag = None
364 364 self.errorCount = None
365 365 self.nCohInt = None
366 366 self.blocksize = None
367 367 self.flagCohInt = False
368 368 self.flagDecodeData = False # asumo q la data no esta decodificada
369 369 self.flagDeflipData = False # asumo q la data no esta sin flip
370 370 self.flagShiftFFT = False
371 371 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
372 372 self.profileIndex = 0
373 373 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
374 374 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
375 375
376 376 def getNoisebyHildebrand(self, channel=None):
377 377 """
378 378 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
379 379
380 380 Return:
381 381 noiselevel
382 382 """
383 383
384 384 if channel != None:
385 385 data = self.data[channel]
386 386 nChannels = 1
387 387 else:
388 388 data = self.data
389 389 nChannels = self.nChannels
390 390
391 391 noise = numpy.zeros(nChannels)
392 392 power = data * numpy.conjugate(data)
393 393
394 394 for thisChannel in range(nChannels):
395 395 if nChannels == 1:
396 396 daux = power[:].real
397 397 else:
398 398 daux = power[thisChannel, :].real
399 399 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
400 400
401 401 return noise
402 402
403 403 def getNoise(self, type=1, channel=None):
404 404
405 405 if type == 1:
406 406 noise = self.getNoisebyHildebrand(channel)
407 407
408 408 return noise
409 409
410 410 def getPower(self, channel=None):
411 411
412 412 if channel != None:
413 413 data = self.data[channel]
414 414 else:
415 415 data = self.data
416 416
417 417 power = data * numpy.conjugate(data)
418 418 powerdB = 10 * numpy.log10(power.real)
419 419 powerdB = numpy.squeeze(powerdB)
420 420
421 421 return powerdB
422 422
423 423 @property
424 424 def timeInterval(self):
425 425
426 426 return self.ippSeconds * self.nCohInt
427 427
428 428 noise = property(getNoise, "I'm the 'nHeights' property.")
429 429
430 430
431 431 class Spectra(JROData):
432 432
433 433 def __init__(self):
434 434 '''
435 435 Constructor
436 436 '''
437 437
438 438 self.data_dc = None
439 439 self.data_spc = None
440 440 self.data_cspc = None
441 441 self.useLocalTime = True
442 442 self.radarControllerHeaderObj = RadarControllerHeader()
443 443 self.systemHeaderObj = SystemHeader()
444 444 self.type = "Spectra"
445 445 self.timeZone = 0
446 446 self.nProfiles = None
447 447 self.heightList = None
448 448 self.channelList = None
449 449 self.pairsList = None
450 450 self.flagNoData = True
451 451 self.flagDiscontinuousBlock = False
452 452 self.utctime = None
453 453 self.nCohInt = None
454 454 self.nIncohInt = None
455 455 self.blocksize = None
456 456 self.nFFTPoints = None
457 457 self.wavelength = None
458 458 self.flagDecodeData = False # asumo q la data no esta decodificada
459 459 self.flagDeflipData = False # asumo q la data no esta sin flip
460 460 self.flagShiftFFT = False
461 461 self.ippFactor = 1
462 462 self.beacon_heiIndexList = []
463 463 self.noise_estimation = None
464 464 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
465 465 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
466 466
467 467 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
468 468 """
469 469 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
470 470
471 471 Return:
472 472 noiselevel
473 473 """
474 474
475 475 noise = numpy.zeros(self.nChannels)
476 476
477 477 for channel in range(self.nChannels):
478 478 daux = self.data_spc[channel,
479 479 xmin_index:xmax_index, ymin_index:ymax_index]
480 480 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
481 481
482 482 return noise
483 483
484 484 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
485 485
486 486 if self.noise_estimation is not None:
487 487 # this was estimated by getNoise Operation defined in jroproc_spectra.py
488 488 return self.noise_estimation
489 489 else:
490 490 noise = self.getNoisebyHildebrand(
491 491 xmin_index, xmax_index, ymin_index, ymax_index)
492 492 return noise
493 493
494 494 def getFreqRangeTimeResponse(self, extrapoints=0):
495 495
496 496 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
497 497 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
498 498
499 499 return freqrange
500 500
501 501 def getAcfRange(self, extrapoints=0):
502 502
503 503 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
504 504 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
505 505
506 506 return freqrange
507 507
508 508 def getFreqRange(self, extrapoints=0):
509 509
510 510 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
511 511 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
512 512
513 513 return freqrange
514 514
515 515 def getVelRange(self, extrapoints=0):
516 516
517 517 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
518 518 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
519 519
520 520 if self.nmodes:
521 521 return velrange/self.nmodes
522 522 else:
523 523 return velrange
524 524
525 525 @property
526 526 def nPairs(self):
527 527
528 528 return len(self.pairsList)
529 529
530 530 @property
531 531 def pairsIndexList(self):
532 532
533 533 return list(range(self.nPairs))
534 534
535 535 @property
536 536 def normFactor(self):
537 537
538 538 pwcode = 1
539 539
540 540 if self.flagDecodeData:
541 541 pwcode = numpy.sum(self.code[0]**2)
542 542 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
543 543 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
544 544
545 545 return normFactor
546 546
547 547 @property
548 548 def flag_cspc(self):
549 549
550 550 if self.data_cspc is None:
551 551 return True
552 552
553 553 return False
554 554
555 555 @property
556 556 def flag_dc(self):
557 557
558 558 if self.data_dc is None:
559 559 return True
560 560
561 561 return False
562 562
563 563 @property
564 564 def timeInterval(self):
565 565
566 566 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
567 567 if self.nmodes:
568 568 return self.nmodes*timeInterval
569 569 else:
570 570 return timeInterval
571 571
572 572 def getPower(self):
573 573
574 574 factor = self.normFactor
575 575 z = self.data_spc / factor
576 576 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
577 577 avg = numpy.average(z, axis=1)
578
579 578 return 10 * numpy.log10(avg)
580 579
581 580 def getCoherence(self, pairsList=None, phase=False):
582 581
583 582 z = []
584 583 if pairsList is None:
585 584 pairsIndexList = self.pairsIndexList
586 585 else:
587 586 pairsIndexList = []
588 587 for pair in pairsList:
589 588 if pair not in self.pairsList:
590 589 raise ValueError("Pair %s is not in dataOut.pairsList" % (
591 590 pair))
592 591 pairsIndexList.append(self.pairsList.index(pair))
593 592 for i in range(len(pairsIndexList)):
594 593 pair = self.pairsList[pairsIndexList[i]]
595 594 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
596 595 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
597 596 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
598 597 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
599 598 if phase:
600 599 data = numpy.arctan2(avgcoherenceComplex.imag,
601 600 avgcoherenceComplex.real) * 180 / numpy.pi
602 601 else:
603 602 data = numpy.abs(avgcoherenceComplex)
604 603
605 604 z.append(data)
606 605
607 606 return numpy.array(z)
608 607
609 608 def setValue(self, value):
610 609
611 610 print("This property should not be initialized")
612 611
613 612 return
614 613
615 614 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
616 615
617 616
618 617 class SpectraHeis(Spectra):
619 618
620 619 def __init__(self):
621 620
622 621 self.radarControllerHeaderObj = RadarControllerHeader()
623 622 self.systemHeaderObj = SystemHeader()
624 623 self.type = "SpectraHeis"
625 624 self.nProfiles = None
626 625 self.heightList = None
627 626 self.channelList = None
628 627 self.flagNoData = True
629 628 self.flagDiscontinuousBlock = False
630 629 self.utctime = None
631 630 self.blocksize = None
632 631 self.profileIndex = 0
633 632 self.nCohInt = 1
634 633 self.nIncohInt = 1
635 634
636 635 @property
637 636 def normFactor(self):
638 637 pwcode = 1
639 638 if self.flagDecodeData:
640 639 pwcode = numpy.sum(self.code[0]**2)
641 640
642 641 normFactor = self.nIncohInt * self.nCohInt * pwcode
643 642
644 643 return normFactor
645 644
646 645 @property
647 646 def timeInterval(self):
648 647
649 648 return self.ippSeconds * self.nCohInt * self.nIncohInt
650 649
651 650
652 651 class Fits(JROData):
653 652
654 653 def __init__(self):
655 654
656 655 self.type = "Fits"
657 656 self.nProfiles = None
658 657 self.heightList = None
659 658 self.channelList = None
660 659 self.flagNoData = True
661 660 self.utctime = None
662 661 self.nCohInt = 1
663 662 self.nIncohInt = 1
664 663 self.useLocalTime = True
665 664 self.profileIndex = 0
666 665 self.timeZone = 0
667 666
668 667 def getTimeRange(self):
669 668
670 669 datatime = []
671 670
672 671 datatime.append(self.ltctime)
673 672 datatime.append(self.ltctime + self.timeInterval)
674 673
675 674 datatime = numpy.array(datatime)
676 675
677 676 return datatime
678 677
679 678 def getChannelIndexList(self):
680 679
681 680 return list(range(self.nChannels))
682 681
683 682 def getNoise(self, type=1):
684 683
685 684
686 685 if type == 1:
687 686 noise = self.getNoisebyHildebrand()
688 687
689 688 if type == 2:
690 689 noise = self.getNoisebySort()
691 690
692 691 if type == 3:
693 692 noise = self.getNoisebyWindow()
694 693
695 694 return noise
696 695
697 696 @property
698 697 def timeInterval(self):
699 698
700 699 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
701 700
702 701 return timeInterval
703 702
704 703 @property
705 704 def ippSeconds(self):
706 705 '''
707 706 '''
708 707 return self.ipp_sec
709 708
710 709 noise = property(getNoise, "I'm the 'nHeights' property.")
711 710
712 711
713 712 class Correlation(JROData):
714 713
715 714 def __init__(self):
716 715 '''
717 716 Constructor
718 717 '''
719 718 self.radarControllerHeaderObj = RadarControllerHeader()
720 719 self.systemHeaderObj = SystemHeader()
721 720 self.type = "Correlation"
722 721 self.data = None
723 722 self.dtype = None
724 723 self.nProfiles = None
725 724 self.heightList = None
726 725 self.channelList = None
727 726 self.flagNoData = True
728 727 self.flagDiscontinuousBlock = False
729 728 self.utctime = None
730 729 self.timeZone = 0
731 730 self.dstFlag = None
732 731 self.errorCount = None
733 732 self.blocksize = None
734 733 self.flagDecodeData = False # asumo q la data no esta decodificada
735 734 self.flagDeflipData = False # asumo q la data no esta sin flip
736 735 self.pairsList = None
737 736 self.nPoints = None
738 737
739 738 def getPairsList(self):
740 739
741 740 return self.pairsList
742 741
743 742 def getNoise(self, mode=2):
744 743
745 744 indR = numpy.where(self.lagR == 0)[0][0]
746 745 indT = numpy.where(self.lagT == 0)[0][0]
747 746
748 747 jspectra0 = self.data_corr[:, :, indR, :]
749 748 jspectra = copy.copy(jspectra0)
750 749
751 750 num_chan = jspectra.shape[0]
752 751 num_hei = jspectra.shape[2]
753 752
754 753 freq_dc = jspectra.shape[1] / 2
755 754 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
756 755
757 756 if ind_vel[0] < 0:
758 757 ind_vel[list(range(0, 1))] = ind_vel[list(
759 758 range(0, 1))] + self.num_prof
760 759
761 760 if mode == 1:
762 761 jspectra[:, freq_dc, :] = (
763 762 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
764 763
765 764 if mode == 2:
766 765
767 766 vel = numpy.array([-2, -1, 1, 2])
768 767 xx = numpy.zeros([4, 4])
769 768
770 769 for fil in range(4):
771 770 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
772 771
773 772 xx_inv = numpy.linalg.inv(xx)
774 773 xx_aux = xx_inv[0, :]
775 774
776 775 for ich in range(num_chan):
777 776 yy = jspectra[ich, ind_vel, :]
778 777 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
779 778
780 779 junkid = jspectra[ich, freq_dc, :] <= 0
781 780 cjunkid = sum(junkid)
782 781
783 782 if cjunkid.any():
784 783 jspectra[ich, freq_dc, junkid.nonzero()] = (
785 784 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
786 785
787 786 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
788 787
789 788 return noise
790 789
791 790 @property
792 791 def timeInterval(self):
793 792
794 793 return self.ippSeconds * self.nCohInt * self.nProfiles
795 794
796 795 def splitFunctions(self):
797 796
798 797 pairsList = self.pairsList
799 798 ccf_pairs = []
800 799 acf_pairs = []
801 800 ccf_ind = []
802 801 acf_ind = []
803 802 for l in range(len(pairsList)):
804 803 chan0 = pairsList[l][0]
805 804 chan1 = pairsList[l][1]
806 805
807 806 # Obteniendo pares de Autocorrelacion
808 807 if chan0 == chan1:
809 808 acf_pairs.append(chan0)
810 809 acf_ind.append(l)
811 810 else:
812 811 ccf_pairs.append(pairsList[l])
813 812 ccf_ind.append(l)
814 813
815 814 data_acf = self.data_cf[acf_ind]
816 815 data_ccf = self.data_cf[ccf_ind]
817 816
818 817 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
819 818
820 819 @property
821 820 def normFactor(self):
822 821 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
823 822 acf_pairs = numpy.array(acf_pairs)
824 823 normFactor = numpy.zeros((self.nPairs, self.nHeights))
825 824
826 825 for p in range(self.nPairs):
827 826 pair = self.pairsList[p]
828 827
829 828 ch0 = pair[0]
830 829 ch1 = pair[1]
831 830
832 831 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
833 832 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
834 833 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
835 834
836 835 return normFactor
837 836
838 837
839 838 class Parameters(Spectra):
840 839
841 840 groupList = None # List of Pairs, Groups, etc
842 841 data_param = None # Parameters obtained
843 842 data_pre = None # Data Pre Parametrization
844 843 data_SNR = None # Signal to Noise Ratio
845 844 abscissaList = None # Abscissa, can be velocities, lags or time
846 845 utctimeInit = None # Initial UTC time
847 846 paramInterval = None # Time interval to calculate Parameters in seconds
848 847 useLocalTime = True
849 848 # Fitting
850 849 data_error = None # Error of the estimation
851 850 constants = None
852 851 library = None
853 852 # Output signal
854 853 outputInterval = None # Time interval to calculate output signal in seconds
855 854 data_output = None # Out signal
856 855 nAvg = None
857 856 noise_estimation = None
858 857 GauSPC = None # Fit gaussian SPC
859 858
860 859 def __init__(self):
861 860 '''
862 861 Constructor
863 862 '''
864 863 self.radarControllerHeaderObj = RadarControllerHeader()
865 864 self.systemHeaderObj = SystemHeader()
866 865 self.type = "Parameters"
867 866 self.timeZone = 0
868 867
869 868 def getTimeRange1(self, interval):
870 869
871 870 datatime = []
872 871
873 872 if self.useLocalTime:
874 873 time1 = self.utctimeInit - self.timeZone * 60
875 874 else:
876 875 time1 = self.utctimeInit
877 876
878 877 datatime.append(time1)
879 878 datatime.append(time1 + interval)
880 879 datatime = numpy.array(datatime)
881 880
882 881 return datatime
883 882
884 883 @property
885 884 def timeInterval(self):
886 885
887 886 if hasattr(self, 'timeInterval1'):
888 887 return self.timeInterval1
889 888 else:
890 889 return self.paramInterval
891 890
892 891 def setValue(self, value):
893 892
894 893 print("This property should not be initialized")
895 894
896 895 return
897 896
898 897 def getNoise(self):
899 898
900 899 return self.spc_noise
901 900
902 901 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
903 902
904 903
905 904 class PlotterData(object):
906 905 '''
907 906 Object to hold data to be plotted
908 907 '''
909 908
910 909 MAXNUMX = 200
911 910 MAXNUMY = 200
912 911
913 912 def __init__(self, code, exp_code, localtime=True):
914 913
915 914 self.key = code
916 915 self.exp_code = exp_code
917 916 self.ready = False
918 917 self.flagNoData = False
919 918 self.localtime = localtime
920 919 self.data = {}
921 920 self.meta = {}
922 921 self.__heights = []
923 922
924 923 def __str__(self):
925 924 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
926 925 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
927 926
928 927 def __len__(self):
929 928 return len(self.data)
930 929
931 930 def __getitem__(self, key):
932 931 if isinstance(key, int):
933 932 return self.data[self.times[key]]
934 933 elif isinstance(key, str):
935 934 ret = numpy.array([self.data[x][key] for x in self.times])
936 935 if ret.ndim > 1:
937 936 ret = numpy.swapaxes(ret, 0, 1)
938 937 return ret
939 938
940 939 def __contains__(self, key):
941 940 return key in self.data[self.min_time]
942 941
943 942 def setup(self):
944 943 '''
945 944 Configure object
946 945 '''
947 946 self.type = ''
948 947 self.ready = False
949 948 del self.data
950 949 self.data = {}
951 950 self.__heights = []
952 951 self.__all_heights = set()
953 952
954 953 def shape(self, key):
955 954 '''
956 955 Get the shape of the one-element data for the given key
957 956 '''
958 957
959 958 if len(self.data[self.min_time][key]):
960 959 return self.data[self.min_time][key].shape
961 960 return (0,)
962 961
963 962 def update(self, data, tm, meta={}):
964 963 '''
965 964 Update data object with new dataOut
966 965 '''
967 966
968 967 self.data[tm] = data
969 968
970 969 for key, value in meta.items():
971 970 setattr(self, key, value)
972 971
973 972 def normalize_heights(self):
974 973 '''
975 974 Ensure same-dimension of the data for different heighList
976 975 '''
977 976
978 977 H = numpy.array(list(self.__all_heights))
979 978 H.sort()
980 979 for key in self.data:
981 980 shape = self.shape(key)[:-1] + H.shape
982 981 for tm, obj in list(self.data[key].items()):
983 982 h = self.__heights[self.times.tolist().index(tm)]
984 983 if H.size == h.size:
985 984 continue
986 985 index = numpy.where(numpy.in1d(H, h))[0]
987 986 dummy = numpy.zeros(shape) + numpy.nan
988 987 if len(shape) == 2:
989 988 dummy[:, index] = obj
990 989 else:
991 990 dummy[index] = obj
992 991 self.data[key][tm] = dummy
993 992
994 993 self.__heights = [H for tm in self.times]
995 994
996 995 def jsonify(self, tm, plot_name, plot_type, decimate=False):
997 996 '''
998 997 Convert data to json
999 998 '''
1000 999
1001 1000 meta = {}
1002 1001 meta['xrange'] = []
1003 1002 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1004 1003 tmp = self.data[tm][self.key]
1005 1004 shape = tmp.shape
1006 1005 if len(shape) == 2:
1007 1006 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1008 1007 elif len(shape) == 3:
1009 1008 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1010 1009 data = self.roundFloats(
1011 1010 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1012 1011 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1013 1012 else:
1014 1013 data = self.roundFloats(self.data[tm][self.key].tolist())
1015 1014
1016 1015 ret = {
1017 1016 'plot': plot_name,
1018 1017 'code': self.exp_code,
1019 1018 'time': float(tm),
1020 1019 'data': data,
1021 1020 }
1022 1021 meta['type'] = plot_type
1023 1022 meta['interval'] = float(self.interval)
1024 1023 meta['localtime'] = self.localtime
1025 1024 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1026 1025 meta.update(self.meta)
1027 1026 ret['metadata'] = meta
1028 1027 return json.dumps(ret)
1029 1028
1030 1029 @property
1031 1030 def times(self):
1032 1031 '''
1033 1032 Return the list of times of the current data
1034 1033 '''
1035 1034
1036 1035 ret = [t for t in self.data]
1037 1036 ret.sort()
1038 1037 return numpy.array(ret)
1039 1038
1040 1039 @property
1041 1040 def min_time(self):
1042 1041 '''
1043 1042 Return the minimun time value
1044 1043 '''
1045 1044
1046 1045 return self.times[0]
1047 1046
1048 1047 @property
1049 1048 def max_time(self):
1050 1049 '''
1051 1050 Return the maximun time value
1052 1051 '''
1053 1052
1054 1053 return self.times[-1]
1055 1054
1056 1055 # @property
1057 1056 # def heights(self):
1058 1057 # '''
1059 1058 # Return the list of heights of the current data
1060 1059 # '''
1061 1060
1062 1061 # return numpy.array(self.__heights[-1])
1063 1062
1064 1063 @staticmethod
1065 1064 def roundFloats(obj):
1066 1065 if isinstance(obj, list):
1067 1066 return list(map(PlotterData.roundFloats, obj))
1068 1067 elif isinstance(obj, float):
1069 1068 return round(obj, 2)
@@ -1,519 +1,509
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
99 98 data = {
100 99 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
101 100 }
102
103 101 return data, {}
104 102
105 103 class SpectralWidthPlot(RTIPlot):
106 104 '''
107 105 Plot for Spectral Width Data (2nd moment)
108 106 '''
109 107
110 108 CODE = 'width'
111 109 colormap = 'jet'
112 110
113 111 def update(self, dataOut):
114 112
115 113 data = {
116 114 'width': dataOut.data_width
117 115 }
118 116
119 117 return data, {}
120 118
121 119 class SkyMapPlot(Plot):
122 120 '''
123 121 Plot for meteors detection data
124 122 '''
125 123
126 124 CODE = 'param'
127 125
128 126 def setup(self):
129 127
130 128 self.ncols = 1
131 129 self.nrows = 1
132 130 self.width = 7.2
133 131 self.height = 7.2
134 132 self.nplots = 1
135 133 self.xlabel = 'Zonal Zenith Angle (deg)'
136 134 self.ylabel = 'Meridional Zenith Angle (deg)'
137 135 self.polar = True
138 136 self.ymin = -180
139 137 self.ymax = 180
140 138 self.colorbar = False
141 139
142 140 def plot(self):
143 141
144 142 arrayParameters = numpy.concatenate(self.data['param'])
145 143 error = arrayParameters[:, -1]
146 144 indValid = numpy.where(error == 0)[0]
147 145 finalMeteor = arrayParameters[indValid, :]
148 146 finalAzimuth = finalMeteor[:, 3]
149 147 finalZenith = finalMeteor[:, 4]
150 148
151 149 x = finalAzimuth * numpy.pi / 180
152 150 y = finalZenith
153 151
154 152 ax = self.axes[0]
155 153
156 154 if ax.firsttime:
157 155 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
158 156 else:
159 157 ax.plot.set_data(x, y)
160 158
161 159 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
162 160 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
163 161 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
164 162 dt2,
165 163 len(x))
166 164 self.titles[0] = title
167 165
168 166
169 167 class GenericRTIPlot(Plot):
170 168 '''
171 169 Plot for data_xxxx object
172 170 '''
173 171
174 172 CODE = 'param'
175 173 colormap = 'viridis'
176 174 plot_type = 'pcolorbuffer'
177 175
178 176 def setup(self):
179 177 self.xaxis = 'time'
180 178 self.ncols = 1
181 179 self.nrows = self.data.shape('param')[0]
182 180 self.nplots = self.nrows
183 181 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
184 182
185 183 if not self.xlabel:
186 184 self.xlabel = 'Time'
187 185
188 186 self.ylabel = 'Range [km]'
189 187 if not self.titles:
190 188 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
191 189
192 190 def update(self, dataOut):
193 191
194 192 data = {
195 193 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
196 194 }
197 195
198 196 meta = {}
199 197
200 198 return data, meta
201 199
202 200 def plot(self):
203 201 # self.data.normalize_heights()
204 202 self.x = self.data.times
205 203 self.y = self.data.yrange
206 204 self.z = self.data['param']
207
208 205 self.z = 10*numpy.log10(self.z)
209
210 206 self.z = numpy.ma.masked_invalid(self.z)
211 207
212 208 if self.decimation is None:
213 209 x, y, z = self.fill_gaps(self.x, self.y, self.z)
214 210 else:
215 211 x, y, z = self.fill_gaps(*self.decimate())
216 212
217 213 for n, ax in enumerate(self.axes):
218 214
219 215 self.zmax = self.zmax if self.zmax is not None else numpy.max(
220 216 self.z[n])
221 217 self.zmin = self.zmin if self.zmin is not None else numpy.min(
222 218 self.z[n])
223 219
224 220 if ax.firsttime:
225 221 if self.zlimits is not None:
226 222 self.zmin, self.zmax = self.zlimits[n]
227 223
228 224 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
229 225 vmin=self.zmin,
230 226 vmax=self.zmax,
231 227 cmap=self.cmaps[n]
232 228 )
233 229 else:
234 230 if self.zlimits is not None:
235 231 self.zmin, self.zmax = self.zlimits[n]
236 232 ax.collections.remove(ax.collections[0])
237 233 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
238 234 vmin=self.zmin,
239 235 vmax=self.zmax,
240 236 cmap=self.cmaps[n]
241 237 )
242 238
243 239
244 240 class PolarMapPlot(Plot):
245 241 '''
246 242 Plot for weather radar
247 243 '''
248 244
249 245 CODE = 'param'
250 246 colormap = 'seismic'
251 247
252 248 def setup(self):
253 249 self.ncols = 1
254 250 self.nrows = 1
255 251 self.width = 9
256 252 self.height = 8
257 253 self.mode = self.data.meta['mode']
258 254 if self.channels is not None:
259 255 self.nplots = len(self.channels)
260 256 self.nrows = len(self.channels)
261 257 else:
262 258 self.nplots = self.data.shape(self.CODE)[0]
263 259 self.nrows = self.nplots
264 260 self.channels = list(range(self.nplots))
265 261 if self.mode == 'E':
266 262 self.xlabel = 'Longitude'
267 263 self.ylabel = 'Latitude'
268 264 else:
269 265 self.xlabel = 'Range (km)'
270 266 self.ylabel = 'Height (km)'
271 267 self.bgcolor = 'white'
272 268 self.cb_labels = self.data.meta['units']
273 269 self.lat = self.data.meta['latitude']
274 270 self.lon = self.data.meta['longitude']
275 271 self.xmin, self.xmax = float(
276 272 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
277 273 self.ymin, self.ymax = float(
278 274 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
279 275 # self.polar = True
280 276
281 277 def plot(self):
282 278
283 279 for n, ax in enumerate(self.axes):
284 280 data = self.data['param'][self.channels[n]]
285 281
286 282 zeniths = numpy.linspace(
287 283 0, self.data.meta['max_range'], data.shape[1])
288 284 if self.mode == 'E':
289 285 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
290 286 r, theta = numpy.meshgrid(zeniths, azimuths)
291 287 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
292 288 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
293 289 x = km2deg(x) + self.lon
294 290 y = km2deg(y) + self.lat
295 291 else:
296 292 azimuths = numpy.radians(self.data.yrange)
297 293 r, theta = numpy.meshgrid(zeniths, azimuths)
298 294 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
299 295 self.y = zeniths
300 296
301 297 if ax.firsttime:
302 298 if self.zlimits is not None:
303 299 self.zmin, self.zmax = self.zlimits[n]
304 300 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
305 301 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
306 302 vmin=self.zmin,
307 303 vmax=self.zmax,
308 304 cmap=self.cmaps[n])
309 305 else:
310 306 if self.zlimits is not None:
311 307 self.zmin, self.zmax = self.zlimits[n]
312 308 ax.collections.remove(ax.collections[0])
313 309 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
314 310 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
315 311 vmin=self.zmin,
316 312 vmax=self.zmax,
317 313 cmap=self.cmaps[n])
318 314
319 315 if self.mode == 'A':
320 316 continue
321 317
322 318 # plot district names
323 319 f = open('/data/workspace/schain_scripts/distrito.csv')
324 320 for line in f:
325 321 label, lon, lat = [s.strip() for s in line.split(',') if s]
326 322 lat = float(lat)
327 323 lon = float(lon)
328 324 # ax.plot(lon, lat, '.b', ms=2)
329 325 ax.text(lon, lat, label.decode('utf8'), ha='center',
330 326 va='bottom', size='8', color='black')
331 327
332 328 # plot limites
333 329 limites = []
334 330 tmp = []
335 331 for line in open('/data/workspace/schain_scripts/lima.csv'):
336 332 if '#' in line:
337 333 if tmp:
338 334 limites.append(tmp)
339 335 tmp = []
340 336 continue
341 337 values = line.strip().split(',')
342 338 tmp.append((float(values[0]), float(values[1])))
343 339 for points in limites:
344 340 ax.add_patch(
345 341 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
346 342
347 343 # plot Cuencas
348 344 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
349 345 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
350 346 values = [line.strip().split(',') for line in f]
351 347 points = [(float(s[0]), float(s[1])) for s in values]
352 348 ax.add_patch(Polygon(points, ec='b', fc='none'))
353 349
354 350 # plot grid
355 351 for r in (15, 30, 45, 60):
356 352 ax.add_artist(plt.Circle((self.lon, self.lat),
357 353 km2deg(r), color='0.6', fill=False, lw=0.2))
358 354 ax.text(
359 355 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
360 356 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
361 357 '{}km'.format(r),
362 358 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
363 359
364 360 if self.mode == 'E':
365 361 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
366 362 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
367 363 else:
368 364 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
369 365 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
370 366
371 367 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
372 368 self.titles = ['{} {}'.format(
373 369 self.data.parameters[x], title) for x in self.channels]
374 370
375 371 class WeatherPlot(Plot):
376 372 CODE = 'weather'
377 373 plot_name = 'weather'
378 374 plot_type = 'ppistyle'
379 375 buffering = False
380 376
381 377 def setup(self):
382 378 self.ncols = 1
383 379 self.nrows = 1
384 380 self.nplots= 1
385 381 self.ylabel= 'Range [Km]'
386 382 self.titles= ['Weather']
387 383 self.colorbar=False
388 384 self.width =8
389 385 self.height =8
390 386 self.ini =0
391 387 self.len_azi =0
392 388 self.buffer_ini = None
393 389 self.buffer_azi = None
394 390 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
395 391 self.flag =0
396 392 self.indicador= 0
397 393
398 394 def update(self, dataOut):
399 395
400 396 data = {}
401 397 meta = {}
402 print("aprox",dataOut.data_360[0])
403 398 data['weather'] = 10*numpy.log10(dataOut.data_360[0]/(250.0))
404 #print(data['weather'])
405 399 data['azi'] = dataOut.data_azi
406 print("UPDATE",data['azi'])
407 400 return data, meta
408 401
409 402 def const_ploteo(self,data_weather,data_azi,step,res):
410 #print("data_weather",data_weather)
411 print("data_azi",data_azi)
412 print("step",step)
413 403 if self.ini==0:
414 404 #------- AZIMUTH
415 405 n = (360/res)-len(data_azi)
416 406 start = data_azi[-1] + res
417 407 end = data_azi[0] - res
418 408 if start>end:
419 409 end = end + 360
420 410 azi_vacia = numpy.linspace(start,end,int(n))
421 411 azi_vacia = numpy.where(azi_vacia>360,azi_vacia-360,azi_vacia)
422 412 data_azi = numpy.hstack((data_azi,azi_vacia))
423 413 # RADAR
424 414 val_mean = numpy.mean(data_weather[:,0])
425 415 data_weather_cmp = numpy.ones([(360-data_weather.shape[0]),data_weather.shape[1]])*val_mean
426 416 data_weather = numpy.vstack((data_weather,data_weather_cmp))
427 417 else:
428 418 # azimuth
429 419 flag=0
430 420 start_azi = self.res_azi[0]
431 421 start = data_azi[0]
432 422 end = data_azi[-1]
433 423 print("start",start)
434 424 print("end",end)
435 425 if start< start_azi:
436 426 start = start +360
437 427 if end <start_azi:
438 428 end = end +360
439 429
440 430 print("start",start)
441 431 print("end",end)
442 432 #### AQUI SERA LA MAGIA
443 433 pos_ini = int((start-start_azi)/res)
444 434 len_azi = len(data_azi)
445 435 if (360-pos_ini)<len_azi:
446 436 if pos_ini+1==360:
447 437 pos_ini=0
448 438 else:
449 439 flag=1
450 440 dif= 360-pos_ini
451 441 comp= len_azi-dif
452 442
453 443 print(pos_ini)
454 444 print(len_azi)
455 445 print("shape",self.res_azi.shape)
456 446 if flag==0:
457 447 # AZIMUTH
458 448 self.res_azi[pos_ini:pos_ini+len_azi] = data_azi
459 449 # RADAR
460 450 self.res_weather[pos_ini:pos_ini+len_azi,:] = data_weather
461 451 else:
462 452 # AZIMUTH
463 453 self.res_azi[pos_ini:pos_ini+dif] = data_azi[0:dif]
464 454 self.res_azi[0:comp] = data_azi[dif:]
465 455 # RADAR
466 456 self.res_weather[pos_ini:pos_ini+dif,:] = data_weather[0:dif,:]
467 457 self.res_weather[0:comp,:] = data_weather[dif:,:]
468 458 flag=0
469 459 data_azi = self.res_azi
470 460 data_weather = self.res_weather
471 461
472 462 return data_weather,data_azi
473 463
474 464 def plot(self):
475 465 print("--------------------------------------",self.ini,"-----------------------------------")
476 466 #numpy.set_printoptions(suppress=True)
477 467 #print(self.data.times)
478 468 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
479 469 data = self.data[-1]
480 470 # ALTURA altura_tmp_h
481 471 altura_h = (data['weather'].shape[1])/10.0
482 472 stoprange = float(altura_h*1.5)#stoprange = float(33*1.5) por ahora 400
483 473 rangestep = float(0.15)
484 474 r = numpy.arange(0, stoprange, rangestep)
485 475 self.y = 2*r
486 476 # RADAR
487 477 #data_weather = data['weather']
488 478 # PEDESTAL
489 479 #data_azi = data['azi']
490 480 res = 1
491 481 # STEP
492 482 step = (360/(res*data['weather'].shape[0]))
493 483 #print("shape wr_data", wr_data.shape)
494 484 #print("shape wr_azi",wr_azi.shape)
495 485 #print("step",step)
496 486 print("Time---->",self.data.times[-1],thisDatetime)
497 487 #print("alturas", len(self.y))
498 488 self.res_weather, self.res_azi = self.const_ploteo(data_weather=data['weather'],data_azi=data['azi'],step=step,res=res)
499 489 #numpy.set_printoptions(suppress=True)
500 490 #print("resultado",self.res_azi)
501 491 ##########################################################
502 492 ################# PLOTEO ###################
503 493 ##########################################################
504 494
505 495 for i,ax in enumerate(self.axes):
506 496 if ax.firsttime:
507 497 plt.clf()
508 498 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)
509 499 else:
510 500 plt.clf()
511 501 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)
512 502 caax = cgax.parasites[0]
513 503 paax = cgax.parasites[1]
514 504 cbar = plt.gcf().colorbar(pm, pad=0.075)
515 505 caax.set_xlabel('x_range [km]')
516 506 caax.set_ylabel('y_range [km]')
517 507 plt.text(1.0, 1.05, 'azimuth '+str(thisDatetime)+"step"+str(self.ini), transform=caax.transAxes, va='bottom',ha='right')
518 508
519 509 self.ini= self.ini+1
@@ -1,745 +1,743
1 1 # Copyright (c) 2012-2021 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13
14 14
15 15 class SpectraPlot(Plot):
16 16 '''
17 17 Plot for Spectra data
18 18 '''
19 19
20 20 CODE = 'spc'
21 21 colormap = 'jet'
22 22 plot_type = 'pcolor'
23 23 buffering = False
24 24
25 25 def setup(self):
26 26 self.nplots = len(self.data.channels)
27 27 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
28 28 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
29 29 self.height = 2.6 * self.nrows
30 30 self.cb_label = 'dB'
31 31 if self.showprofile:
32 32 self.width = 4 * self.ncols
33 33 else:
34 34 self.width = 3.5 * self.ncols
35 35 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
36 36 self.ylabel = 'Range [km]'
37 37
38 38 def update(self, dataOut):
39 39
40 40 data = {}
41 41 meta = {}
42 42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 43 data['spc'] = spc
44 44 data['rti'] = dataOut.getPower()
45 45 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 46 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
47 47
48 48 if self.CODE == 'spc_moments':
49 49 data['moments'] = dataOut.moments
50 50 # data['spc'] = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
51 51 if self.CODE == 'gaussian_fit':
52 52 # data['moments'] = dataOut.moments
53 53 data['gaussfit'] = dataOut.DGauFitParams
54 54 # data['spc'] = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
55 55
56 56 return data, meta
57 57
58 58 def plot(self):
59 59 if self.xaxis == "frequency":
60 60 x = self.data.xrange[0]
61 61 self.xlabel = "Frequency (kHz)"
62 62 elif self.xaxis == "time":
63 63 x = self.data.xrange[1]
64 64 self.xlabel = "Time (ms)"
65 65 else:
66 66 x = self.data.xrange[2]
67 67 self.xlabel = "Velocity (m/s)"
68 68
69 69 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
70 70 x = self.data.xrange[2]
71 71 self.xlabel = "Velocity (m/s)"
72 72
73 73 self.titles = []
74 74
75 75 y = self.data.yrange
76 76 self.y = y
77 77
78 78 data = self.data[-1]
79 79 z = data['spc']
80 80
81 81 for n, ax in enumerate(self.axes):
82 82 noise = data['noise'][n]
83 83 if self.CODE == 'spc_moments':
84 84 mean = data['moments'][n, 1]
85 85 if self.CODE == 'gaussian_fit':
86 86 # mean = data['moments'][n, 1]
87 87 gau0 = data['gaussfit'][n][2,:,0]
88 88 gau1 = data['gaussfit'][n][2,:,1]
89 89 if ax.firsttime:
90 90 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
91 91 self.xmin = self.xmin if self.xmin else -self.xmax
92 92 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
93 93 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
94 94 ax.plt = ax.pcolormesh(x, y, z[n].T,
95 95 vmin=self.zmin,
96 96 vmax=self.zmax,
97 97 cmap=plt.get_cmap(self.colormap)
98 98 )
99 99
100 100 if self.showprofile:
101 101 ax.plt_profile = self.pf_axes[n].plot(
102 102 data['rti'][n], y)[0]
103 103 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
104 104 color="k", linestyle="dashed", lw=1)[0]
105 105 if self.CODE == 'spc_moments':
106 106 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
107 107 if self.CODE == 'gaussian_fit':
108 108 # ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
109 109 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
110 110 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
111 111 else:
112 112 ax.plt.set_array(z[n].T.ravel())
113 113 if self.showprofile:
114 114 ax.plt_profile.set_data(data['rti'][n], y)
115 115 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
116 116 if self.CODE == 'spc_moments':
117 117 ax.plt_mean.set_data(mean, y)
118 118 if self.CODE == 'gaussian_fit':
119 119 # ax.plt_mean.set_data(mean, y)
120 120 ax.plt_gau0.set_data(gau0, y)
121 121 ax.plt_gau1.set_data(gau1, y)
122 122 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
123 123
124 124
125 125 class CrossSpectraPlot(Plot):
126 126
127 127 CODE = 'cspc'
128 128 colormap = 'jet'
129 129 plot_type = 'pcolor'
130 130 zmin_coh = None
131 131 zmax_coh = None
132 132 zmin_phase = None
133 133 zmax_phase = None
134 134
135 135 def setup(self):
136 136
137 137 self.ncols = 4
138 138 self.nplots = len(self.data.pairs) * 2
139 139 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
140 140 self.width = 3.1 * self.ncols
141 141 self.height = 2.6 * self.nrows
142 142 self.ylabel = 'Range [km]'
143 143 self.showprofile = False
144 144 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
145 145
146 146 def update(self, dataOut):
147 147
148 148 data = {}
149 149 meta = {}
150 150
151 151 spc = dataOut.data_spc
152 152 cspc = dataOut.data_cspc
153 153 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
154 154 meta['pairs'] = dataOut.pairsList
155 155
156 156 tmp = []
157 157
158 158 for n, pair in enumerate(meta['pairs']):
159 159 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
160 160 coh = numpy.abs(out)
161 161 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
162 162 tmp.append(coh)
163 163 tmp.append(phase)
164 164
165 165 data['cspc'] = numpy.array(tmp)
166 166
167 167 return data, meta
168 168
169 169 def plot(self):
170 170
171 171 if self.xaxis == "frequency":
172 172 x = self.data.xrange[0]
173 173 self.xlabel = "Frequency (kHz)"
174 174 elif self.xaxis == "time":
175 175 x = self.data.xrange[1]
176 176 self.xlabel = "Time (ms)"
177 177 else:
178 178 x = self.data.xrange[2]
179 179 self.xlabel = "Velocity (m/s)"
180 180
181 181 self.titles = []
182 182
183 183 y = self.data.yrange
184 184 self.y = y
185 185
186 186 data = self.data[-1]
187 187 cspc = data['cspc']
188 188
189 189 for n in range(len(self.data.pairs)):
190 190 pair = self.data.pairs[n]
191 191 coh = cspc[n*2]
192 192 phase = cspc[n*2+1]
193 193 ax = self.axes[2 * n]
194 194 if ax.firsttime:
195 195 ax.plt = ax.pcolormesh(x, y, coh.T,
196 196 vmin=0,
197 197 vmax=1,
198 198 cmap=plt.get_cmap(self.colormap_coh)
199 199 )
200 200 else:
201 201 ax.plt.set_array(coh.T.ravel())
202 202 self.titles.append(
203 203 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
204 204
205 205 ax = self.axes[2 * n + 1]
206 206 if ax.firsttime:
207 207 ax.plt = ax.pcolormesh(x, y, phase.T,
208 208 vmin=-180,
209 209 vmax=180,
210 210 cmap=plt.get_cmap(self.colormap_phase)
211 211 )
212 212 else:
213 213 ax.plt.set_array(phase.T.ravel())
214 214 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
215 215
216 216
217 217 class RTIPlot(Plot):
218 218 '''
219 219 Plot for RTI data
220 220 '''
221 221
222 222 CODE = 'rti'
223 223 colormap = 'jet'
224 224 plot_type = 'pcolorbuffer'
225 225
226 226 def setup(self):
227 227 self.xaxis = 'time'
228 228 self.ncols = 1
229 print("ch",self.data.channels)
230 229 self.nrows = len(self.data.channels)
231 230 self.nplots = len(self.data.channels)
232 231 self.ylabel = 'Range [km]'
233 232 self.xlabel = 'Time'
234 233 self.cb_label = 'dB'
235 234 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
236 235 self.titles = ['{} Channel {}'.format(
237 236 self.CODE.upper(), x) for x in range(self.nrows)]
238 237
239 238 def update(self, dataOut):
240 239
241 240 data = {}
242 241 meta = {}
243 242 data['rti'] = dataOut.getPower()
244 243 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
245 244
246 245 return data, meta
247 246
248 247 def plot(self):
249 248 self.x = self.data.times
250 249 self.y = self.data.yrange
251 250 self.z = self.data[self.CODE]
252 251 self.z = numpy.ma.masked_invalid(self.z)
253 252
254 253 if self.decimation is None:
255 254 x, y, z = self.fill_gaps(self.x, self.y, self.z)
256 255 else:
257 256 x, y, z = self.fill_gaps(*self.decimate())
258 257
259 258 for n, ax in enumerate(self.axes):
260 259 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
261 260 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
262 261 data = self.data[-1]
263 262 if ax.firsttime:
264 263 ax.plt = ax.pcolormesh(x, y, z[n].T,
265 264 vmin=self.zmin,
266 265 vmax=self.zmax,
267 266 cmap=plt.get_cmap(self.colormap)
268 267 )
269 268 if self.showprofile:
270 print("test-------------------------------------1")
271 269 ax.plot_profile = self.pf_axes[n].plot(
272 270 data['rti'][n], self.y)[0]
273 271 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
274 272 color="k", linestyle="dashed", lw=1)[0]
275 273 else:
276 274 ax.collections.remove(ax.collections[0])
277 275 ax.plt = ax.pcolormesh(x, y, z[n].T,
278 276 vmin=self.zmin,
279 277 vmax=self.zmax,
280 278 cmap=plt.get_cmap(self.colormap)
281 279 )
282 280 if self.showprofile:
283 281 ax.plot_profile.set_data(data['rti'][n], self.y)
284 282 ax.plot_noise.set_data(numpy.repeat(
285 283 data['noise'][n], len(self.y)), self.y)
286 284
287 285
288 286 class CoherencePlot(RTIPlot):
289 287 '''
290 288 Plot for Coherence data
291 289 '''
292 290
293 291 CODE = 'coh'
294 292
295 293 def setup(self):
296 294 self.xaxis = 'time'
297 295 self.ncols = 1
298 296 self.nrows = len(self.data.pairs)
299 297 self.nplots = len(self.data.pairs)
300 298 self.ylabel = 'Range [km]'
301 299 self.xlabel = 'Time'
302 300 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
303 301 if self.CODE == 'coh':
304 302 self.cb_label = ''
305 303 self.titles = [
306 304 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
307 305 else:
308 306 self.cb_label = 'Degrees'
309 307 self.titles = [
310 308 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
311 309
312 310 def update(self, dataOut):
313 311
314 312 data = {}
315 313 meta = {}
316 314 data['coh'] = dataOut.getCoherence()
317 315 meta['pairs'] = dataOut.pairsList
318 316
319 317 return data, meta
320 318
321 319 class PhasePlot(CoherencePlot):
322 320 '''
323 321 Plot for Phase map data
324 322 '''
325 323
326 324 CODE = 'phase'
327 325 colormap = 'seismic'
328 326
329 327 def update(self, dataOut):
330 328
331 329 data = {}
332 330 meta = {}
333 331 data['phase'] = dataOut.getCoherence(phase=True)
334 332 meta['pairs'] = dataOut.pairsList
335 333
336 334 return data, meta
337 335
338 336 class NoisePlot(Plot):
339 337 '''
340 338 Plot for noise
341 339 '''
342 340
343 341 CODE = 'noise'
344 342 plot_type = 'scatterbuffer'
345 343
346 344 def setup(self):
347 345 self.xaxis = 'time'
348 346 self.ncols = 1
349 347 self.nrows = 1
350 348 self.nplots = 1
351 349 self.ylabel = 'Intensity [dB]'
352 350 self.xlabel = 'Time'
353 351 self.titles = ['Noise']
354 352 self.colorbar = False
355 353 self.plots_adjust.update({'right': 0.85 })
356 354
357 355 def update(self, dataOut):
358 356
359 357 data = {}
360 358 meta = {}
361 359 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
362 360 meta['yrange'] = numpy.array([])
363 361
364 362 return data, meta
365 363
366 364 def plot(self):
367 365
368 366 x = self.data.times
369 367 xmin = self.data.min_time
370 368 xmax = xmin + self.xrange * 60 * 60
371 369 Y = self.data['noise']
372 370
373 371 if self.axes[0].firsttime:
374 372 self.ymin = numpy.nanmin(Y) - 5
375 373 self.ymax = numpy.nanmax(Y) + 5
376 374 for ch in self.data.channels:
377 375 y = Y[ch]
378 376 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
379 377 plt.legend(bbox_to_anchor=(1.18, 1.0))
380 378 else:
381 379 for ch in self.data.channels:
382 380 y = Y[ch]
383 381 self.axes[0].lines[ch].set_data(x, y)
384 382
385 383
386 384 class PowerProfilePlot(Plot):
387 385
388 386 CODE = 'pow_profile'
389 387 plot_type = 'scatter'
390 388
391 389 def setup(self):
392 390
393 391 self.ncols = 1
394 392 self.nrows = 1
395 393 self.nplots = 1
396 394 self.height = 4
397 395 self.width = 3
398 396 self.ylabel = 'Range [km]'
399 397 self.xlabel = 'Intensity [dB]'
400 398 self.titles = ['Power Profile']
401 399 self.colorbar = False
402 400
403 401 def update(self, dataOut):
404 402
405 403 data = {}
406 404 meta = {}
407 405 data[self.CODE] = dataOut.getPower()
408 406
409 407 return data, meta
410 408
411 409 def plot(self):
412 410
413 411 y = self.data.yrange
414 412 self.y = y
415 413
416 414 x = self.data[-1][self.CODE]
417 415
418 416 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
419 417 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
420 418
421 419 if self.axes[0].firsttime:
422 420 for ch in self.data.channels:
423 421 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
424 422 plt.legend()
425 423 else:
426 424 for ch in self.data.channels:
427 425 self.axes[0].lines[ch].set_data(x[ch], y)
428 426
429 427
430 428 class SpectraCutPlot(Plot):
431 429
432 430 CODE = 'spc_cut'
433 431 plot_type = 'scatter'
434 432 buffering = False
435 433
436 434 def setup(self):
437 435
438 436 self.nplots = len(self.data.channels)
439 437 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
440 438 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
441 439 self.width = 3.4 * self.ncols + 1.5
442 440 self.height = 3 * self.nrows
443 441 self.ylabel = 'Power [dB]'
444 442 self.colorbar = False
445 443 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
446 444
447 445 def update(self, dataOut):
448 446
449 447 data = {}
450 448 meta = {}
451 449 spc = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
452 450 data['spc'] = spc
453 451 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
454 452 if self.CODE == 'cut_gaussian_fit':
455 453 data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
456 454 data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
457 455 return data, meta
458 456
459 457 def plot(self):
460 458 if self.xaxis == "frequency":
461 459 x = self.data.xrange[0][1:]
462 460 self.xlabel = "Frequency (kHz)"
463 461 elif self.xaxis == "time":
464 462 x = self.data.xrange[1]
465 463 self.xlabel = "Time (ms)"
466 464 else:
467 465 x = self.data.xrange[2][:-1]
468 466 self.xlabel = "Velocity (m/s)"
469 467
470 468 if self.CODE == 'cut_gaussian_fit':
471 469 x = self.data.xrange[2][:-1]
472 470 self.xlabel = "Velocity (m/s)"
473 471
474 472 self.titles = []
475 473
476 474 y = self.data.yrange
477 475 data = self.data[-1]
478 476 z = data['spc']
479 477
480 478 if self.height_index:
481 479 index = numpy.array(self.height_index)
482 480 else:
483 481 index = numpy.arange(0, len(y), int((len(y))/9))
484 482
485 483 for n, ax in enumerate(self.axes):
486 484 if self.CODE == 'cut_gaussian_fit':
487 485 gau0 = data['gauss_fit0']
488 486 gau1 = data['gauss_fit1']
489 487 if ax.firsttime:
490 488 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
491 489 self.xmin = self.xmin if self.xmin else -self.xmax
492 490 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
493 491 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
494 492 ax.plt = ax.plot(x, z[n, :, index].T, lw=0.25)
495 493 if self.CODE == 'cut_gaussian_fit':
496 494 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
497 495 for i, line in enumerate(ax.plt_gau0):
498 496 line.set_color(ax.plt[i].get_color())
499 497 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
500 498 for i, line in enumerate(ax.plt_gau1):
501 499 line.set_color(ax.plt[i].get_color())
502 500 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
503 501 self.figures[0].legend(ax.plt, labels, loc='center right')
504 502 else:
505 503 for i, line in enumerate(ax.plt):
506 504 line.set_data(x, z[n, :, index[i]].T)
507 505 for i, line in enumerate(ax.plt_gau0):
508 506 line.set_data(x, gau0[n, :, index[i]].T)
509 507 line.set_color(ax.plt[i].get_color())
510 508 for i, line in enumerate(ax.plt_gau1):
511 509 line.set_data(x, gau1[n, :, index[i]].T)
512 510 line.set_color(ax.plt[i].get_color())
513 511 self.titles.append('CH {}'.format(n))
514 512
515 513
516 514 class BeaconPhase(Plot):
517 515
518 516 __isConfig = None
519 517 __nsubplots = None
520 518
521 519 PREFIX = 'beacon_phase'
522 520
523 521 def __init__(self):
524 522 Plot.__init__(self)
525 523 self.timerange = 24*60*60
526 524 self.isConfig = False
527 525 self.__nsubplots = 1
528 526 self.counter_imagwr = 0
529 527 self.WIDTH = 800
530 528 self.HEIGHT = 400
531 529 self.WIDTHPROF = 120
532 530 self.HEIGHTPROF = 0
533 531 self.xdata = None
534 532 self.ydata = None
535 533
536 534 self.PLOT_CODE = BEACON_CODE
537 535
538 536 self.FTP_WEI = None
539 537 self.EXP_CODE = None
540 538 self.SUB_EXP_CODE = None
541 539 self.PLOT_POS = None
542 540
543 541 self.filename_phase = None
544 542
545 543 self.figfile = None
546 544
547 545 self.xmin = None
548 546 self.xmax = None
549 547
550 548 def getSubplots(self):
551 549
552 550 ncol = 1
553 551 nrow = 1
554 552
555 553 return nrow, ncol
556 554
557 555 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
558 556
559 557 self.__showprofile = showprofile
560 558 self.nplots = nplots
561 559
562 560 ncolspan = 7
563 561 colspan = 6
564 562 self.__nsubplots = 2
565 563
566 564 self.createFigure(id = id,
567 565 wintitle = wintitle,
568 566 widthplot = self.WIDTH+self.WIDTHPROF,
569 567 heightplot = self.HEIGHT+self.HEIGHTPROF,
570 568 show=show)
571 569
572 570 nrow, ncol = self.getSubplots()
573 571
574 572 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
575 573
576 574 def save_phase(self, filename_phase):
577 575 f = open(filename_phase,'w+')
578 576 f.write('\n\n')
579 577 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
580 578 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
581 579 f.close()
582 580
583 581 def save_data(self, filename_phase, data, data_datetime):
584 582 f=open(filename_phase,'a')
585 583 timetuple_data = data_datetime.timetuple()
586 584 day = str(timetuple_data.tm_mday)
587 585 month = str(timetuple_data.tm_mon)
588 586 year = str(timetuple_data.tm_year)
589 587 hour = str(timetuple_data.tm_hour)
590 588 minute = str(timetuple_data.tm_min)
591 589 second = str(timetuple_data.tm_sec)
592 590 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
593 591 f.close()
594 592
595 593 def plot(self):
596 594 log.warning('TODO: Not yet implemented...')
597 595
598 596 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
599 597 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
600 598 timerange=None,
601 599 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
602 600 server=None, folder=None, username=None, password=None,
603 601 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
604 602
605 603 if dataOut.flagNoData:
606 604 return dataOut
607 605
608 606 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
609 607 return
610 608
611 609 if pairsList == None:
612 610 pairsIndexList = dataOut.pairsIndexList[:10]
613 611 else:
614 612 pairsIndexList = []
615 613 for pair in pairsList:
616 614 if pair not in dataOut.pairsList:
617 615 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
618 616 pairsIndexList.append(dataOut.pairsList.index(pair))
619 617
620 618 if pairsIndexList == []:
621 619 return
622 620
623 621 # if len(pairsIndexList) > 4:
624 622 # pairsIndexList = pairsIndexList[0:4]
625 623
626 624 hmin_index = None
627 625 hmax_index = None
628 626
629 627 if hmin != None and hmax != None:
630 628 indexes = numpy.arange(dataOut.nHeights)
631 629 hmin_list = indexes[dataOut.heightList >= hmin]
632 630 hmax_list = indexes[dataOut.heightList <= hmax]
633 631
634 632 if hmin_list.any():
635 633 hmin_index = hmin_list[0]
636 634
637 635 if hmax_list.any():
638 636 hmax_index = hmax_list[-1]+1
639 637
640 638 x = dataOut.getTimeRange()
641 639
642 640 thisDatetime = dataOut.datatime
643 641
644 642 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
645 643 xlabel = "Local Time"
646 644 ylabel = "Phase (degrees)"
647 645
648 646 update_figfile = False
649 647
650 648 nplots = len(pairsIndexList)
651 649 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
652 650 phase_beacon = numpy.zeros(len(pairsIndexList))
653 651 for i in range(nplots):
654 652 pair = dataOut.pairsList[pairsIndexList[i]]
655 653 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
656 654 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
657 655 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
658 656 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
659 657 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
660 658
661 659 if dataOut.beacon_heiIndexList:
662 660 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
663 661 else:
664 662 phase_beacon[i] = numpy.average(phase)
665 663
666 664 if not self.isConfig:
667 665
668 666 nplots = len(pairsIndexList)
669 667
670 668 self.setup(id=id,
671 669 nplots=nplots,
672 670 wintitle=wintitle,
673 671 showprofile=showprofile,
674 672 show=show)
675 673
676 674 if timerange != None:
677 675 self.timerange = timerange
678 676
679 677 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
680 678
681 679 if ymin == None: ymin = 0
682 680 if ymax == None: ymax = 360
683 681
684 682 self.FTP_WEI = ftp_wei
685 683 self.EXP_CODE = exp_code
686 684 self.SUB_EXP_CODE = sub_exp_code
687 685 self.PLOT_POS = plot_pos
688 686
689 687 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
690 688 self.isConfig = True
691 689 self.figfile = figfile
692 690 self.xdata = numpy.array([])
693 691 self.ydata = numpy.array([])
694 692
695 693 update_figfile = True
696 694
697 695 #open file beacon phase
698 696 path = '%s%03d' %(self.PREFIX, self.id)
699 697 beacon_file = os.path.join(path,'%s.txt'%self.name)
700 698 self.filename_phase = os.path.join(figpath,beacon_file)
701 699 #self.save_phase(self.filename_phase)
702 700
703 701
704 702 #store data beacon phase
705 703 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
706 704
707 705 self.setWinTitle(title)
708 706
709 707
710 708 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
711 709
712 710 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
713 711
714 712 axes = self.axesList[0]
715 713
716 714 self.xdata = numpy.hstack((self.xdata, x[0:1]))
717 715
718 716 if len(self.ydata)==0:
719 717 self.ydata = phase_beacon.reshape(-1,1)
720 718 else:
721 719 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
722 720
723 721
724 722 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
725 723 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
726 724 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
727 725 XAxisAsTime=True, grid='both'
728 726 )
729 727
730 728 self.draw()
731 729
732 730 if dataOut.ltctime >= self.xmax:
733 731 self.counter_imagwr = wr_period
734 732 self.isConfig = False
735 733 update_figfile = True
736 734
737 735 self.save(figpath=figpath,
738 736 figfile=figfile,
739 737 save=save,
740 738 ftp=ftp,
741 739 wr_period=wr_period,
742 740 thisDatetime=thisDatetime,
743 741 update_figfile=update_figfile)
744 742
745 743 return dataOut
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
@@ -1,1627 +1,1628
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62 62 for channel in channelList:
63 63 if channel not in self.dataOut.channelList:
64 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
65 65
66 66 index = self.dataOut.channelList.index(channel)
67 67 channelIndexList.append(index)
68 68 self.selectChannelsByIndex(channelIndexList)
69 69 return self.dataOut
70 70
71 71 def selectChannelsByIndex(self, channelIndexList):
72 72 """
73 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
74 74
75 75 Input:
76 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
77 77
78 78 Affected:
79 79 self.dataOut.data
80 80 self.dataOut.channelIndexList
81 81 self.dataOut.nChannels
82 82 self.dataOut.m_ProcessingHeader.totalSpectra
83 83 self.dataOut.systemHeaderObj.numChannels
84 84 self.dataOut.m_ProcessingHeader.blockSize
85 85
86 86 Return:
87 87 None
88 88 """
89 89
90 90 for channelIndex in channelIndexList:
91 91 if channelIndex not in self.dataOut.channelIndexList:
92 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
93 93
94 94 if self.dataOut.type == 'Voltage':
95 95 if self.dataOut.flagDataAsBlock:
96 96 """
97 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
98 98 """
99 99 data = self.dataOut.data[channelIndexList,:,:]
100 100 else:
101 101 data = self.dataOut.data[channelIndexList,:]
102 102
103 103 self.dataOut.data = data
104 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
105 105 self.dataOut.channelList = range(len(channelIndexList))
106 106
107 107 elif self.dataOut.type == 'Spectra':
108 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
109 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
110 110
111 111 self.dataOut.data_spc = data_spc
112 112 self.dataOut.data_dc = data_dc
113 113
114 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 self.dataOut.channelList = range(len(channelIndexList))
116 116 self.__selectPairsByChannel(channelIndexList)
117 117
118 118 return 1
119 119
120 120 def __selectPairsByChannel(self, channelList=None):
121 121
122 122 if channelList == None:
123 123 return
124 124
125 125 pairsIndexListSelected = []
126 126 for pairIndex in self.dataOut.pairsIndexList:
127 127 # First pair
128 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
129 129 continue
130 130 # Second pair
131 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
132 132 continue
133 133
134 134 pairsIndexListSelected.append(pairIndex)
135 135
136 136 if not pairsIndexListSelected:
137 137 self.dataOut.data_cspc = None
138 138 self.dataOut.pairsList = []
139 139 return
140 140
141 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
142 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
143 143 for i in pairsIndexListSelected]
144 144
145 145 return
146 146
147 147 class selectHeights(Operation):
148 148
149 149 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
150 150 """
151 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
152 152 minHei <= height <= maxHei
153 153
154 154 Input:
155 155 minHei : valor minimo de altura a considerar
156 156 maxHei : valor maximo de altura a considerar
157 157
158 158 Affected:
159 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
160 160
161 161 Return:
162 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
163 163 """
164 164
165 165 self.dataOut = dataOut
166 166
167 167 if minHei and maxHei:
168 168
169 169 if (minHei < self.dataOut.heightList[0]):
170 170 minHei = self.dataOut.heightList[0]
171 171
172 172 if (maxHei > self.dataOut.heightList[-1]):
173 173 maxHei = self.dataOut.heightList[-1]
174 174
175 175 minIndex = 0
176 176 maxIndex = 0
177 177 heights = self.dataOut.heightList
178 178
179 179 inda = numpy.where(heights >= minHei)
180 180 indb = numpy.where(heights <= maxHei)
181 181
182 182 try:
183 183 minIndex = inda[0][0]
184 184 except:
185 185 minIndex = 0
186 186
187 187 try:
188 188 maxIndex = indb[0][-1]
189 189 except:
190 190 maxIndex = len(heights)
191 191
192 192 self.selectHeightsByIndex(minIndex, maxIndex)
193 193
194 194 return self.dataOut
195 195
196 196 def selectHeightsByIndex(self, minIndex, maxIndex):
197 197 """
198 198 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
199 199 minIndex <= index <= maxIndex
200 200
201 201 Input:
202 202 minIndex : valor de indice minimo de altura a considerar
203 203 maxIndex : valor de indice maximo de altura a considerar
204 204
205 205 Affected:
206 206 self.dataOut.data
207 207 self.dataOut.heightList
208 208
209 209 Return:
210 210 1 si el metodo se ejecuto con exito caso contrario devuelve 0
211 211 """
212 212
213 213 if self.dataOut.type == 'Voltage':
214 214 if (minIndex < 0) or (minIndex > maxIndex):
215 215 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
216 216
217 217 if (maxIndex >= self.dataOut.nHeights):
218 218 maxIndex = self.dataOut.nHeights
219 219
220 220 #voltage
221 221 if self.dataOut.flagDataAsBlock:
222 222 """
223 223 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
224 224 """
225 225 data = self.dataOut.data[:,:, minIndex:maxIndex]
226 226 else:
227 227 data = self.dataOut.data[:, minIndex:maxIndex]
228 228
229 229 # firstHeight = self.dataOut.heightList[minIndex]
230 230
231 231 self.dataOut.data = data
232 232 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
233 233
234 234 if self.dataOut.nHeights <= 1:
235 235 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
236 236 elif self.dataOut.type == 'Spectra':
237 237 if (minIndex < 0) or (minIndex > maxIndex):
238 238 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
239 239 minIndex, maxIndex))
240 240
241 241 if (maxIndex >= self.dataOut.nHeights):
242 242 maxIndex = self.dataOut.nHeights - 1
243 243
244 244 # Spectra
245 245 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
246 246
247 247 data_cspc = None
248 248 if self.dataOut.data_cspc is not None:
249 249 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
250 250
251 251 data_dc = None
252 252 if self.dataOut.data_dc is not None:
253 253 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
254 254
255 255 self.dataOut.data_spc = data_spc
256 256 self.dataOut.data_cspc = data_cspc
257 257 self.dataOut.data_dc = data_dc
258 258
259 259 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
260 260
261 261 return 1
262 262
263 263
264 264 class filterByHeights(Operation):
265 265
266 266 def run(self, dataOut, window):
267 267
268 268 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
269 269
270 270 if window == None:
271 271 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
272 272
273 273 newdelta = deltaHeight * window
274 274 r = dataOut.nHeights % window
275 275 newheights = (dataOut.nHeights-r)/window
276 276
277 277 if newheights <= 1:
278 278 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
279 279
280 280 if dataOut.flagDataAsBlock:
281 281 """
282 282 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
283 283 """
284 284 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
285 285 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
286 286 buffer = numpy.sum(buffer,3)
287 287
288 288 else:
289 289 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
290 290 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
291 291 buffer = numpy.sum(buffer,2)
292 292
293 293 dataOut.data = buffer
294 294 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
295 295 dataOut.windowOfFilter = window
296 296
297 297 return dataOut
298 298
299 299
300 300 class setH0(Operation):
301 301
302 302 def run(self, dataOut, h0, deltaHeight = None):
303 303
304 304 if not deltaHeight:
305 305 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
306 306
307 307 nHeights = dataOut.nHeights
308 308
309 309 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
310 310
311 311 dataOut.heightList = newHeiRange
312 312
313 313 return dataOut
314 314
315 315
316 316 class deFlip(Operation):
317 317
318 318 def run(self, dataOut, channelList = []):
319 319
320 320 data = dataOut.data.copy()
321 321
322 322 if dataOut.flagDataAsBlock:
323 323 flip = self.flip
324 324 profileList = list(range(dataOut.nProfiles))
325 325
326 326 if not channelList:
327 327 for thisProfile in profileList:
328 328 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
329 329 flip *= -1.0
330 330 else:
331 331 for thisChannel in channelList:
332 332 if thisChannel not in dataOut.channelList:
333 333 continue
334 334
335 335 for thisProfile in profileList:
336 336 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
337 337 flip *= -1.0
338 338
339 339 self.flip = flip
340 340
341 341 else:
342 342 if not channelList:
343 343 data[:,:] = data[:,:]*self.flip
344 344 else:
345 345 for thisChannel in channelList:
346 346 if thisChannel not in dataOut.channelList:
347 347 continue
348 348
349 349 data[thisChannel,:] = data[thisChannel,:]*self.flip
350 350
351 351 self.flip *= -1.
352 352
353 353 dataOut.data = data
354 354
355 355 return dataOut
356 356
357 357
358 358 class setAttribute(Operation):
359 359 '''
360 360 Set an arbitrary attribute(s) to dataOut
361 361 '''
362 362
363 363 def __init__(self):
364 364
365 365 Operation.__init__(self)
366 366 self._ready = False
367 367
368 368 def run(self, dataOut, **kwargs):
369 369
370 370 for key, value in kwargs.items():
371 371 setattr(dataOut, key, value)
372 372
373 373 return dataOut
374 374
375 375
376 376 @MPDecorator
377 377 class printAttribute(Operation):
378 378 '''
379 379 Print an arbitrary attribute of dataOut
380 380 '''
381 381
382 382 def __init__(self):
383 383
384 384 Operation.__init__(self)
385 385
386 386 def run(self, dataOut, attributes):
387 387
388 388 if isinstance(attributes, str):
389 389 attributes = [attributes]
390 390 for attr in attributes:
391 391 if hasattr(dataOut, attr):
392 392 log.log(getattr(dataOut, attr), attr)
393 393
394 394
395 395 class interpolateHeights(Operation):
396 396
397 397 def run(self, dataOut, topLim, botLim):
398 398 #69 al 72 para julia
399 399 #82-84 para meteoros
400 400 if len(numpy.shape(dataOut.data))==2:
401 401 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
402 402 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
403 403 #dataOut.data[:,botLim:limSup+1] = sampInterp
404 404 dataOut.data[:,botLim:topLim+1] = sampInterp
405 405 else:
406 406 nHeights = dataOut.data.shape[2]
407 407 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
408 408 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
409 409 f = interpolate.interp1d(x, y, axis = 2)
410 410 xnew = numpy.arange(botLim,topLim+1)
411 411 ynew = f(xnew)
412 412 dataOut.data[:,:,botLim:topLim+1] = ynew
413 413
414 414 return dataOut
415 415
416 416
417 417 class CohInt(Operation):
418 418
419 419 isConfig = False
420 420 __profIndex = 0
421 421 __byTime = False
422 422 __initime = None
423 423 __lastdatatime = None
424 424 __integrationtime = None
425 425 __buffer = None
426 426 __bufferStride = []
427 427 __dataReady = False
428 428 __profIndexStride = 0
429 429 __dataToPutStride = False
430 430 n = None
431 431
432 432 def __init__(self, **kwargs):
433 433
434 434 Operation.__init__(self, **kwargs)
435 435
436 436 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
437 437 """
438 438 Set the parameters of the integration class.
439 439
440 440 Inputs:
441 441
442 442 n : Number of coherent integrations
443 443 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
444 444 overlapping :
445 445 """
446 446
447 447 self.__initime = None
448 448 self.__lastdatatime = 0
449 449 self.__buffer = None
450 450 self.__dataReady = False
451 451 self.byblock = byblock
452 452 self.stride = stride
453 453
454 454 if n == None and timeInterval == None:
455 455 raise ValueError("n or timeInterval should be specified ...")
456 456
457 457 if n != None:
458 458 self.n = n
459 459 self.__byTime = False
460 460 else:
461 461 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
462 462 self.n = 9999
463 463 self.__byTime = True
464 464
465 465 if overlapping:
466 466 self.__withOverlapping = True
467 467 self.__buffer = None
468 468 else:
469 469 self.__withOverlapping = False
470 470 self.__buffer = 0
471 471
472 472 self.__profIndex = 0
473 473
474 474 def putData(self, data):
475 475
476 476 """
477 477 Add a profile to the __buffer and increase in one the __profileIndex
478 478
479 479 """
480 480
481 481 if not self.__withOverlapping:
482 482 self.__buffer += data.copy()
483 483 self.__profIndex += 1
484 484 return
485 485
486 486 #Overlapping data
487 487 nChannels, nHeis = data.shape
488 488 data = numpy.reshape(data, (1, nChannels, nHeis))
489 489
490 490 #If the buffer is empty then it takes the data value
491 491 if self.__buffer is None:
492 492 self.__buffer = data
493 493 self.__profIndex += 1
494 494 return
495 495
496 496 #If the buffer length is lower than n then stakcing the data value
497 497 if self.__profIndex < self.n:
498 498 self.__buffer = numpy.vstack((self.__buffer, data))
499 499 self.__profIndex += 1
500 500 return
501 501
502 502 #If the buffer length is equal to n then replacing the last buffer value with the data value
503 503 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
504 504 self.__buffer[self.n-1] = data
505 505 self.__profIndex = self.n
506 506 return
507 507
508 508
509 509 def pushData(self):
510 510 """
511 511 Return the sum of the last profiles and the profiles used in the sum.
512 512
513 513 Affected:
514 514
515 515 self.__profileIndex
516 516
517 517 """
518 518
519 519 if not self.__withOverlapping:
520 520 data = self.__buffer
521 521 n = self.__profIndex
522 522
523 523 self.__buffer = 0
524 524 self.__profIndex = 0
525 525
526 526 return data, n
527 527
528 528 #Integration with Overlapping
529 529 data = numpy.sum(self.__buffer, axis=0)
530 530 # print data
531 531 # raise
532 532 n = self.__profIndex
533 533
534 534 return data, n
535 535
536 536 def byProfiles(self, data):
537 537
538 538 self.__dataReady = False
539 539 avgdata = None
540 540 # n = None
541 541 # print data
542 542 # raise
543 543 self.putData(data)
544 544
545 545 if self.__profIndex == self.n:
546 546 avgdata, n = self.pushData()
547 547 self.__dataReady = True
548 548
549 549 return avgdata
550 550
551 551 def byTime(self, data, datatime):
552 552
553 553 self.__dataReady = False
554 554 avgdata = None
555 555 n = None
556 556
557 557 self.putData(data)
558 558
559 559 if (datatime - self.__initime) >= self.__integrationtime:
560 560 avgdata, n = self.pushData()
561 561 self.n = n
562 562 self.__dataReady = True
563 563
564 564 return avgdata
565 565
566 566 def integrateByStride(self, data, datatime):
567 567 # print data
568 568 if self.__profIndex == 0:
569 569 self.__buffer = [[data.copy(), datatime]]
570 570 else:
571 571 self.__buffer.append([data.copy(),datatime])
572 572 self.__profIndex += 1
573 573 self.__dataReady = False
574 574
575 575 if self.__profIndex == self.n * self.stride :
576 576 self.__dataToPutStride = True
577 577 self.__profIndexStride = 0
578 578 self.__profIndex = 0
579 579 self.__bufferStride = []
580 580 for i in range(self.stride):
581 581 current = self.__buffer[i::self.stride]
582 582 data = numpy.sum([t[0] for t in current], axis=0)
583 583 avgdatatime = numpy.average([t[1] for t in current])
584 584 # print data
585 585 self.__bufferStride.append((data, avgdatatime))
586 586
587 587 if self.__dataToPutStride:
588 588 self.__dataReady = True
589 589 self.__profIndexStride += 1
590 590 if self.__profIndexStride == self.stride:
591 591 self.__dataToPutStride = False
592 592 # print self.__bufferStride[self.__profIndexStride - 1]
593 593 # raise
594 594 return self.__bufferStride[self.__profIndexStride - 1]
595 595
596 596
597 597 return None, None
598 598
599 599 def integrate(self, data, datatime=None):
600 600
601 601 if self.__initime == None:
602 602 self.__initime = datatime
603 603
604 604 if self.__byTime:
605 605 avgdata = self.byTime(data, datatime)
606 606 else:
607 607 avgdata = self.byProfiles(data)
608 608
609 609
610 610 self.__lastdatatime = datatime
611 611
612 612 if avgdata is None:
613 613 return None, None
614 614
615 615 avgdatatime = self.__initime
616 616
617 617 deltatime = datatime - self.__lastdatatime
618 618
619 619 if not self.__withOverlapping:
620 620 self.__initime = datatime
621 621 else:
622 622 self.__initime += deltatime
623 623
624 624 return avgdata, avgdatatime
625 625
626 626 def integrateByBlock(self, dataOut):
627 627
628 628 times = int(dataOut.data.shape[1]/self.n)
629 629 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
630 630
631 631 id_min = 0
632 632 id_max = self.n
633 633
634 634 for i in range(times):
635 635 junk = dataOut.data[:,id_min:id_max,:]
636 636 avgdata[:,i,:] = junk.sum(axis=1)
637 637 id_min += self.n
638 638 id_max += self.n
639 639
640 640 timeInterval = dataOut.ippSeconds*self.n
641 641 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
642 642 self.__dataReady = True
643 643 return avgdata, avgdatatime
644 644
645 645 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
646 646
647 647 if not self.isConfig:
648 648 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
649 649 self.isConfig = True
650 650
651 651 if dataOut.flagDataAsBlock:
652 652 """
653 653 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
654 654 """
655 655 avgdata, avgdatatime = self.integrateByBlock(dataOut)
656 656 dataOut.nProfiles /= self.n
657 657 else:
658 658 if stride is None:
659 659 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
660 660 else:
661 661 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
662 662
663 663
664 664 # dataOut.timeInterval *= n
665 665 dataOut.flagNoData = True
666 666
667 667 if self.__dataReady:
668 668 dataOut.data = avgdata
669 669 if not dataOut.flagCohInt:
670 670 dataOut.nCohInt *= self.n
671 671 dataOut.flagCohInt = True
672 672 dataOut.utctime = avgdatatime
673 673 # print avgdata, avgdatatime
674 674 # raise
675 675 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
676 676 dataOut.flagNoData = False
677 677 return dataOut
678 678
679 679 class Decoder(Operation):
680 680
681 681 isConfig = False
682 682 __profIndex = 0
683 683
684 684 code = None
685 685
686 686 nCode = None
687 687 nBaud = None
688 688
689 689 def __init__(self, **kwargs):
690 690
691 691 Operation.__init__(self, **kwargs)
692 692
693 693 self.times = None
694 694 self.osamp = None
695 695 # self.__setValues = False
696 696 self.isConfig = False
697 697 self.setupReq = False
698 698 def setup(self, code, osamp, dataOut):
699 699
700 700 self.__profIndex = 0
701 701
702 702 self.code = code
703 703
704 704 self.nCode = len(code)
705 705 self.nBaud = len(code[0])
706 706
707 707 if (osamp != None) and (osamp >1):
708 708 self.osamp = osamp
709 709 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
710 710 self.nBaud = self.nBaud*self.osamp
711 711
712 712 self.__nChannels = dataOut.nChannels
713 713 self.__nProfiles = dataOut.nProfiles
714 714 self.__nHeis = dataOut.nHeights
715 715
716 716 if self.__nHeis < self.nBaud:
717 717 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
718 718
719 719 #Frequency
720 720 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
721 721
722 722 __codeBuffer[:,0:self.nBaud] = self.code
723 723
724 724 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
725 725
726 726 if dataOut.flagDataAsBlock:
727 727
728 728 self.ndatadec = self.__nHeis #- self.nBaud + 1
729 729
730 730 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
731 731
732 732 else:
733 733
734 734 #Time
735 735 self.ndatadec = self.__nHeis #- self.nBaud + 1
736 736
737 737 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
738 738
739 739 def __convolutionInFreq(self, data):
740 740
741 741 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
742 742
743 743 fft_data = numpy.fft.fft(data, axis=1)
744 744
745 745 conv = fft_data*fft_code
746 746
747 747 data = numpy.fft.ifft(conv,axis=1)
748 748
749 749 return data
750 750
751 751 def __convolutionInFreqOpt(self, data):
752 752
753 753 raise NotImplementedError
754 754
755 755 def __convolutionInTime(self, data):
756 756
757 757 code = self.code[self.__profIndex]
758 758 for i in range(self.__nChannels):
759 759 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
760 760
761 761 return self.datadecTime
762 762
763 763 def __convolutionByBlockInTime(self, data):
764 764
765 765 repetitions = int(self.__nProfiles / self.nCode)
766 766 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
767 767 junk = junk.flatten()
768 768 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
769 769 profilesList = range(self.__nProfiles)
770 770
771 771 for i in range(self.__nChannels):
772 772 for j in profilesList:
773 773 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
774 774 return self.datadecTime
775 775
776 776 def __convolutionByBlockInFreq(self, data):
777 777
778 778 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
779 779
780 780
781 781 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
782 782
783 783 fft_data = numpy.fft.fft(data, axis=2)
784 784
785 785 conv = fft_data*fft_code
786 786
787 787 data = numpy.fft.ifft(conv,axis=2)
788 788
789 789 return data
790 790
791 791
792 792 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
793 793
794 794 if dataOut.flagDecodeData:
795 795 print("This data is already decoded, recoding again ...")
796 796
797 797 if not self.isConfig:
798 798
799 799 if code is None:
800 800 if dataOut.code is None:
801 801 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
802 802
803 803 code = dataOut.code
804 804 else:
805 805 code = numpy.array(code).reshape(nCode,nBaud)
806 806 self.setup(code, osamp, dataOut)
807 807
808 808 self.isConfig = True
809 809
810 810 if mode == 3:
811 811 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
812 812
813 813 if times != None:
814 814 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
815 815
816 816 if self.code is None:
817 817 print("Fail decoding: Code is not defined.")
818 818 return
819 819
820 820 self.__nProfiles = dataOut.nProfiles
821 821 datadec = None
822 822
823 823 if mode == 3:
824 824 mode = 0
825 825
826 826 if dataOut.flagDataAsBlock:
827 827 """
828 828 Decoding when data have been read as block,
829 829 """
830 830
831 831 if mode == 0:
832 832 datadec = self.__convolutionByBlockInTime(dataOut.data)
833 833 if mode == 1:
834 834 datadec = self.__convolutionByBlockInFreq(dataOut.data)
835 835 else:
836 836 """
837 837 Decoding when data have been read profile by profile
838 838 """
839 839 if mode == 0:
840 840 datadec = self.__convolutionInTime(dataOut.data)
841 841
842 842 if mode == 1:
843 843 datadec = self.__convolutionInFreq(dataOut.data)
844 844
845 845 if mode == 2:
846 846 datadec = self.__convolutionInFreqOpt(dataOut.data)
847 847
848 848 if datadec is None:
849 849 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
850 850
851 851 dataOut.code = self.code
852 852 dataOut.nCode = self.nCode
853 853 dataOut.nBaud = self.nBaud
854 854
855 855 dataOut.data = datadec
856 856
857 857 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
858 858
859 859 dataOut.flagDecodeData = True #asumo q la data esta decodificada
860 860
861 861 if self.__profIndex == self.nCode-1:
862 862 self.__profIndex = 0
863 863 return dataOut
864 864
865 865 self.__profIndex += 1
866 866
867 867 return dataOut
868 868 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
869 869
870 870
871 871 class ProfileConcat(Operation):
872 872
873 873 isConfig = False
874 874 buffer = None
875 875
876 876 def __init__(self, **kwargs):
877 877
878 878 Operation.__init__(self, **kwargs)
879 879 self.profileIndex = 0
880 880
881 881 def reset(self):
882 882 self.buffer = numpy.zeros_like(self.buffer)
883 883 self.start_index = 0
884 884 self.times = 1
885 885
886 886 def setup(self, data, m, n=1):
887 887 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
888 888 self.nHeights = data.shape[1]#.nHeights
889 889 self.start_index = 0
890 890 self.times = 1
891 891
892 892 def concat(self, data):
893 893
894 894 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
895 895 self.start_index = self.start_index + self.nHeights
896 896
897 897 def run(self, dataOut, m):
898 898 dataOut.flagNoData = True
899 899
900 900 if not self.isConfig:
901 901 self.setup(dataOut.data, m, 1)
902 902 self.isConfig = True
903 903
904 904 if dataOut.flagDataAsBlock:
905 905 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
906 906
907 907 else:
908 908 self.concat(dataOut.data)
909 909 self.times += 1
910 910 if self.times > m:
911 911 dataOut.data = self.buffer
912 912 self.reset()
913 913 dataOut.flagNoData = False
914 914 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
915 915 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
916 916 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
917 917 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
918 918 dataOut.ippSeconds *= m
919 919 return dataOut
920 920
921 921 class ProfileSelector(Operation):
922 922
923 923 profileIndex = None
924 924 # Tamanho total de los perfiles
925 925 nProfiles = None
926 926
927 927 def __init__(self, **kwargs):
928 928
929 929 Operation.__init__(self, **kwargs)
930 930 self.profileIndex = 0
931 931
932 932 def incProfileIndex(self):
933 933
934 934 self.profileIndex += 1
935 935
936 936 if self.profileIndex >= self.nProfiles:
937 937 self.profileIndex = 0
938 938
939 939 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
940 940
941 941 if profileIndex < minIndex:
942 942 return False
943 943
944 944 if profileIndex > maxIndex:
945 945 return False
946 946
947 947 return True
948 948
949 949 def isThisProfileInList(self, profileIndex, profileList):
950 950
951 951 if profileIndex not in profileList:
952 952 return False
953 953
954 954 return True
955 955
956 956 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
957 957
958 958 """
959 959 ProfileSelector:
960 960
961 961 Inputs:
962 962 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
963 963
964 964 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
965 965
966 966 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
967 967
968 968 """
969 969
970 970 if rangeList is not None:
971 971 if type(rangeList[0]) not in (tuple, list):
972 972 rangeList = [rangeList]
973 973
974 974 dataOut.flagNoData = True
975 975
976 976 if dataOut.flagDataAsBlock:
977 977 """
978 978 data dimension = [nChannels, nProfiles, nHeis]
979 979 """
980 980 if profileList != None:
981 981 dataOut.data = dataOut.data[:,profileList,:]
982 982
983 983 if profileRangeList != None:
984 984 minIndex = profileRangeList[0]
985 985 maxIndex = profileRangeList[1]
986 986 profileList = list(range(minIndex, maxIndex+1))
987 987
988 988 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
989 989
990 990 if rangeList != None:
991 991
992 992 profileList = []
993 993
994 994 for thisRange in rangeList:
995 995 minIndex = thisRange[0]
996 996 maxIndex = thisRange[1]
997 997
998 998 profileList.extend(list(range(minIndex, maxIndex+1)))
999 999
1000 1000 dataOut.data = dataOut.data[:,profileList,:]
1001 1001
1002 1002 dataOut.nProfiles = len(profileList)
1003 1003 dataOut.profileIndex = dataOut.nProfiles - 1
1004 1004 dataOut.flagNoData = False
1005 1005
1006 1006 return dataOut
1007 1007
1008 1008 """
1009 1009 data dimension = [nChannels, nHeis]
1010 1010 """
1011 1011
1012 1012 if profileList != None:
1013 1013
1014 1014 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1015 1015
1016 1016 self.nProfiles = len(profileList)
1017 1017 dataOut.nProfiles = self.nProfiles
1018 1018 dataOut.profileIndex = self.profileIndex
1019 1019 dataOut.flagNoData = False
1020 1020
1021 1021 self.incProfileIndex()
1022 1022 return dataOut
1023 1023
1024 1024 if profileRangeList != None:
1025 1025
1026 1026 minIndex = profileRangeList[0]
1027 1027 maxIndex = profileRangeList[1]
1028 1028
1029 1029 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1030 1030
1031 1031 self.nProfiles = maxIndex - minIndex + 1
1032 1032 dataOut.nProfiles = self.nProfiles
1033 1033 dataOut.profileIndex = self.profileIndex
1034 1034 dataOut.flagNoData = False
1035 1035
1036 1036 self.incProfileIndex()
1037 1037 return dataOut
1038 1038
1039 1039 if rangeList != None:
1040 1040
1041 1041 nProfiles = 0
1042 1042
1043 1043 for thisRange in rangeList:
1044 1044 minIndex = thisRange[0]
1045 1045 maxIndex = thisRange[1]
1046 1046
1047 1047 nProfiles += maxIndex - minIndex + 1
1048 1048
1049 1049 for thisRange in rangeList:
1050 1050
1051 1051 minIndex = thisRange[0]
1052 1052 maxIndex = thisRange[1]
1053 1053
1054 1054 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1055 1055
1056 1056 self.nProfiles = nProfiles
1057 1057 dataOut.nProfiles = self.nProfiles
1058 1058 dataOut.profileIndex = self.profileIndex
1059 1059 dataOut.flagNoData = False
1060 1060
1061 1061 self.incProfileIndex()
1062 1062
1063 1063 break
1064 1064
1065 1065 return dataOut
1066 1066
1067 1067
1068 1068 if beam != None: #beam is only for AMISR data
1069 1069 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1070 1070 dataOut.flagNoData = False
1071 1071 dataOut.profileIndex = self.profileIndex
1072 1072
1073 1073 self.incProfileIndex()
1074 1074
1075 1075 return dataOut
1076 1076
1077 1077 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1078 1078
1079 1079
1080 1080 class Reshaper(Operation):
1081 1081
1082 1082 def __init__(self, **kwargs):
1083 1083
1084 1084 Operation.__init__(self, **kwargs)
1085 1085
1086 1086 self.__buffer = None
1087 1087 self.__nitems = 0
1088 1088
1089 1089 def __appendProfile(self, dataOut, nTxs):
1090 1090
1091 1091 if self.__buffer is None:
1092 1092 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1093 1093 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1094 1094
1095 1095 ini = dataOut.nHeights * self.__nitems
1096 1096 end = ini + dataOut.nHeights
1097 1097
1098 1098 self.__buffer[:, ini:end] = dataOut.data
1099 1099
1100 1100 self.__nitems += 1
1101 1101
1102 1102 return int(self.__nitems*nTxs)
1103 1103
1104 1104 def __getBuffer(self):
1105 1105
1106 1106 if self.__nitems == int(1./self.__nTxs):
1107 1107
1108 1108 self.__nitems = 0
1109 1109
1110 1110 return self.__buffer.copy()
1111 1111
1112 1112 return None
1113 1113
1114 1114 def __checkInputs(self, dataOut, shape, nTxs):
1115 1115
1116 1116 if shape is None and nTxs is None:
1117 1117 raise ValueError("Reshaper: shape of factor should be defined")
1118 1118
1119 1119 if nTxs:
1120 1120 if nTxs < 0:
1121 1121 raise ValueError("nTxs should be greater than 0")
1122 1122
1123 1123 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1124 1124 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1125 1125
1126 1126 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1127 1127
1128 1128 return shape, nTxs
1129 1129
1130 1130 if len(shape) != 2 and len(shape) != 3:
1131 1131 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1132 1132
1133 1133 if len(shape) == 2:
1134 1134 shape_tuple = [dataOut.nChannels]
1135 1135 shape_tuple.extend(shape)
1136 1136 else:
1137 1137 shape_tuple = list(shape)
1138 1138
1139 1139 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1140 1140
1141 1141 return shape_tuple, nTxs
1142 1142
1143 1143 def run(self, dataOut, shape=None, nTxs=None):
1144 1144
1145 1145 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1146 1146
1147 1147 dataOut.flagNoData = True
1148 1148 profileIndex = None
1149 1149
1150 1150 if dataOut.flagDataAsBlock:
1151 1151
1152 1152 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1153 1153 dataOut.flagNoData = False
1154 1154
1155 1155 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1156 1156
1157 1157 else:
1158 1158
1159 1159 if self.__nTxs < 1:
1160 1160
1161 1161 self.__appendProfile(dataOut, self.__nTxs)
1162 1162 new_data = self.__getBuffer()
1163 1163
1164 1164 if new_data is not None:
1165 1165 dataOut.data = new_data
1166 1166 dataOut.flagNoData = False
1167 1167
1168 1168 profileIndex = dataOut.profileIndex*nTxs
1169 1169
1170 1170 else:
1171 1171 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1172 1172
1173 1173 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1174 1174
1175 1175 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1176 1176
1177 1177 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1178 1178
1179 1179 dataOut.profileIndex = profileIndex
1180 1180
1181 1181 dataOut.ippSeconds /= self.__nTxs
1182 1182
1183 1183 return dataOut
1184 1184
1185 1185 class SplitProfiles(Operation):
1186 1186
1187 1187 def __init__(self, **kwargs):
1188 1188
1189 1189 Operation.__init__(self, **kwargs)
1190 1190
1191 1191 def run(self, dataOut, n):
1192 1192
1193 1193 dataOut.flagNoData = True
1194 1194 profileIndex = None
1195 1195
1196 1196 if dataOut.flagDataAsBlock:
1197 1197
1198 1198 #nchannels, nprofiles, nsamples
1199 1199 shape = dataOut.data.shape
1200 1200
1201 1201 if shape[2] % n != 0:
1202 1202 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1203 1203
1204 1204 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1205 1205
1206 1206 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1207 1207 dataOut.flagNoData = False
1208 1208
1209 1209 profileIndex = int(dataOut.nProfiles/n) - 1
1210 1210
1211 1211 else:
1212 1212
1213 1213 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1214 1214
1215 1215 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1216 1216
1217 1217 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1218 1218
1219 1219 dataOut.nProfiles = int(dataOut.nProfiles*n)
1220 1220
1221 1221 dataOut.profileIndex = profileIndex
1222 1222
1223 1223 dataOut.ippSeconds /= n
1224 1224
1225 1225 return dataOut
1226 1226
1227 1227 class CombineProfiles(Operation):
1228 1228 def __init__(self, **kwargs):
1229 1229
1230 1230 Operation.__init__(self, **kwargs)
1231 1231
1232 1232 self.__remData = None
1233 1233 self.__profileIndex = 0
1234 1234
1235 1235 def run(self, dataOut, n):
1236 1236
1237 1237 dataOut.flagNoData = True
1238 1238 profileIndex = None
1239 1239
1240 1240 if dataOut.flagDataAsBlock:
1241 1241
1242 1242 #nchannels, nprofiles, nsamples
1243 1243 shape = dataOut.data.shape
1244 1244 new_shape = shape[0], shape[1]/n, shape[2]*n
1245 1245
1246 1246 if shape[1] % n != 0:
1247 1247 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1248 1248
1249 1249 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1250 1250 dataOut.flagNoData = False
1251 1251
1252 1252 profileIndex = int(dataOut.nProfiles*n) - 1
1253 1253
1254 1254 else:
1255 1255
1256 1256 #nchannels, nsamples
1257 1257 if self.__remData is None:
1258 1258 newData = dataOut.data
1259 1259 else:
1260 1260 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1261 1261
1262 1262 self.__profileIndex += 1
1263 1263
1264 1264 if self.__profileIndex < n:
1265 1265 self.__remData = newData
1266 1266 #continue
1267 1267 return
1268 1268
1269 1269 self.__profileIndex = 0
1270 1270 self.__remData = None
1271 1271
1272 1272 dataOut.data = newData
1273 1273 dataOut.flagNoData = False
1274 1274
1275 1275 profileIndex = dataOut.profileIndex/n
1276 1276
1277 1277
1278 1278 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1279 1279
1280 1280 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1281 1281
1282 1282 dataOut.nProfiles = int(dataOut.nProfiles/n)
1283 1283
1284 1284 dataOut.profileIndex = profileIndex
1285 1285
1286 1286 dataOut.ippSeconds *= n
1287 1287
1288 1288 return dataOut
1289 1289
1290 1290 class PulsePair(Operation):
1291 1291 '''
1292 1292 Function PulsePair(Signal Power, Velocity)
1293 1293 The real component of Lag[0] provides Intensity Information
1294 1294 The imag component of Lag[1] Phase provides Velocity Information
1295 1295
1296 1296 Configuration Parameters:
1297 1297 nPRF = Number of Several PRF
1298 1298 theta = Degree Azimuth angel Boundaries
1299 1299
1300 1300 Input:
1301 1301 self.dataOut
1302 1302 lag[N]
1303 1303 Affected:
1304 1304 self.dataOut.spc
1305 1305 '''
1306 1306 isConfig = False
1307 1307 __profIndex = 0
1308 1308 __initime = None
1309 1309 __lastdatatime = None
1310 1310 __buffer = None
1311 1311 noise = None
1312 1312 __dataReady = False
1313 1313 n = None
1314 1314 __nch = 0
1315 1315 __nHeis = 0
1316 1316 removeDC = False
1317 1317 ipp = None
1318 1318 lambda_ = 0
1319 1319
1320 1320 def __init__(self,**kwargs):
1321 1321 Operation.__init__(self,**kwargs)
1322 1322
1323 1323 def setup(self, dataOut, n = None, removeDC=False):
1324 1324 '''
1325 1325 n= Numero de PRF's de entrada
1326 1326 '''
1327 1327 print("[INICIO]-setup del METODO PULSE PAIR")
1328 1328 self.__initime = None
1329 1329 self.__lastdatatime = 0
1330 1330 self.__dataReady = False
1331 1331 self.__buffer = 0
1332 1332 self.__profIndex = 0
1333 1333 self.noise = None
1334 1334 self.__nch = dataOut.nChannels
1335 1335 self.__nHeis = dataOut.nHeights
1336 1336 self.removeDC = removeDC
1337 1337 self.lambda_ = 3.0e8/(9345.0e6)
1338 1338 self.ippSec = dataOut.ippSeconds
1339 1339 self.nCohInt = dataOut.nCohInt
1340 1340 print("IPPseconds",dataOut.ippSeconds)
1341 1341
1342 1342 print("ELVALOR DE n es:", n)
1343 1343 if n == None:
1344 1344 raise ValueError("n should be specified.")
1345 1345
1346 1346 if n != None:
1347 1347 if n<2:
1348 1348 raise ValueError("n should be greater than 2")
1349 1349
1350 1350 self.n = n
1351 1351 self.__nProf = n
1352 1352
1353 1353 self.__buffer = numpy.zeros((dataOut.nChannels,
1354 1354 n,
1355 1355 dataOut.nHeights),
1356 1356 dtype='complex')
1357 1357
1358 1358 def putData(self,data):
1359 1359 '''
1360 1360 Add a profile to he __buffer and increase in one the __profiel Index
1361 1361 '''
1362 1362 self.__buffer[:,self.__profIndex,:]= data
1363 1363 self.__profIndex += 1
1364 1364 return
1365 1365
1366 1366 def pushData(self,dataOut):
1367 1367 '''
1368 1368 Return the PULSEPAIR and the profiles used in the operation
1369 1369 Affected : self.__profileIndex
1370 1370 '''
1371 1371 #----------------- Remove DC-----------------------------------
1372 1372 if self.removeDC==True:
1373 1373 mean = numpy.mean(self.__buffer,1)
1374 1374 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1375 1375 dc= numpy.tile(tmp,[1,self.__nProf,1])
1376 1376 self.__buffer = self.__buffer - dc
1377 1377 #------------------Calculo de Potencia ------------------------
1378 1378 pair0 = self.__buffer*numpy.conj(self.__buffer)
1379 1379 pair0 = pair0.real
1380 1380 lag_0 = numpy.sum(pair0,1)
1381 1381 #------------------Calculo de Ruido x canal--------------------
1382 1382 self.noise = numpy.zeros(self.__nch)
1383 1383 for i in range(self.__nch):
1384 1384 daux = numpy.sort(pair0[i,:,:],axis= None)
1385 1385 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1386 1386
1387 1387 self.noise = self.noise.reshape(self.__nch,1)
1388 1388 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1389 1389 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1390 1390 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1391 1391 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1392 1392 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1393 1393 #-------------------- Power --------------------------------------------------
1394 1394 data_power = lag_0/(self.n*self.nCohInt)
1395 1395 #------------------ Senal ---------------------------------------------------
1396 1396 data_intensity = pair0 - noise_buffer
1397 1397 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1398 1398 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1399 1399 for i in range(self.__nch):
1400 1400 for j in range(self.__nHeis):
1401 1401 if data_intensity[i][j] < 0:
1402 1402 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1403 1403
1404 1404 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1405 1405 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1406 1406 lag_1 = numpy.sum(pair1,1)
1407 1407 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1408 1408 data_velocity = (self.lambda_/2.0)*data_freq
1409 1409
1410 1410 #---------------- Potencia promedio estimada de la Senal-----------
1411 1411 lag_0 = lag_0/self.n
1412 1412 S = lag_0-self.noise
1413 1413
1414 1414 #---------------- Frecuencia Doppler promedio ---------------------
1415 1415 lag_1 = lag_1/(self.n-1)
1416 1416 R1 = numpy.abs(lag_1)
1417 1417
1418 1418 #---------------- Calculo del SNR----------------------------------
1419 1419 data_snrPP = S/self.noise
1420 1420 for i in range(self.__nch):
1421 1421 for j in range(self.__nHeis):
1422 1422 if data_snrPP[i][j] < 1.e-20:
1423 1423 data_snrPP[i][j] = 1.e-20
1424 1424
1425 1425 #----------------- Calculo del ancho espectral ----------------------
1426 1426 L = S/R1
1427 1427 L = numpy.where(L<0,1,L)
1428 1428 L = numpy.log(L)
1429 1429 tmp = numpy.sqrt(numpy.absolute(L))
1430 1430 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1431 1431 n = self.__profIndex
1432 1432
1433 1433 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1434 1434 self.__profIndex = 0
1435 1435 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1436 1436
1437 1437
1438 1438 def pulsePairbyProfiles(self,dataOut):
1439 1439
1440 1440 self.__dataReady = False
1441 1441 data_power = None
1442 1442 data_intensity = None
1443 1443 data_velocity = None
1444 1444 data_specwidth = None
1445 1445 data_snrPP = None
1446 1446 self.putData(data=dataOut.data)
1447 1447 if self.__profIndex == self.n:
1448 1448 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1449 1449 self.__dataReady = True
1450 1450
1451 1451 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1452 1452
1453 1453
1454 1454 def pulsePairOp(self, dataOut, datatime= None):
1455 1455
1456 1456 if self.__initime == None:
1457 1457 self.__initime = datatime
1458 1458 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1459 1459 self.__lastdatatime = datatime
1460 1460
1461 1461 if data_power is None:
1462 1462 return None, None, None,None,None,None
1463 1463
1464 1464 avgdatatime = self.__initime
1465 1465 deltatime = datatime - self.__lastdatatime
1466 1466 self.__initime = datatime
1467 1467
1468 1468 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1469 1469
1470 1470 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1471 1471
1472 1472 if not self.isConfig:
1473 1473 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1474 1474 self.isConfig = True
1475 1475 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1476 1476 dataOut.flagNoData = True
1477 1477
1478 1478 if self.__dataReady:
1479 1479 dataOut.nCohInt *= self.n
1480 1480 dataOut.dataPP_POW = data_intensity # S
1481 print("help",data_power)
1481 1482 dataOut.dataPP_POWER = data_power # P
1482 1483 dataOut.dataPP_DOP = data_velocity
1483 1484 dataOut.dataPP_SNR = data_snrPP
1484 1485 dataOut.dataPP_WIDTH = data_specwidth
1485 1486 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1486 1487 dataOut.nProfiles = int(dataOut.nProfiles/n)
1487 1488 dataOut.utctime = avgdatatime
1488 1489 dataOut.flagNoData = False
1489 1490 return dataOut
1490 1491
1491 1492
1492 1493
1493 1494 # import collections
1494 1495 # from scipy.stats import mode
1495 1496 #
1496 1497 # class Synchronize(Operation):
1497 1498 #
1498 1499 # isConfig = False
1499 1500 # __profIndex = 0
1500 1501 #
1501 1502 # def __init__(self, **kwargs):
1502 1503 #
1503 1504 # Operation.__init__(self, **kwargs)
1504 1505 # # self.isConfig = False
1505 1506 # self.__powBuffer = None
1506 1507 # self.__startIndex = 0
1507 1508 # self.__pulseFound = False
1508 1509 #
1509 1510 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1510 1511 #
1511 1512 # #Read data
1512 1513 #
1513 1514 # powerdB = dataOut.getPower(channel = channel)
1514 1515 # noisedB = dataOut.getNoise(channel = channel)[0]
1515 1516 #
1516 1517 # self.__powBuffer.extend(powerdB.flatten())
1517 1518 #
1518 1519 # dataArray = numpy.array(self.__powBuffer)
1519 1520 #
1520 1521 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1521 1522 #
1522 1523 # maxValue = numpy.nanmax(filteredPower)
1523 1524 #
1524 1525 # if maxValue < noisedB + 10:
1525 1526 # #No se encuentra ningun pulso de transmision
1526 1527 # return None
1527 1528 #
1528 1529 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1529 1530 #
1530 1531 # if len(maxValuesIndex) < 2:
1531 1532 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1532 1533 # return None
1533 1534 #
1534 1535 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1535 1536 #
1536 1537 # #Seleccionar solo valores con un espaciamiento de nSamples
1537 1538 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1538 1539 #
1539 1540 # if len(pulseIndex) < 2:
1540 1541 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1541 1542 # return None
1542 1543 #
1543 1544 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1544 1545 #
1545 1546 # #remover senales que se distancien menos de 10 unidades o muestras
1546 1547 # #(No deberian existir IPP menor a 10 unidades)
1547 1548 #
1548 1549 # realIndex = numpy.where(spacing > 10 )[0]
1549 1550 #
1550 1551 # if len(realIndex) < 2:
1551 1552 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1552 1553 # return None
1553 1554 #
1554 1555 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1555 1556 # realPulseIndex = pulseIndex[realIndex]
1556 1557 #
1557 1558 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1558 1559 #
1559 1560 # print "IPP = %d samples" %period
1560 1561 #
1561 1562 # self.__newNSamples = dataOut.nHeights #int(period)
1562 1563 # self.__startIndex = int(realPulseIndex[0])
1563 1564 #
1564 1565 # return 1
1565 1566 #
1566 1567 #
1567 1568 # def setup(self, nSamples, nChannels, buffer_size = 4):
1568 1569 #
1569 1570 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1570 1571 # maxlen = buffer_size*nSamples)
1571 1572 #
1572 1573 # bufferList = []
1573 1574 #
1574 1575 # for i in range(nChannels):
1575 1576 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1576 1577 # maxlen = buffer_size*nSamples)
1577 1578 #
1578 1579 # bufferList.append(bufferByChannel)
1579 1580 #
1580 1581 # self.__nSamples = nSamples
1581 1582 # self.__nChannels = nChannels
1582 1583 # self.__bufferList = bufferList
1583 1584 #
1584 1585 # def run(self, dataOut, channel = 0):
1585 1586 #
1586 1587 # if not self.isConfig:
1587 1588 # nSamples = dataOut.nHeights
1588 1589 # nChannels = dataOut.nChannels
1589 1590 # self.setup(nSamples, nChannels)
1590 1591 # self.isConfig = True
1591 1592 #
1592 1593 # #Append new data to internal buffer
1593 1594 # for thisChannel in range(self.__nChannels):
1594 1595 # bufferByChannel = self.__bufferList[thisChannel]
1595 1596 # bufferByChannel.extend(dataOut.data[thisChannel])
1596 1597 #
1597 1598 # if self.__pulseFound:
1598 1599 # self.__startIndex -= self.__nSamples
1599 1600 #
1600 1601 # #Finding Tx Pulse
1601 1602 # if not self.__pulseFound:
1602 1603 # indexFound = self.__findTxPulse(dataOut, channel)
1603 1604 #
1604 1605 # if indexFound == None:
1605 1606 # dataOut.flagNoData = True
1606 1607 # return
1607 1608 #
1608 1609 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1609 1610 # self.__pulseFound = True
1610 1611 # self.__startIndex = indexFound
1611 1612 #
1612 1613 # #If pulse was found ...
1613 1614 # for thisChannel in range(self.__nChannels):
1614 1615 # bufferByChannel = self.__bufferList[thisChannel]
1615 1616 # #print self.__startIndex
1616 1617 # x = numpy.array(bufferByChannel)
1617 1618 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1618 1619 #
1619 1620 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1620 1621 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1621 1622 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1622 1623 #
1623 1624 # dataOut.data = self.__arrayBuffer
1624 1625 #
1625 1626 # self.__startIndex += self.__newNSamples
1626 1627 #
1627 1628 # return
@@ -1,157 +1,159
1 1 #!python
2 2 '''
3 3 '''
4 4
5 5 import os, sys
6 6 import datetime
7 7 import time
8 8
9 9 #path = os.path.dirname(os.getcwd())
10 10 #path = os.path.dirname(path)
11 11 #sys.path.insert(0, path)
12 12
13 13 from schainpy.controller import Project
14 14
15 15 desc = "USRP_test"
16 16 filename = "USRP_processing.xml"
17 17 controllerObj = Project()
18 18 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
19 19
20 20 ############## USED TO PLOT IQ VOLTAGE, POWER AND SPECTRA #############
21 21
22 22 #######################################################################
23 23 ######PATH DE LECTURA, ESCRITURA, GRAFICOS Y ENVIO WEB#################
24 24 #######################################################################
25 25 #path = '/media/data/data/vientos/57.2063km/echoes/NCO_Woodman'
26 26 #path = '/DATA_RM/TEST_INTEGRACION'
27 27 #path = '/DATA_RM/TEST_ONLINE'
28 28 #path_pp = '/DATA_RM/TEST_HDF5'
29 29
30 30 #figpath = '/home/soporte/Pictures/TEST_INTEGRACION_IMG'
31 31 ###path = '/DATA_RM/TEST_INTEGRACION/ADQ_OFFLINE/'
32 32 ###path_pp = '/DATA_RM/TEST_HDF5_SPEC'
33 33
34 34 #path = '/DATA_RM/USRP_22'
35 path = '/DATA_RM/23/6v'
35 ###path = '/DATA_RM/23/6v'
36 path = '/DATA_RM/TEST_19OCTUBRE/10MHZ'
36 37 #path_pp = '/DATA_RM/TEST_HDF5'
38 path_pp = '/DATA_RM/TEST_HDF5_19OCT'
37 39 # UTIMO TEST 22 DE SEPTIEMBRE
38 40 #path_pp = '/DATA_RM/TEST_HDF5_SPEC_22'
39 41 #path_pp = '/DATA_RM/TEST_HDF5_SPEC_3v'
40 path_pp = '/DATA_RM/TEST_HDF5_SPEC_23/6v'
42 ###path_pp = '/DATA_RM/TEST_HDF5_SPEC_23/6v'
41 43
42 44
43 45 #remotefolder = "/home/wmaster/graficos"
44 46 #######################################################################
45 47 ################# RANGO DE PLOTEO######################################
46 48 #######################################################################
47 49 dBmin = '-5'
48 50 dBmax = '20'
49 51 xmin = '0'
50 52 xmax ='24'
51 53 ymin = '0'
52 54 ymax = '600'
53 55 #######################################################################
54 56 ########################FECHA##########################################
55 57 #######################################################################
56 58 str = datetime.date.today()
57 59 today = str.strftime("%Y/%m/%d")
58 60 str2 = str - datetime.timedelta(days=1)
59 61 yesterday = str2.strftime("%Y/%m/%d")
60 62 #######################################################################
61 63 ######################## UNIDAD DE LECTURA#############################
62 64 #######################################################################
63 65 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
64 66 path=path,
65 67 startDate="2021/01/01",#today,
66 68 endDate="2021/12/30",#today,
67 69 startTime='00:00:00',
68 70 endTime='23:59:59',
69 71 delay=0,
70 72 #set=0,
71 73 online=0,
72 74 walk=1,
73 75 ippKm = 60)
74 76
75 77 opObj11 = readUnitConfObj.addOperation(name='printInfo')
76 78 #opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
77 79 #######################################################################
78 80 ################ OPERACIONES DOMINIO DEL TIEMPO########################
79 81 #######################################################################
80 82
81 83
82 V=6
84 V=10
83 85 IPP=400*1e-6
84 86 n= int(1/(V*IPP))
85 87 print("n numero de Perfiles a procesar con nFFTPoints ", n)
86 88
87 89 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
88 90
89 91 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
90 92 procUnitConfObjB.addParameter(name='nFFTPoints', value=n, format='int')
91 93 procUnitConfObjB.addParameter(name='nProfiles' , value=n, format='int')
92 94
93 95
94 96
95 97 #
96 98 # 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,'+\
97 99 # '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'
98 100
99 101 #opObj11 = procUnitConfObjA.addOperation(name='setRadarFrequency')
100 102 #opObj11.addParameter(name='frequency', value='70312500')
101 103 #opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
102 104 #opObj11.addParameter(name='n', value='625', format='int')#10
103 105 #opObj11.addParameter(name='removeDC', value=1, format='int')
104 106
105 107 # Ploteo TEST
106 108 '''
107 109 opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')
108 110 opObj11 = procUnitConfObjA.addOperation(name='PulsepairSignalPlot', optype='other')
109 111 opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')
110 112 #opObj11.addParameter(name='xmax', value=8)
111 113 opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')
112 114 '''
113 115 # OJO SCOPE
114 116 #opObj10 = procUnitConfObjA.addOperation(name='ScopePlot', optype='external')
115 117 #opObj10.addParameter(name='buffer_sizeid', value='10', format='int')
116 118 ##opObj10.addParameter(name='xmin', value='0', format='int')
117 119 ##opObj10.addParameter(name='xmax', value='50', format='int')
118 120 #opObj10.addParameter(name='type', value='iq')
119 121 ##opObj10.addParameter(name='ymin', value='-5000', format='int')
120 122 ##opObj10.addParameter(name='ymax', value='8500', format='int')
121 123 #opObj11.addParameter(name='save', value=figpath, format='str')
122 124 #opObj11.addParameter(name='save_period', value=10, format='int')
123 125
124 126 #opObj10 = procUnitConfObjA.addOperation(name='setH0')
125 127 #opObj10.addParameter(name='h0', value='-5000', format='float')
126 128
127 129 #opObj11 = procUnitConfObjA.addOperation(name='filterByHeights')
128 130 #opObj11.addParameter(name='window', value='1', format='int')
129 131
130 132 #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'
131 133 #opObj11 = procUnitConfObjSousy.addOperation(name='Decoder', optype='other')
132 134 #opObj11.addParameter(name='code', value=codigo, formatyesterday='floatlist')
133 135 #opObj11.addParameter(name='nCode', value='1', format='int')
134 136 #opObj11.addParameter(name='nBaud', value='28', format='int')
135 137
136 138 #opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other')
137 139 #opObj11.addParameter(name='n', value='100', format='int')
138 140
139 141 #######################################################################
140 142 ########## OPERACIONES ParametersProc########################
141 143 #######################################################################
142 144
143 145 procUnitConfObjC= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjB.getId())
144 146
145 147 procUnitConfObjC.addOperation(name='SpectralMoments')
146 148
147 149
148 150 opObj10 = procUnitConfObjC.addOperation(name='HDFWriter')
149 151 opObj10.addParameter(name='path',value=path_pp)
150 152 #opObj10.addParameter(name='mode',value=0)
151 153 opObj10.addParameter(name='blocksPerFile',value='100',format='int')
152 154 #opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
153 155 opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
154 156
155 157 opObj10.addParameter(name='dataList',value='data_pow,data_dop,utctime',format='list')#,format='list'
156 158
157 159 controllerObj.start()
@@ -1,213 +1,216
1 1 # Ing. AVP
2 2 # 06/10/2021
3 3 # ARCHIVO DE LECTURA
4 4 import os, sys
5 5 import datetime
6 6 import time
7 7 from schainpy.controller import Project
8 8 #### NOTA###########################################
9 9 # INPUT :
10 10 # VELOCIDAD PARAMETRO : V = 2Β°/seg
11 11 # MODO PULSE PAIR O MOMENTOS: 0 : Pulse Pair ,1 : Momentos
12 12 ######################################################
13 13 ##### PROCESAMIENTO ##################################
14 14 ##### OJO TENER EN CUENTA EL n= para el Pulse Pair ##
15 15 ##### O EL n= nFFTPoints ###
16 16 ######################################################
17 17 ######## BUSCAMOS EL numero de IPP equivalente 1Β°#####
18 18 ######## Sea V la velocidad del Pedestal en Β°/seg#####
19 19 ######## 1Β° sera Recorrido en un tiempo de 1/V ######
20 20 ######## IPP del Radar 400 useg --> 60 Km ############
21 21 ######## n = 1/(V(Β°/seg)*IPP(Km)) , NUMERO DE IPP ##
22 22 ######## n = 1/(V*IPP) #############################
23 23 ######## VELOCIDAD DEL PEDESTAL ######################
24 24 print("SETUP- RADAR METEOROLOGICO")
25 25 V = 10
26 26 mode = 1
27 27 #path = '/DATA_RM/23/6v'
28 path = '/DATA_RM/TEST_INTEGRACION_2M'
29 path_ped='/DATA_RM/TEST_PEDESTAL/P20211012-082745'
28 ####path = '/DATA_RM/TEST_INTEGRACION_2M'
29 #path = '/DATA_RM/TEST_19OCTUBRE/10MHZ'
30 path = '/DATA_RM/WR_20_OCT'
31 #### path_ped='/DATA_RM/TEST_PEDESTAL/P20211012-082745'
32 #### path_ped='/DATA_RM/TEST_PEDESTAL/P20211019-192244'
30 33 figpath_pp = "/home/soporte/Pictures/TEST_PP"
31 figpath_mom = "/home/soporte/Pictures/TEST_MOM"
32 plot = 0
33 integration = 1
34 figpath_spec = "/home/soporte/Pictures/TEST_MOM"
35 plot = 1
36 integration = 0
34 37 save = 0
35 38 if save == 1:
36 39 if mode==0:
37 40 path_save = '/DATA_RM/TEST_HDF5_PP_23/6v'
38 41 path_save = '/DATA_RM/TEST_HDF5_PP'
39 42 path_save = '/DATA_RM/TEST_HDF5_PP_100'
40 43 else:
41 44 path_save = '/DATA_RM/TEST_HDF5_SPEC_23_V2/6v'
42 45
43 46 print("* PATH data ADQ :", path)
44 47 print("* Velocidad Pedestal :",V,"Β°/seg")
45 48 ############################ NRO Perfiles PROCESAMIENTO ###################
46 49 V=V
47 50 IPP=400*1e-6
48 51 n= int(1/(V*IPP))
49 52 print("* n - NRO Perfiles Proc:", n )
50 53 ################################## MODE ###################################
51 54 print("* Modo de Operacion :",mode)
52 55 if mode ==0:
53 56 print("* Met. Seleccionado : Pulse Pair")
54 57 else:
55 58 print("* Met. Momentos : Momentos")
56 59
57 60 ################################## MODE ###################################
58 61 print("* Grabado de datos :",save)
59 62 if save ==1:
60 63 if mode==0:
61 64 ope= "Pulse Pair"
62 65 else:
63 66 ope= "Momentos"
64 67 print("* Path-Save Data -", ope , path_save)
65 68
66 69 print("* Integracion de datos :",integration)
67 70
68 71 time.sleep(15)
69 72 #remotefolder = "/home/wmaster/graficos"
70 73 #######################################################################
71 74 ################# RANGO DE PLOTEO######################################
72 75 dBmin = '1'
73 dBmax = '85'
74 xmin = '15'
75 xmax = '15.25'
76 dBmax = '65'
77 xmin = '13.2'
78 xmax = '13.5'
76 79 ymin = '0'
77 ymax = '600'
80 ymax = '60'
78 81 #######################################################################
79 82 ########################FECHA##########################################
80 83 str = datetime.date.today()
81 84 today = str.strftime("%Y/%m/%d")
82 85 str2 = str - datetime.timedelta(days=1)
83 86 yesterday = str2.strftime("%Y/%m/%d")
84 87 #######################################################################
85 88 ########################SIGNAL CHAIN ##################################
86 89 #######################################################################
87 90 desc = "USRP_test"
88 91 filename = "USRP_processing.xml"
89 92 controllerObj = Project()
90 93 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
91 94 #######################################################################
92 95 ######################## UNIDAD DE LECTURA#############################
93 96 #######################################################################
94 97 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
95 98 path=path,
96 99 startDate="2021/01/01",#today,
97 100 endDate="2021/12/30",#today,
98 101 startTime='00:00:00',
99 102 endTime='23:59:59',
100 103 delay=0,
101 104 #set=0,
102 105 online=0,
103 106 walk=1,
104 107 ippKm = 60)
105 108
106 109 opObj11 = readUnitConfObj.addOperation(name='printInfo')
107 110
108 111 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
109 112
110 113 if mode ==0:
111 114 ####################### METODO PULSE PAIR ######################################################################
112 115 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
113 116 opObj11.addParameter(name='n', value=int(n), format='int')#10 VOY A USAR 250 DADO QUE LA VELOCIDAD ES 10 GRADOS
114 117 #opObj11.addParameter(name='removeDC', value=1, format='int')
115 118 ####################### METODO Parametros ######################################################################
116 119 procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
117 120 if plot==1:
118 121 opObj11 = procUnitConfObjB.addOperation(name='GenericRTIPlot',optype='external')
119 opObj11.addParameter(name='attr_data', value='dataPP_POW')
122 opObj11.addParameter(name='attr_data', value='dataPP_POWER')
120 123 opObj11.addParameter(name='colormap', value='jet')
121 124 opObj11.addParameter(name='xmin', value=xmin)
122 125 opObj11.addParameter(name='xmax', value=xmax)
123 126 opObj11.addParameter(name='zmin', value=dBmin)
124 127 opObj11.addParameter(name='zmax', value=dBmax)
125 128 opObj11.addParameter(name='save', value=figpath_pp)
126 129 opObj11.addParameter(name='showprofile', value=0)
127 opObj11.addParameter(name='save_period', value=50)
130 opObj11.addParameter(name='save_period', value=10)
128 131
129 132 ####################### METODO ESCRITURA #######################################################################
130 133 if save==1:
131 134 opObj10 = procUnitConfObjB.addOperation(name='HDFWriter')
132 135 opObj10.addParameter(name='path',value=path_save)
133 136 #opObj10.addParameter(name='mode',value=0)
134 137 opObj10.addParameter(name='blocksPerFile',value='100',format='int')
135 138 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'
139 opObj10.addParameter(name='dataList',value='dataPP_POWER,dataPP_DOP,utctime',format='list')#,format='list'
137 140 if integration==1:
138 141 V=10
139 142 blocksPerfile=360
140 143 print("* Velocidad del Pedestal:",V)
141 144 tmp_blocksPerfile = 100
142 145 f_a_p= int(tmp_blocksPerfile/V)
143 146
144 147 opObj11 = procUnitConfObjB.addOperation(name='PedestalInformation')
145 148 opObj11.addParameter(name='path_ped', value=path_ped)
146 149 #opObj11.addParameter(name='path_adq', value=path_adq)
147 150 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
148 151 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
149 152 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
150 153 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
151 154 opObj11.addParameter(name='online', value='0', format='int')
152 155
153 156 opObj11 = procUnitConfObjB.addOperation(name='Block360')
154 157 opObj11.addParameter(name='n', value='10', format='int')
155 158 opObj11.addParameter(name='mode', value=mode, format='int')
156 159
157 160 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
158 161
159 162 opObj11= procUnitConfObjB.addOperation(name='WeatherPlot',optype='other')
160 163
161 164
162 165 else:
163 166 ####################### METODO SPECTROS ######################################################################
164 167 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
165 168 procUnitConfObjB.addParameter(name='nFFTPoints', value=n, format='int')
166 169 procUnitConfObjB.addParameter(name='nProfiles' , value=n, format='int')
167 170
168 171 procUnitConfObjC = controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjB.getId())
169 172 procUnitConfObjC.addOperation(name='SpectralMoments')
170 173 if plot==1:
171 174 dBmin = '1'
172 175 dBmax = '65'
173 176 opObj11 = procUnitConfObjC.addOperation(name='PowerPlot',optype='external')
174 177 opObj11.addParameter(name='xmin', value=xmin)
175 178 opObj11.addParameter(name='xmax', value=xmax)
176 179 opObj11.addParameter(name='zmin', value=dBmin)
177 180 opObj11.addParameter(name='zmax', value=dBmax)
178 opObj11.addParameter(name='save', value=figpath_mom)
181 opObj11.addParameter(name='save', value=figpath_spec)
179 182 opObj11.addParameter(name='showprofile', value=0)
180 opObj11.addParameter(name='save_period', value=100)
183 opObj11.addParameter(name='save_period', value=10)
181 184
182 185 if save==1:
183 186 opObj10 = procUnitConfObjC.addOperation(name='HDFWriter')
184 187 opObj10.addParameter(name='path',value=path_save)
185 188 #opObj10.addParameter(name='mode',value=0)
186 189 opObj10.addParameter(name='blocksPerFile',value='360',format='int')
187 190 #opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
188 191 opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
189 192 opObj10.addParameter(name='dataList',value='data_pow,data_dop,utctime',format='list')#,format='list'
190 193
191 194 if integration==1:
192 195 V=10
193 196 blocksPerfile=360
194 197 print("* Velocidad del Pedestal:",V)
195 198 tmp_blocksPerfile = 100
196 199 f_a_p= int(tmp_blocksPerfile/V)
197 200
198 201 opObj11 = procUnitConfObjC.addOperation(name='PedestalInformation')
199 202 opObj11.addParameter(name='path_ped', value=path_ped)
200 203 #opObj11.addParameter(name='path_adq', value=path_adq)
201 204 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
202 205 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
203 206 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
204 207 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
205 208 opObj11.addParameter(name='online', value='0', format='int')
206 209
207 210 opObj11 = procUnitConfObjC.addOperation(name='Block360')
208 211 opObj11.addParameter(name='n', value='10', format='int')
209 212 opObj11.addParameter(name='mode', value=mode, format='int')
210 213
211 214 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
212 215 opObj11= procUnitConfObjC.addOperation(name='WeatherPlot',optype='other')
213 216 controllerObj.start()
@@ -1,112 +1,119
1 1 # Ing-AlexanderValdez
2 2 # Monitoreo de Pedestal
3 3
4 4 ############## IMPORTA LIBRERIAS ###################
5 5 import os,numpy,h5py
6 6 import sys,time
7 7 import matplotlib.pyplot as plt
8 8 ####################################################
9 path_ped = '/DATA_RM/TEST_PEDESTAL/P20211012-082745'
9 #################################################################
10 # LA FECHA 21-10-20 CORRESPONDE A LAS PRUEBAS DEL DIA MIERCOLES
11 # 1:15:51 pm hasta 3:49:32 pm
12 #################################################################
13
14 #path_ped = '/DATA_RM/TEST_PEDESTAL/P20211012-082745'
15 path_ped = '/DATA_RM/TEST_PEDESTAL/P20211020-131248'
10 16 # Metodo para verificar numero
11 17 def isNumber(str):
12 18 try:
13 19 float(str)
14 20 return True
15 21 except:
16 22 return False
17 23 # Metodo para extraer el arreglo
18 24 def getDatavaluefromDirFilename(path,file,value):
19 25 dir_file= path+"/"+file
20 26 fp = h5py.File(dir_file,'r')
21 27 array = fp['Data'].get(value)[()]
22 28 fp.close()
23 29 return array
24 30
25 31 # LISTA COMPLETA DE ARCHIVOS HDF5 Pedestal
26 32 LIST= sorted(os.listdir(path_ped))
27 33 m=len(LIST)
28 34 print("TOTAL DE ARCHIVOS DE PEDESTAL:",m)
29 35 # Contadores temporales
30 36 k= 0
31 37 l= 0
32 38 t= 0
33 39 # Marca de tiempo temporal
34 40 time_ = numpy.zeros([m])
35 41 # creacion de
36 42 for i in range(m):
43 print("order:",i)
37 44 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_pos")
38 45 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="ele_pos")
39 46 tmp_azi_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_vel")
40 47 tmp_ele_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_vel")# nuevo :D
41 48
42 49 time_[i] = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="utc")
43 50
44 51 k=k +tmp_azi_pos.shape[0]
45 52 l=l +tmp_ele_pos.shape[0]
46 53 t=t +tmp_azi_vel.shape[0]
47 54
48 55 print("TOTAL DE MUESTRAS, ARCHIVOS X100:",k)
49 56 time.sleep(5)
50 57 ######CREACION DE ARREGLOS CANTIDAD DE VALORES POR MUESTRA#################
51 58 azi_pos = numpy.zeros([k])
52 59 ele_pos = numpy.zeros([l])
53 60 time_azi_pos= numpy.zeros([k])
54 61 # Contadores temporales
55 62 p=0
56 63 r=0
57 64 z=0
58 65 # VARIABLES TMP para almacenar azimuth, elevacion y tiempo
59 66
60 67 #for filename in sorted(os.listdir(path_ped)):
61 68 # CONDICION POR LEER EN TIEMPO REAL NO OFFLINE
62 69
63 70 for filename in LIST:
64 71 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="azi_pos")
65 72 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="ele_pos")
66 73 # CONDICION POR LEER EN TIEMPO REAL NO OFFLINE
67 74
68 75 if z==(m-1):
69 76 tmp_azi_time=numpy.arange(time_[z],time_[z]+1,1/(tmp_azi_pos.shape[0]))
70 77 else:
71 78 tmp_azi_time=numpy.arange(time_[z],time_[z+1],(time_[z+1]-time_[z])/(tmp_azi_pos.shape[0]))
72 79
73 80 print(filename,time_[z])
74 81 print(z,tmp_azi_pos.shape[0])
75 82
76 83 i=0
77 84 for i in range(tmp_azi_pos.shape[0]):
78 85 index=p+i
79 86 azi_pos[index]=tmp_azi_pos[i]
80 87 time_azi_pos[index]=tmp_azi_time[i]
81 88 p=p+tmp_azi_pos.shape[0]
82 89 i=0
83 90 for i in range(tmp_ele_pos.shape[0]):
84 91 index=r+i
85 92 ele_pos[index]=tmp_ele_pos[i]
86 93 r=r+tmp_ele_pos.shape[0]
87 94
88 95
89 96 z+=1
90 97
91 98
92 99 ######## GRAFIQUEMOS Y VEAMOS LOS DATOS DEL Pedestal
93 100 fig, ax = plt.subplots(figsize=(16,8))
94 101 print(time_azi_pos.shape)
95 102 print(azi_pos.shape)
96 103 t=numpy.arange(time_azi_pos.shape[0])*0.01/(60.0)
97 104 plt.plot(t,azi_pos,label='AZIMUTH_POS',color='blue')
98 105
99 106 # AQUI ESTOY ADICIONANDO LA POSICION EN elevaciont=numpy.arange(len(ele_pos))*0.01/60.0
100 107 t=numpy.arange(len(ele_pos))*0.01/60.0
101 108 plt.plot(t,ele_pos,label='ELEVATION_POS',color='red')#*10
102 109
103 110 #ax.set_xlim(0, 9)
104 111 ax.set_ylim(-5, 400)
105 112 plt.ylabel("Azimuth Position")
106 113 plt.xlabel("Muestra")
107 114 plt.title('Azimuth Position vs Muestra ', fontsize=20)
108 115 axes = plt.gca()
109 116 axes.yaxis.grid()
110 117 plt.xticks(fontsize=16)
111 118 plt.yticks(fontsize=16)
112 119 plt.show()
@@ -1,78 +1,79
1 1 import os,sys,json
2 2 import datetime
3 3 import time
4 4 from schainpy.controller import Project
5 5 '''
6 6 NOTA:
7 7 Este script de prueba.
8 8 - Unidad del lectura 'HDFReader'.
9 9 - Unidad de procesamiento ParametersProc
10 10 - Operacion SpectralMomentsPlot
11 11
12 12 '''
13 13
14 14 #######################################################################
15 15 ################# RANGO DE PLOTEO######################################
16 16 #######################################################################
17 17 dBmin = '1'
18 18 dBmax = '65'
19 19 xmin = '0'
20 20 xmax ='24'
21 21 #tmmin = 16.2
22 22 #tmmax = 16.25
23 23 tmmin =15
24 24 tmmax =15.5
25 25 ymin = '0'
26 26 ymax = '600'
27 27 #######################################################################
28 28 #######################################################################
29 29 #######################################################################
30 30 #path = '/DATA_RM/TEST_HDF5_SPEC'
31 path = '/DATA_RM/TEST_HDF5_SPEC_23/6v/'
31 #path = '/DATA_RM/TEST_HDF5_SPEC_23/6v/'
32 path = '/DATA_RM/TEST_HDF5_19OCT'
32 33 figpath = '/home/soporte/Downloads/23/6v'
33 34 desc = "Simulator Test"
34 35 desc_data = {
35 36 'Data': {
36 37 'data_pow': ['Data/data_pow/channel00','Data/data_pow/channel01'],
37 38 'data_dop': ['Data/data_dop/channel00','Data/data_dop/channel01'],
38 39 'utctime':'Data/utctime'
39 40 },
40 41 'Metadata': {
41 42 'heightList':'Metadata/heightList',
42 43 'nIncohInt' :'Metadata/nIncohInt',
43 44 'nCohInt' :'Metadata/nCohInt',
44 45 'nProfiles' :'Metadata/nProfiles',
45 46 'channelList' :'Metadata/channelList',
46 47 'utctimeInit' :'Metadata/utctimeInit'
47 48
48 49 }
49 50 }
50 51
51 52 controllerObj = Project()
52 53
53 54 controllerObj.setup(id='10',name='Test Simulator',description=desc)
54 55
55 56 readUnitConfObj = controllerObj.addReadUnit(datatype='HDFReader',
56 57 path=path,
57 58 startDate="2021/01/01", #"2020/01/01",#today,
58 59 endDate= "2021/12/01", #"2020/12/30",#today,
59 60 startTime='00:00:00',
60 61 endTime='23:59:59',
61 62 delay=0,
62 63 #set=0,
63 64 online=0,
64 65 walk=1,
65 66 description= json.dumps(desc_data))#1
66 67
67 68 procUnitConfObjA = controllerObj.addProcUnit(datatype='ParametersProc',inputId=readUnitConfObj.getId())
68 69
69 70 opObj11 = procUnitConfObjA.addOperation(name='PowerPlot',optype='external')
70 71 opObj11.addParameter(name='xmin', value=tmmin)
71 72 opObj11.addParameter(name='xmax', value=tmmax)
72 73 opObj11.addParameter(name='zmin', value=dBmin)
73 74 opObj11.addParameter(name='zmax', value=dBmax)
74 75 opObj11.addParameter(name='save', value=figpath)
75 76 opObj11.addParameter(name='showprofile', value=0)
76 77 opObj11.addParameter(name='save_period', value=10)
77 78
78 79 controllerObj.start()
General Comments 0
You need to be logged in to leave comments. Login now