##// END OF EJS Templates
schain
RSS
LAST UPDATE INTEGRACION
avaldez -
r1388:ec0a56a3ee09
parent child
Show More
Side by Side Unified
Context file:
r1388:ec0a56a3ee09 Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,213
1 # Ing. AVP
2 # 06/10/2021
3 # ARCHIVO DE LECTURA
4 import os, sys
5 import datetime
6 import time
7 from schainpy.controller import Project
8 #### NOTA###########################################
9 # INPUT :
10 # VELOCIDAD PARAMETRO : V = 2Β°/seg
11 # MODO PULSE PAIR O MOMENTOS: 0 : Pulse Pair ,1 : Momentos
12 ######################################################
13 ##### PROCESAMIENTO ##################################
14 ##### OJO TENER EN CUENTA EL n= para el Pulse Pair ##
15 ##### O EL n= nFFTPoints ###
16 ######################################################
17 ######## BUSCAMOS EL numero de IPP equivalente 1Β°#####
18 ######## Sea V la velocidad del Pedestal en Β°/seg#####
19 ######## 1Β° sera Recorrido en un tiempo de 1/V ######
20 ######## IPP del Radar 400 useg --> 60 Km ############
21 ######## n = 1/(V(Β°/seg)*IPP(Km)) , NUMERO DE IPP ##
22 ######## n = 1/(V*IPP) #############################
23 ######## VELOCIDAD DEL PEDESTAL ######################
24 print("SETUP- RADAR METEOROLOGICO")
25 V = 10
26 mode = 1
27 #path = '/DATA_RM/23/6v'
28 path = '/DATA_RM/TEST_INTEGRACION_2M'
29 path_ped='/DATA_RM/TEST_PEDESTAL/P20211012-082745'
30 figpath_pp = "/home/soporte/Pictures/TEST_PP"
31 figpath_mom = "/home/soporte/Pictures/TEST_MOM"
32 plot = 0
33 integration = 1
34 save = 0
35 if save == 1:
36 if mode==0:
37 path_save = '/DATA_RM/TEST_HDF5_PP_23/6v'
38 path_save = '/DATA_RM/TEST_HDF5_PP'
39 path_save = '/DATA_RM/TEST_HDF5_PP_100'
40 else:
41 path_save = '/DATA_RM/TEST_HDF5_SPEC_23_V2/6v'
42
43 print("* PATH data ADQ :", path)
44 print("* Velocidad Pedestal :",V,"Β°/seg")
45 ############################ NRO Perfiles PROCESAMIENTO ###################
46 V=V
47 IPP=400*1e-6
48 n= int(1/(V*IPP))
49 print("* n - NRO Perfiles Proc:", n )
50 ################################## MODE ###################################
51 print("* Modo de Operacion :",mode)
52 if mode ==0:
53 print("* Met. Seleccionado : Pulse Pair")
54 else:
55 print("* Met. Momentos : Momentos")
56
57 ################################## MODE ###################################
58 print("* Grabado de datos :",save)
59 if save ==1:
60 if mode==0:
61 ope= "Pulse Pair"
62 else:
63 ope= "Momentos"
64 print("* Path-Save Data -", ope , path_save)
65
66 print("* Integracion de datos :",integration)
67
68 time.sleep(15)
69 #remotefolder = "/home/wmaster/graficos"
70 #######################################################################
71 ################# RANGO DE PLOTEO######################################
72 dBmin = '1'
73 dBmax = '85'
74 xmin = '15'
75 xmax = '15.25'
76 ymin = '0'
77 ymax = '600'
78 #######################################################################
79 ########################FECHA##########################################
80 str = datetime.date.today()
81 today = str.strftime("%Y/%m/%d")
82 str2 = str - datetime.timedelta(days=1)
83 yesterday = str2.strftime("%Y/%m/%d")
84 #######################################################################
85 ########################SIGNAL CHAIN ##################################
86 #######################################################################
87 desc = "USRP_test"
88 filename = "USRP_processing.xml"
89 controllerObj = Project()
90 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
91 #######################################################################
92 ######################## UNIDAD DE LECTURA#############################
93 #######################################################################
94 readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
95 path=path,
96 startDate="2021/01/01",#today,
97 endDate="2021/12/30",#today,
98 startTime='00:00:00',
99 endTime='23:59:59',
100 delay=0,
101 #set=0,
102 online=0,
103 walk=1,
104 ippKm = 60)
105
106 opObj11 = readUnitConfObj.addOperation(name='printInfo')
107
108 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
109
110 if mode ==0:
111 ####################### METODO PULSE PAIR ######################################################################
112 opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
113 opObj11.addParameter(name='n', value=int(n), format='int')#10 VOY A USAR 250 DADO QUE LA VELOCIDAD ES 10 GRADOS
114 #opObj11.addParameter(name='removeDC', value=1, format='int')
115 ####################### METODO Parametros ######################################################################
116 procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
117 if plot==1:
118 opObj11 = procUnitConfObjB.addOperation(name='GenericRTIPlot',optype='external')
119 opObj11.addParameter(name='attr_data', value='dataPP_POW')
120 opObj11.addParameter(name='colormap', value='jet')
121 opObj11.addParameter(name='xmin', value=xmin)
122 opObj11.addParameter(name='xmax', value=xmax)
123 opObj11.addParameter(name='zmin', value=dBmin)
124 opObj11.addParameter(name='zmax', value=dBmax)
125 opObj11.addParameter(name='save', value=figpath_pp)
126 opObj11.addParameter(name='showprofile', value=0)
127 opObj11.addParameter(name='save_period', value=50)
128
129 ####################### METODO ESCRITURA #######################################################################
130 if save==1:
131 opObj10 = procUnitConfObjB.addOperation(name='HDFWriter')
132 opObj10.addParameter(name='path',value=path_save)
133 #opObj10.addParameter(name='mode',value=0)
134 opObj10.addParameter(name='blocksPerFile',value='100',format='int')
135 opObj10.addParameter(name='metadataList',value='utctimeInit,timeZone,paramInterval,profileIndex,channelList,heightList,flagDataAsBlock',format='list')
136 opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,utctime',format='list')#,format='list'
137 if integration==1:
138 V=10
139 blocksPerfile=360
140 print("* Velocidad del Pedestal:",V)
141 tmp_blocksPerfile = 100
142 f_a_p= int(tmp_blocksPerfile/V)
143
144 opObj11 = procUnitConfObjB.addOperation(name='PedestalInformation')
145 opObj11.addParameter(name='path_ped', value=path_ped)
146 #opObj11.addParameter(name='path_adq', value=path_adq)
147 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
148 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
149 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
150 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
151 opObj11.addParameter(name='online', value='0', format='int')
152
153 opObj11 = procUnitConfObjB.addOperation(name='Block360')
154 opObj11.addParameter(name='n', value='10', format='int')
155 opObj11.addParameter(name='mode', value=mode, format='int')
156
157 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
158
159 opObj11= procUnitConfObjB.addOperation(name='WeatherPlot',optype='other')
160
161
162 else:
163 ####################### METODO SPECTROS ######################################################################
164 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
165 procUnitConfObjB.addParameter(name='nFFTPoints', value=n, format='int')
166 procUnitConfObjB.addParameter(name='nProfiles' , value=n, format='int')
167
168 procUnitConfObjC = controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjB.getId())
169 procUnitConfObjC.addOperation(name='SpectralMoments')
170 if plot==1:
171 dBmin = '1'
172 dBmax = '65'
173 opObj11 = procUnitConfObjC.addOperation(name='PowerPlot',optype='external')
174 opObj11.addParameter(name='xmin', value=xmin)
175 opObj11.addParameter(name='xmax', value=xmax)
176 opObj11.addParameter(name='zmin', value=dBmin)
177 opObj11.addParameter(name='zmax', value=dBmax)
178 opObj11.addParameter(name='save', value=figpath_mom)
179 opObj11.addParameter(name='showprofile', value=0)
180 opObj11.addParameter(name='save_period', value=100)
181
182 if save==1:
183 opObj10 = procUnitConfObjC.addOperation(name='HDFWriter')
184 opObj10.addParameter(name='path',value=path_save)
185 #opObj10.addParameter(name='mode',value=0)
186 opObj10.addParameter(name='blocksPerFile',value='360',format='int')
187 #opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
188 opObj10.addParameter(name='metadataList',value='utctimeInit,heightList,nIncohInt,nCohInt,nProfiles,channelList',format='list')#profileIndex
189 opObj10.addParameter(name='dataList',value='data_pow,data_dop,utctime',format='list')#,format='list'
190
191 if integration==1:
192 V=10
193 blocksPerfile=360
194 print("* Velocidad del Pedestal:",V)
195 tmp_blocksPerfile = 100
196 f_a_p= int(tmp_blocksPerfile/V)
197
198 opObj11 = procUnitConfObjC.addOperation(name='PedestalInformation')
199 opObj11.addParameter(name='path_ped', value=path_ped)
200 #opObj11.addParameter(name='path_adq', value=path_adq)
201 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
202 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
203 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
204 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
205 opObj11.addParameter(name='online', value='0', format='int')
206
207 opObj11 = procUnitConfObjC.addOperation(name='Block360')
208 opObj11.addParameter(name='n', value='10', format='int')
209 opObj11.addParameter(name='mode', value=mode, format='int')
210
211 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
212 opObj11= procUnitConfObjC.addOperation(name='WeatherPlot',optype='other')
213 controllerObj.start()
Show file before | Show file after | Hide comments
@@ -0,0 +1,112
1 # Ing-AlexanderValdez
2 # Monitoreo de Pedestal
3
4 ############## IMPORTA LIBRERIAS ###################
5 import os,numpy,h5py
6 import sys,time
7 import matplotlib.pyplot as plt
8 ####################################################
9 path_ped = '/DATA_RM/TEST_PEDESTAL/P20211012-082745'
10 # Metodo para verificar numero
11 def isNumber(str):
12 try:
13 float(str)
14 return True
15 except:
16 return False
17 # Metodo para extraer el arreglo
18 def getDatavaluefromDirFilename(path,file,value):
19 dir_file= path+"/"+file
20 fp = h5py.File(dir_file,'r')
21 array = fp['Data'].get(value)[()]
22 fp.close()
23 return array
24
25 # LISTA COMPLETA DE ARCHIVOS HDF5 Pedestal
26 LIST= sorted(os.listdir(path_ped))
27 m=len(LIST)
28 print("TOTAL DE ARCHIVOS DE PEDESTAL:",m)
29 # Contadores temporales
30 k= 0
31 l= 0
32 t= 0
33 # Marca de tiempo temporal
34 time_ = numpy.zeros([m])
35 # creacion de
36 for i in range(m):
37 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_pos")
38 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="ele_pos")
39 tmp_azi_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_vel")
40 tmp_ele_vel = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="azi_vel")# nuevo :D
41
42 time_[i] = getDatavaluefromDirFilename(path=path_ped,file=LIST[i],value="utc")
43
44 k=k +tmp_azi_pos.shape[0]
45 l=l +tmp_ele_pos.shape[0]
46 t=t +tmp_azi_vel.shape[0]
47
48 print("TOTAL DE MUESTRAS, ARCHIVOS X100:",k)
49 time.sleep(5)
50 ######CREACION DE ARREGLOS CANTIDAD DE VALORES POR MUESTRA#################
51 azi_pos = numpy.zeros([k])
52 ele_pos = numpy.zeros([l])
53 time_azi_pos= numpy.zeros([k])
54 # Contadores temporales
55 p=0
56 r=0
57 z=0
58 # VARIABLES TMP para almacenar azimuth, elevacion y tiempo
59
60 #for filename in sorted(os.listdir(path_ped)):
61 # CONDICION POR LEER EN TIEMPO REAL NO OFFLINE
62
63 for filename in LIST:
64 tmp_azi_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="azi_pos")
65 tmp_ele_pos = getDatavaluefromDirFilename(path=path_ped,file=filename,value="ele_pos")
66 # CONDICION POR LEER EN TIEMPO REAL NO OFFLINE
67
68 if z==(m-1):
69 tmp_azi_time=numpy.arange(time_[z],time_[z]+1,1/(tmp_azi_pos.shape[0]))
70 else:
71 tmp_azi_time=numpy.arange(time_[z],time_[z+1],(time_[z+1]-time_[z])/(tmp_azi_pos.shape[0]))
72
73 print(filename,time_[z])
74 print(z,tmp_azi_pos.shape[0])
75
76 i=0
77 for i in range(tmp_azi_pos.shape[0]):
78 index=p+i
79 azi_pos[index]=tmp_azi_pos[i]
80 time_azi_pos[index]=tmp_azi_time[i]
81 p=p+tmp_azi_pos.shape[0]
82 i=0
83 for i in range(tmp_ele_pos.shape[0]):
84 index=r+i
85 ele_pos[index]=tmp_ele_pos[i]
86 r=r+tmp_ele_pos.shape[0]
87
88
89 z+=1
90
91
92 ######## GRAFIQUEMOS Y VEAMOS LOS DATOS DEL Pedestal
93 fig, ax = plt.subplots(figsize=(16,8))
94 print(time_azi_pos.shape)
95 print(azi_pos.shape)
96 t=numpy.arange(time_azi_pos.shape[0])*0.01/(60.0)
97 plt.plot(t,azi_pos,label='AZIMUTH_POS',color='blue')
98
99 # AQUI ESTOY ADICIONANDO LA POSICION EN elevaciont=numpy.arange(len(ele_pos))*0.01/60.0
100 t=numpy.arange(len(ele_pos))*0.01/60.0
101 plt.plot(t,ele_pos,label='ELEVATION_POS',color='red')#*10
102
103 #ax.set_xlim(0, 9)
104 ax.set_ylim(-5, 400)
105 plt.ylabel("Azimuth Position")
106 plt.xlabel("Muestra")
107 plt.title('Azimuth Position vs Muestra ', fontsize=20)
108 axes = plt.gca()
109 axes.yaxis.grid()
110 plt.xticks(fontsize=16)
111 plt.yticks(fontsize=16)
112 plt.show()
Show file before | Show file after | Hide comments
@@ -0,0 +1,90
1 import os,sys,json
2 import datetime
3 import time
4 from schainpy.controller import Project
5 '''
6 NOTA:
7 Este script de prueba.
8 - Unidad del lectura 'HDFReader'.
9 - Unidad de procesamiento ParametersProc
10 - Operacion SpectralMomentsPlot
11
12 '''
13
14 #######################################################################
15 ################# RANGO DE PLOTEO######################################
16 #######################################################################
17 dBmin = '1'
18 dBmax = '65'
19 xmin = '0'
20 xmax ='24'
21 #tmmin = 16.2
22 #tmmax = 16.25
23 tmmin =15
24 tmmax =15.5
25 ymin = '0'
26 ymax = '600'
27 #######################################################################
28 #######################################################################
29 #######################################################################
30 #path = '/DATA_RM/TEST_HDF5_SPEC'
31 path = '/DATA_RM/TEST_HDF5_SPEC_23/6v/'
32 figpath = '/home/soporte/Downloads/23/6v'
33 path="/home/soporte/Downloads/params-20211015T174046Z-001/params"
34 desc = "Simulator Test"
35 desc_data = {
36 'Data': {
37 'data_spc': ['Data/data_spc/channel00','Data/data_spc/channel01'\
38 ,'Data/data_spc/channel02','Data/data_spc/channel03'\
39 ,'Data/data_spc/channel04','Data/data_spc/channel05'\
40 ,'Data/data_spc/channel06','Data/data_spc/channel07'\
41 ,'Data/data_spc/channel08','Data/data_spc/channel09'],
42 'utctime':'Data/utctime'
43 },
44 'Metadata': {
45 'type' :'Metadata/type',
46 'channelList' :'Metadata/channelList',
47 'heightList' :'Metadata/heightList',
48 'ippSeconds' :'Metadata/ippSeconds',
49 'nProfiles' :'Metadata/nProfiles',
50 'codeList' :'Metadata/codeList',
51 'timeZone' :'Metadata/timeZone',
52 'azimuthList' :'Metadata/azimuthList',
53 'elevationList' :'Metadata/elevationList',
54 'nCohInt' :'Metadata/nCohInt',
55 'nIncohInt' :'Metadata/nIncohInt',
56 'nFFTPoints' :'Metadata/nFFTPoints'
57
58 }
59 }
60
61 controllerObj = Project()
62
63 controllerObj.setup(id='10',name='Test Simulator',description=desc)
64
65 readUnitConfObj = controllerObj.addReadUnit(datatype='HDFReader',
66 path=path,
67 startDate="2021/01/01", #"2020/01/01",#today,
68 endDate= "2021/12/01", #"2020/12/30",#today,
69 startTime='00:00:00',
70 endTime='23:59:59',
71 delay=0,
72 #set=0,
73 online=0,
74 walk=0,
75 description= json.dumps(desc_data))#1
76
77 procUnitConfObjA = controllerObj.addProcUnit(datatype='ParametersProc',inputId=readUnitConfObj.getId())
78 procUnitConfObjA.addOperation(name='SpectralMoments')
79
80 '''
81 opObj11 = readUnitConfObj.addOperation(name='SpectraPlot',optype='external')
82 opObj11.addParameter(name='xmin', value=tmmin)
83 opObj11.addParameter(name='xmax', value=tmmax)
84 opObj11.addParameter(name='zmin', value=dBmin)
85 opObj11.addParameter(name='zmax', value=dBmax)
86 opObj11.addParameter(name='save', value=figpath)
87 opObj11.addParameter(name='showprofile', value=0)
88 opObj11.addParameter(name='save_period', value=10)
89 '''
90 controllerObj.start()
Show file before | Show file after | Hide comments
@@ -1,518 +1,519
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
7 7 from schainpy.utils import log
8 8 # libreria wradlib
9 9 import wradlib as wrl
10 10
11 11 EARTH_RADIUS = 6.3710e3
12 12
13 13
14 14 def ll2xy(lat1, lon1, lat2, lon2):
15 15
16 16 p = 0.017453292519943295
17 17 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
18 18 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
19 19 r = 12742 * numpy.arcsin(numpy.sqrt(a))
20 20 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
21 21 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
22 22 theta = -theta + numpy.pi/2
23 23 return r*numpy.cos(theta), r*numpy.sin(theta)
24 24
25 25
26 26 def km2deg(km):
27 27 '''
28 28 Convert distance in km to degrees
29 29 '''
30 30
31 31 return numpy.rad2deg(km/EARTH_RADIUS)
32 32
33 33
34 34
35 35 class SpectralMomentsPlot(SpectraPlot):
36 36 '''
37 37 Plot for Spectral Moments
38 38 '''
39 39 CODE = 'spc_moments'
40 40 # colormap = 'jet'
41 41 # plot_type = 'pcolor'
42 42
43 43 class DobleGaussianPlot(SpectraPlot):
44 44 '''
45 45 Plot for Double Gaussian Plot
46 46 '''
47 47 CODE = 'gaussian_fit'
48 48 # colormap = 'jet'
49 49 # plot_type = 'pcolor'
50 50
51 51 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
52 52 '''
53 53 Plot SpectraCut with Double Gaussian Fit
54 54 '''
55 55 CODE = 'cut_gaussian_fit'
56 56
57 57 class SnrPlot(RTIPlot):
58 58 '''
59 59 Plot for SNR Data
60 60 '''
61 61
62 62 CODE = 'snr'
63 63 colormap = 'jet'
64 64
65 65 def update(self, dataOut):
66 66
67 67 data = {
68 68 'snr': 10*numpy.log10(dataOut.data_snr)
69 69 }
70 70
71 71 return data, {}
72 72
73 73 class DopplerPlot(RTIPlot):
74 74 '''
75 75 Plot for DOPPLER Data (1st moment)
76 76 '''
77 77
78 78 CODE = 'dop'
79 79 colormap = 'jet'
80 80
81 81 def update(self, dataOut):
82 82
83 83 data = {
84 84 'dop': 10*numpy.log10(dataOut.data_dop)
85 85 }
86 86
87 87 return data, {}
88 88
89 89 class PowerPlot(RTIPlot):
90 90 '''
91 91 Plot for Power Data (0 moment)
92 92 '''
93 93
94 94 CODE = 'pow'
95 95 colormap = 'jet'
96 96
97 97 def update(self, dataOut):
98 98
99 99 data = {
100 100 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
101 101 }
102 102
103 103 return data, {}
104 104
105 105 class SpectralWidthPlot(RTIPlot):
106 106 '''
107 107 Plot for Spectral Width Data (2nd moment)
108 108 '''
109 109
110 110 CODE = 'width'
111 111 colormap = 'jet'
112 112
113 113 def update(self, dataOut):
114 114
115 115 data = {
116 116 'width': dataOut.data_width
117 117 }
118 118
119 119 return data, {}
120 120
121 121 class SkyMapPlot(Plot):
122 122 '''
123 123 Plot for meteors detection data
124 124 '''
125 125
126 126 CODE = 'param'
127 127
128 128 def setup(self):
129 129
130 130 self.ncols = 1
131 131 self.nrows = 1
132 132 self.width = 7.2
133 133 self.height = 7.2
134 134 self.nplots = 1
135 135 self.xlabel = 'Zonal Zenith Angle (deg)'
136 136 self.ylabel = 'Meridional Zenith Angle (deg)'
137 137 self.polar = True
138 138 self.ymin = -180
139 139 self.ymax = 180
140 140 self.colorbar = False
141 141
142 142 def plot(self):
143 143
144 144 arrayParameters = numpy.concatenate(self.data['param'])
145 145 error = arrayParameters[:, -1]
146 146 indValid = numpy.where(error == 0)[0]
147 147 finalMeteor = arrayParameters[indValid, :]
148 148 finalAzimuth = finalMeteor[:, 3]
149 149 finalZenith = finalMeteor[:, 4]
150 150
151 151 x = finalAzimuth * numpy.pi / 180
152 152 y = finalZenith
153 153
154 154 ax = self.axes[0]
155 155
156 156 if ax.firsttime:
157 157 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
158 158 else:
159 159 ax.plot.set_data(x, y)
160 160
161 161 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
162 162 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
163 163 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
164 164 dt2,
165 165 len(x))
166 166 self.titles[0] = title
167 167
168 168
169 169 class GenericRTIPlot(Plot):
170 170 '''
171 171 Plot for data_xxxx object
172 172 '''
173 173
174 174 CODE = 'param'
175 175 colormap = 'viridis'
176 176 plot_type = 'pcolorbuffer'
177 177
178 178 def setup(self):
179 179 self.xaxis = 'time'
180 180 self.ncols = 1
181 181 self.nrows = self.data.shape('param')[0]
182 182 self.nplots = self.nrows
183 183 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
184 184
185 185 if not self.xlabel:
186 186 self.xlabel = 'Time'
187 187
188 188 self.ylabel = 'Range [km]'
189 189 if not self.titles:
190 190 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
191 191
192 192 def update(self, dataOut):
193 193
194 194 data = {
195 195 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
196 196 }
197 197
198 198 meta = {}
199 199
200 200 return data, meta
201 201
202 202 def plot(self):
203 203 # self.data.normalize_heights()
204 204 self.x = self.data.times
205 205 self.y = self.data.yrange
206 206 self.z = self.data['param']
207 207
208 208 self.z = 10*numpy.log10(self.z)
209 209
210 210 self.z = numpy.ma.masked_invalid(self.z)
211 211
212 212 if self.decimation is None:
213 213 x, y, z = self.fill_gaps(self.x, self.y, self.z)
214 214 else:
215 215 x, y, z = self.fill_gaps(*self.decimate())
216 216
217 217 for n, ax in enumerate(self.axes):
218 218
219 219 self.zmax = self.zmax if self.zmax is not None else numpy.max(
220 220 self.z[n])
221 221 self.zmin = self.zmin if self.zmin is not None else numpy.min(
222 222 self.z[n])
223 223
224 224 if ax.firsttime:
225 225 if self.zlimits is not None:
226 226 self.zmin, self.zmax = self.zlimits[n]
227 227
228 228 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
229 229 vmin=self.zmin,
230 230 vmax=self.zmax,
231 231 cmap=self.cmaps[n]
232 232 )
233 233 else:
234 234 if self.zlimits is not None:
235 235 self.zmin, self.zmax = self.zlimits[n]
236 236 ax.collections.remove(ax.collections[0])
237 237 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
238 238 vmin=self.zmin,
239 239 vmax=self.zmax,
240 240 cmap=self.cmaps[n]
241 241 )
242 242
243 243
244 244 class PolarMapPlot(Plot):
245 245 '''
246 246 Plot for weather radar
247 247 '''
248 248
249 249 CODE = 'param'
250 250 colormap = 'seismic'
251 251
252 252 def setup(self):
253 253 self.ncols = 1
254 254 self.nrows = 1
255 255 self.width = 9
256 256 self.height = 8
257 257 self.mode = self.data.meta['mode']
258 258 if self.channels is not None:
259 259 self.nplots = len(self.channels)
260 260 self.nrows = len(self.channels)
261 261 else:
262 262 self.nplots = self.data.shape(self.CODE)[0]
263 263 self.nrows = self.nplots
264 264 self.channels = list(range(self.nplots))
265 265 if self.mode == 'E':
266 266 self.xlabel = 'Longitude'
267 267 self.ylabel = 'Latitude'
268 268 else:
269 269 self.xlabel = 'Range (km)'
270 270 self.ylabel = 'Height (km)'
271 271 self.bgcolor = 'white'
272 272 self.cb_labels = self.data.meta['units']
273 273 self.lat = self.data.meta['latitude']
274 274 self.lon = self.data.meta['longitude']
275 275 self.xmin, self.xmax = float(
276 276 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
277 277 self.ymin, self.ymax = float(
278 278 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
279 279 # self.polar = True
280 280
281 281 def plot(self):
282 282
283 283 for n, ax in enumerate(self.axes):
284 284 data = self.data['param'][self.channels[n]]
285 285
286 286 zeniths = numpy.linspace(
287 287 0, self.data.meta['max_range'], data.shape[1])
288 288 if self.mode == 'E':
289 289 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
290 290 r, theta = numpy.meshgrid(zeniths, azimuths)
291 291 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
292 292 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
293 293 x = km2deg(x) + self.lon
294 294 y = km2deg(y) + self.lat
295 295 else:
296 296 azimuths = numpy.radians(self.data.yrange)
297 297 r, theta = numpy.meshgrid(zeniths, azimuths)
298 298 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
299 299 self.y = zeniths
300 300
301 301 if ax.firsttime:
302 302 if self.zlimits is not None:
303 303 self.zmin, self.zmax = self.zlimits[n]
304 304 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
305 305 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
306 306 vmin=self.zmin,
307 307 vmax=self.zmax,
308 308 cmap=self.cmaps[n])
309 309 else:
310 310 if self.zlimits is not None:
311 311 self.zmin, self.zmax = self.zlimits[n]
312 312 ax.collections.remove(ax.collections[0])
313 313 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
314 314 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
315 315 vmin=self.zmin,
316 316 vmax=self.zmax,
317 317 cmap=self.cmaps[n])
318 318
319 319 if self.mode == 'A':
320 320 continue
321 321
322 322 # plot district names
323 323 f = open('/data/workspace/schain_scripts/distrito.csv')
324 324 for line in f:
325 325 label, lon, lat = [s.strip() for s in line.split(',') if s]
326 326 lat = float(lat)
327 327 lon = float(lon)
328 328 # ax.plot(lon, lat, '.b', ms=2)
329 329 ax.text(lon, lat, label.decode('utf8'), ha='center',
330 330 va='bottom', size='8', color='black')
331 331
332 332 # plot limites
333 333 limites = []
334 334 tmp = []
335 335 for line in open('/data/workspace/schain_scripts/lima.csv'):
336 336 if '#' in line:
337 337 if tmp:
338 338 limites.append(tmp)
339 339 tmp = []
340 340 continue
341 341 values = line.strip().split(',')
342 342 tmp.append((float(values[0]), float(values[1])))
343 343 for points in limites:
344 344 ax.add_patch(
345 345 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
346 346
347 347 # plot Cuencas
348 348 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
349 349 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
350 350 values = [line.strip().split(',') for line in f]
351 351 points = [(float(s[0]), float(s[1])) for s in values]
352 352 ax.add_patch(Polygon(points, ec='b', fc='none'))
353 353
354 354 # plot grid
355 355 for r in (15, 30, 45, 60):
356 356 ax.add_artist(plt.Circle((self.lon, self.lat),
357 357 km2deg(r), color='0.6', fill=False, lw=0.2))
358 358 ax.text(
359 359 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
360 360 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
361 361 '{}km'.format(r),
362 362 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
363 363
364 364 if self.mode == 'E':
365 365 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
366 366 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
367 367 else:
368 368 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
369 369 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
370 370
371 371 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
372 372 self.titles = ['{} {}'.format(
373 373 self.data.parameters[x], title) for x in self.channels]
374 374
375 375 class WeatherPlot(Plot):
376 376 CODE = 'weather'
377 377 plot_name = 'weather'
378 378 plot_type = 'ppistyle'
379 379 buffering = False
380 380
381 381 def setup(self):
382 382 self.ncols = 1
383 383 self.nrows = 1
384 384 self.nplots= 1
385 385 self.ylabel= 'Range [Km]'
386 386 self.titles= ['Weather']
387 387 self.colorbar=False
388 388 self.width =8
389 389 self.height =8
390 390 self.ini =0
391 391 self.len_azi =0
392 392 self.buffer_ini = None
393 393 self.buffer_azi = None
394 394 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
395 395 self.flag =0
396 396 self.indicador= 0
397 397
398 398 def update(self, dataOut):
399 399
400 400 data = {}
401 401 meta = {}
402 data['weather'] = 10*numpy.log10(dataOut.data_360[0]/(250**2))
403 print(data['weather'])
402 print("aprox",dataOut.data_360[0])
403 data['weather'] = 10*numpy.log10(dataOut.data_360[0]/(250.0))
404 #print(data['weather'])
404 405 data['azi'] = dataOut.data_azi
405 406 print("UPDATE",data['azi'])
406 407 return data, meta
407 408
408 409 def const_ploteo(self,data_weather,data_azi,step,res):
409 410 #print("data_weather",data_weather)
410 411 print("data_azi",data_azi)
411 412 print("step",step)
412 413 if self.ini==0:
413 414 #------- AZIMUTH
414 415 n = (360/res)-len(data_azi)
415 416 start = data_azi[-1] + res
416 417 end = data_azi[0] - res
417 418 if start>end:
418 419 end = end + 360
419 420 azi_vacia = numpy.linspace(start,end,int(n))
420 421 azi_vacia = numpy.where(azi_vacia>360,azi_vacia-360,azi_vacia)
421 422 data_azi = numpy.hstack((data_azi,azi_vacia))
422 423 # RADAR
423 424 val_mean = numpy.mean(data_weather[:,0])
424 425 data_weather_cmp = numpy.ones([(360-data_weather.shape[0]),data_weather.shape[1]])*val_mean
425 426 data_weather = numpy.vstack((data_weather,data_weather_cmp))
426 427 else:
427 428 # azimuth
428 429 flag=0
429 430 start_azi = self.res_azi[0]
430 431 start = data_azi[0]
431 432 end = data_azi[-1]
432 433 print("start",start)
433 434 print("end",end)
434 435 if start< start_azi:
435 436 start = start +360
436 437 if end <start_azi:
437 438 end = end +360
438 439
439 440 print("start",start)
440 441 print("end",end)
441 442 #### AQUI SERA LA MAGIA
442 443 pos_ini = int((start-start_azi)/res)
443 444 len_azi = len(data_azi)
444 445 if (360-pos_ini)<len_azi:
445 446 if pos_ini+1==360:
446 447 pos_ini=0
447 448 else:
448 449 flag=1
449 450 dif= 360-pos_ini
450 451 comp= len_azi-dif
451 452
452 453 print(pos_ini)
453 454 print(len_azi)
454 455 print("shape",self.res_azi.shape)
455 456 if flag==0:
456 457 # AZIMUTH
457 458 self.res_azi[pos_ini:pos_ini+len_azi] = data_azi
458 459 # RADAR
459 460 self.res_weather[pos_ini:pos_ini+len_azi,:] = data_weather
460 461 else:
461 462 # AZIMUTH
462 463 self.res_azi[pos_ini:pos_ini+dif] = data_azi[0:dif]
463 464 self.res_azi[0:comp] = data_azi[dif:]
464 465 # RADAR
465 466 self.res_weather[pos_ini:pos_ini+dif,:] = data_weather[0:dif,:]
466 467 self.res_weather[0:comp,:] = data_weather[dif:,:]
467 468 flag=0
468 469 data_azi = self.res_azi
469 470 data_weather = self.res_weather
470 471
471 472 return data_weather,data_azi
472 473
473 474 def plot(self):
474 475 print("--------------------------------------",self.ini,"-----------------------------------")
475 476 #numpy.set_printoptions(suppress=True)
476 477 #print(self.data.times)
477 478 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
478 479 data = self.data[-1]
479 480 # ALTURA altura_tmp_h
480 481 altura_h = (data['weather'].shape[1])/10.0
481 482 stoprange = float(altura_h*1.5)#stoprange = float(33*1.5) por ahora 400
482 483 rangestep = float(0.15)
483 484 r = numpy.arange(0, stoprange, rangestep)
484 485 self.y = 2*r
485 486 # RADAR
486 487 #data_weather = data['weather']
487 488 # PEDESTAL
488 489 #data_azi = data['azi']
489 490 res = 1
490 491 # STEP
491 492 step = (360/(res*data['weather'].shape[0]))
492 493 #print("shape wr_data", wr_data.shape)
493 494 #print("shape wr_azi",wr_azi.shape)
494 495 #print("step",step)
495 496 print("Time---->",self.data.times[-1],thisDatetime)
496 497 #print("alturas", len(self.y))
497 498 self.res_weather, self.res_azi = self.const_ploteo(data_weather=data['weather'],data_azi=data['azi'],step=step,res=res)
498 499 #numpy.set_printoptions(suppress=True)
499 500 #print("resultado",self.res_azi)
500 501 ##########################################################
501 502 ################# PLOTEO ###################
502 503 ##########################################################
503 504
504 505 for i,ax in enumerate(self.axes):
505 506 if ax.firsttime:
506 507 plt.clf()
507 508 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=1, vmax=60)
508 509 else:
509 510 plt.clf()
510 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=1, vmax=60)
511 cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=0, vmax=60)
511 512 caax = cgax.parasites[0]
512 513 paax = cgax.parasites[1]
513 514 cbar = plt.gcf().colorbar(pm, pad=0.075)
514 515 caax.set_xlabel('x_range [km]')
515 516 caax.set_ylabel('y_range [km]')
516 517 plt.text(1.0, 1.05, 'azimuth '+str(thisDatetime)+"step"+str(self.ini), transform=caax.transAxes, va='bottom',ha='right')
517 518
518 519 self.ini= self.ini+1
Show file before | Show file after | Hide comments
@@ -1,4450 +1,4475
1 1 import numpy,os,h5py
2 2 import math
3 3 from scipy import optimize, interpolate, signal, stats, ndimage
4 4 import scipy
5 5 import re
6 6 import datetime
7 7 import copy
8 8 import sys
9 9 import importlib
10 10 import itertools
11 11 from multiprocessing import Pool, TimeoutError
12 12 from multiprocessing.pool import ThreadPool
13 13 import time
14 14
15 15 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
16 16 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
17 17 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
18 18 from scipy import asarray as ar,exp
19 19 from scipy.optimize import curve_fit
20 20 from schainpy.utils import log
21 21 import warnings
22 22 from numpy import NaN
23 23 from scipy.optimize.optimize import OptimizeWarning
24 24 warnings.filterwarnings('ignore')
25 25
26 26 import matplotlib.pyplot as plt
27 27
28 28 SPEED_OF_LIGHT = 299792458
29 29
30 30 '''solving pickling issue'''
31 31
32 32 def _pickle_method(method):
33 33 func_name = method.__func__.__name__
34 34 obj = method.__self__
35 35 cls = method.__self__.__class__
36 36 return _unpickle_method, (func_name, obj, cls)
37 37
38 38 def _unpickle_method(func_name, obj, cls):
39 39 for cls in cls.mro():
40 40 try:
41 41 func = cls.__dict__[func_name]
42 42 except KeyError:
43 43 pass
44 44 else:
45 45 break
46 46 return func.__get__(obj, cls)
47 47
48 48 def isNumber(str):
49 49 try:
50 50 float(str)
51 51 return True
52 52 except:
53 53 return False
54 54
55 55 class ParametersProc(ProcessingUnit):
56 56
57 57 METHODS = {}
58 58 nSeconds = None
59 59
60 60 def __init__(self):
61 61 ProcessingUnit.__init__(self)
62 62
63 63 # self.objectDict = {}
64 64 self.buffer = None
65 65 self.firstdatatime = None
66 66 self.profIndex = 0
67 67 self.dataOut = Parameters()
68 68 self.setupReq = False #Agregar a todas las unidades de proc
69 69
70 70 def __updateObjFromInput(self):
71 71
72 72 self.dataOut.inputUnit = self.dataIn.type
73 73
74 74 self.dataOut.timeZone = self.dataIn.timeZone
75 75 self.dataOut.dstFlag = self.dataIn.dstFlag
76 76 self.dataOut.errorCount = self.dataIn.errorCount
77 77 self.dataOut.useLocalTime = self.dataIn.useLocalTime
78 78
79 79 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
80 80 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
81 81 self.dataOut.channelList = self.dataIn.channelList
82 82 self.dataOut.heightList = self.dataIn.heightList
83 83 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
84 84 # self.dataOut.nHeights = self.dataIn.nHeights
85 85 # self.dataOut.nChannels = self.dataIn.nChannels
86 86 # self.dataOut.nBaud = self.dataIn.nBaud
87 87 # self.dataOut.nCode = self.dataIn.nCode
88 88 # self.dataOut.code = self.dataIn.code
89 89 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
90 90 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
91 91 # self.dataOut.utctime = self.firstdatatime
92 92 self.dataOut.utctime = self.dataIn.utctime
93 93 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
94 94 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
95 95 self.dataOut.nCohInt = self.dataIn.nCohInt
96 96 # self.dataOut.nIncohInt = 1
97 97 # self.dataOut.ippSeconds = self.dataIn.ippSeconds
98 98 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
99 99 self.dataOut.timeInterval1 = self.dataIn.timeInterval
100 100 self.dataOut.heightList = self.dataIn.heightList
101 101 self.dataOut.frequency = self.dataIn.frequency
102 102 # self.dataOut.noise = self.dataIn.noise
103 103
104 104 def run(self):
105 105
106 106
107 107 #print("HOLA MUNDO SOY YO")
108 108 #---------------------- Voltage Data ---------------------------
109 109
110 110 if self.dataIn.type == "Voltage":
111 111
112 112 self.__updateObjFromInput()
113 113 self.dataOut.data_pre = self.dataIn.data.copy()
114 114 self.dataOut.flagNoData = False
115 115 self.dataOut.utctimeInit = self.dataIn.utctime
116 116 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
117 117
118 118 if hasattr(self.dataIn, 'flagDataAsBlock'):
119 119 self.dataOut.flagDataAsBlock = self.dataIn.flagDataAsBlock
120 120
121 121 if hasattr(self.dataIn, 'profileIndex'):
122 122 self.dataOut.profileIndex = self.dataIn.profileIndex
123 123
124 124 if hasattr(self.dataIn, 'dataPP_POW'):
125 125 self.dataOut.dataPP_POW = self.dataIn.dataPP_POW
126 126
127 127 if hasattr(self.dataIn, 'dataPP_POWER'):
128 128 self.dataOut.dataPP_POWER = self.dataIn.dataPP_POWER
129 129
130 130 if hasattr(self.dataIn, 'dataPP_DOP'):
131 131 self.dataOut.dataPP_DOP = self.dataIn.dataPP_DOP
132 132
133 133 if hasattr(self.dataIn, 'dataPP_SNR'):
134 134 self.dataOut.dataPP_SNR = self.dataIn.dataPP_SNR
135 135
136 136 if hasattr(self.dataIn, 'dataPP_WIDTH'):
137 137 self.dataOut.dataPP_WIDTH = self.dataIn.dataPP_WIDTH
138 138 return
139 139
140 140 #---------------------- Spectra Data ---------------------------
141 141
142 142 if self.dataIn.type == "Spectra":
143 143 #print("que paso en spectra")
144 144 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
145 145 self.dataOut.data_spc = self.dataIn.data_spc
146 146 self.dataOut.data_cspc = self.dataIn.data_cspc
147 147 self.dataOut.nProfiles = self.dataIn.nProfiles
148 148 self.dataOut.nIncohInt = self.dataIn.nIncohInt
149 149 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
150 150 self.dataOut.ippFactor = self.dataIn.ippFactor
151 151 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
152 152 self.dataOut.spc_noise = self.dataIn.getNoise()
153 153 self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
154 154 # self.dataOut.normFactor = self.dataIn.normFactor
155 155 self.dataOut.pairsList = self.dataIn.pairsList
156 156 self.dataOut.groupList = self.dataIn.pairsList
157 157 self.dataOut.flagNoData = False
158 158
159 159 if hasattr(self.dataIn, 'flagDataAsBlock'):
160 160 self.dataOut.flagDataAsBlock = self.dataIn.flagDataAsBlock
161 161
162 162 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
163 163 self.dataOut.ChanDist = self.dataIn.ChanDist
164 164 else: self.dataOut.ChanDist = None
165 165
166 166 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
167 167 # self.dataOut.VelRange = self.dataIn.VelRange
168 168 #else: self.dataOut.VelRange = None
169 169
170 170 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
171 171 self.dataOut.RadarConst = self.dataIn.RadarConst
172 172
173 173 if hasattr(self.dataIn, 'NPW'): #NPW
174 174 self.dataOut.NPW = self.dataIn.NPW
175 175
176 176 if hasattr(self.dataIn, 'COFA'): #COFA
177 177 self.dataOut.COFA = self.dataIn.COFA
178 178
179 179
180 180
181 181 #---------------------- Correlation Data ---------------------------
182 182
183 183 if self.dataIn.type == "Correlation":
184 184 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
185 185
186 186 self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
187 187 self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
188 188 self.dataOut.groupList = (acf_pairs, ccf_pairs)
189 189
190 190 self.dataOut.abscissaList = self.dataIn.lagRange
191 191 self.dataOut.noise = self.dataIn.noise
192 192 self.dataOut.data_snr = self.dataIn.SNR
193 193 self.dataOut.flagNoData = False
194 194 self.dataOut.nAvg = self.dataIn.nAvg
195 195
196 196 #---------------------- Parameters Data ---------------------------
197 197
198 198 if self.dataIn.type == "Parameters":
199 199 self.dataOut.copy(self.dataIn)
200 200 self.dataOut.flagNoData = False
201 #print("yo si entre")
201 202
202 203 return True
203 204
204 205 self.__updateObjFromInput()
206 #print("yo si entre2")
205 207
206 208 self.dataOut.utctimeInit = self.dataIn.utctime
207 209 self.dataOut.paramInterval = self.dataIn.timeInterval
208 210 #print("soy spectra ",self.dataOut.utctimeInit)
209 211 return
210 212
211 213
212 214 def target(tups):
213 215
214 216 obj, args = tups
215 217
216 218 return obj.FitGau(args)
217 219
218 220 class RemoveWideGC(Operation):
219 221 ''' This class remove the wide clutter and replace it with a simple interpolation points
220 222 This mainly applies to CLAIRE radar
221 223
222 224 ClutterWidth : Width to look for the clutter peak
223 225
224 226 Input:
225 227
226 228 self.dataOut.data_pre : SPC and CSPC
227 229 self.dataOut.spc_range : To select wind and rainfall velocities
228 230
229 231 Affected:
230 232
231 233 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
232 234
233 235 Written by D. ScipiΓ³n 25.02.2021
234 236 '''
235 237 def __init__(self):
236 238 Operation.__init__(self)
237 239 self.i = 0
238 240 self.ich = 0
239 241 self.ir = 0
240 242
241 243 def run(self, dataOut, ClutterWidth=2.5):
242 244 # print ('Entering RemoveWideGC ... ')
243 245
244 246 self.spc = dataOut.data_pre[0].copy()
245 247 self.spc_out = dataOut.data_pre[0].copy()
246 248 self.Num_Chn = self.spc.shape[0]
247 249 self.Num_Hei = self.spc.shape[2]
248 250 VelRange = dataOut.spc_range[2][:-1]
249 251 dv = VelRange[1]-VelRange[0]
250 252
251 253 # Find the velocities that corresponds to zero
252 254 gc_values = numpy.squeeze(numpy.where(numpy.abs(VelRange) <= ClutterWidth))
253 255
254 256 # Removing novalid data from the spectra
255 257 for ich in range(self.Num_Chn) :
256 258 for ir in range(self.Num_Hei) :
257 259 # Estimate the noise at each range
258 260 HSn = hildebrand_sekhon(self.spc[ich,:,ir],dataOut.nIncohInt)
259 261
260 262 # Removing the noise floor at each range
261 263 novalid = numpy.where(self.spc[ich,:,ir] < HSn)
262 264 self.spc[ich,novalid,ir] = HSn
263 265
264 266 junk = numpy.append(numpy.insert(numpy.squeeze(self.spc[ich,gc_values,ir]),0,HSn),HSn)
265 267 j1index = numpy.squeeze(numpy.where(numpy.diff(junk)>0))
266 268 j2index = numpy.squeeze(numpy.where(numpy.diff(junk)<0))
267 269 if ((numpy.size(j1index)<=1) | (numpy.size(j2index)<=1)) :
268 270 continue
269 271 junk3 = numpy.squeeze(numpy.diff(j1index))
270 272 junk4 = numpy.squeeze(numpy.diff(j2index))
271 273
272 274 valleyindex = j2index[numpy.where(junk4>1)]
273 275 peakindex = j1index[numpy.where(junk3>1)]
274 276
275 277 isvalid = numpy.squeeze(numpy.where(numpy.abs(VelRange[gc_values[peakindex]]) <= 2.5*dv))
276 278 if numpy.size(isvalid) == 0 :
277 279 continue
278 280 if numpy.size(isvalid) >1 :
279 281 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
280 282 isvalid = isvalid[vindex]
281 283
282 284 # clutter peak
283 285 gcpeak = peakindex[isvalid]
284 286 vl = numpy.where(valleyindex < gcpeak)
285 287 if numpy.size(vl) == 0:
286 288 continue
287 289 gcvl = valleyindex[vl[0][-1]]
288 290 vr = numpy.where(valleyindex > gcpeak)
289 291 if numpy.size(vr) == 0:
290 292 continue
291 293 gcvr = valleyindex[vr[0][0]]
292 294
293 295 # Removing the clutter
294 296 interpindex = numpy.array([gc_values[gcvl], gc_values[gcvr]])
295 297 gcindex = gc_values[gcvl+1:gcvr-1]
296 298 self.spc_out[ich,gcindex,ir] = numpy.interp(VelRange[gcindex],VelRange[interpindex],self.spc[ich,interpindex,ir])
297 299
298 300 dataOut.data_pre[0] = self.spc_out
299 301 #print ('Leaving RemoveWideGC ... ')
300 302 return dataOut
301 303
302 304 class SpectralFilters(Operation):
303 305 ''' This class allows to replace the novalid values with noise for each channel
304 306 This applies to CLAIRE RADAR
305 307
306 308 PositiveLimit : RightLimit of novalid data
307 309 NegativeLimit : LeftLimit of novalid data
308 310
309 311 Input:
310 312
311 313 self.dataOut.data_pre : SPC and CSPC
312 314 self.dataOut.spc_range : To select wind and rainfall velocities
313 315
314 316 Affected:
315 317
316 318 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
317 319
318 320 Written by D. ScipiΓ³n 29.01.2021
319 321 '''
320 322 def __init__(self):
321 323 Operation.__init__(self)
322 324 self.i = 0
323 325
324 326 def run(self, dataOut, ):
325 327
326 328 self.spc = dataOut.data_pre[0].copy()
327 329 self.Num_Chn = self.spc.shape[0]
328 330 VelRange = dataOut.spc_range[2]
329 331
330 332 # novalid corresponds to data within the Negative and PositiveLimit
331 333
332 334
333 335 # Removing novalid data from the spectra
334 336 for i in range(self.Num_Chn):
335 337 self.spc[i,novalid,:] = dataOut.noise[i]
336 338 dataOut.data_pre[0] = self.spc
337 339 return dataOut
338 340
339 341 class GaussianFit(Operation):
340 342
341 343 '''
342 344 Function that fit of one and two generalized gaussians (gg) based
343 345 on the PSD shape across an "power band" identified from a cumsum of
344 346 the measured spectrum - noise.
345 347
346 348 Input:
347 349 self.dataOut.data_pre : SelfSpectra
348 350
349 351 Output:
350 352 self.dataOut.SPCparam : SPC_ch1, SPC_ch2
351 353
352 354 '''
353 355 def __init__(self):
354 356 Operation.__init__(self)
355 357 self.i=0
356 358
357 359
358 360 # def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
359 361 def run(self, dataOut, SNRdBlimit=-9, method='generalized'):
360 362 """This routine will find a couple of generalized Gaussians to a power spectrum
361 363 methods: generalized, squared
362 364 input: spc
363 365 output:
364 366 noise, amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1
365 367 """
366 368 print ('Entering ',method,' double Gaussian fit')
367 369 self.spc = dataOut.data_pre[0].copy()
368 370 self.Num_Hei = self.spc.shape[2]
369 371 self.Num_Bin = self.spc.shape[1]
370 372 self.Num_Chn = self.spc.shape[0]
371 373
372 374 start_time = time.time()
373 375
374 376 pool = Pool(processes=self.Num_Chn)
375 377 args = [(dataOut.spc_range[2], ich, dataOut.spc_noise[ich], dataOut.nIncohInt, SNRdBlimit) for ich in range(self.Num_Chn)]
376 378 objs = [self for __ in range(self.Num_Chn)]
377 379 attrs = list(zip(objs, args))
378 380 DGauFitParam = pool.map(target, attrs)
379 381 # Parameters:
380 382 # 0. Noise, 1. Amplitude, 2. Shift, 3. Width 4. Power
381 383 dataOut.DGauFitParams = numpy.asarray(DGauFitParam)
382 384
383 385 # Double Gaussian Curves
384 386 gau0 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
385 387 gau0[:] = numpy.NaN
386 388 gau1 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
387 389 gau1[:] = numpy.NaN
388 390 x_mtr = numpy.transpose(numpy.tile(dataOut.getVelRange(1)[:-1], (self.Num_Hei,1)))
389 391 for iCh in range(self.Num_Chn):
390 392 N0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,0]] * self.Num_Bin))
391 393 N1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,1]] * self.Num_Bin))
392 394 A0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,0]] * self.Num_Bin))
393 395 A1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,1]] * self.Num_Bin))
394 396 v0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,0]] * self.Num_Bin))
395 397 v1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,1]] * self.Num_Bin))
396 398 s0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,0]] * self.Num_Bin))
397 399 s1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,1]] * self.Num_Bin))
398 400 if method == 'genealized':
399 401 p0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,0]] * self.Num_Bin))
400 402 p1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,1]] * self.Num_Bin))
401 403 elif method == 'squared':
402 404 p0 = 2.
403 405 p1 = 2.
404 406 gau0[iCh] = A0*numpy.exp(-0.5*numpy.abs((x_mtr-v0)/s0)**p0)+N0
405 407 gau1[iCh] = A1*numpy.exp(-0.5*numpy.abs((x_mtr-v1)/s1)**p1)+N1
406 408 dataOut.GaussFit0 = gau0
407 409 dataOut.GaussFit1 = gau1
408 410
409 411 print('Leaving ',method ,' double Gaussian fit')
410 412 return dataOut
411 413
412 414 def FitGau(self, X):
413 415 # print('Entering FitGau')
414 416 # Assigning the variables
415 417 Vrange, ch, wnoise, num_intg, SNRlimit = X
416 418 # Noise Limits
417 419 noisebl = wnoise * 0.9
418 420 noisebh = wnoise * 1.1
419 421 # Radar Velocity
420 422 Va = max(Vrange)
421 423 deltav = Vrange[1] - Vrange[0]
422 424 x = numpy.arange(self.Num_Bin)
423 425
424 426 # print ('stop 0')
425 427
426 428 # 5 parameters, 2 Gaussians
427 429 DGauFitParam = numpy.zeros([5, self.Num_Hei,2])
428 430 DGauFitParam[:] = numpy.NaN
429 431
430 432 # SPCparam = []
431 433 # SPC_ch1 = numpy.zeros([self.Num_Bin,self.Num_Hei])
432 434 # SPC_ch2 = numpy.zeros([self.Num_Bin,self.Num_Hei])
433 435 # SPC_ch1[:] = 0 #numpy.NaN
434 436 # SPC_ch2[:] = 0 #numpy.NaN
435 437 # print ('stop 1')
436 438 for ht in range(self.Num_Hei):
437 439 # print (ht)
438 440 # print ('stop 2')
439 441 # Spectra at each range
440 442 spc = numpy.asarray(self.spc)[ch,:,ht]
441 443 snr = ( spc.mean() - wnoise ) / wnoise
442 444 snrdB = 10.*numpy.log10(snr)
443 445
444 446 #print ('stop 3')
445 447 if snrdB < SNRlimit :
446 448 # snr = numpy.NaN
447 449 # SPC_ch1[:,ht] = 0#numpy.NaN
448 450 # SPC_ch1[:,ht] = 0#numpy.NaN
449 451 # SPCparam = (SPC_ch1,SPC_ch2)
450 452 # print ('SNR less than SNRth')
451 453 continue
452 454 # wnoise = hildebrand_sekhon(spc,num_intg)
453 455 # print ('stop 2.01')
454 456 #############################################
455 457 # normalizing spc and noise
456 458 # This part differs from gg1
457 459 # spc_norm_max = max(spc) #commented by D. ScipiΓ³n 19.03.2021
458 460 #spc = spc / spc_norm_max
459 461 # pnoise = pnoise #/ spc_norm_max #commented by D. ScipiΓ³n 19.03.2021
460 462 #############################################
461 463
462 464 # print ('stop 2.1')
463 465 fatspectra=1.0
464 466 # noise per channel.... we might want to use the noise at each range
465 467
466 468 # wnoise = noise_ #/ spc_norm_max #commented by D. ScipiΓ³n 19.03.2021
467 469 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
468 470 #if wnoise>1.1*pnoise: # to be tested later
469 471 # wnoise=pnoise
470 472 # noisebl = wnoise*0.9
471 473 # noisebh = wnoise*1.1
472 474 spc = spc - wnoise # signal
473 475
474 476 # print ('stop 2.2')
475 477 minx = numpy.argmin(spc)
476 478 #spcs=spc.copy()
477 479 spcs = numpy.roll(spc,-minx)
478 480 cum = numpy.cumsum(spcs)
479 481 # tot_noise = wnoise * self.Num_Bin #64;
480 482
481 483 # print ('stop 2.3')
482 484 # snr = sum(spcs) / tot_noise
483 485 # snrdB = 10.*numpy.log10(snr)
484 486 #print ('stop 3')
485 487 # if snrdB < SNRlimit :
486 488 # snr = numpy.NaN
487 489 # SPC_ch1[:,ht] = 0#numpy.NaN
488 490 # SPC_ch1[:,ht] = 0#numpy.NaN
489 491 # SPCparam = (SPC_ch1,SPC_ch2)
490 492 # print ('SNR less than SNRth')
491 493 # continue
492 494
493 495
494 496 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
495 497 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
496 498 # print ('stop 4')
497 499 cummax = max(cum)
498 500 epsi = 0.08 * fatspectra # cumsum to narrow down the energy region
499 501 cumlo = cummax * epsi
500 502 cumhi = cummax * (1-epsi)
501 503 powerindex = numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
502 504
503 505 # print ('stop 5')
504 506 if len(powerindex) < 1:# case for powerindex 0
505 507 # print ('powerindex < 1')
506 508 continue
507 509 powerlo = powerindex[0]
508 510 powerhi = powerindex[-1]
509 511 powerwidth = powerhi-powerlo
510 512 if powerwidth <= 1:
511 513 # print('powerwidth <= 1')
512 514 continue
513 515
514 516 # print ('stop 6')
515 517 firstpeak = powerlo + powerwidth/10.# first gaussian energy location
516 518 secondpeak = powerhi - powerwidth/10. #second gaussian energy location
517 519 midpeak = (firstpeak + secondpeak)/2.
518 520 firstamp = spcs[int(firstpeak)]
519 521 secondamp = spcs[int(secondpeak)]
520 522 midamp = spcs[int(midpeak)]
521 523
522 524 y_data = spc + wnoise
523 525
524 526 ''' single Gaussian '''
525 527 shift0 = numpy.mod(midpeak+minx, self.Num_Bin )
526 528 width0 = powerwidth/4.#Initialization entire power of spectrum divided by 4
527 529 power0 = 2.
528 530 amplitude0 = midamp
529 531 state0 = [shift0,width0,amplitude0,power0,wnoise]
530 532 bnds = ((0,self.Num_Bin-1),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
531 533 lsq1 = fmin_l_bfgs_b(self.misfit1, state0, args=(y_data,x,num_intg), bounds=bnds, approx_grad=True)
532 534 # print ('stop 7.1')
533 535 # print (bnds)
534 536
535 537 chiSq1=lsq1[1]
536 538
537 539 # print ('stop 8')
538 540 if fatspectra<1.0 and powerwidth<4:
539 541 choice=0
540 542 Amplitude0=lsq1[0][2]
541 543 shift0=lsq1[0][0]
542 544 width0=lsq1[0][1]
543 545 p0=lsq1[0][3]
544 546 Amplitude1=0.
545 547 shift1=0.
546 548 width1=0.
547 549 p1=0.
548 550 noise=lsq1[0][4]
549 551 #return (numpy.array([shift0,width0,Amplitude0,p0]),
550 552 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
551 553
552 554 # print ('stop 9')
553 555 ''' two Gaussians '''
554 556 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
555 557 shift0 = numpy.mod(firstpeak+minx, self.Num_Bin )
556 558 shift1 = numpy.mod(secondpeak+minx, self.Num_Bin )
557 559 width0 = powerwidth/6.
558 560 width1 = width0
559 561 power0 = 2.
560 562 power1 = power0
561 563 amplitude0 = firstamp
562 564 amplitude1 = secondamp
563 565 state0 = [shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
564 566 #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
565 567 bnds=((0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
566 568 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
567 569
568 570 # print ('stop 10')
569 571 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
570 572
571 573 # print ('stop 11')
572 574 chiSq2 = lsq2[1]
573 575
574 576 # print ('stop 12')
575 577
576 578 oneG = (chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
577 579
578 580 # print ('stop 13')
579 581 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
580 582 if oneG:
581 583 choice = 0
582 584 else:
583 585 w1 = lsq2[0][1]; w2 = lsq2[0][5]
584 586 a1 = lsq2[0][2]; a2 = lsq2[0][6]
585 587 p1 = lsq2[0][3]; p2 = lsq2[0][7]
586 588 s1 = (2**(1+1./p1))*scipy.special.gamma(1./p1)/p1
587 589 s2 = (2**(1+1./p2))*scipy.special.gamma(1./p2)/p2
588 590 gp1 = a1*w1*s1; gp2 = a2*w2*s2 # power content of each ggaussian with proper p scaling
589 591
590 592 if gp1>gp2:
591 593 if a1>0.7*a2:
592 594 choice = 1
593 595 else:
594 596 choice = 2
595 597 elif gp2>gp1:
596 598 if a2>0.7*a1:
597 599 choice = 2
598 600 else:
599 601 choice = 1
600 602 else:
601 603 choice = numpy.argmax([a1,a2])+1
602 604 #else:
603 605 #choice=argmin([std2a,std2b])+1
604 606
605 607 else: # with low SNR go to the most energetic peak
606 608 choice = numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
607 609
608 610 # print ('stop 14')
609 611 shift0 = lsq2[0][0]
610 612 vel0 = Vrange[0] + shift0 * deltav
611 613 shift1 = lsq2[0][4]
612 614 # vel1=Vrange[0] + shift1 * deltav
613 615
614 616 # max_vel = 1.0
615 617 # Va = max(Vrange)
616 618 # deltav = Vrange[1]-Vrange[0]
617 619 # print ('stop 15')
618 620 #first peak will be 0, second peak will be 1
619 621 # if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range # Commented by D.ScipiΓ³n 19.03.2021
620 622 if vel0 > -Va and vel0 < Va : #first peak is in the correct range
621 623 shift0 = lsq2[0][0]
622 624 width0 = lsq2[0][1]
623 625 Amplitude0 = lsq2[0][2]
624 626 p0 = lsq2[0][3]
625 627
626 628 shift1 = lsq2[0][4]
627 629 width1 = lsq2[0][5]
628 630 Amplitude1 = lsq2[0][6]
629 631 p1 = lsq2[0][7]
630 632 noise = lsq2[0][8]
631 633 else:
632 634 shift1 = lsq2[0][0]
633 635 width1 = lsq2[0][1]
634 636 Amplitude1 = lsq2[0][2]
635 637 p1 = lsq2[0][3]
636 638
637 639 shift0 = lsq2[0][4]
638 640 width0 = lsq2[0][5]
639 641 Amplitude0 = lsq2[0][6]
640 642 p0 = lsq2[0][7]
641 643 noise = lsq2[0][8]
642 644
643 645 if Amplitude0<0.05: # in case the peak is noise
644 646 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
645 647 if Amplitude1<0.05:
646 648 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
647 649
648 650 # print ('stop 16 ')
649 651 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0)/width0)**p0)
650 652 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1)/width1)**p1)
651 653 # SPCparam = (SPC_ch1,SPC_ch2)
652 654
653 655 DGauFitParam[0,ht,0] = noise
654 656 DGauFitParam[0,ht,1] = noise
655 657 DGauFitParam[1,ht,0] = Amplitude0
656 658 DGauFitParam[1,ht,1] = Amplitude1
657 659 DGauFitParam[2,ht,0] = Vrange[0] + shift0 * deltav
658 660 DGauFitParam[2,ht,1] = Vrange[0] + shift1 * deltav
659 661 DGauFitParam[3,ht,0] = width0 * deltav
660 662 DGauFitParam[3,ht,1] = width1 * deltav
661 663 DGauFitParam[4,ht,0] = p0
662 664 DGauFitParam[4,ht,1] = p1
663 665
664 666 # print (DGauFitParam.shape)
665 667 # print ('Leaving FitGau')
666 668 return DGauFitParam
667 669 # return SPCparam
668 670 # return GauSPC
669 671
670 672 def y_model1(self,x,state):
671 673 shift0, width0, amplitude0, power0, noise = state
672 674 model0 = amplitude0*numpy.exp(-0.5*abs((x - shift0)/width0)**power0)
673 675 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
674 676 model0d = amplitude0*numpy.exp(-0.5*abs((x - shift0 + self.Num_Bin)/width0)**power0)
675 677 return model0 + model0u + model0d + noise
676 678
677 679 def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
678 680 shift0, width0, amplitude0, power0, shift1, width1, amplitude1, power1, noise = state
679 681 model0 = amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
680 682 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
681 683 model0d = amplitude0*numpy.exp(-0.5*abs((x - shift0 + self.Num_Bin)/width0)**power0)
682 684
683 685 model1 = amplitude1*numpy.exp(-0.5*abs((x - shift1)/width1)**power1)
684 686 model1u = amplitude1*numpy.exp(-0.5*abs((x - shift1 - self.Num_Bin)/width1)**power1)
685 687 model1d = amplitude1*numpy.exp(-0.5*abs((x - shift1 + self.Num_Bin)/width1)**power1)
686 688 return model0 + model0u + model0d + model1 + model1u + model1d + noise
687 689
688 690 def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
689 691
690 692 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))**2)#/(64-5.) # /(64-5.) can be commented
691 693
692 694 def misfit2(self,state,y_data,x,num_intg):
693 695 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
694 696
695 697
696 698
697 699 class PrecipitationProc(Operation):
698 700
699 701 '''
700 702 Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R)
701 703
702 704 Input:
703 705 self.dataOut.data_pre : SelfSpectra
704 706
705 707 Output:
706 708
707 709 self.dataOut.data_output : Reflectivity factor, rainfall Rate
708 710
709 711
710 712 Parameters affected:
711 713 '''
712 714
713 715 def __init__(self):
714 716 Operation.__init__(self)
715 717 self.i=0
716 718
717 719 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
718 720 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350,SNRdBlimit=-30):
719 721
720 722 # print ('Entering PrecepitationProc ... ')
721 723
722 724 if radar == "MIRA35C" :
723 725
724 726 self.spc = dataOut.data_pre[0].copy()
725 727 self.Num_Hei = self.spc.shape[2]
726 728 self.Num_Bin = self.spc.shape[1]
727 729 self.Num_Chn = self.spc.shape[0]
728 730 Ze = self.dBZeMODE2(dataOut)
729 731
730 732 else:
731 733
732 734 self.spc = dataOut.data_pre[0].copy()
733 735
734 736 #NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX
735 737 self.spc[:,:,0:7]= numpy.NaN
736 738
737 739 self.Num_Hei = self.spc.shape[2]
738 740 self.Num_Bin = self.spc.shape[1]
739 741 self.Num_Chn = self.spc.shape[0]
740 742
741 743 VelRange = dataOut.spc_range[2]
742 744
743 745 ''' Se obtiene la constante del RADAR '''
744 746
745 747 self.Pt = Pt
746 748 self.Gt = Gt
747 749 self.Gr = Gr
748 750 self.Lambda = Lambda
749 751 self.aL = aL
750 752 self.tauW = tauW
751 753 self.ThetaT = ThetaT
752 754 self.ThetaR = ThetaR
753 755 self.GSys = 10**(36.63/10) # Ganancia de los LNA 36.63 dB
754 756 self.lt = 10**(1.67/10) # Perdida en cables Tx 1.67 dB
755 757 self.lr = 10**(5.73/10) # Perdida en cables Rx 5.73 dB
756 758
757 759 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
758 760 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
759 761 RadarConstant = 10e-26 * Numerator / Denominator #
760 762 ExpConstant = 10**(40/10) #Constante Experimental
761 763
762 764 SignalPower = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
763 765 for i in range(self.Num_Chn):
764 766 SignalPower[i,:,:] = self.spc[i,:,:] - dataOut.noise[i]
765 767 SignalPower[numpy.where(SignalPower < 0)] = 1e-20
766 768
767 769 SPCmean = numpy.mean(SignalPower, 0)
768 770 Pr = SPCmean[:,:]/dataOut.normFactor
769 771
770 772 # Declaring auxiliary variables
771 773 Range = dataOut.heightList*1000. #Range in m
772 774 # replicate the heightlist to obtain a matrix [Num_Bin,Num_Hei]
773 775 rMtrx = numpy.transpose(numpy.transpose([dataOut.heightList*1000.] * self.Num_Bin))
774 776 zMtrx = rMtrx+Altitude
775 777 # replicate the VelRange to obtain a matrix [Num_Bin,Num_Hei]
776 778 VelMtrx = numpy.transpose(numpy.tile(VelRange[:-1], (self.Num_Hei,1)))
777 779
778 780 # height dependence to air density Foote and Du Toit (1969)
779 781 delv_z = 1 + 3.68e-5 * zMtrx + 1.71e-9 * zMtrx**2
780 782 VMtrx = VelMtrx / delv_z #Normalized velocity
781 783 VMtrx[numpy.where(VMtrx> 9.6)] = numpy.NaN
782 784 # Diameter is related to the fall speed of falling drops
783 785 D_Vz = -1.667 * numpy.log( 0.9369 - 0.097087 * VMtrx ) # D in [mm]
784 786 # Only valid for D>= 0.16 mm
785 787 D_Vz[numpy.where(D_Vz < 0.16)] = numpy.NaN
786 788
787 789 #Calculate Radar Reflectivity ETAn
788 790 ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
789 791 ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
790 792 # Radar Cross Section
791 793 sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
792 794 # Drop Size Distribution
793 795 DSD = ETAn / sigmaD
794 796 # Equivalente Reflectivy
795 797 Ze_eqn = numpy.nansum( DSD * D_Vz**6 ,axis=0)
796 798 Ze_org = numpy.nansum(ETAn * Lambda**4, axis=0) / (1e-18*numpy.pi**5 * Km2) # [mm^6 /m^3]
797 799 # RainFall Rate
798 800 RR = 0.0006*numpy.pi * numpy.nansum( D_Vz**3 * DSD * VelMtrx ,0) #mm/hr
799 801
800 802 # Censoring the data
801 803 # Removing data with SNRth < 0dB se debe considerar el SNR por canal
802 804 SNRth = 10**(SNRdBlimit/10) #-30dB
803 805 novalid = numpy.where((dataOut.data_snr[0,:] <SNRth) | (dataOut.data_snr[1,:] <SNRth) | (dataOut.data_snr[2,:] <SNRth)) # AND condition. Maybe OR condition better
804 806 W = numpy.nanmean(dataOut.data_dop,0)
805 807 W[novalid] = numpy.NaN
806 808 Ze_org[novalid] = numpy.NaN
807 809 RR[novalid] = numpy.NaN
808 810
809 811 dataOut.data_output = RR[8]
810 812 dataOut.data_param = numpy.ones([3,self.Num_Hei])
811 813 dataOut.channelList = [0,1,2]
812 814
813 815 dataOut.data_param[0]=10*numpy.log10(Ze_org)
814 816 dataOut.data_param[1]=-W
815 817 dataOut.data_param[2]=RR
816 818
817 819 # print ('Leaving PrecepitationProc ... ')
818 820 return dataOut
819 821
820 822 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
821 823
822 824 NPW = dataOut.NPW
823 825 COFA = dataOut.COFA
824 826
825 827 SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]])
826 828 RadarConst = dataOut.RadarConst
827 829 #frequency = 34.85*10**9
828 830
829 831 ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei]))
830 832 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
831 833
832 834 ETA = numpy.sum(SNR,1)
833 835
834 836 ETA = numpy.where(ETA != 0. , ETA, numpy.NaN)
835 837
836 838 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
837 839
838 840 for r in range(self.Num_Hei):
839 841
840 842 Ze[0,r] = ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2)
841 843 #Ze[1,r] = ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2)
842 844
843 845 return Ze
844 846
845 847 # def GetRadarConstant(self):
846 848 #
847 849 # """
848 850 # Constants:
849 851 #
850 852 # Pt: Transmission Power dB 5kW 5000
851 853 # Gt: Transmission Gain dB 24.7 dB 295.1209
852 854 # Gr: Reception Gain dB 18.5 dB 70.7945
853 855 # Lambda: Wavelenght m 0.6741 m 0.6741
854 856 # aL: Attenuation loses dB 4dB 2.5118
855 857 # tauW: Width of transmission pulse s 4us 4e-6
856 858 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317
857 859 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087
858 860 #
859 861 # """
860 862 #
861 863 # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
862 864 # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
863 865 # RadarConstant = Numerator / Denominator
864 866 #
865 867 # return RadarConstant
866 868
867 869
868 870
869 871 class FullSpectralAnalysis(Operation):
870 872
871 873 """
872 874 Function that implements Full Spectral Analysis technique.
873 875
874 876 Input:
875 877 self.dataOut.data_pre : SelfSpectra and CrossSpectra data
876 878 self.dataOut.groupList : Pairlist of channels
877 879 self.dataOut.ChanDist : Physical distance between receivers
878 880
879 881
880 882 Output:
881 883
882 884 self.dataOut.data_output : Zonal wind, Meridional wind, and Vertical wind
883 885
884 886
885 887 Parameters affected: Winds, height range, SNR
886 888
887 889 """
888 890 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRdBlimit=-30,
889 891 minheight=None, maxheight=None, NegativeLimit=None, PositiveLimit=None):
890 892
891 893 spc = dataOut.data_pre[0].copy()
892 894 cspc = dataOut.data_pre[1]
893 895 nHeights = spc.shape[2]
894 896
895 897 # first_height = 0.75 #km (ref: data header 20170822)
896 898 # resolution_height = 0.075 #km
897 899 '''
898 900 finding height range. check this when radar parameters are changed!
899 901 '''
900 902 if maxheight is not None:
901 903 # range_max = math.ceil((maxheight - first_height) / resolution_height) # theoretical
902 904 range_max = math.ceil(13.26 * maxheight - 3) # empirical, works better
903 905 else:
904 906 range_max = nHeights
905 907 if minheight is not None:
906 908 # range_min = int((minheight - first_height) / resolution_height) # theoretical
907 909 range_min = int(13.26 * minheight - 5) # empirical, works better
908 910 if range_min < 0:
909 911 range_min = 0
910 912 else:
911 913 range_min = 0
912 914
913 915 pairsList = dataOut.groupList
914 916 if dataOut.ChanDist is not None :
915 917 ChanDist = dataOut.ChanDist
916 918 else:
917 919 ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
918 920
919 921 # 4 variables: zonal, meridional, vertical, and average SNR
920 922 data_param = numpy.zeros([4,nHeights]) * numpy.NaN
921 923 velocityX = numpy.zeros([nHeights]) * numpy.NaN
922 924 velocityY = numpy.zeros([nHeights]) * numpy.NaN
923 925 velocityZ = numpy.zeros([nHeights]) * numpy.NaN
924 926
925 927 dbSNR = 10*numpy.log10(numpy.average(dataOut.data_snr,0))
926 928
927 929 '''***********************************************WIND ESTIMATION**************************************'''
928 930 for Height in range(nHeights):
929 931
930 932 if Height >= range_min and Height < range_max:
931 933 # error_code will be useful in future analysis
932 934 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList,
933 935 ChanDist, Height, dataOut.noise, dataOut.spc_range, dbSNR[Height], SNRdBlimit, NegativeLimit, PositiveLimit,dataOut.frequency)
934 936
935 937 if abs(Vzon) < 100. and abs(Vmer) < 100.:
936 938 velocityX[Height] = Vzon
937 939 velocityY[Height] = -Vmer
938 940 velocityZ[Height] = Vver
939 941
940 942 # Censoring data with SNR threshold
941 943 dbSNR [dbSNR < SNRdBlimit] = numpy.NaN
942 944
943 945 data_param[0] = velocityX
944 946 data_param[1] = velocityY
945 947 data_param[2] = velocityZ
946 948 data_param[3] = dbSNR
947 949 dataOut.data_param = data_param
948 950 return dataOut
949 951
950 952 def moving_average(self,x, N=2):
951 953 """ convolution for smoothenig data. note that last N-1 values are convolution with zeroes """
952 954 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
953 955
954 956 def gaus(self,xSamples,Amp,Mu,Sigma):
955 957 return Amp * numpy.exp(-0.5*((xSamples - Mu)/Sigma)**2)
956 958
957 959 def Moments(self, ySamples, xSamples):
958 960 Power = numpy.nanmean(ySamples) # Power, 0th Moment
959 961 yNorm = ySamples / numpy.nansum(ySamples)
960 962 RadVel = numpy.nansum(xSamples * yNorm) # Radial Velocity, 1st Moment
961 963 Sigma2 = numpy.nansum(yNorm * (xSamples - RadVel)**2) # Spectral Width, 2nd Moment
962 964 StdDev = numpy.sqrt(numpy.abs(Sigma2)) # Desv. Estandar, Ancho espectral
963 965 return numpy.array([Power,RadVel,StdDev])
964 966
965 967 def StopWindEstimation(self, error_code):
966 968 Vzon = numpy.NaN
967 969 Vmer = numpy.NaN
968 970 Vver = numpy.NaN
969 971 return Vzon, Vmer, Vver, error_code
970 972
971 973 def AntiAliasing(self, interval, maxstep):
972 974 """
973 975 function to prevent errors from aliased values when computing phaseslope
974 976 """
975 977 antialiased = numpy.zeros(len(interval))
976 978 copyinterval = interval.copy()
977 979
978 980 antialiased[0] = copyinterval[0]
979 981
980 982 for i in range(1,len(antialiased)):
981 983 step = interval[i] - interval[i-1]
982 984 if step > maxstep:
983 985 copyinterval -= 2*numpy.pi
984 986 antialiased[i] = copyinterval[i]
985 987 elif step < maxstep*(-1):
986 988 copyinterval += 2*numpy.pi
987 989 antialiased[i] = copyinterval[i]
988 990 else:
989 991 antialiased[i] = copyinterval[i].copy()
990 992
991 993 return antialiased
992 994
993 995 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit, NegativeLimit, PositiveLimit, radfreq):
994 996 """
995 997 Function that Calculates Zonal, Meridional and Vertical wind velocities.
996 998 Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
997 999
998 1000 Input:
999 1001 spc, cspc : self spectra and cross spectra data. In Briggs notation something like S_i*(S_i)_conj, (S_j)_conj respectively.
1000 1002 pairsList : Pairlist of channels
1001 1003 ChanDist : array of xi_ij and eta_ij
1002 1004 Height : height at which data is processed
1003 1005 noise : noise in [channels] format for specific height
1004 1006 Abbsisarange : range of the frequencies or velocities
1005 1007 dbSNR, SNRlimit : signal to noise ratio in db, lower limit
1006 1008
1007 1009 Output:
1008 1010 Vzon, Vmer, Vver : wind velocities
1009 1011 error_code : int that states where code is terminated
1010 1012
1011 1013 0 : no error detected
1012 1014 1 : Gaussian of mean spc exceeds widthlimit
1013 1015 2 : no Gaussian of mean spc found
1014 1016 3 : SNR to low or velocity to high -> prec. e.g.
1015 1017 4 : at least one Gaussian of cspc exceeds widthlimit
1016 1018 5 : zero out of three cspc Gaussian fits converged
1017 1019 6 : phase slope fit could not be found
1018 1020 7 : arrays used to fit phase have different length
1019 1021 8 : frequency range is either too short (len <= 5) or very long (> 30% of cspc)
1020 1022
1021 1023 """
1022 1024
1023 1025 error_code = 0
1024 1026
1025 1027 nChan = spc.shape[0]
1026 1028 nProf = spc.shape[1]
1027 1029 nPair = cspc.shape[0]
1028 1030
1029 1031 SPC_Samples = numpy.zeros([nChan, nProf]) # for normalized spc values for one height
1030 1032 CSPC_Samples = numpy.zeros([nPair, nProf], dtype=numpy.complex_) # for normalized cspc values
1031 1033 phase = numpy.zeros([nPair, nProf]) # phase between channels
1032 1034 PhaseSlope = numpy.zeros(nPair) # slope of the phases, channelwise
1033 1035 PhaseInter = numpy.zeros(nPair) # intercept to the slope of the phases, channelwise
1034 1036 xFrec = AbbsisaRange[0][:-1] # frequency range
1035 1037 xVel = AbbsisaRange[2][:-1] # velocity range
1036 1038 xSamples = xFrec # the frequency range is taken
1037 1039 delta_x = xSamples[1] - xSamples[0] # delta_f or delta_x
1038 1040
1039 1041 # only consider velocities with in NegativeLimit and PositiveLimit
1040 1042 if (NegativeLimit is None):
1041 1043 NegativeLimit = numpy.min(xVel)
1042 1044 if (PositiveLimit is None):
1043 1045 PositiveLimit = numpy.max(xVel)
1044 1046 xvalid = numpy.where((xVel > NegativeLimit) & (xVel < PositiveLimit))
1045 1047 xSamples_zoom = xSamples[xvalid]
1046 1048
1047 1049 '''Getting Eij and Nij'''
1048 1050 Xi01, Xi02, Xi12 = ChanDist[:,0]
1049 1051 Eta01, Eta02, Eta12 = ChanDist[:,1]
1050 1052
1051 1053 # spwd limit - updated by D. ScipiΓ³n 30.03.2021
1052 1054 widthlimit = 10
1053 1055 '''************************* SPC is normalized ********************************'''
1054 1056 spc_norm = spc.copy()
1055 1057 # For each channel
1056 1058 for i in range(nChan):
1057 1059 spc_sub = spc_norm[i,:] - noise[i] # only the signal power
1058 1060 SPC_Samples[i] = spc_sub / (numpy.nansum(spc_sub) * delta_x)
1059 1061
1060 1062 '''********************** FITTING MEAN SPC GAUSSIAN **********************'''
1061 1063
1062 1064 """ the gaussian of the mean: first subtract noise, then normalize. this is legal because
1063 1065 you only fit the curve and don't need the absolute value of height for calculation,
1064 1066 only for estimation of width. for normalization of cross spectra, you need initial,
1065 1067 unnormalized self-spectra With noise.
1066 1068
1067 1069 Technically, you don't even need to normalize the self-spectra, as you only need the
1068 1070 width of the peak. However, it was left this way. Note that the normalization has a flaw:
1069 1071 due to subtraction of the noise, some values are below zero. Raw "spc" values should be
1070 1072 >= 0, as it is the modulus squared of the signals (complex * it's conjugate)
1071 1073 """
1072 1074 # initial conditions
1073 1075 popt = [1e-10,0,1e-10]
1074 1076 # Spectra average
1075 1077 SPCMean = numpy.average(SPC_Samples,0)
1076 1078 # Moments in frequency
1077 1079 SPCMoments = self.Moments(SPCMean[xvalid], xSamples_zoom)
1078 1080
1079 1081 # Gauss Fit SPC in frequency domain
1080 1082 if dbSNR > SNRlimit: # only if SNR > SNRth
1081 1083 try:
1082 1084 popt,pcov = curve_fit(self.gaus,xSamples_zoom,SPCMean[xvalid],p0=SPCMoments)
1083 1085 if popt[2] <= 0 or popt[2] > widthlimit: # CONDITION
1084 1086 return self.StopWindEstimation(error_code = 1)
1085 1087 FitGauss = self.gaus(xSamples_zoom,*popt)
1086 1088 except :#RuntimeError:
1087 1089 return self.StopWindEstimation(error_code = 2)
1088 1090 else:
1089 1091 return self.StopWindEstimation(error_code = 3)
1090 1092
1091 1093 '''***************************** CSPC Normalization *************************
1092 1094 The Spc spectra are used to normalize the crossspectra. Peaks from precipitation
1093 1095 influence the norm which is not desired. First, a range is identified where the
1094 1096 wind peak is estimated -> sum_wind is sum of those frequencies. Next, the area
1095 1097 around it gets cut off and values replaced by mean determined by the boundary
1096 1098 data -> sum_noise (spc is not normalized here, thats why the noise is important)
1097 1099
1098 1100 The sums are then added and multiplied by range/datapoints, because you need
1099 1101 an integral and not a sum for normalization.
1100 1102
1101 1103 A norm is found according to Briggs 92.
1102 1104 '''
1103 1105 # for each pair
1104 1106 for i in range(nPair):
1105 1107 cspc_norm = cspc[i,:].copy()
1106 1108 chan_index0 = pairsList[i][0]
1107 1109 chan_index1 = pairsList[i][1]
1108 1110 CSPC_Samples[i] = cspc_norm / (numpy.sqrt(numpy.nansum(spc_norm[chan_index0])*numpy.nansum(spc_norm[chan_index1])) * delta_x)
1109 1111 phase[i] = numpy.arctan2(CSPC_Samples[i].imag, CSPC_Samples[i].real)
1110 1112
1111 1113 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0,xvalid]), xSamples_zoom),
1112 1114 self.Moments(numpy.abs(CSPC_Samples[1,xvalid]), xSamples_zoom),
1113 1115 self.Moments(numpy.abs(CSPC_Samples[2,xvalid]), xSamples_zoom)])
1114 1116
1115 1117 popt01, popt02, popt12 = [1e-10,0,1e-10], [1e-10,0,1e-10] ,[1e-10,0,1e-10]
1116 1118 FitGauss01, FitGauss02, FitGauss12 = numpy.zeros(len(xSamples)), numpy.zeros(len(xSamples)), numpy.zeros(len(xSamples))
1117 1119
1118 1120 '''*******************************FIT GAUSS CSPC************************************'''
1119 1121 try:
1120 1122 popt01,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[0][xvalid]),p0=CSPCmoments[0])
1121 1123 if popt01[2] > widthlimit: # CONDITION
1122 1124 return self.StopWindEstimation(error_code = 4)
1123 1125 popt02,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[1][xvalid]),p0=CSPCmoments[1])
1124 1126 if popt02[2] > widthlimit: # CONDITION
1125 1127 return self.StopWindEstimation(error_code = 4)
1126 1128 popt12,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[2][xvalid]),p0=CSPCmoments[2])
1127 1129 if popt12[2] > widthlimit: # CONDITION
1128 1130 return self.StopWindEstimation(error_code = 4)
1129 1131
1130 1132 FitGauss01 = self.gaus(xSamples_zoom, *popt01)
1131 1133 FitGauss02 = self.gaus(xSamples_zoom, *popt02)
1132 1134 FitGauss12 = self.gaus(xSamples_zoom, *popt12)
1133 1135 except:
1134 1136 return self.StopWindEstimation(error_code = 5)
1135 1137
1136 1138
1137 1139 '''************* Getting Fij ***************'''
1138 1140 # x-axis point of the gaussian where the center is located from GaussFit of spectra
1139 1141 GaussCenter = popt[1]
1140 1142 ClosestCenter = xSamples_zoom[numpy.abs(xSamples_zoom-GaussCenter).argmin()]
1141 1143 PointGauCenter = numpy.where(xSamples_zoom==ClosestCenter)[0][0]
1142 1144
1143 1145 # Point where e^-1 is located in the gaussian
1144 1146 PeMinus1 = numpy.max(FitGauss) * numpy.exp(-1)
1145 1147 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # The closest point to"Peminus1" in "FitGauss"
1146 1148 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1147 1149 Fij = numpy.abs(xSamples_zoom[PointFij] - xSamples_zoom[PointGauCenter])
1148 1150
1149 1151 '''********** Taking frequency ranges from mean SPCs **********'''
1150 1152 GauWidth = popt[2] * 3/2 # Bandwidth of Gau01
1151 1153 Range = numpy.empty(2)
1152 1154 Range[0] = GaussCenter - GauWidth
1153 1155 Range[1] = GaussCenter + GauWidth
1154 1156 # Point in x-axis where the bandwidth is located (min:max)
1155 1157 ClosRangeMin = xSamples_zoom[numpy.abs(xSamples_zoom-Range[0]).argmin()]
1156 1158 ClosRangeMax = xSamples_zoom[numpy.abs(xSamples_zoom-Range[1]).argmin()]
1157 1159 PointRangeMin = numpy.where(xSamples_zoom==ClosRangeMin)[0][0]
1158 1160 PointRangeMax = numpy.where(xSamples_zoom==ClosRangeMax)[0][0]
1159 1161 Range = numpy.array([ PointRangeMin, PointRangeMax ])
1160 1162 FrecRange = xSamples_zoom[ Range[0] : Range[1] ]
1161 1163
1162 1164 '''************************** Getting Phase Slope ***************************'''
1163 1165 for i in range(nPair):
1164 1166 if len(FrecRange) > 5:
1165 1167 PhaseRange = phase[i, xvalid[0][Range[0]:Range[1]]].copy()
1166 1168 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1167 1169 if len(FrecRange) == len(PhaseRange):
1168 1170 try:
1169 1171 slope, intercept, _, _, _ = stats.linregress(FrecRange[mask], self.AntiAliasing(PhaseRange[mask], 4.5))
1170 1172 PhaseSlope[i] = slope
1171 1173 PhaseInter[i] = intercept
1172 1174 except:
1173 1175 return self.StopWindEstimation(error_code = 6)
1174 1176 else:
1175 1177 return self.StopWindEstimation(error_code = 7)
1176 1178 else:
1177 1179 return self.StopWindEstimation(error_code = 8)
1178 1180
1179 1181 '''*** Constants A-H correspond to the convention as in Briggs and Vincent 1992 ***'''
1180 1182
1181 1183 '''Getting constant C'''
1182 1184 cC=(Fij*numpy.pi)**2
1183 1185
1184 1186 '''****** Getting constants F and G ******'''
1185 1187 MijEijNij = numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1186 1188 # MijEijNij = numpy.array([[Xi01,Eta01], [Xi02,Eta02], [Xi12,Eta12]])
1187 1189 # MijResult0 = (-PhaseSlope[0] * cC) / (2*numpy.pi)
1188 1190 MijResult1 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
1189 1191 MijResult2 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
1190 1192 # MijResults = numpy.array([MijResult0, MijResult1, MijResult2])
1191 1193 MijResults = numpy.array([MijResult1, MijResult2])
1192 1194 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1193 1195
1194 1196 '''****** Getting constants A, B and H ******'''
1195 1197 W01 = numpy.nanmax( FitGauss01 )
1196 1198 W02 = numpy.nanmax( FitGauss02 )
1197 1199 W12 = numpy.nanmax( FitGauss12 )
1198 1200
1199 1201 WijResult01 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
1200 1202 WijResult02 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC))
1201 1203 WijResult12 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC))
1202 1204 WijResults = numpy.array([WijResult01, WijResult02, WijResult12])
1203 1205
1204 1206 WijEijNij = numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1205 1207 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1206 1208
1207 1209 VxVy = numpy.array([[cA,cH],[cH,cB]])
1208 1210 VxVyResults = numpy.array([-cF,-cG])
1209 1211 (Vmer,Vzon) = numpy.linalg.solve(VxVy, VxVyResults)
1210 1212 Vver = -SPCMoments[1]*SPEED_OF_LIGHT/(2*radfreq)
1211 1213 error_code = 0
1212 1214
1213 1215 return Vzon, Vmer, Vver, error_code
1214 1216
1215 1217 class SpectralMoments(Operation):
1216 1218
1217 1219 '''
1218 1220 Function SpectralMoments()
1219 1221
1220 1222 Calculates moments (power, mean, standard deviation) and SNR of the signal
1221 1223
1222 1224 Type of dataIn: Spectra
1223 1225
1224 1226 Configuration Parameters:
1225 1227
1226 1228 dirCosx : Cosine director in X axis
1227 1229 dirCosy : Cosine director in Y axis
1228 1230
1229 1231 elevation :
1230 1232 azimuth :
1231 1233
1232 1234 Input:
1233 1235 channelList : simple channel list to select e.g. [2,3,7]
1234 1236 self.dataOut.data_pre : Spectral data
1235 1237 self.dataOut.abscissaList : List of frequencies
1236 1238 self.dataOut.noise : Noise level per channel
1237 1239
1238 1240 Affected:
1239 1241 self.dataOut.moments : Parameters per channel
1240 1242 self.dataOut.data_snr : SNR per channel
1241 1243
1242 1244 '''
1243 1245
1244 1246 def run(self, dataOut):
1245 1247
1246 1248 data = dataOut.data_pre[0]
1247 1249 absc = dataOut.abscissaList[:-1]
1248 1250 noise = dataOut.noise
1249 1251 nChannel = data.shape[0]
1250 1252 data_param = numpy.zeros((nChannel, 4, data.shape[2]))
1251 1253
1252 1254 for ind in range(nChannel):
1253 1255 data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
1254 1256
1255 1257 dataOut.moments = data_param[:,1:,:]
1256 1258 dataOut.data_snr = data_param[:,0]
1257 1259 dataOut.data_pow = data_param[:,1]
1258 1260 dataOut.data_dop = data_param[:,2]
1259 1261 dataOut.data_width = data_param[:,3]
1260 1262
1261 1263 return dataOut
1262 1264
1263 1265 def __calculateMoments(self, oldspec, oldfreq, n0,
1264 1266 nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
1265 1267
1266 1268 if (nicoh is None): nicoh = 1
1267 1269 if (graph is None): graph = 0
1268 1270 if (smooth is None): smooth = 0
1269 1271 elif (self.smooth < 3): smooth = 0
1270 1272
1271 1273 if (type1 is None): type1 = 0
1272 1274 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
1273 1275 if (snrth is None): snrth = -3
1274 1276 if (dc is None): dc = 0
1275 1277 if (aliasing is None): aliasing = 0
1276 1278 if (oldfd is None): oldfd = 0
1277 1279 if (wwauto is None): wwauto = 0
1278 1280
1279 1281 if (n0 < 1.e-20): n0 = 1.e-20
1280 1282
1281 1283 freq = oldfreq
1282 1284 vec_power = numpy.zeros(oldspec.shape[1])
1283 1285 vec_fd = numpy.zeros(oldspec.shape[1])
1284 1286 vec_w = numpy.zeros(oldspec.shape[1])
1285 1287 vec_snr = numpy.zeros(oldspec.shape[1])
1286 1288
1287 1289 # oldspec = numpy.ma.masked_invalid(oldspec)
1288 1290
1289 1291 for ind in range(oldspec.shape[1]):
1290 1292
1291 1293 spec = oldspec[:,ind]
1292 1294 aux = spec*fwindow
1293 1295 max_spec = aux.max()
1294 1296 m = aux.tolist().index(max_spec)
1295 1297
1296 1298 # Smooth
1297 1299 if (smooth == 0):
1298 1300 spec2 = spec
1299 1301 else:
1300 1302 spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
1301 1303
1302 1304 # Moments Estimation
1303 1305 bb = spec2[numpy.arange(m,spec2.size)]
1304 1306 bb = (bb<n0).nonzero()
1305 1307 bb = bb[0]
1306 1308
1307 1309 ss = spec2[numpy.arange(0,m + 1)]
1308 1310 ss = (ss<n0).nonzero()
1309 1311 ss = ss[0]
1310 1312
1311 1313 if (bb.size == 0):
1312 1314 bb0 = spec.size - 1 - m
1313 1315 else:
1314 1316 bb0 = bb[0] - 1
1315 1317 if (bb0 < 0):
1316 1318 bb0 = 0
1317 1319
1318 1320 if (ss.size == 0):
1319 1321 ss1 = 1
1320 1322 else:
1321 1323 ss1 = max(ss) + 1
1322 1324
1323 1325 if (ss1 > m):
1324 1326 ss1 = m
1325 1327
1326 1328 valid = numpy.arange(int(m + bb0 - ss1 + 1)) + ss1
1327 1329
1328 1330 signal_power = ((spec2[valid] - n0) * fwindow[valid]).mean() # D. ScipiΓ³n added with correct definition
1329 1331 total_power = (spec2[valid] * fwindow[valid]).mean() # D. ScipiΓ³n added with correct definition
1330 1332 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1331 1333 fd = ((spec2[valid]- n0)*freq[valid] * fwindow[valid]).sum() / power
1332 1334 w = numpy.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum() / power)
1333 1335 snr = (spec2.mean()-n0)/n0
1334 1336 if (snr < 1.e-20) :
1335 1337 snr = 1.e-20
1336 1338
1337 1339 # vec_power[ind] = power #D. ScipiΓ³n replaced with the line below
1338 1340 vec_power[ind] = total_power
1339 1341 vec_fd[ind] = fd
1340 1342 vec_w[ind] = w
1341 1343 vec_snr[ind] = snr
1342 1344
1343 1345 return numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1344 1346
1345 1347 #------------------ Get SA Parameters --------------------------
1346 1348
1347 1349 def GetSAParameters(self):
1348 1350 #SA en frecuencia
1349 1351 pairslist = self.dataOut.groupList
1350 1352 num_pairs = len(pairslist)
1351 1353
1352 1354 vel = self.dataOut.abscissaList
1353 1355 spectra = self.dataOut.data_pre
1354 1356 cspectra = self.dataIn.data_cspc
1355 1357 delta_v = vel[1] - vel[0]
1356 1358
1357 1359 #Calculating the power spectrum
1358 1360 spc_pow = numpy.sum(spectra, 3)*delta_v
1359 1361 #Normalizing Spectra
1360 1362 norm_spectra = spectra/spc_pow
1361 1363 #Calculating the norm_spectra at peak
1362 1364 max_spectra = numpy.max(norm_spectra, 3)
1363 1365
1364 1366 #Normalizing Cross Spectra
1365 1367 norm_cspectra = numpy.zeros(cspectra.shape)
1366 1368
1367 1369 for i in range(num_chan):
1368 1370 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
1369 1371
1370 1372 max_cspectra = numpy.max(norm_cspectra,2)
1371 1373 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
1372 1374
1373 1375 for i in range(num_pairs):
1374 1376 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
1375 1377 #------------------- Get Lags ----------------------------------
1376 1378
1377 1379 class SALags(Operation):
1378 1380 '''
1379 1381 Function GetMoments()
1380 1382
1381 1383 Input:
1382 1384 self.dataOut.data_pre
1383 1385 self.dataOut.abscissaList
1384 1386 self.dataOut.noise
1385 1387 self.dataOut.normFactor
1386 1388 self.dataOut.data_snr
1387 1389 self.dataOut.groupList
1388 1390 self.dataOut.nChannels
1389 1391
1390 1392 Affected:
1391 1393 self.dataOut.data_param
1392 1394
1393 1395 '''
1394 1396 def run(self, dataOut):
1395 1397 data_acf = dataOut.data_pre[0]
1396 1398 data_ccf = dataOut.data_pre[1]
1397 1399 normFactor_acf = dataOut.normFactor[0]
1398 1400 normFactor_ccf = dataOut.normFactor[1]
1399 1401 pairs_acf = dataOut.groupList[0]
1400 1402 pairs_ccf = dataOut.groupList[1]
1401 1403
1402 1404 nHeights = dataOut.nHeights
1403 1405 absc = dataOut.abscissaList
1404 1406 noise = dataOut.noise
1405 1407 SNR = dataOut.data_snr
1406 1408 nChannels = dataOut.nChannels
1407 1409 # pairsList = dataOut.groupList
1408 1410 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1409 1411
1410 1412 for l in range(len(pairs_acf)):
1411 1413 data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
1412 1414
1413 1415 for l in range(len(pairs_ccf)):
1414 1416 data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
1415 1417
1416 1418 dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
1417 1419 dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
1418 1420 dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
1419 1421 return
1420 1422
1421 1423 # def __getPairsAutoCorr(self, pairsList, nChannels):
1422 1424 #
1423 1425 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1424 1426 #
1425 1427 # for l in range(len(pairsList)):
1426 1428 # firstChannel = pairsList[l][0]
1427 1429 # secondChannel = pairsList[l][1]
1428 1430 #
1429 1431 # #Obteniendo pares de Autocorrelacion
1430 1432 # if firstChannel == secondChannel:
1431 1433 # pairsAutoCorr[firstChannel] = int(l)
1432 1434 #
1433 1435 # pairsAutoCorr = pairsAutoCorr.astype(int)
1434 1436 #
1435 1437 # pairsCrossCorr = range(len(pairsList))
1436 1438 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1437 1439 #
1438 1440 # return pairsAutoCorr, pairsCrossCorr
1439 1441
1440 1442 def __calculateTaus(self, data_acf, data_ccf, lagRange):
1441 1443
1442 1444 lag0 = data_acf.shape[1]/2
1443 1445 #Funcion de Autocorrelacion
1444 1446 mean_acf = stats.nanmean(data_acf, axis = 0)
1445 1447
1446 1448 #Obtencion Indice de TauCross
1447 1449 ind_ccf = data_ccf.argmax(axis = 1)
1448 1450 #Obtencion Indice de TauAuto
1449 1451 ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
1450 1452 ccf_lag0 = data_ccf[:,lag0,:]
1451 1453
1452 1454 for i in range(ccf_lag0.shape[0]):
1453 1455 ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
1454 1456
1455 1457 #Obtencion de TauCross y TauAuto
1456 1458 tau_ccf = lagRange[ind_ccf]
1457 1459 tau_acf = lagRange[ind_acf]
1458 1460
1459 1461 Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
1460 1462
1461 1463 tau_ccf[Nan1,Nan2] = numpy.nan
1462 1464 tau_acf[Nan1,Nan2] = numpy.nan
1463 1465 tau = numpy.vstack((tau_ccf,tau_acf))
1464 1466
1465 1467 return tau
1466 1468
1467 1469 def __calculateLag1Phase(self, data, lagTRange):
1468 1470 data1 = stats.nanmean(data, axis = 0)
1469 1471 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
1470 1472
1471 1473 phase = numpy.angle(data1[lag1,:])
1472 1474
1473 1475 return phase
1474 1476
1475 1477 class SpectralFitting(Operation):
1476 1478 '''
1477 1479 Function GetMoments()
1478 1480
1479 1481 Input:
1480 1482 Output:
1481 1483 Variables modified:
1482 1484 '''
1483 1485
1484 1486 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
1485 1487
1486 1488
1487 1489 if path != None:
1488 1490 sys.path.append(path)
1489 1491 self.dataOut.library = importlib.import_module(file)
1490 1492
1491 1493 #To be inserted as a parameter
1492 1494 groupArray = numpy.array(groupList)
1493 1495 # groupArray = numpy.array([[0,1],[2,3]])
1494 1496 self.dataOut.groupList = groupArray
1495 1497
1496 1498 nGroups = groupArray.shape[0]
1497 1499 nChannels = self.dataIn.nChannels
1498 1500 nHeights=self.dataIn.heightList.size
1499 1501
1500 1502 #Parameters Array
1501 1503 self.dataOut.data_param = None
1502 1504
1503 1505 #Set constants
1504 1506 constants = self.dataOut.library.setConstants(self.dataIn)
1505 1507 self.dataOut.constants = constants
1506 1508 M = self.dataIn.normFactor
1507 1509 N = self.dataIn.nFFTPoints
1508 1510 ippSeconds = self.dataIn.ippSeconds
1509 1511 K = self.dataIn.nIncohInt
1510 1512 pairsArray = numpy.array(self.dataIn.pairsList)
1511 1513
1512 1514 #List of possible combinations
1513 1515 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1514 1516 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1515 1517
1516 1518 if getSNR:
1517 1519 listChannels = groupArray.reshape((groupArray.size))
1518 1520 listChannels.sort()
1519 1521 noise = self.dataIn.getNoise()
1520 1522 self.dataOut.data_snr = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1521 1523
1522 1524 for i in range(nGroups):
1523 1525 coord = groupArray[i,:]
1524 1526
1525 1527 #Input data array
1526 1528 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1527 1529 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1528 1530
1529 1531 #Cross Spectra data array for Covariance Matrixes
1530 1532 ind = 0
1531 1533 for pairs in listComb:
1532 1534 pairsSel = numpy.array([coord[x],coord[y]])
1533 1535 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1534 1536 ind += 1
1535 1537 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1536 1538 dataCross = dataCross**2/K
1537 1539
1538 1540 for h in range(nHeights):
1539 1541
1540 1542 #Input
1541 1543 d = data[:,h]
1542 1544
1543 1545 #Covariance Matrix
1544 1546 D = numpy.diag(d**2/K)
1545 1547 ind = 0
1546 1548 for pairs in listComb:
1547 1549 #Coordinates in Covariance Matrix
1548 1550 x = pairs[0]
1549 1551 y = pairs[1]
1550 1552 #Channel Index
1551 1553 S12 = dataCross[ind,:,h]
1552 1554 D12 = numpy.diag(S12)
1553 1555 #Completing Covariance Matrix with Cross Spectras
1554 1556 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1555 1557 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1556 1558 ind += 1
1557 1559 Dinv=numpy.linalg.inv(D)
1558 1560 L=numpy.linalg.cholesky(Dinv)
1559 1561 LT=L.T
1560 1562
1561 1563 dp = numpy.dot(LT,d)
1562 1564
1563 1565 #Initial values
1564 1566 data_spc = self.dataIn.data_spc[coord,:,h]
1565 1567
1566 1568 if (h>0)and(error1[3]<5):
1567 1569 p0 = self.dataOut.data_param[i,:,h-1]
1568 1570 else:
1569 1571 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1570 1572
1571 1573 try:
1572 1574 #Least Squares
1573 1575 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1574 1576 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1575 1577 #Chi square error
1576 1578 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1577 1579 #Error with Jacobian
1578 1580 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1579 1581 except:
1580 1582 minp = p0*numpy.nan
1581 1583 error0 = numpy.nan
1582 1584 error1 = p0*numpy.nan
1583 1585
1584 1586 #Save
1585 1587 if self.dataOut.data_param is None:
1586 1588 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1587 1589 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1588 1590
1589 1591 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1590 1592 self.dataOut.data_param[i,:,h] = minp
1591 1593 return
1592 1594
1593 1595 def __residFunction(self, p, dp, LT, constants):
1594 1596
1595 1597 fm = self.dataOut.library.modelFunction(p, constants)
1596 1598 fmp=numpy.dot(LT,fm)
1597 1599
1598 1600 return dp-fmp
1599 1601
1600 1602 def __getSNR(self, z, noise):
1601 1603
1602 1604 avg = numpy.average(z, axis=1)
1603 1605 SNR = (avg.T-noise)/noise
1604 1606 SNR = SNR.T
1605 1607 return SNR
1606 1608
1607 1609 def __chisq(p,chindex,hindex):
1608 1610 #similar to Resid but calculates CHI**2
1609 1611 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1610 1612 dp=numpy.dot(LT,d)
1611 1613 fmp=numpy.dot(LT,fm)
1612 1614 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1613 1615 return chisq
1614 1616
1615 1617 class WindProfiler(Operation):
1616 1618
1617 1619 __isConfig = False
1618 1620
1619 1621 __initime = None
1620 1622 __lastdatatime = None
1621 1623 __integrationtime = None
1622 1624
1623 1625 __buffer = None
1624 1626
1625 1627 __dataReady = False
1626 1628
1627 1629 __firstdata = None
1628 1630
1629 1631 n = None
1630 1632
1631 1633 def __init__(self):
1632 1634 Operation.__init__(self)
1633 1635
1634 1636 def __calculateCosDir(self, elev, azim):
1635 1637 zen = (90 - elev)*numpy.pi/180
1636 1638 azim = azim*numpy.pi/180
1637 1639 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1638 1640 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1639 1641
1640 1642 signX = numpy.sign(numpy.cos(azim))
1641 1643 signY = numpy.sign(numpy.sin(azim))
1642 1644
1643 1645 cosDirX = numpy.copysign(cosDirX, signX)
1644 1646 cosDirY = numpy.copysign(cosDirY, signY)
1645 1647 return cosDirX, cosDirY
1646 1648
1647 1649 def __calculateAngles(self, theta_x, theta_y, azimuth):
1648 1650
1649 1651 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1650 1652 zenith_arr = numpy.arccos(dir_cosw)
1651 1653 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1652 1654
1653 1655 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1654 1656 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1655 1657
1656 1658 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1657 1659
1658 1660 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1659 1661
1660 1662 #
1661 1663 if horOnly:
1662 1664 A = numpy.c_[dir_cosu,dir_cosv]
1663 1665 else:
1664 1666 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1665 1667 A = numpy.asmatrix(A)
1666 1668 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1667 1669
1668 1670 return A1
1669 1671
1670 1672 def __correctValues(self, heiRang, phi, velRadial, SNR):
1671 1673 listPhi = phi.tolist()
1672 1674 maxid = listPhi.index(max(listPhi))
1673 1675 minid = listPhi.index(min(listPhi))
1674 1676
1675 1677 rango = list(range(len(phi)))
1676 1678 # rango = numpy.delete(rango,maxid)
1677 1679
1678 1680 heiRang1 = heiRang*math.cos(phi[maxid])
1679 1681 heiRangAux = heiRang*math.cos(phi[minid])
1680 1682 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1681 1683 heiRang1 = numpy.delete(heiRang1,indOut)
1682 1684
1683 1685 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1684 1686 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1685 1687
1686 1688 for i in rango:
1687 1689 x = heiRang*math.cos(phi[i])
1688 1690 y1 = velRadial[i,:]
1689 1691 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1690 1692
1691 1693 x1 = heiRang1
1692 1694 y11 = f1(x1)
1693 1695
1694 1696 y2 = SNR[i,:]
1695 1697 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1696 1698 y21 = f2(x1)
1697 1699
1698 1700 velRadial1[i,:] = y11
1699 1701 SNR1[i,:] = y21
1700 1702
1701 1703 return heiRang1, velRadial1, SNR1
1702 1704
1703 1705 def __calculateVelUVW(self, A, velRadial):
1704 1706
1705 1707 #Operacion Matricial
1706 1708 # velUVW = numpy.zeros((velRadial.shape[1],3))
1707 1709 # for ind in range(velRadial.shape[1]):
1708 1710 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1709 1711 # velUVW = velUVW.transpose()
1710 1712 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1711 1713 velUVW[:,:] = numpy.dot(A,velRadial)
1712 1714
1713 1715
1714 1716 return velUVW
1715 1717
1716 1718 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1717 1719
1718 1720 def techniqueDBS(self, kwargs):
1719 1721 """
1720 1722 Function that implements Doppler Beam Swinging (DBS) technique.
1721 1723
1722 1724 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1723 1725 Direction correction (if necessary), Ranges and SNR
1724 1726
1725 1727 Output: Winds estimation (Zonal, Meridional and Vertical)
1726 1728
1727 1729 Parameters affected: Winds, height range, SNR
1728 1730 """
1729 1731 velRadial0 = kwargs['velRadial']
1730 1732 heiRang = kwargs['heightList']
1731 1733 SNR0 = kwargs['SNR']
1732 1734
1733 1735 if 'dirCosx' in kwargs and 'dirCosy' in kwargs:
1734 1736 theta_x = numpy.array(kwargs['dirCosx'])
1735 1737 theta_y = numpy.array(kwargs['dirCosy'])
1736 1738 else:
1737 1739 elev = numpy.array(kwargs['elevation'])
1738 1740 azim = numpy.array(kwargs['azimuth'])
1739 1741 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1740 1742 azimuth = kwargs['correctAzimuth']
1741 1743 if 'horizontalOnly' in kwargs:
1742 1744 horizontalOnly = kwargs['horizontalOnly']
1743 1745 else: horizontalOnly = False
1744 1746 if 'correctFactor' in kwargs:
1745 1747 correctFactor = kwargs['correctFactor']
1746 1748 else: correctFactor = 1
1747 1749 if 'channelList' in kwargs:
1748 1750 channelList = kwargs['channelList']
1749 1751 if len(channelList) == 2:
1750 1752 horizontalOnly = True
1751 1753 arrayChannel = numpy.array(channelList)
1752 1754 param = param[arrayChannel,:,:]
1753 1755 theta_x = theta_x[arrayChannel]
1754 1756 theta_y = theta_y[arrayChannel]
1755 1757
1756 1758 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
1757 1759 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
1758 1760 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1759 1761
1760 1762 #Calculo de Componentes de la velocidad con DBS
1761 1763 winds = self.__calculateVelUVW(A,velRadial1)
1762 1764
1763 1765 return winds, heiRang1, SNR1
1764 1766
1765 1767 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
1766 1768
1767 1769 nPairs = len(pairs_ccf)
1768 1770 posx = numpy.asarray(posx)
1769 1771 posy = numpy.asarray(posy)
1770 1772
1771 1773 #Rotacion Inversa para alinear con el azimuth
1772 1774 if azimuth!= None:
1773 1775 azimuth = azimuth*math.pi/180
1774 1776 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1775 1777 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1776 1778 else:
1777 1779 posx1 = posx
1778 1780 posy1 = posy
1779 1781
1780 1782 #Calculo de Distancias
1781 1783 distx = numpy.zeros(nPairs)
1782 1784 disty = numpy.zeros(nPairs)
1783 1785 dist = numpy.zeros(nPairs)
1784 1786 ang = numpy.zeros(nPairs)
1785 1787
1786 1788 for i in range(nPairs):
1787 1789 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
1788 1790 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
1789 1791 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1790 1792 ang[i] = numpy.arctan2(disty[i],distx[i])
1791 1793
1792 1794 return distx, disty, dist, ang
1793 1795 #Calculo de Matrices
1794 1796 # nPairs = len(pairs)
1795 1797 # ang1 = numpy.zeros((nPairs, 2, 1))
1796 1798 # dist1 = numpy.zeros((nPairs, 2, 1))
1797 1799 #
1798 1800 # for j in range(nPairs):
1799 1801 # dist1[j,0,0] = dist[pairs[j][0]]
1800 1802 # dist1[j,1,0] = dist[pairs[j][1]]
1801 1803 # ang1[j,0,0] = ang[pairs[j][0]]
1802 1804 # ang1[j,1,0] = ang[pairs[j][1]]
1803 1805 #
1804 1806 # return distx,disty, dist1,ang1
1805 1807
1806 1808
1807 1809 def __calculateVelVer(self, phase, lagTRange, _lambda):
1808 1810
1809 1811 Ts = lagTRange[1] - lagTRange[0]
1810 1812 velW = -_lambda*phase/(4*math.pi*Ts)
1811 1813
1812 1814 return velW
1813 1815
1814 1816 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1815 1817 nPairs = tau1.shape[0]
1816 1818 nHeights = tau1.shape[1]
1817 1819 vel = numpy.zeros((nPairs,3,nHeights))
1818 1820 dist1 = numpy.reshape(dist, (dist.size,1))
1819 1821
1820 1822 angCos = numpy.cos(ang)
1821 1823 angSin = numpy.sin(ang)
1822 1824
1823 1825 vel0 = dist1*tau1/(2*tau2**2)
1824 1826 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1825 1827 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1826 1828
1827 1829 ind = numpy.where(numpy.isinf(vel))
1828 1830 vel[ind] = numpy.nan
1829 1831
1830 1832 return vel
1831 1833
1832 1834 # def __getPairsAutoCorr(self, pairsList, nChannels):
1833 1835 #
1834 1836 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1835 1837 #
1836 1838 # for l in range(len(pairsList)):
1837 1839 # firstChannel = pairsList[l][0]
1838 1840 # secondChannel = pairsList[l][1]
1839 1841 #
1840 1842 # #Obteniendo pares de Autocorrelacion
1841 1843 # if firstChannel == secondChannel:
1842 1844 # pairsAutoCorr[firstChannel] = int(l)
1843 1845 #
1844 1846 # pairsAutoCorr = pairsAutoCorr.astype(int)
1845 1847 #
1846 1848 # pairsCrossCorr = range(len(pairsList))
1847 1849 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1848 1850 #
1849 1851 # return pairsAutoCorr, pairsCrossCorr
1850 1852
1851 1853 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1852 1854 def techniqueSA(self, kwargs):
1853 1855
1854 1856 """
1855 1857 Function that implements Spaced Antenna (SA) technique.
1856 1858
1857 1859 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1858 1860 Direction correction (if necessary), Ranges and SNR
1859 1861
1860 1862 Output: Winds estimation (Zonal, Meridional and Vertical)
1861 1863
1862 1864 Parameters affected: Winds
1863 1865 """
1864 1866 position_x = kwargs['positionX']
1865 1867 position_y = kwargs['positionY']
1866 1868 azimuth = kwargs['azimuth']
1867 1869
1868 1870 if 'correctFactor' in kwargs:
1869 1871 correctFactor = kwargs['correctFactor']
1870 1872 else:
1871 1873 correctFactor = 1
1872 1874
1873 1875 groupList = kwargs['groupList']
1874 1876 pairs_ccf = groupList[1]
1875 1877 tau = kwargs['tau']
1876 1878 _lambda = kwargs['_lambda']
1877 1879
1878 1880 #Cross Correlation pairs obtained
1879 1881 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
1880 1882 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1881 1883 # pairsSelArray = numpy.array(pairsSelected)
1882 1884 # pairs = []
1883 1885 #
1884 1886 # #Wind estimation pairs obtained
1885 1887 # for i in range(pairsSelArray.shape[0]/2):
1886 1888 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1887 1889 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1888 1890 # pairs.append((ind1,ind2))
1889 1891
1890 1892 indtau = tau.shape[0]/2
1891 1893 tau1 = tau[:indtau,:]
1892 1894 tau2 = tau[indtau:-1,:]
1893 1895 # tau1 = tau1[pairs,:]
1894 1896 # tau2 = tau2[pairs,:]
1895 1897 phase1 = tau[-1,:]
1896 1898
1897 1899 #---------------------------------------------------------------------
1898 1900 #Metodo Directo
1899 1901 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
1900 1902 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1901 1903 winds = stats.nanmean(winds, axis=0)
1902 1904 #---------------------------------------------------------------------
1903 1905 #Metodo General
1904 1906 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1905 1907 # #Calculo Coeficientes de Funcion de Correlacion
1906 1908 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1907 1909 # #Calculo de Velocidades
1908 1910 # winds = self.calculateVelUV(F,G,A,B,H)
1909 1911
1910 1912 #---------------------------------------------------------------------
1911 1913 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1912 1914 winds = correctFactor*winds
1913 1915 return winds
1914 1916
1915 1917 def __checkTime(self, currentTime, paramInterval, outputInterval):
1916 1918
1917 1919 dataTime = currentTime + paramInterval
1918 1920 deltaTime = dataTime - self.__initime
1919 1921
1920 1922 if deltaTime >= outputInterval or deltaTime < 0:
1921 1923 self.__dataReady = True
1922 1924 return
1923 1925
1924 1926 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1925 1927 '''
1926 1928 Function that implements winds estimation technique with detected meteors.
1927 1929
1928 1930 Input: Detected meteors, Minimum meteor quantity to wind estimation
1929 1931
1930 1932 Output: Winds estimation (Zonal and Meridional)
1931 1933
1932 1934 Parameters affected: Winds
1933 1935 '''
1934 1936 #Settings
1935 1937 nInt = (heightMax - heightMin)/2
1936 1938 nInt = int(nInt)
1937 1939 winds = numpy.zeros((2,nInt))*numpy.nan
1938 1940
1939 1941 #Filter errors
1940 1942 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1941 1943 finalMeteor = arrayMeteor[error,:]
1942 1944
1943 1945 #Meteor Histogram
1944 1946 finalHeights = finalMeteor[:,2]
1945 1947 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1946 1948 nMeteorsPerI = hist[0]
1947 1949 heightPerI = hist[1]
1948 1950
1949 1951 #Sort of meteors
1950 1952 indSort = finalHeights.argsort()
1951 1953 finalMeteor2 = finalMeteor[indSort,:]
1952 1954
1953 1955 # Calculating winds
1954 1956 ind1 = 0
1955 1957 ind2 = 0
1956 1958
1957 1959 for i in range(nInt):
1958 1960 nMet = nMeteorsPerI[i]
1959 1961 ind1 = ind2
1960 1962 ind2 = ind1 + nMet
1961 1963
1962 1964 meteorAux = finalMeteor2[ind1:ind2,:]
1963 1965
1964 1966 if meteorAux.shape[0] >= meteorThresh:
1965 1967 vel = meteorAux[:, 6]
1966 1968 zen = meteorAux[:, 4]*numpy.pi/180
1967 1969 azim = meteorAux[:, 3]*numpy.pi/180
1968 1970
1969 1971 n = numpy.cos(zen)
1970 1972 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1971 1973 # l = m*numpy.tan(azim)
1972 1974 l = numpy.sin(zen)*numpy.sin(azim)
1973 1975 m = numpy.sin(zen)*numpy.cos(azim)
1974 1976
1975 1977 A = numpy.vstack((l, m)).transpose()
1976 1978 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1977 1979 windsAux = numpy.dot(A1, vel)
1978 1980
1979 1981 winds[0,i] = windsAux[0]
1980 1982 winds[1,i] = windsAux[1]
1981 1983
1982 1984 return winds, heightPerI[:-1]
1983 1985
1984 1986 def techniqueNSM_SA(self, **kwargs):
1985 1987 metArray = kwargs['metArray']
1986 1988 heightList = kwargs['heightList']
1987 1989 timeList = kwargs['timeList']
1988 1990
1989 1991 rx_location = kwargs['rx_location']
1990 1992 groupList = kwargs['groupList']
1991 1993 azimuth = kwargs['azimuth']
1992 1994 dfactor = kwargs['dfactor']
1993 1995 k = kwargs['k']
1994 1996
1995 1997 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
1996 1998 d = dist*dfactor
1997 1999 #Phase calculation
1998 2000 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
1999 2001
2000 2002 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
2001 2003
2002 2004 velEst = numpy.zeros((heightList.size,2))*numpy.nan
2003 2005 azimuth1 = azimuth1*numpy.pi/180
2004 2006
2005 2007 for i in range(heightList.size):
2006 2008 h = heightList[i]
2007 2009 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
2008 2010 metHeight = metArray1[indH,:]
2009 2011 if metHeight.shape[0] >= 2:
2010 2012 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
2011 2013 iazim = metHeight[:,1].astype(int)
2012 2014 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
2013 2015 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
2014 2016 A = numpy.asmatrix(A)
2015 2017 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
2016 2018 velHor = numpy.dot(A1,velAux)
2017 2019
2018 2020 velEst[i,:] = numpy.squeeze(velHor)
2019 2021 return velEst
2020 2022
2021 2023 def __getPhaseSlope(self, metArray, heightList, timeList):
2022 2024 meteorList = []
2023 2025 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
2024 2026 #Putting back together the meteor matrix
2025 2027 utctime = metArray[:,0]
2026 2028 uniqueTime = numpy.unique(utctime)
2027 2029
2028 2030 phaseDerThresh = 0.5
2029 2031 ippSeconds = timeList[1] - timeList[0]
2030 2032 sec = numpy.where(timeList>1)[0][0]
2031 2033 nPairs = metArray.shape[1] - 6
2032 2034 nHeights = len(heightList)
2033 2035
2034 2036 for t in uniqueTime:
2035 2037 metArray1 = metArray[utctime==t,:]
2036 2038 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
2037 2039 tmet = metArray1[:,1].astype(int)
2038 2040 hmet = metArray1[:,2].astype(int)
2039 2041
2040 2042 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
2041 2043 metPhase[:,:] = numpy.nan
2042 2044 metPhase[:,hmet,tmet] = metArray1[:,6:].T
2043 2045
2044 2046 #Delete short trails
2045 2047 metBool = ~numpy.isnan(metPhase[0,:,:])
2046 2048 heightVect = numpy.sum(metBool, axis = 1)
2047 2049 metBool[heightVect<sec,:] = False
2048 2050 metPhase[:,heightVect<sec,:] = numpy.nan
2049 2051
2050 2052 #Derivative
2051 2053 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
2052 2054 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
2053 2055 metPhase[phDerAux] = numpy.nan
2054 2056
2055 2057 #--------------------------METEOR DETECTION -----------------------------------------
2056 2058 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
2057 2059
2058 2060 for p in numpy.arange(nPairs):
2059 2061 phase = metPhase[p,:,:]
2060 2062 phDer = metDer[p,:,:]
2061 2063
2062 2064 for h in indMet:
2063 2065 height = heightList[h]
2064 2066 phase1 = phase[h,:] #82
2065 2067 phDer1 = phDer[h,:]
2066 2068
2067 2069 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
2068 2070
2069 2071 indValid = numpy.where(~numpy.isnan(phase1))[0]
2070 2072 initMet = indValid[0]
2071 2073 endMet = 0
2072 2074
2073 2075 for i in range(len(indValid)-1):
2074 2076
2075 2077 #Time difference
2076 2078 inow = indValid[i]
2077 2079 inext = indValid[i+1]
2078 2080 idiff = inext - inow
2079 2081 #Phase difference
2080 2082 phDiff = numpy.abs(phase1[inext] - phase1[inow])
2081 2083
2082 2084 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
2083 2085 sizeTrail = inow - initMet + 1
2084 2086 if sizeTrail>3*sec: #Too short meteors
2085 2087 x = numpy.arange(initMet,inow+1)*ippSeconds
2086 2088 y = phase1[initMet:inow+1]
2087 2089 ynnan = ~numpy.isnan(y)
2088 2090 x = x[ynnan]
2089 2091 y = y[ynnan]
2090 2092 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
2091 2093 ylin = x*slope + intercept
2092 2094 rsq = r_value**2
2093 2095 if rsq > 0.5:
2094 2096 vel = slope#*height*1000/(k*d)
2095 2097 estAux = numpy.array([utctime,p,height, vel, rsq])
2096 2098 meteorList.append(estAux)
2097 2099 initMet = inext
2098 2100 metArray2 = numpy.array(meteorList)
2099 2101
2100 2102 return metArray2
2101 2103
2102 2104 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
2103 2105
2104 2106 azimuth1 = numpy.zeros(len(pairslist))
2105 2107 dist = numpy.zeros(len(pairslist))
2106 2108
2107 2109 for i in range(len(rx_location)):
2108 2110 ch0 = pairslist[i][0]
2109 2111 ch1 = pairslist[i][1]
2110 2112
2111 2113 diffX = rx_location[ch0][0] - rx_location[ch1][0]
2112 2114 diffY = rx_location[ch0][1] - rx_location[ch1][1]
2113 2115 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
2114 2116 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
2115 2117
2116 2118 azimuth1 -= azimuth0
2117 2119 return azimuth1, dist
2118 2120
2119 2121 def techniqueNSM_DBS(self, **kwargs):
2120 2122 metArray = kwargs['metArray']
2121 2123 heightList = kwargs['heightList']
2122 2124 timeList = kwargs['timeList']
2123 2125 azimuth = kwargs['azimuth']
2124 2126 theta_x = numpy.array(kwargs['theta_x'])
2125 2127 theta_y = numpy.array(kwargs['theta_y'])
2126 2128
2127 2129 utctime = metArray[:,0]
2128 2130 cmet = metArray[:,1].astype(int)
2129 2131 hmet = metArray[:,3].astype(int)
2130 2132 SNRmet = metArray[:,4]
2131 2133 vmet = metArray[:,5]
2132 2134 spcmet = metArray[:,6]
2133 2135
2134 2136 nChan = numpy.max(cmet) + 1
2135 2137 nHeights = len(heightList)
2136 2138
2137 2139 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
2138 2140 hmet = heightList[hmet]
2139 2141 h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights
2140 2142
2141 2143 velEst = numpy.zeros((heightList.size,2))*numpy.nan
2142 2144
2143 2145 for i in range(nHeights - 1):
2144 2146 hmin = heightList[i]
2145 2147 hmax = heightList[i + 1]
2146 2148
2147 2149 thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10)
2148 2150 indthisH = numpy.where(thisH)
2149 2151
2150 2152 if numpy.size(indthisH) > 3:
2151 2153
2152 2154 vel_aux = vmet[thisH]
2153 2155 chan_aux = cmet[thisH]
2154 2156 cosu_aux = dir_cosu[chan_aux]
2155 2157 cosv_aux = dir_cosv[chan_aux]
2156 2158 cosw_aux = dir_cosw[chan_aux]
2157 2159
2158 2160 nch = numpy.size(numpy.unique(chan_aux))
2159 2161 if nch > 1:
2160 2162 A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True)
2161 2163 velEst[i,:] = numpy.dot(A,vel_aux)
2162 2164
2163 2165 return velEst
2164 2166
2165 2167 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
2166 2168
2167 2169 param = dataOut.data_param
2168 2170 if dataOut.abscissaList != None:
2169 2171 absc = dataOut.abscissaList[:-1]
2170 2172 # noise = dataOut.noise
2171 2173 heightList = dataOut.heightList
2172 2174 SNR = dataOut.data_snr
2173 2175
2174 2176 if technique == 'DBS':
2175 2177
2176 2178 kwargs['velRadial'] = param[:,1,:] #Radial velocity
2177 2179 kwargs['heightList'] = heightList
2178 2180 kwargs['SNR'] = SNR
2179 2181
2180 2182 dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function
2181 2183 dataOut.utctimeInit = dataOut.utctime
2182 2184 dataOut.outputInterval = dataOut.paramInterval
2183 2185
2184 2186 elif technique == 'SA':
2185 2187
2186 2188 #Parameters
2187 2189 # position_x = kwargs['positionX']
2188 2190 # position_y = kwargs['positionY']
2189 2191 # azimuth = kwargs['azimuth']
2190 2192 #
2191 2193 # if kwargs.has_key('crosspairsList'):
2192 2194 # pairs = kwargs['crosspairsList']
2193 2195 # else:
2194 2196 # pairs = None
2195 2197 #
2196 2198 # if kwargs.has_key('correctFactor'):
2197 2199 # correctFactor = kwargs['correctFactor']
2198 2200 # else:
2199 2201 # correctFactor = 1
2200 2202
2201 2203 # tau = dataOut.data_param
2202 2204 # _lambda = dataOut.C/dataOut.frequency
2203 2205 # pairsList = dataOut.groupList
2204 2206 # nChannels = dataOut.nChannels
2205 2207
2206 2208 kwargs['groupList'] = dataOut.groupList
2207 2209 kwargs['tau'] = dataOut.data_param
2208 2210 kwargs['_lambda'] = dataOut.C/dataOut.frequency
2209 2211 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
2210 2212 dataOut.data_output = self.techniqueSA(kwargs)
2211 2213 dataOut.utctimeInit = dataOut.utctime
2212 2214 dataOut.outputInterval = dataOut.timeInterval
2213 2215
2214 2216 elif technique == 'Meteors':
2215 2217 dataOut.flagNoData = True
2216 2218 self.__dataReady = False
2217 2219
2218 2220 if 'nHours' in kwargs:
2219 2221 nHours = kwargs['nHours']
2220 2222 else:
2221 2223 nHours = 1
2222 2224
2223 2225 if 'meteorsPerBin' in kwargs:
2224 2226 meteorThresh = kwargs['meteorsPerBin']
2225 2227 else:
2226 2228 meteorThresh = 6
2227 2229
2228 2230 if 'hmin' in kwargs:
2229 2231 hmin = kwargs['hmin']
2230 2232 else: hmin = 70
2231 2233 if 'hmax' in kwargs:
2232 2234 hmax = kwargs['hmax']
2233 2235 else: hmax = 110
2234 2236
2235 2237 dataOut.outputInterval = nHours*3600
2236 2238
2237 2239 if self.__isConfig == False:
2238 2240 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2239 2241 #Get Initial LTC time
2240 2242 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2241 2243 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2242 2244
2243 2245 self.__isConfig = True
2244 2246
2245 2247 if self.__buffer is None:
2246 2248 self.__buffer = dataOut.data_param
2247 2249 self.__firstdata = copy.copy(dataOut)
2248 2250
2249 2251 else:
2250 2252 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2251 2253
2252 2254 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2253 2255
2254 2256 if self.__dataReady:
2255 2257 dataOut.utctimeInit = self.__initime
2256 2258
2257 2259 self.__initime += dataOut.outputInterval #to erase time offset
2258 2260
2259 2261 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
2260 2262 dataOut.flagNoData = False
2261 2263 self.__buffer = None
2262 2264
2263 2265 elif technique == 'Meteors1':
2264 2266 dataOut.flagNoData = True
2265 2267 self.__dataReady = False
2266 2268
2267 2269 if 'nMins' in kwargs:
2268 2270 nMins = kwargs['nMins']
2269 2271 else: nMins = 20
2270 2272 if 'rx_location' in kwargs:
2271 2273 rx_location = kwargs['rx_location']
2272 2274 else: rx_location = [(0,1),(1,1),(1,0)]
2273 2275 if 'azimuth' in kwargs:
2274 2276 azimuth = kwargs['azimuth']
2275 2277 else: azimuth = 51.06
2276 2278 if 'dfactor' in kwargs:
2277 2279 dfactor = kwargs['dfactor']
2278 2280 if 'mode' in kwargs:
2279 2281 mode = kwargs['mode']
2280 2282 if 'theta_x' in kwargs:
2281 2283 theta_x = kwargs['theta_x']
2282 2284 if 'theta_y' in kwargs:
2283 2285 theta_y = kwargs['theta_y']
2284 2286 else: mode = 'SA'
2285 2287
2286 2288 #Borrar luego esto
2287 2289 if dataOut.groupList is None:
2288 2290 dataOut.groupList = [(0,1),(0,2),(1,2)]
2289 2291 groupList = dataOut.groupList
2290 2292 C = 3e8
2291 2293 freq = 50e6
2292 2294 lamb = C/freq
2293 2295 k = 2*numpy.pi/lamb
2294 2296
2295 2297 timeList = dataOut.abscissaList
2296 2298 heightList = dataOut.heightList
2297 2299
2298 2300 if self.__isConfig == False:
2299 2301 dataOut.outputInterval = nMins*60
2300 2302 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2301 2303 #Get Initial LTC time
2302 2304 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2303 2305 minuteAux = initime.minute
2304 2306 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
2305 2307 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2306 2308
2307 2309 self.__isConfig = True
2308 2310
2309 2311 if self.__buffer is None:
2310 2312 self.__buffer = dataOut.data_param
2311 2313 self.__firstdata = copy.copy(dataOut)
2312 2314
2313 2315 else:
2314 2316 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2315 2317
2316 2318 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2317 2319
2318 2320 if self.__dataReady:
2319 2321 dataOut.utctimeInit = self.__initime
2320 2322 self.__initime += dataOut.outputInterval #to erase time offset
2321 2323
2322 2324 metArray = self.__buffer
2323 2325 if mode == 'SA':
2324 2326 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
2325 2327 elif mode == 'DBS':
2326 2328 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y)
2327 2329 dataOut.data_output = dataOut.data_output.T
2328 2330 dataOut.flagNoData = False
2329 2331 self.__buffer = None
2330 2332
2331 2333 return
2332 2334
2333 2335 class EWDriftsEstimation(Operation):
2334 2336
2335 2337 def __init__(self):
2336 2338 Operation.__init__(self)
2337 2339
2338 2340 def __correctValues(self, heiRang, phi, velRadial, SNR):
2339 2341 listPhi = phi.tolist()
2340 2342 maxid = listPhi.index(max(listPhi))
2341 2343 minid = listPhi.index(min(listPhi))
2342 2344
2343 2345 rango = list(range(len(phi)))
2344 2346 # rango = numpy.delete(rango,maxid)
2345 2347
2346 2348 heiRang1 = heiRang*math.cos(phi[maxid])
2347 2349 heiRangAux = heiRang*math.cos(phi[minid])
2348 2350 indOut = (heiRang1 < heiRangAux[0]).nonzero()
2349 2351 heiRang1 = numpy.delete(heiRang1,indOut)
2350 2352
2351 2353 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
2352 2354 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
2353 2355
2354 2356 for i in rango:
2355 2357 x = heiRang*math.cos(phi[i])
2356 2358 y1 = velRadial[i,:]
2357 2359 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
2358 2360
2359 2361 x1 = heiRang1
2360 2362 y11 = f1(x1)
2361 2363
2362 2364 y2 = SNR[i,:]
2363 2365 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
2364 2366 y21 = f2(x1)
2365 2367
2366 2368 velRadial1[i,:] = y11
2367 2369 SNR1[i,:] = y21
2368 2370
2369 2371 return heiRang1, velRadial1, SNR1
2370 2372
2371 2373 def run(self, dataOut, zenith, zenithCorrection):
2372 2374 heiRang = dataOut.heightList
2373 2375 velRadial = dataOut.data_param[:,3,:]
2374 2376 SNR = dataOut.data_snr
2375 2377
2376 2378 zenith = numpy.array(zenith)
2377 2379 zenith -= zenithCorrection
2378 2380 zenith *= numpy.pi/180
2379 2381
2380 2382 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
2381 2383
2382 2384 alp = zenith[0]
2383 2385 bet = zenith[1]
2384 2386
2385 2387 w_w = velRadial1[0,:]
2386 2388 w_e = velRadial1[1,:]
2387 2389
2388 2390 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
2389 2391 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
2390 2392
2391 2393 winds = numpy.vstack((u,w))
2392 2394
2393 2395 dataOut.heightList = heiRang1
2394 2396 dataOut.data_output = winds
2395 2397 dataOut.data_snr = SNR1
2396 2398
2397 2399 dataOut.utctimeInit = dataOut.utctime
2398 2400 dataOut.outputInterval = dataOut.timeInterval
2399 2401 return
2400 2402
2401 2403 #--------------- Non Specular Meteor ----------------
2402 2404
2403 2405 class NonSpecularMeteorDetection(Operation):
2404 2406
2405 2407 def run(self, dataOut, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
2406 2408 data_acf = dataOut.data_pre[0]
2407 2409 data_ccf = dataOut.data_pre[1]
2408 2410 pairsList = dataOut.groupList[1]
2409 2411
2410 2412 lamb = dataOut.C/dataOut.frequency
2411 2413 tSamp = dataOut.ippSeconds*dataOut.nCohInt
2412 2414 paramInterval = dataOut.paramInterval
2413 2415
2414 2416 nChannels = data_acf.shape[0]
2415 2417 nLags = data_acf.shape[1]
2416 2418 nProfiles = data_acf.shape[2]
2417 2419 nHeights = dataOut.nHeights
2418 2420 nCohInt = dataOut.nCohInt
2419 2421 sec = numpy.round(nProfiles/dataOut.paramInterval)
2420 2422 heightList = dataOut.heightList
2421 2423 ippSeconds = dataOut.ippSeconds*dataOut.nCohInt*dataOut.nAvg
2422 2424 utctime = dataOut.utctime
2423 2425
2424 2426 dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
2425 2427
2426 2428 #------------------------ SNR --------------------------------------
2427 2429 power = data_acf[:,0,:,:].real
2428 2430 noise = numpy.zeros(nChannels)
2429 2431 SNR = numpy.zeros(power.shape)
2430 2432 for i in range(nChannels):
2431 2433 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
2432 2434 SNR[i] = (power[i]-noise[i])/noise[i]
2433 2435 SNRm = numpy.nanmean(SNR, axis = 0)
2434 2436 SNRdB = 10*numpy.log10(SNR)
2435 2437
2436 2438 if mode == 'SA':
2437 2439 dataOut.groupList = dataOut.groupList[1]
2438 2440 nPairs = data_ccf.shape[0]
2439 2441 #---------------------- Coherence and Phase --------------------------
2440 2442 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
2441 2443 # phase1 = numpy.copy(phase)
2442 2444 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
2443 2445
2444 2446 for p in range(nPairs):
2445 2447 ch0 = pairsList[p][0]
2446 2448 ch1 = pairsList[p][1]
2447 2449 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
2448 2450 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
2449 2451 # phase1[p,:,:] = numpy.angle(ccf) #median filter
2450 2452 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
2451 2453 # coh1[p,:,:] = numpy.abs(ccf) #median filter
2452 2454 coh = numpy.nanmax(coh1, axis = 0)
2453 2455 # struc = numpy.ones((5,1))
2454 2456 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
2455 2457 #---------------------- Radial Velocity ----------------------------
2456 2458 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
2457 2459 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
2458 2460
2459 2461 if allData:
2460 2462 boolMetFin = ~numpy.isnan(SNRm)
2461 2463 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
2462 2464 else:
2463 2465 #------------------------ Meteor mask ---------------------------------
2464 2466 # #SNR mask
2465 2467 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
2466 2468 #
2467 2469 # #Erase small objects
2468 2470 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
2469 2471 #
2470 2472 # auxEEJ = numpy.sum(boolMet1,axis=0)
2471 2473 # indOver = auxEEJ>nProfiles*0.8 #Use this later
2472 2474 # indEEJ = numpy.where(indOver)[0]
2473 2475 # indNEEJ = numpy.where(~indOver)[0]
2474 2476 #
2475 2477 # boolMetFin = boolMet1
2476 2478 #
2477 2479 # if indEEJ.size > 0:
2478 2480 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
2479 2481 #
2480 2482 # boolMet2 = coh > cohThresh
2481 2483 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
2482 2484 #
2483 2485 # #Final Meteor mask
2484 2486 # boolMetFin = boolMet1|boolMet2
2485 2487
2486 2488 #Coherence mask
2487 2489 boolMet1 = coh > 0.75
2488 2490 struc = numpy.ones((30,1))
2489 2491 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
2490 2492
2491 2493 #Derivative mask
2492 2494 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
2493 2495 boolMet2 = derPhase < 0.2
2494 2496 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
2495 2497 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
2496 2498 boolMet2 = ndimage.median_filter(boolMet2,size=5)
2497 2499 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
2498 2500 # #Final mask
2499 2501 # boolMetFin = boolMet2
2500 2502 boolMetFin = boolMet1&boolMet2
2501 2503 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
2502 2504 #Creating data_param
2503 2505 coordMet = numpy.where(boolMetFin)
2504 2506
2505 2507 tmet = coordMet[0]
2506 2508 hmet = coordMet[1]
2507 2509
2508 2510 data_param = numpy.zeros((tmet.size, 6 + nPairs))
2509 2511 data_param[:,0] = utctime
2510 2512 data_param[:,1] = tmet
2511 2513 data_param[:,2] = hmet
2512 2514 data_param[:,3] = SNRm[tmet,hmet]
2513 2515 data_param[:,4] = velRad[tmet,hmet]
2514 2516 data_param[:,5] = coh[tmet,hmet]
2515 2517 data_param[:,6:] = phase[:,tmet,hmet].T
2516 2518
2517 2519 elif mode == 'DBS':
2518 2520 dataOut.groupList = numpy.arange(nChannels)
2519 2521
2520 2522 #Radial Velocities
2521 2523 phase = numpy.angle(data_acf[:,1,:,:])
2522 2524 # phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
2523 2525 velRad = phase*lamb/(4*numpy.pi*tSamp)
2524 2526
2525 2527 #Spectral width
2526 2528 # acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
2527 2529 # acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
2528 2530 acf1 = data_acf[:,1,:,:]
2529 2531 acf2 = data_acf[:,2,:,:]
2530 2532
2531 2533 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
2532 2534 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
2533 2535 if allData:
2534 2536 boolMetFin = ~numpy.isnan(SNRdB)
2535 2537 else:
2536 2538 #SNR
2537 2539 boolMet1 = (SNRdB>SNRthresh) #SNR mask
2538 2540 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
2539 2541
2540 2542 #Radial velocity
2541 2543 boolMet2 = numpy.abs(velRad) < 20
2542 2544 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
2543 2545
2544 2546 #Spectral Width
2545 2547 boolMet3 = spcWidth < 30
2546 2548 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
2547 2549 # boolMetFin = self.__erase_small(boolMet1, 10,5)
2548 2550 boolMetFin = boolMet1&boolMet2&boolMet3
2549 2551
2550 2552 #Creating data_param
2551 2553 coordMet = numpy.where(boolMetFin)
2552 2554
2553 2555 cmet = coordMet[0]
2554 2556 tmet = coordMet[1]
2555 2557 hmet = coordMet[2]
2556 2558
2557 2559 data_param = numpy.zeros((tmet.size, 7))
2558 2560 data_param[:,0] = utctime
2559 2561 data_param[:,1] = cmet
2560 2562 data_param[:,2] = tmet
2561 2563 data_param[:,3] = hmet
2562 2564 data_param[:,4] = SNR[cmet,tmet,hmet].T
2563 2565 data_param[:,5] = velRad[cmet,tmet,hmet].T
2564 2566 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
2565 2567
2566 2568 # self.dataOut.data_param = data_int
2567 2569 if len(data_param) == 0:
2568 2570 dataOut.flagNoData = True
2569 2571 else:
2570 2572 dataOut.data_param = data_param
2571 2573
2572 2574 def __erase_small(self, binArray, threshX, threshY):
2573 2575 labarray, numfeat = ndimage.measurements.label(binArray)
2574 2576 binArray1 = numpy.copy(binArray)
2575 2577
2576 2578 for i in range(1,numfeat + 1):
2577 2579 auxBin = (labarray==i)
2578 2580 auxSize = auxBin.sum()
2579 2581
2580 2582 x,y = numpy.where(auxBin)
2581 2583 widthX = x.max() - x.min()
2582 2584 widthY = y.max() - y.min()
2583 2585
2584 2586 #width X: 3 seg -> 12.5*3
2585 2587 #width Y:
2586 2588
2587 2589 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
2588 2590 binArray1[auxBin] = False
2589 2591
2590 2592 return binArray1
2591 2593
2592 2594 #--------------- Specular Meteor ----------------
2593 2595
2594 2596 class SMDetection(Operation):
2595 2597 '''
2596 2598 Function DetectMeteors()
2597 2599 Project developed with paper:
2598 2600 HOLDSWORTH ET AL. 2004
2599 2601
2600 2602 Input:
2601 2603 self.dataOut.data_pre
2602 2604
2603 2605 centerReceiverIndex: From the channels, which is the center receiver
2604 2606
2605 2607 hei_ref: Height reference for the Beacon signal extraction
2606 2608 tauindex:
2607 2609 predefinedPhaseShifts: Predefined phase offset for the voltge signals
2608 2610
2609 2611 cohDetection: Whether to user Coherent detection or not
2610 2612 cohDet_timeStep: Coherent Detection calculation time step
2611 2613 cohDet_thresh: Coherent Detection phase threshold to correct phases
2612 2614
2613 2615 noise_timeStep: Noise calculation time step
2614 2616 noise_multiple: Noise multiple to define signal threshold
2615 2617
2616 2618 multDet_timeLimit: Multiple Detection Removal time limit in seconds
2617 2619 multDet_rangeLimit: Multiple Detection Removal range limit in km
2618 2620
2619 2621 phaseThresh: Maximum phase difference between receiver to be consider a meteor
2620 2622 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
2621 2623
2622 2624 hmin: Minimum Height of the meteor to use it in the further wind estimations
2623 2625 hmax: Maximum Height of the meteor to use it in the further wind estimations
2624 2626 azimuth: Azimuth angle correction
2625 2627
2626 2628 Affected:
2627 2629 self.dataOut.data_param
2628 2630
2629 2631 Rejection Criteria (Errors):
2630 2632 0: No error; analysis OK
2631 2633 1: SNR < SNR threshold
2632 2634 2: angle of arrival (AOA) ambiguously determined
2633 2635 3: AOA estimate not feasible
2634 2636 4: Large difference in AOAs obtained from different antenna baselines
2635 2637 5: echo at start or end of time series
2636 2638 6: echo less than 5 examples long; too short for analysis
2637 2639 7: echo rise exceeds 0.3s
2638 2640 8: echo decay time less than twice rise time
2639 2641 9: large power level before echo
2640 2642 10: large power level after echo
2641 2643 11: poor fit to amplitude for estimation of decay time
2642 2644 12: poor fit to CCF phase variation for estimation of radial drift velocity
2643 2645 13: height unresolvable echo: not valid height within 70 to 110 km
2644 2646 14: height ambiguous echo: more then one possible height within 70 to 110 km
2645 2647 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
2646 2648 16: oscilatory echo, indicating event most likely not an underdense echo
2647 2649
2648 2650 17: phase difference in meteor Reestimation
2649 2651
2650 2652 Data Storage:
2651 2653 Meteors for Wind Estimation (8):
2652 2654 Utc Time | Range Height
2653 2655 Azimuth Zenith errorCosDir
2654 2656 VelRad errorVelRad
2655 2657 Phase0 Phase1 Phase2 Phase3
2656 2658 TypeError
2657 2659
2658 2660 '''
2659 2661
2660 2662 def run(self, dataOut, hei_ref = None, tauindex = 0,
2661 2663 phaseOffsets = None,
2662 2664 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
2663 2665 noise_timeStep = 4, noise_multiple = 4,
2664 2666 multDet_timeLimit = 1, multDet_rangeLimit = 3,
2665 2667 phaseThresh = 20, SNRThresh = 5,
2666 2668 hmin = 50, hmax=150, azimuth = 0,
2667 2669 channelPositions = None) :
2668 2670
2669 2671
2670 2672 #Getting Pairslist
2671 2673 if channelPositions is None:
2672 2674 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2673 2675 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2674 2676 meteorOps = SMOperations()
2675 2677 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2676 2678 heiRang = dataOut.heightList
2677 2679 #Get Beacon signal - No Beacon signal anymore
2678 2680 # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
2679 2681 #
2680 2682 # if hei_ref != None:
2681 2683 # newheis = numpy.where(self.dataOut.heightList>hei_ref)
2682 2684 #
2683 2685
2684 2686
2685 2687 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
2686 2688 # see if the user put in pre defined phase shifts
2687 2689 voltsPShift = dataOut.data_pre.copy()
2688 2690
2689 2691 # if predefinedPhaseShifts != None:
2690 2692 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
2691 2693 #
2692 2694 # # elif beaconPhaseShifts:
2693 2695 # # #get hardware phase shifts using beacon signal
2694 2696 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
2695 2697 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
2696 2698 #
2697 2699 # else:
2698 2700 # hardwarePhaseShifts = numpy.zeros(5)
2699 2701 #
2700 2702 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
2701 2703 # for i in range(self.dataOut.data_pre.shape[0]):
2702 2704 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
2703 2705
2704 2706 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
2705 2707
2706 2708 #Remove DC
2707 2709 voltsDC = numpy.mean(voltsPShift,1)
2708 2710 voltsDC = numpy.mean(voltsDC,1)
2709 2711 for i in range(voltsDC.shape[0]):
2710 2712 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
2711 2713
2712 2714 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
2713 2715 # voltsPShift = voltsPShift[:,:,:newheis[0][0]]
2714 2716
2715 2717 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
2716 2718 #Coherent Detection
2717 2719 if cohDetection:
2718 2720 #use coherent detection to get the net power
2719 2721 cohDet_thresh = cohDet_thresh*numpy.pi/180
2720 2722 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
2721 2723
2722 2724 #Non-coherent detection!
2723 2725 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
2724 2726 #********** END OF COH/NON-COH POWER CALCULATION**********************
2725 2727
2726 2728 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
2727 2729 #Get noise
2728 2730 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
2729 2731 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
2730 2732 #Get signal threshold
2731 2733 signalThresh = noise_multiple*noise
2732 2734 #Meteor echoes detection
2733 2735 listMeteors = self.__findMeteors(powerNet, signalThresh)
2734 2736 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
2735 2737
2736 2738 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
2737 2739 #Parameters
2738 2740 heiRange = dataOut.heightList
2739 2741 rangeInterval = heiRange[1] - heiRange[0]
2740 2742 rangeLimit = multDet_rangeLimit/rangeInterval
2741 2743 timeLimit = multDet_timeLimit/dataOut.timeInterval
2742 2744 #Multiple detection removals
2743 2745 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
2744 2746 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
2745 2747
2746 2748 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
2747 2749 #Parameters
2748 2750 phaseThresh = phaseThresh*numpy.pi/180
2749 2751 thresh = [phaseThresh, noise_multiple, SNRThresh]
2750 2752 #Meteor reestimation (Errors N 1, 6, 12, 17)
2751 2753 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
2752 2754 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
2753 2755 #Estimation of decay times (Errors N 7, 8, 11)
2754 2756 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
2755 2757 #******************* END OF METEOR REESTIMATION *******************
2756 2758
2757 2759 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
2758 2760 #Calculating Radial Velocity (Error N 15)
2759 2761 radialStdThresh = 10
2760 2762 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
2761 2763
2762 2764 if len(listMeteors4) > 0:
2763 2765 #Setting New Array
2764 2766 date = dataOut.utctime
2765 2767 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
2766 2768
2767 2769 #Correcting phase offset
2768 2770 if phaseOffsets != None:
2769 2771 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2770 2772 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2771 2773
2772 2774 #Second Pairslist
2773 2775 pairsList = []
2774 2776 pairx = (0,1)
2775 2777 pairy = (2,3)
2776 2778 pairsList.append(pairx)
2777 2779 pairsList.append(pairy)
2778 2780
2779 2781 jph = numpy.array([0,0,0,0])
2780 2782 h = (hmin,hmax)
2781 2783 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2782 2784
2783 2785 # #Calculate AOA (Error N 3, 4)
2784 2786 # #JONES ET AL. 1998
2785 2787 # error = arrayParameters[:,-1]
2786 2788 # AOAthresh = numpy.pi/8
2787 2789 # phases = -arrayParameters[:,9:13]
2788 2790 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
2789 2791 #
2790 2792 # #Calculate Heights (Error N 13 and 14)
2791 2793 # error = arrayParameters[:,-1]
2792 2794 # Ranges = arrayParameters[:,2]
2793 2795 # zenith = arrayParameters[:,5]
2794 2796 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
2795 2797 # error = arrayParameters[:,-1]
2796 2798 #********************* END OF PARAMETERS CALCULATION **************************
2797 2799
2798 2800 #***************************+ PASS DATA TO NEXT STEP **********************
2799 2801 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
2800 2802 dataOut.data_param = arrayParameters
2801 2803
2802 2804 if arrayParameters is None:
2803 2805 dataOut.flagNoData = True
2804 2806 else:
2805 2807 dataOut.flagNoData = True
2806 2808
2807 2809 return
2808 2810
2809 2811 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
2810 2812
2811 2813 minIndex = min(newheis[0])
2812 2814 maxIndex = max(newheis[0])
2813 2815
2814 2816 voltage = voltage0[:,:,minIndex:maxIndex+1]
2815 2817 nLength = voltage.shape[1]/n
2816 2818 nMin = 0
2817 2819 nMax = 0
2818 2820 phaseOffset = numpy.zeros((len(pairslist),n))
2819 2821
2820 2822 for i in range(n):
2821 2823 nMax += nLength
2822 2824 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
2823 2825 phaseCCF = numpy.mean(phaseCCF, axis = 2)
2824 2826 phaseOffset[:,i] = phaseCCF.transpose()
2825 2827 nMin = nMax
2826 2828 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
2827 2829
2828 2830 #Remove Outliers
2829 2831 factor = 2
2830 2832 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
2831 2833 dw = numpy.std(wt,axis = 1)
2832 2834 dw = dw.reshape((dw.size,1))
2833 2835 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
2834 2836 phaseOffset[ind] = numpy.nan
2835 2837 phaseOffset = stats.nanmean(phaseOffset, axis=1)
2836 2838
2837 2839 return phaseOffset
2838 2840
2839 2841 def __shiftPhase(self, data, phaseShift):
2840 2842 #this will shift the phase of a complex number
2841 2843 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
2842 2844 return dataShifted
2843 2845
2844 2846 def __estimatePhaseDifference(self, array, pairslist):
2845 2847 nChannel = array.shape[0]
2846 2848 nHeights = array.shape[2]
2847 2849 numPairs = len(pairslist)
2848 2850 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
2849 2851 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
2850 2852
2851 2853 #Correct phases
2852 2854 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
2853 2855 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2854 2856
2855 2857 if indDer[0].shape[0] > 0:
2856 2858 for i in range(indDer[0].shape[0]):
2857 2859 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
2858 2860 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
2859 2861
2860 2862 # for j in range(numSides):
2861 2863 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
2862 2864 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
2863 2865 #
2864 2866 #Linear
2865 2867 phaseInt = numpy.zeros((numPairs,1))
2866 2868 angAllCCF = phaseCCF[:,[0,1,3,4],0]
2867 2869 for j in range(numPairs):
2868 2870 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
2869 2871 phaseInt[j] = fit[1]
2870 2872 #Phase Differences
2871 2873 phaseDiff = phaseInt - phaseCCF[:,2,:]
2872 2874 phaseArrival = phaseInt.reshape(phaseInt.size)
2873 2875
2874 2876 #Dealias
2875 2877 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
2876 2878 # indAlias = numpy.where(phaseArrival > numpy.pi)
2877 2879 # phaseArrival[indAlias] -= 2*numpy.pi
2878 2880 # indAlias = numpy.where(phaseArrival < -numpy.pi)
2879 2881 # phaseArrival[indAlias] += 2*numpy.pi
2880 2882
2881 2883 return phaseDiff, phaseArrival
2882 2884
2883 2885 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
2884 2886 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
2885 2887 #find the phase shifts of each channel over 1 second intervals
2886 2888 #only look at ranges below the beacon signal
2887 2889 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
2888 2890 numBlocks = int(volts.shape[1]/numProfPerBlock)
2889 2891 numHeights = volts.shape[2]
2890 2892 nChannel = volts.shape[0]
2891 2893 voltsCohDet = volts.copy()
2892 2894
2893 2895 pairsarray = numpy.array(pairslist)
2894 2896 indSides = pairsarray[:,1]
2895 2897 # indSides = numpy.array(range(nChannel))
2896 2898 # indSides = numpy.delete(indSides, indCenter)
2897 2899 #
2898 2900 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
2899 2901 listBlocks = numpy.array_split(volts, numBlocks, 1)
2900 2902
2901 2903 startInd = 0
2902 2904 endInd = 0
2903 2905
2904 2906 for i in range(numBlocks):
2905 2907 startInd = endInd
2906 2908 endInd = endInd + listBlocks[i].shape[1]
2907 2909
2908 2910 arrayBlock = listBlocks[i]
2909 2911 # arrayBlockCenter = listCenter[i]
2910 2912
2911 2913 #Estimate the Phase Difference
2912 2914 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
2913 2915 #Phase Difference RMS
2914 2916 arrayPhaseRMS = numpy.abs(phaseDiff)
2915 2917 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
2916 2918 indPhase = numpy.where(phaseRMSaux==4)
2917 2919 #Shifting
2918 2920 if indPhase[0].shape[0] > 0:
2919 2921 for j in range(indSides.size):
2920 2922 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
2921 2923 voltsCohDet[:,startInd:endInd,:] = arrayBlock
2922 2924
2923 2925 return voltsCohDet
2924 2926
2925 2927 def __calculateCCF(self, volts, pairslist ,laglist):
2926 2928
2927 2929 nHeights = volts.shape[2]
2928 2930 nPoints = volts.shape[1]
2929 2931 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
2930 2932
2931 2933 for i in range(len(pairslist)):
2932 2934 volts1 = volts[pairslist[i][0]]
2933 2935 volts2 = volts[pairslist[i][1]]
2934 2936
2935 2937 for t in range(len(laglist)):
2936 2938 idxT = laglist[t]
2937 2939 if idxT >= 0:
2938 2940 vStacked = numpy.vstack((volts2[idxT:,:],
2939 2941 numpy.zeros((idxT, nHeights),dtype='complex')))
2940 2942 else:
2941 2943 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
2942 2944 volts2[:(nPoints + idxT),:]))
2943 2945 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
2944 2946
2945 2947 vStacked = None
2946 2948 return voltsCCF
2947 2949
2948 2950 def __getNoise(self, power, timeSegment, timeInterval):
2949 2951 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
2950 2952 numBlocks = int(power.shape[0]/numProfPerBlock)
2951 2953 numHeights = power.shape[1]
2952 2954
2953 2955 listPower = numpy.array_split(power, numBlocks, 0)
2954 2956 noise = numpy.zeros((power.shape[0], power.shape[1]))
2955 2957 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
2956 2958
2957 2959 startInd = 0
2958 2960 endInd = 0
2959 2961
2960 2962 for i in range(numBlocks): #split por canal
2961 2963 startInd = endInd
2962 2964 endInd = endInd + listPower[i].shape[0]
2963 2965
2964 2966 arrayBlock = listPower[i]
2965 2967 noiseAux = numpy.mean(arrayBlock, 0)
2966 2968 # noiseAux = numpy.median(noiseAux)
2967 2969 # noiseAux = numpy.mean(arrayBlock)
2968 2970 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
2969 2971
2970 2972 noiseAux1 = numpy.mean(arrayBlock)
2971 2973 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
2972 2974
2973 2975 return noise, noise1
2974 2976
2975 2977 def __findMeteors(self, power, thresh):
2976 2978 nProf = power.shape[0]
2977 2979 nHeights = power.shape[1]
2978 2980 listMeteors = []
2979 2981
2980 2982 for i in range(nHeights):
2981 2983 powerAux = power[:,i]
2982 2984 threshAux = thresh[:,i]
2983 2985
2984 2986 indUPthresh = numpy.where(powerAux > threshAux)[0]
2985 2987 indDNthresh = numpy.where(powerAux <= threshAux)[0]
2986 2988
2987 2989 j = 0
2988 2990
2989 2991 while (j < indUPthresh.size - 2):
2990 2992 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
2991 2993 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
2992 2994 indDNthresh = indDNthresh[indDNAux]
2993 2995
2994 2996 if (indDNthresh.size > 0):
2995 2997 indEnd = indDNthresh[0] - 1
2996 2998 indInit = indUPthresh[j]
2997 2999
2998 3000 meteor = powerAux[indInit:indEnd + 1]
2999 3001 indPeak = meteor.argmax() + indInit
3000 3002 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
3001 3003
3002 3004 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
3003 3005 j = numpy.where(indUPthresh == indEnd)[0] + 1
3004 3006 else: j+=1
3005 3007 else: j+=1
3006 3008
3007 3009 return listMeteors
3008 3010
3009 3011 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
3010 3012
3011 3013 arrayMeteors = numpy.asarray(listMeteors)
3012 3014 listMeteors1 = []
3013 3015
3014 3016 while arrayMeteors.shape[0] > 0:
3015 3017 FLAs = arrayMeteors[:,4]
3016 3018 maxFLA = FLAs.argmax()
3017 3019 listMeteors1.append(arrayMeteors[maxFLA,:])
3018 3020
3019 3021 MeteorInitTime = arrayMeteors[maxFLA,1]
3020 3022 MeteorEndTime = arrayMeteors[maxFLA,3]
3021 3023 MeteorHeight = arrayMeteors[maxFLA,0]
3022 3024
3023 3025 #Check neighborhood
3024 3026 maxHeightIndex = MeteorHeight + rangeLimit
3025 3027 minHeightIndex = MeteorHeight - rangeLimit
3026 3028 minTimeIndex = MeteorInitTime - timeLimit
3027 3029 maxTimeIndex = MeteorEndTime + timeLimit
3028 3030
3029 3031 #Check Heights
3030 3032 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
3031 3033 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
3032 3034 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
3033 3035
3034 3036 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
3035 3037
3036 3038 return listMeteors1
3037 3039
3038 3040 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
3039 3041 numHeights = volts.shape[2]
3040 3042 nChannel = volts.shape[0]
3041 3043
3042 3044 thresholdPhase = thresh[0]
3043 3045 thresholdNoise = thresh[1]
3044 3046 thresholdDB = float(thresh[2])
3045 3047
3046 3048 thresholdDB1 = 10**(thresholdDB/10)
3047 3049 pairsarray = numpy.array(pairslist)
3048 3050 indSides = pairsarray[:,1]
3049 3051
3050 3052 pairslist1 = list(pairslist)
3051 3053 pairslist1.append((0,1))
3052 3054 pairslist1.append((3,4))
3053 3055
3054 3056 listMeteors1 = []
3055 3057 listPowerSeries = []
3056 3058 listVoltageSeries = []
3057 3059 #volts has the war data
3058 3060
3059 3061 if frequency == 30e6:
3060 3062 timeLag = 45*10**-3
3061 3063 else:
3062 3064 timeLag = 15*10**-3
3063 3065 lag = numpy.ceil(timeLag/timeInterval)
3064 3066
3065 3067 for i in range(len(listMeteors)):
3066 3068
3067 3069 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
3068 3070 meteorAux = numpy.zeros(16)
3069 3071
3070 3072 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
3071 3073 mHeight = listMeteors[i][0]
3072 3074 mStart = listMeteors[i][1]
3073 3075 mPeak = listMeteors[i][2]
3074 3076 mEnd = listMeteors[i][3]
3075 3077
3076 3078 #get the volt data between the start and end times of the meteor
3077 3079 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
3078 3080 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
3079 3081
3080 3082 #3.6. Phase Difference estimation
3081 3083 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
3082 3084
3083 3085 #3.7. Phase difference removal & meteor start, peak and end times reestimated
3084 3086 #meteorVolts0.- all Channels, all Profiles
3085 3087 meteorVolts0 = volts[:,:,mHeight]
3086 3088 meteorThresh = noise[:,mHeight]*thresholdNoise
3087 3089 meteorNoise = noise[:,mHeight]
3088 3090 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
3089 3091 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
3090 3092
3091 3093 #Times reestimation
3092 3094 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
3093 3095 if mStart1.size > 0:
3094 3096 mStart1 = mStart1[-1] + 1
3095 3097
3096 3098 else:
3097 3099 mStart1 = mPeak
3098 3100
3099 3101 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
3100 3102 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
3101 3103 if mEndDecayTime1.size == 0:
3102 3104 mEndDecayTime1 = powerNet0.size
3103 3105 else:
3104 3106 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
3105 3107 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
3106 3108
3107 3109 #meteorVolts1.- all Channels, from start to end
3108 3110 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
3109 3111 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
3110 3112 if meteorVolts2.shape[1] == 0:
3111 3113 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
3112 3114 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
3113 3115 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
3114 3116 ##################### END PARAMETERS REESTIMATION #########################
3115 3117
3116 3118 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
3117 3119 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
3118 3120 if meteorVolts2.shape[1] > 0:
3119 3121 #Phase Difference re-estimation
3120 3122 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
3121 3123 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
3122 3124 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
3123 3125 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
3124 3126 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
3125 3127
3126 3128 #Phase Difference RMS
3127 3129 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
3128 3130 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
3129 3131 #Data from Meteor
3130 3132 mPeak1 = powerNet1.argmax() + mStart1
3131 3133 mPeakPower1 = powerNet1.max()
3132 3134 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
3133 3135 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
3134 3136 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
3135 3137 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
3136 3138 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
3137 3139 #Vectorize
3138 3140 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
3139 3141 meteorAux[7:11] = phaseDiffint[0:4]
3140 3142
3141 3143 #Rejection Criterions
3142 3144 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
3143 3145 meteorAux[-1] = 17
3144 3146 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
3145 3147 meteorAux[-1] = 1
3146 3148
3147 3149
3148 3150 else:
3149 3151 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
3150 3152 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
3151 3153 PowerSeries = 0
3152 3154
3153 3155 listMeteors1.append(meteorAux)
3154 3156 listPowerSeries.append(PowerSeries)
3155 3157 listVoltageSeries.append(meteorVolts1)
3156 3158
3157 3159 return listMeteors1, listPowerSeries, listVoltageSeries
3158 3160
3159 3161 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
3160 3162
3161 3163 threshError = 10
3162 3164 #Depending if it is 30 or 50 MHz
3163 3165 if frequency == 30e6:
3164 3166 timeLag = 45*10**-3
3165 3167 else:
3166 3168 timeLag = 15*10**-3
3167 3169 lag = numpy.ceil(timeLag/timeInterval)
3168 3170
3169 3171 listMeteors1 = []
3170 3172
3171 3173 for i in range(len(listMeteors)):
3172 3174 meteorPower = listPower[i]
3173 3175 meteorAux = listMeteors[i]
3174 3176
3175 3177 if meteorAux[-1] == 0:
3176 3178
3177 3179 try:
3178 3180 indmax = meteorPower.argmax()
3179 3181 indlag = indmax + lag
3180 3182
3181 3183 y = meteorPower[indlag:]
3182 3184 x = numpy.arange(0, y.size)*timeLag
3183 3185
3184 3186 #first guess
3185 3187 a = y[0]
3186 3188 tau = timeLag
3187 3189 #exponential fit
3188 3190 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
3189 3191 y1 = self.__exponential_function(x, *popt)
3190 3192 #error estimation
3191 3193 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
3192 3194
3193 3195 decayTime = popt[1]
3194 3196 riseTime = indmax*timeInterval
3195 3197 meteorAux[11:13] = [decayTime, error]
3196 3198
3197 3199 #Table items 7, 8 and 11
3198 3200 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
3199 3201 meteorAux[-1] = 7
3200 3202 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
3201 3203 meteorAux[-1] = 8
3202 3204 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
3203 3205 meteorAux[-1] = 11
3204 3206
3205 3207
3206 3208 except:
3207 3209 meteorAux[-1] = 11
3208 3210
3209 3211
3210 3212 listMeteors1.append(meteorAux)
3211 3213
3212 3214 return listMeteors1
3213 3215
3214 3216 #Exponential Function
3215 3217
3216 3218 def __exponential_function(self, x, a, tau):
3217 3219 y = a*numpy.exp(-x/tau)
3218 3220 return y
3219 3221
3220 3222 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
3221 3223
3222 3224 pairslist1 = list(pairslist)
3223 3225 pairslist1.append((0,1))
3224 3226 pairslist1.append((3,4))
3225 3227 numPairs = len(pairslist1)
3226 3228 #Time Lag
3227 3229 timeLag = 45*10**-3
3228 3230 c = 3e8
3229 3231 lag = numpy.ceil(timeLag/timeInterval)
3230 3232 freq = 30e6
3231 3233
3232 3234 listMeteors1 = []
3233 3235
3234 3236 for i in range(len(listMeteors)):
3235 3237 meteorAux = listMeteors[i]
3236 3238 if meteorAux[-1] == 0:
3237 3239 mStart = listMeteors[i][1]
3238 3240 mPeak = listMeteors[i][2]
3239 3241 mLag = mPeak - mStart + lag
3240 3242
3241 3243 #get the volt data between the start and end times of the meteor
3242 3244 meteorVolts = listVolts[i]
3243 3245 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
3244 3246
3245 3247 #Get CCF
3246 3248 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
3247 3249
3248 3250 #Method 2
3249 3251 slopes = numpy.zeros(numPairs)
3250 3252 time = numpy.array([-2,-1,1,2])*timeInterval
3251 3253 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
3252 3254
3253 3255 #Correct phases
3254 3256 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
3255 3257 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
3256 3258
3257 3259 if indDer[0].shape[0] > 0:
3258 3260 for i in range(indDer[0].shape[0]):
3259 3261 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
3260 3262 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
3261 3263
3262 3264 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
3263 3265 for j in range(numPairs):
3264 3266 fit = stats.linregress(time, angAllCCF[j,:])
3265 3267 slopes[j] = fit[0]
3266 3268
3267 3269 #Remove Outlier
3268 3270 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
3269 3271 # slopes = numpy.delete(slopes,indOut)
3270 3272 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
3271 3273 # slopes = numpy.delete(slopes,indOut)
3272 3274
3273 3275 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
3274 3276 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
3275 3277 meteorAux[-2] = radialError
3276 3278 meteorAux[-3] = radialVelocity
3277 3279
3278 3280 #Setting Error
3279 3281 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
3280 3282 if numpy.abs(radialVelocity) > 200:
3281 3283 meteorAux[-1] = 15
3282 3284 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
3283 3285 elif radialError > radialStdThresh:
3284 3286 meteorAux[-1] = 12
3285 3287
3286 3288 listMeteors1.append(meteorAux)
3287 3289 return listMeteors1
3288 3290
3289 3291 def __setNewArrays(self, listMeteors, date, heiRang):
3290 3292
3291 3293 #New arrays
3292 3294 arrayMeteors = numpy.array(listMeteors)
3293 3295 arrayParameters = numpy.zeros((len(listMeteors), 13))
3294 3296
3295 3297 #Date inclusion
3296 3298 # date = re.findall(r'\((.*?)\)', date)
3297 3299 # date = date[0].split(',')
3298 3300 # date = map(int, date)
3299 3301 #
3300 3302 # if len(date)<6:
3301 3303 # date.append(0)
3302 3304 #
3303 3305 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
3304 3306 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
3305 3307 arrayDate = numpy.tile(date, (len(listMeteors)))
3306 3308
3307 3309 #Meteor array
3308 3310 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
3309 3311 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
3310 3312
3311 3313 #Parameters Array
3312 3314 arrayParameters[:,0] = arrayDate #Date
3313 3315 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
3314 3316 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
3315 3317 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
3316 3318 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
3317 3319
3318 3320
3319 3321 return arrayParameters
3320 3322
3321 3323 class CorrectSMPhases(Operation):
3322 3324
3323 3325 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
3324 3326
3325 3327 arrayParameters = dataOut.data_param
3326 3328 pairsList = []
3327 3329 pairx = (0,1)
3328 3330 pairy = (2,3)
3329 3331 pairsList.append(pairx)
3330 3332 pairsList.append(pairy)
3331 3333 jph = numpy.zeros(4)
3332 3334
3333 3335 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
3334 3336 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
3335 3337 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
3336 3338
3337 3339 meteorOps = SMOperations()
3338 3340 if channelPositions is None:
3339 3341 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3340 3342 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3341 3343
3342 3344 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
3343 3345 h = (hmin,hmax)
3344 3346
3345 3347 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
3346 3348
3347 3349 dataOut.data_param = arrayParameters
3348 3350 return
3349 3351
3350 3352 class SMPhaseCalibration(Operation):
3351 3353
3352 3354 __buffer = None
3353 3355
3354 3356 __initime = None
3355 3357
3356 3358 __dataReady = False
3357 3359
3358 3360 __isConfig = False
3359 3361
3360 3362 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
3361 3363
3362 3364 dataTime = currentTime + paramInterval
3363 3365 deltaTime = dataTime - initTime
3364 3366
3365 3367 if deltaTime >= outputInterval or deltaTime < 0:
3366 3368 return True
3367 3369
3368 3370 return False
3369 3371
3370 3372 def __getGammas(self, pairs, d, phases):
3371 3373 gammas = numpy.zeros(2)
3372 3374
3373 3375 for i in range(len(pairs)):
3374 3376
3375 3377 pairi = pairs[i]
3376 3378
3377 3379 phip3 = phases[:,pairi[0]]
3378 3380 d3 = d[pairi[0]]
3379 3381 phip2 = phases[:,pairi[1]]
3380 3382 d2 = d[pairi[1]]
3381 3383 #Calculating gamma
3382 3384 # jdcos = alp1/(k*d1)
3383 3385 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
3384 3386 jgamma = -phip2*d3/d2 - phip3
3385 3387 jgamma = numpy.angle(numpy.exp(1j*jgamma))
3386 3388 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
3387 3389 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
3388 3390
3389 3391 #Revised distribution
3390 3392 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
3391 3393
3392 3394 #Histogram
3393 3395 nBins = 64
3394 3396 rmin = -0.5*numpy.pi
3395 3397 rmax = 0.5*numpy.pi
3396 3398 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
3397 3399
3398 3400 meteorsY = phaseHisto[0]
3399 3401 phasesX = phaseHisto[1][:-1]
3400 3402 width = phasesX[1] - phasesX[0]
3401 3403 phasesX += width/2
3402 3404
3403 3405 #Gaussian aproximation
3404 3406 bpeak = meteorsY.argmax()
3405 3407 peak = meteorsY.max()
3406 3408 jmin = bpeak - 5
3407 3409 jmax = bpeak + 5 + 1
3408 3410
3409 3411 if jmin<0:
3410 3412 jmin = 0
3411 3413 jmax = 6
3412 3414 elif jmax > meteorsY.size:
3413 3415 jmin = meteorsY.size - 6
3414 3416 jmax = meteorsY.size
3415 3417
3416 3418 x0 = numpy.array([peak,bpeak,50])
3417 3419 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
3418 3420
3419 3421 #Gammas
3420 3422 gammas[i] = coeff[0][1]
3421 3423
3422 3424 return gammas
3423 3425
3424 3426 def __residualFunction(self, coeffs, y, t):
3425 3427
3426 3428 return y - self.__gauss_function(t, coeffs)
3427 3429
3428 3430 def __gauss_function(self, t, coeffs):
3429 3431
3430 3432 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
3431 3433
3432 3434 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
3433 3435 meteorOps = SMOperations()
3434 3436 nchan = 4
3435 3437 pairx = pairsList[0] #x es 0
3436 3438 pairy = pairsList[1] #y es 1
3437 3439 center_xangle = 0
3438 3440 center_yangle = 0
3439 3441 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
3440 3442 ntimes = len(range_angle)
3441 3443
3442 3444 nstepsx = 20
3443 3445 nstepsy = 20
3444 3446
3445 3447 for iz in range(ntimes):
3446 3448 min_xangle = -range_angle[iz]/2 + center_xangle
3447 3449 max_xangle = range_angle[iz]/2 + center_xangle
3448 3450 min_yangle = -range_angle[iz]/2 + center_yangle
3449 3451 max_yangle = range_angle[iz]/2 + center_yangle
3450 3452
3451 3453 inc_x = (max_xangle-min_xangle)/nstepsx
3452 3454 inc_y = (max_yangle-min_yangle)/nstepsy
3453 3455
3454 3456 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
3455 3457 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
3456 3458 penalty = numpy.zeros((nstepsx,nstepsy))
3457 3459 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
3458 3460 jph = numpy.zeros(nchan)
3459 3461
3460 3462 # Iterations looking for the offset
3461 3463 for iy in range(int(nstepsy)):
3462 3464 for ix in range(int(nstepsx)):
3463 3465 d3 = d[pairsList[1][0]]
3464 3466 d2 = d[pairsList[1][1]]
3465 3467 d5 = d[pairsList[0][0]]
3466 3468 d4 = d[pairsList[0][1]]
3467 3469
3468 3470 alp2 = alpha_y[iy] #gamma 1
3469 3471 alp4 = alpha_x[ix] #gamma 0
3470 3472
3471 3473 alp3 = -alp2*d3/d2 - gammas[1]
3472 3474 alp5 = -alp4*d5/d4 - gammas[0]
3473 3475 # jph[pairy[1]] = alpha_y[iy]
3474 3476 # jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
3475 3477
3476 3478 # jph[pairx[1]] = alpha_x[ix]
3477 3479 # jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
3478 3480 jph[pairsList[0][1]] = alp4
3479 3481 jph[pairsList[0][0]] = alp5
3480 3482 jph[pairsList[1][0]] = alp3
3481 3483 jph[pairsList[1][1]] = alp2
3482 3484 jph_array[:,ix,iy] = jph
3483 3485 # d = [2.0,2.5,2.5,2.0]
3484 3486 #falta chequear si va a leer bien los meteoros
3485 3487 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
3486 3488 error = meteorsArray1[:,-1]
3487 3489 ind1 = numpy.where(error==0)[0]
3488 3490 penalty[ix,iy] = ind1.size
3489 3491
3490 3492 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
3491 3493 phOffset = jph_array[:,i,j]
3492 3494
3493 3495 center_xangle = phOffset[pairx[1]]
3494 3496 center_yangle = phOffset[pairy[1]]
3495 3497
3496 3498 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
3497 3499 phOffset = phOffset*180/numpy.pi
3498 3500 return phOffset
3499 3501
3500 3502
3501 3503 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
3502 3504
3503 3505 dataOut.flagNoData = True
3504 3506 self.__dataReady = False
3505 3507 dataOut.outputInterval = nHours*3600
3506 3508
3507 3509 if self.__isConfig == False:
3508 3510 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
3509 3511 #Get Initial LTC time
3510 3512 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
3511 3513 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
3512 3514
3513 3515 self.__isConfig = True
3514 3516
3515 3517 if self.__buffer is None:
3516 3518 self.__buffer = dataOut.data_param.copy()
3517 3519
3518 3520 else:
3519 3521 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
3520 3522
3521 3523 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
3522 3524
3523 3525 if self.__dataReady:
3524 3526 dataOut.utctimeInit = self.__initime
3525 3527 self.__initime += dataOut.outputInterval #to erase time offset
3526 3528
3527 3529 freq = dataOut.frequency
3528 3530 c = dataOut.C #m/s
3529 3531 lamb = c/freq
3530 3532 k = 2*numpy.pi/lamb
3531 3533 azimuth = 0
3532 3534 h = (hmin, hmax)
3533 3535 # pairs = ((0,1),(2,3)) #Estrella
3534 3536 # pairs = ((1,0),(2,3)) #T
3535 3537
3536 3538 if channelPositions is None:
3537 3539 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3538 3540 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3539 3541 meteorOps = SMOperations()
3540 3542 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
3541 3543
3542 3544 #Checking correct order of pairs
3543 3545 pairs = []
3544 3546 if distances[1] > distances[0]:
3545 3547 pairs.append((1,0))
3546 3548 else:
3547 3549 pairs.append((0,1))
3548 3550
3549 3551 if distances[3] > distances[2]:
3550 3552 pairs.append((3,2))
3551 3553 else:
3552 3554 pairs.append((2,3))
3553 3555 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
3554 3556
3555 3557 meteorsArray = self.__buffer
3556 3558 error = meteorsArray[:,-1]
3557 3559 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
3558 3560 ind1 = numpy.where(boolError)[0]
3559 3561 meteorsArray = meteorsArray[ind1,:]
3560 3562 meteorsArray[:,-1] = 0
3561 3563 phases = meteorsArray[:,8:12]
3562 3564
3563 3565 #Calculate Gammas
3564 3566 gammas = self.__getGammas(pairs, distances, phases)
3565 3567 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
3566 3568 #Calculate Phases
3567 3569 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
3568 3570 phasesOff = phasesOff.reshape((1,phasesOff.size))
3569 3571 dataOut.data_output = -phasesOff
3570 3572 dataOut.flagNoData = False
3571 3573 self.__buffer = None
3572 3574
3573 3575
3574 3576 return
3575 3577
3576 3578 class SMOperations():
3577 3579
3578 3580 def __init__(self):
3579 3581
3580 3582 return
3581 3583
3582 3584 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
3583 3585
3584 3586 arrayParameters = arrayParameters0.copy()
3585 3587 hmin = h[0]
3586 3588 hmax = h[1]
3587 3589
3588 3590 #Calculate AOA (Error N 3, 4)
3589 3591 #JONES ET AL. 1998
3590 3592 AOAthresh = numpy.pi/8
3591 3593 error = arrayParameters[:,-1]
3592 3594 phases = -arrayParameters[:,8:12] + jph
3593 3595 # phases = numpy.unwrap(phases)
3594 3596 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
3595 3597
3596 3598 #Calculate Heights (Error N 13 and 14)
3597 3599 error = arrayParameters[:,-1]
3598 3600 Ranges = arrayParameters[:,1]
3599 3601 zenith = arrayParameters[:,4]
3600 3602 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
3601 3603
3602 3604 #----------------------- Get Final data ------------------------------------
3603 3605 # error = arrayParameters[:,-1]
3604 3606 # ind1 = numpy.where(error==0)[0]
3605 3607 # arrayParameters = arrayParameters[ind1,:]
3606 3608
3607 3609 return arrayParameters
3608 3610
3609 3611 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
3610 3612
3611 3613 arrayAOA = numpy.zeros((phases.shape[0],3))
3612 3614 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
3613 3615
3614 3616 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
3615 3617 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
3616 3618 arrayAOA[:,2] = cosDirError
3617 3619
3618 3620 azimuthAngle = arrayAOA[:,0]
3619 3621 zenithAngle = arrayAOA[:,1]
3620 3622
3621 3623 #Setting Error
3622 3624 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
3623 3625 error[indError] = 0
3624 3626 #Number 3: AOA not fesible
3625 3627 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
3626 3628 error[indInvalid] = 3
3627 3629 #Number 4: Large difference in AOAs obtained from different antenna baselines
3628 3630 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
3629 3631 error[indInvalid] = 4
3630 3632 return arrayAOA, error
3631 3633
3632 3634 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
3633 3635
3634 3636 #Initializing some variables
3635 3637 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
3636 3638 ang_aux = ang_aux.reshape(1,ang_aux.size)
3637 3639
3638 3640 cosdir = numpy.zeros((arrayPhase.shape[0],2))
3639 3641 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
3640 3642
3641 3643
3642 3644 for i in range(2):
3643 3645 ph0 = arrayPhase[:,pairsList[i][0]]
3644 3646 ph1 = arrayPhase[:,pairsList[i][1]]
3645 3647 d0 = distances[pairsList[i][0]]
3646 3648 d1 = distances[pairsList[i][1]]
3647 3649
3648 3650 ph0_aux = ph0 + ph1
3649 3651 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
3650 3652 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
3651 3653 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
3652 3654 #First Estimation
3653 3655 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
3654 3656
3655 3657 #Most-Accurate Second Estimation
3656 3658 phi1_aux = ph0 - ph1
3657 3659 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
3658 3660 #Direction Cosine 1
3659 3661 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
3660 3662
3661 3663 #Searching the correct Direction Cosine
3662 3664 cosdir0_aux = cosdir0[:,i]
3663 3665 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
3664 3666 #Minimum Distance
3665 3667 cosDiff = (cosdir1 - cosdir0_aux)**2
3666 3668 indcos = cosDiff.argmin(axis = 1)
3667 3669 #Saving Value obtained
3668 3670 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
3669 3671
3670 3672 return cosdir0, cosdir
3671 3673
3672 3674 def __calculateAOA(self, cosdir, azimuth):
3673 3675 cosdirX = cosdir[:,0]
3674 3676 cosdirY = cosdir[:,1]
3675 3677
3676 3678 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
3677 3679 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
3678 3680 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
3679 3681
3680 3682 return angles
3681 3683
3682 3684 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
3683 3685
3684 3686 Ramb = 375 #Ramb = c/(2*PRF)
3685 3687 Re = 6371 #Earth Radius
3686 3688 heights = numpy.zeros(Ranges.shape)
3687 3689
3688 3690 R_aux = numpy.array([0,1,2])*Ramb
3689 3691 R_aux = R_aux.reshape(1,R_aux.size)
3690 3692
3691 3693 Ranges = Ranges.reshape(Ranges.size,1)
3692 3694
3693 3695 Ri = Ranges + R_aux
3694 3696 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
3695 3697
3696 3698 #Check if there is a height between 70 and 110 km
3697 3699 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
3698 3700 ind_h = numpy.where(h_bool == 1)[0]
3699 3701
3700 3702 hCorr = hi[ind_h, :]
3701 3703 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
3702 3704
3703 3705 hCorr = hi[ind_hCorr][:len(ind_h)]
3704 3706 heights[ind_h] = hCorr
3705 3707
3706 3708 #Setting Error
3707 3709 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
3708 3710 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
3709 3711 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
3710 3712 error[indError] = 0
3711 3713 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
3712 3714 error[indInvalid2] = 14
3713 3715 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
3714 3716 error[indInvalid1] = 13
3715 3717
3716 3718 return heights, error
3717 3719
3718 3720 def getPhasePairs(self, channelPositions):
3719 3721 chanPos = numpy.array(channelPositions)
3720 3722 listOper = list(itertools.combinations(list(range(5)),2))
3721 3723
3722 3724 distances = numpy.zeros(4)
3723 3725 axisX = []
3724 3726 axisY = []
3725 3727 distX = numpy.zeros(3)
3726 3728 distY = numpy.zeros(3)
3727 3729 ix = 0
3728 3730 iy = 0
3729 3731
3730 3732 pairX = numpy.zeros((2,2))
3731 3733 pairY = numpy.zeros((2,2))
3732 3734
3733 3735 for i in range(len(listOper)):
3734 3736 pairi = listOper[i]
3735 3737
3736 3738 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
3737 3739
3738 3740 if posDif[0] == 0:
3739 3741 axisY.append(pairi)
3740 3742 distY[iy] = posDif[1]
3741 3743 iy += 1
3742 3744 elif posDif[1] == 0:
3743 3745 axisX.append(pairi)
3744 3746 distX[ix] = posDif[0]
3745 3747 ix += 1
3746 3748
3747 3749 for i in range(2):
3748 3750 if i==0:
3749 3751 dist0 = distX
3750 3752 axis0 = axisX
3751 3753 else:
3752 3754 dist0 = distY
3753 3755 axis0 = axisY
3754 3756
3755 3757 side = numpy.argsort(dist0)[:-1]
3756 3758 axis0 = numpy.array(axis0)[side,:]
3757 3759 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
3758 3760 axis1 = numpy.unique(numpy.reshape(axis0,4))
3759 3761 side = axis1[axis1 != chanC]
3760 3762 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
3761 3763 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
3762 3764 if diff1<0:
3763 3765 chan2 = side[0]
3764 3766 d2 = numpy.abs(diff1)
3765 3767 chan1 = side[1]
3766 3768 d1 = numpy.abs(diff2)
3767 3769 else:
3768 3770 chan2 = side[1]
3769 3771 d2 = numpy.abs(diff2)
3770 3772 chan1 = side[0]
3771 3773 d1 = numpy.abs(diff1)
3772 3774
3773 3775 if i==0:
3774 3776 chanCX = chanC
3775 3777 chan1X = chan1
3776 3778 chan2X = chan2
3777 3779 distances[0:2] = numpy.array([d1,d2])
3778 3780 else:
3779 3781 chanCY = chanC
3780 3782 chan1Y = chan1
3781 3783 chan2Y = chan2
3782 3784 distances[2:4] = numpy.array([d1,d2])
3783 3785 # axisXsides = numpy.reshape(axisX[ix,:],4)
3784 3786 #
3785 3787 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
3786 3788 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
3787 3789 #
3788 3790 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
3789 3791 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
3790 3792 # channel25X = int(pairX[0,ind25X])
3791 3793 # channel20X = int(pairX[1,ind20X])
3792 3794 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
3793 3795 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
3794 3796 # channel25Y = int(pairY[0,ind25Y])
3795 3797 # channel20Y = int(pairY[1,ind20Y])
3796 3798
3797 3799 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
3798 3800 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
3799 3801
3800 3802 return pairslist, distances
3801 3803 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
3802 3804 #
3803 3805 # arrayAOA = numpy.zeros((phases.shape[0],3))
3804 3806 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
3805 3807 #
3806 3808 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
3807 3809 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
3808 3810 # arrayAOA[:,2] = cosDirError
3809 3811 #
3810 3812 # azimuthAngle = arrayAOA[:,0]
3811 3813 # zenithAngle = arrayAOA[:,1]
3812 3814 #
3813 3815 # #Setting Error
3814 3816 # #Number 3: AOA not fesible
3815 3817 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
3816 3818 # error[indInvalid] = 3
3817 3819 # #Number 4: Large difference in AOAs obtained from different antenna baselines
3818 3820 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
3819 3821 # error[indInvalid] = 4
3820 3822 # return arrayAOA, error
3821 3823 #
3822 3824 # def __getDirectionCosines(self, arrayPhase, pairsList):
3823 3825 #
3824 3826 # #Initializing some variables
3825 3827 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
3826 3828 # ang_aux = ang_aux.reshape(1,ang_aux.size)
3827 3829 #
3828 3830 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
3829 3831 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
3830 3832 #
3831 3833 #
3832 3834 # for i in range(2):
3833 3835 # #First Estimation
3834 3836 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
3835 3837 # #Dealias
3836 3838 # indcsi = numpy.where(phi0_aux > numpy.pi)
3837 3839 # phi0_aux[indcsi] -= 2*numpy.pi
3838 3840 # indcsi = numpy.where(phi0_aux < -numpy.pi)
3839 3841 # phi0_aux[indcsi] += 2*numpy.pi
3840 3842 # #Direction Cosine 0
3841 3843 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
3842 3844 #
3843 3845 # #Most-Accurate Second Estimation
3844 3846 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
3845 3847 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
3846 3848 # #Direction Cosine 1
3847 3849 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
3848 3850 #
3849 3851 # #Searching the correct Direction Cosine
3850 3852 # cosdir0_aux = cosdir0[:,i]
3851 3853 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
3852 3854 # #Minimum Distance
3853 3855 # cosDiff = (cosdir1 - cosdir0_aux)**2
3854 3856 # indcos = cosDiff.argmin(axis = 1)
3855 3857 # #Saving Value obtained
3856 3858 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
3857 3859 #
3858 3860 # return cosdir0, cosdir
3859 3861 #
3860 3862 # def __calculateAOA(self, cosdir, azimuth):
3861 3863 # cosdirX = cosdir[:,0]
3862 3864 # cosdirY = cosdir[:,1]
3863 3865 #
3864 3866 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
3865 3867 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
3866 3868 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
3867 3869 #
3868 3870 # return angles
3869 3871 #
3870 3872 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
3871 3873 #
3872 3874 # Ramb = 375 #Ramb = c/(2*PRF)
3873 3875 # Re = 6371 #Earth Radius
3874 3876 # heights = numpy.zeros(Ranges.shape)
3875 3877 #
3876 3878 # R_aux = numpy.array([0,1,2])*Ramb
3877 3879 # R_aux = R_aux.reshape(1,R_aux.size)
3878 3880 #
3879 3881 # Ranges = Ranges.reshape(Ranges.size,1)
3880 3882 #
3881 3883 # Ri = Ranges + R_aux
3882 3884 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
3883 3885 #
3884 3886 # #Check if there is a height between 70 and 110 km
3885 3887 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
3886 3888 # ind_h = numpy.where(h_bool == 1)[0]
3887 3889 #
3888 3890 # hCorr = hi[ind_h, :]
3889 3891 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
3890 3892 #
3891 3893 # hCorr = hi[ind_hCorr]
3892 3894 # heights[ind_h] = hCorr
3893 3895 #
3894 3896 # #Setting Error
3895 3897 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
3896 3898 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
3897 3899 #
3898 3900 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
3899 3901 # error[indInvalid2] = 14
3900 3902 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
3901 3903 # error[indInvalid1] = 13
3902 3904 #
3903 3905 # return heights, error
3904 3906
3905 3907
3906 3908 class WeatherRadar(Operation):
3907 3909 '''
3908 3910 Function tat implements Weather Radar operations-
3909 3911 Input:
3910 3912 Output:
3911 3913 Parameters affected:
3912 3914 '''
3913 3915 isConfig = False
3914 3916
3915 3917 def __init__(self):
3916 3918 Operation.__init__(self)
3917 3919
3918 3920 def setup(self,dataOut,Pt=0,Gt=0,Gr=0,lambda_=0, aL=0,
3919 3921 tauW= 0,thetaT=0,thetaR=0,Km =0):
3920 3922 self.nCh = dataOut.nChannels
3921 3923 self.nHeis = dataOut.nHeights
3922 3924 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
3923 3925 self.Range = numpy.arange(dataOut.nHeights)*deltaHeight + dataOut.heightList[0]
3924 3926 self.Range = self.Range.reshape(1,self.nHeis)
3925 3927 self.Range = numpy.tile(self.Range,[self.nCh,1])
3926 3928 '''-----------1 Constante del Radar----------'''
3927 3929 self.Pt = Pt
3928 3930 self.Gt = Gt
3929 3931 self.Gr = Gr
3930 3932 self.lambda_ = lambda_
3931 3933 self.aL = aL
3932 3934 self.tauW = tauW
3933 3935 self.thetaT = thetaT
3934 3936 self.thetaR = thetaR
3935 3937 self.Km = Km
3936 3938 Numerator = ((4*numpy.pi)**3 * aL**2 * 16 *numpy.log(2))
3937 3939 Denominator = (Pt * Gt * Gr * lambda_**2 * SPEED_OF_LIGHT * tauW * numpy.pi*thetaT*thetaR)
3938 3940 self.RadarConstant = Numerator/Denominator
3939 3941 '''-----------2 Reflectividad del Radar y Factor de Reflectividad------'''
3940 3942 self.n_radar = numpy.zeros((self.nCh,self.nHeis))
3941 3943 self.Z_radar = numpy.zeros((self.nCh,self.nHeis))
3942 3944
3943 3945 def setMoments(self,dataOut,i):
3944 3946
3945 3947 type = dataOut.inputUnit
3946 3948 nCh = dataOut.nChannels
3947 3949 nHeis= dataOut.nHeights
3948 3950 data_param = numpy.zeros((nCh,4,nHeis))
3949 3951 if type == "Voltage":
3950 3952 data_param[:,0,:] = dataOut.dataPP_POW/(dataOut.nCohInt**2)
3951 3953 data_param[:,1,:] = dataOut.dataPP_DOP
3952 3954 data_param[:,2,:] = dataOut.dataPP_WIDTH
3953 3955 data_param[:,3,:] = dataOut.dataPP_SNR
3954 3956 if type == "Spectra":
3955 3957 data_param[:,0,:] = dataOut.data_POW
3956 3958 data_param[:,1,:] = dataOut.data_DOP
3957 3959 data_param[:,2,:] = dataOut.data_WIDTH
3958 3960 def setMoments(self,dataOut,i):
3959 3961 data_param[:,3,:] = dataOut.data_SNR
3960 3962
3961 3963 return data_param[:,i,:]
3962 3964
3963 3965
3964 3966 def run(self,dataOut,Pt=25,Gt=200.0,Gr=50.0,lambda_=0.32, aL=2.5118,
3965 3967 tauW= 4.0e-6,thetaT=0.165,thetaR=0.367,Km =0.93):
3966 3968
3967 3969 if not self.isConfig:
3968 3970 self.setup(dataOut= dataOut,Pt=25,Gt=200.0,Gr=50.0,lambda_=0.32, aL=2.5118,
3969 3971 tauW= 4.0e-6,thetaT=0.165,thetaR=0.367,Km =0.93)
3970 3972 self.isConfig = True
3971 3973 '''-----------------------------Potencia de Radar -Signal S-----------------------------'''
3972 3974 Pr = self.setMoments(dataOut,0)
3973 3975
3974 3976 for R in range(self.nHeis):
3975 3977 self.n_radar[:,R] = self.RadarConstant*Pr[:,R]* (self.Range[:,R])**2
3976 3978
3977 3979 self.Z_radar[:,R] = self.n_radar[:,R]* self.lambda_**4/( numpy.pi**5 * self.Km**2)
3978 3980
3979 3981 '''----------- Factor de Reflectividad Equivalente lamda_ < 10 cm , lamda_= 3.2cm-------'''
3980 3982 Zeh = self.Z_radar
3981 3983 dBZeh = 10*numpy.log10(Zeh)
3982 3984 dataOut.factor_Zeh= dBZeh
3983 3985 self.n_radar = numpy.zeros((self.nCh,self.nHeis))
3984 3986 self.Z_radar = numpy.zeros((self.nCh,self.nHeis))
3985 3987
3986 3988 return dataOut
3987 3989
3988 3990 class PedestalInformation(Operation):
3989 3991 path_ped = None
3990 3992 path_adq = None
3991 3993 t_Interval_p = None
3992 3994 n_Muestras_p = None
3993 3995 isConfig = False
3994 3996 blocksPerfile= None
3995 3997 f_a_p = None
3996 3998 online = None
3997 3999 angulo_adq = None
3998 4000 nro_file = None
3999 4001 nro_key_p = None
4000 4002 tmp = None
4001 4003
4002 4004
4003 4005 def __init__(self):
4004 4006 Operation.__init__(self)
4005 4007
4006 4008 def getfirstFilefromPath(self,path,meta,ext):
4007 4009 validFilelist = []
4008 4010 #print("SEARH",path)
4009 4011 try:
4010 4012 fileList = os.listdir(path)
4011 4013 except:
4012 4014 print("check path - fileList")
4013 4015 if len(fileList)<1:
4014 4016 return None
4015 4017 # meta 1234 567 8-18 BCDE
4016 4018 # H,D,PE YYYY DDD EPOC .ext
4017 4019
4018 4020 for thisFile in fileList:
4019 4021 #print("HI",thisFile)
4020 4022 if meta =="PE":
4021 4023 try:
4022 4024 number= int(thisFile[len(meta)+7:len(meta)+17])
4023 4025 except:
4024 4026 print("There is a file or folder with different format")
4025 4027 if meta == "D":
4026 4028 try:
4027 4029 number= int(thisFile[8:11])
4028 4030 except:
4029 4031 print("There is a file or folder with different format")
4030 4032
4031 4033 if not isNumber(str=number):
4032 4034 continue
4033 4035 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
4034 4036 continue
4035 4037 validFilelist.sort()
4036 4038 validFilelist.append(thisFile)
4037 4039 if len(validFilelist)>0:
4038 4040 validFilelist = sorted(validFilelist,key=str.lower)
4039 4041 return validFilelist
4040 4042 return None
4041 4043
4042 4044 def gettimeutcfromDirFilename(self,path,file):
4043 4045 dir_file= path+"/"+file
4044 4046 fp = h5py.File(dir_file,'r')
4045 4047 #epoc = fp['Metadata'].get('utctimeInit')[()]
4046 4048 epoc = fp['Data'].get('utc')[()]
4047 4049 fp.close()
4048 4050 return epoc
4049 4051
4050 4052 def gettimeutcadqfromDirFilename(self,path,file):
4051 4053 dir_file= path+"/"+file
4052 4054 fp = h5py.File(dir_file,'r')
4053 4055 epoc = fp['Metadata'].get('utctimeInit')[()]
4054 4056 #epoc = fp['Data'].get('utc')[()]
4055 4057 fp.close()
4056 4058 return epoc
4057 4059
4058 4060 def getDatavaluefromDirFilename(self,path,file,value):
4059 4061 dir_file= path+"/"+file
4060 4062 fp = h5py.File(dir_file,'r')
4061 4063 array = fp['Data'].get(value)[()]
4062 4064 fp.close()
4063 4065 return array
4064 4066
4065 4067 def getFile_KeyP(self,list_pedestal,list_adq):
4066 4068 print(list_pedestal)
4067 4069 print(list_adq)
4068 4070
4069 4071 def getNROFile(self,utc_adq,utc_ped_list):
4070 4072 c=0
4071 4073 print("insidegetNROFile")
4072 4074 print(utc_adq)
4073 4075 print(len(utc_ped_list))
4074 4076 for i in range(len(utc_ped_list)):
4075 4077 if utc_adq>utc_ped_list[i]:
4076 4078 #print("mayor")
4077 4079 #print("utc_ped_list",utc_ped_list[i])
4078 4080 c +=1
4079 4081
4080 4082 return c-1,utc_ped_list[c-1],utc_ped_list[c]
4081 4083
4082 4084 def verificarNROFILE(self,dataOut,utc_ped,f_a_p,n_Muestras_p):
4083 4085 var =int(f_a_p/n_Muestras_p)
4084 4086 flag=0
4085 4087 for i in range(var):
4086 4088 if dataOut.utctime+i==utc_ped:
4087 4089 flag==1
4088 4090 break
4089 4091 return flag
4090 4092
4091 def setup_offline(self,dataOut,list_pedestal,list_adq):
4093 #def setup_offline(self,dataOut,list_pedestal,list_adq):
4094 def setup_offline(self,dataOut,list_pedestal):
4095
4092 4096 print("SETUP OFFLINE")
4093 4097 print(self.path_ped)
4094 print(self.path_adq)
4098 #print(self.path_adq)
4095 4099 print(len(self.list_pedestal))
4096 print(len(self.list_adq))
4100 #print(len(self.list_adq))
4097 4101 utc_ped_list=[]
4098 4102 for i in range(len(self.list_pedestal)):
4103 print(i)
4099 4104 utc_ped_list.append(self.gettimeutcfromDirFilename(path=self.path_ped,file=self.list_pedestal[i]))
4100 4105
4101 4106 #utc_ped_list= utc_ped_list
4102 utc_adq = self.gettimeutcadqfromDirFilename(path=self.path_adq,file=self.list_adq[0])
4107 ###utc_adq = self.gettimeutcadqfromDirFilename(path=self.path_adq,file=self.list_adq[0])
4103 4108 print("dios existe donde esta")
4109
4104 4110 #print("utc_ped_list",utc_ped_list)
4105 print("utc_adq",utc_adq)
4111 ###print("utc_adq",utc_adq)
4106 4112 # utc_adq_dataOut
4107 4113 utc_adq_dataOut =dataOut.utctime
4108 4114 print("Offline-utc_adq_dataout",utc_adq_dataOut)
4109 4115
4110 nro_file,utc_ped,utc_ped_1 = self.getNROFile(utc_adq=utc_adq, utc_ped_list= utc_ped_list)
4116 nro_file,utc_ped,utc_ped_1 = self.getNROFile(utc_adq=utc_adq_dataOut, utc_ped_list= utc_ped_list)
4111 4117
4112 4118 print("nro_file",nro_file,"utc_ped",utc_ped)
4113 4119 print("nro_file",i)
4114 nro_key_p = int((utc_adq-utc_ped)/self.t_Interval_p)-1 # ojito al -1 estimado alex
4120 nro_key_p = int((utc_adq_dataOut-utc_ped)/self.t_Interval_p)-1 # ojito al -1 estimado alex
4115 4121 print("nro_key_p",nro_key_p)
4116 4122
4117 4123 ff_pedestal = self.list_pedestal[nro_file]
4118 4124 #angulo = self.getDatavaluefromDirFilename(path=self.path_ped,file=ff_pedestal,value="azimuth")
4119 4125 angulo = self.getDatavaluefromDirFilename(path=self.path_ped,file=ff_pedestal,value="azi_pos")
4120 4126
4121 4127 print("utc_pedestal_init :",utc_ped+nro_key_p*self.t_Interval_p)
4122 4128 print("angulo_array :",angulo[nro_key_p])
4123 4129 self.nro_file = nro_file
4124 4130 self.nro_key_p = nro_key_p
4125 4131
4126 4132 def setup_online(self,dataOut):
4127 4133 utc_adq =dataOut.utctime
4128 4134 print("Online-utc_adq",utc_adq)
4129 4135 print(len(self.list_pedestal))
4130 4136 utc_ped_list=[]
4131 4137 for i in range(len(self.list_pedestal)):
4132 4138 utc_ped_list.append(self.gettimeutcfromDirFilename(path=self.path_ped,file=self.list_pedestal[i]))
4133 4139 print(utc_ped_list[:20])
4134 4140 #print(utc_ped_list[488:498])
4135 4141 print("ultimo UTC-PEDESTAL",utc_ped_list[-1])
4136 4142 nro_file,utc_ped,utc_ped_1 = self.getNROFile(utc_adq=utc_adq, utc_ped_list= utc_ped_list)
4137 4143 print("nro_file",nro_file,"utc_ped",utc_ped,"utc_ped_1",utc_ped_1)
4138 4144 print("name_PEDESTAL",self.list_pedestal[nro_file])
4139 4145 nro_key_p = int((utc_adq-utc_ped)/self.t_Interval_p)-1
4140 4146 print("nro_key_p",nro_key_p)
4141 4147 ff_pedestal = self.list_pedestal[nro_file]
4142 4148 #angulo = self.getDatavaluefromDirFilename(path=self.path_ped,file=ff_pedestal,value="azimuth")
4143 4149 angulo = self.getDatavaluefromDirFilename(path=self.path_ped,file=ff_pedestal,value="azi_pos")
4144 4150
4145 4151 print("utc_pedestal_init :",utc_ped+nro_key_p*self.t_Interval_p)
4146 4152 print("angulo_array :",angulo[nro_key_p])
4147 4153 self.nro_file = nro_file
4148 4154 self.nro_key_p = nro_key_p
4149 4155
4150 def setup(self,dataOut,path_ped,path_adq,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online):
4156 #def setup(self,dataOut,path_ped,path_adq,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online):
4157 def setup(self,dataOut,path_ped,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online):
4158 print("SETUP PEDESTAL")
4151 4159 self.__dataReady = False
4152 4160 self.path_ped = path_ped
4153 self.path_adq = path_adq
4161 #self.path_adq = path_adq
4154 4162 self.t_Interval_p = t_Interval_p
4155 4163 self.n_Muestras_p = n_Muestras_p
4156 4164 self.blocksPerfile= blocksPerfile
4157 4165 self.f_a_p = f_a_p
4158 4166 self.online = online
4159 4167 self.angulo_adq = numpy.zeros(self.blocksPerfile)
4160 4168 self.__profIndex = 0
4161 4169 self.tmp = 0
4162 4170 self.c_ped = 0
4163 4171 print(self.path_ped)
4164 print(self.path_adq)
4172 #print(self.path_adq)
4165 4173 self.list_pedestal = self.getfirstFilefromPath(path=self.path_ped,meta="PE",ext=".hdf5")
4166 4174 print("LIST NEW", self.list_pedestal[:20])
4167 self.list_adq = self.getfirstFilefromPath(path=self.path_adq,meta="D",ext=".hdf5")
4175 #self.list_adq = self.getfirstFilefromPath(path=self.path_adq,meta="D",ext=".hdf5")
4168 4176 print("*************Longitud list pedestal****************",len(self.list_pedestal))
4169 4177
4170 4178 if self.online:
4171 4179 print("Enable Online")
4172 4180 self.setup_online(dataOut)
4173 4181 else:
4174 self.setup_offline(dataOut,list_pedestal=self.list_pedestal,list_adq=self.list_adq)
4182 #self.setup_offline(dataOut,list_pedestal=self.list_pedestal,list_adq=self.list_adq)
4183 self.setup_offline(dataOut,list_pedestal=self.list_pedestal)
4184
4175 4185
4176 4186 def setNextFileP(self,dataOut):
4177 4187 if self.online:
4178 4188 data_pedestal = self.setNextFileonline()
4179 4189 else:
4180 4190 data_pedestal = self.setNextFileoffline(dataOut)
4181 4191
4182 4192 return data_pedestal
4183 4193
4184 4194
4185 4195 def setNextFileoffline(self,dataOut):
4186 4196 ##tmp=0
4187 4197 for j in range(self.blocksPerfile):
4188 4198 ###print("NUMERO DEL BLOQUE---->",j)
4189 4199 ###print("nro_key_p",self.nro_key_p)
4190 4200
4191 4201 #iterador = self.nro_key_p +self.f_a_p*(j-tmp)
4192 4202 iterador = self.nro_key_p +self.f_a_p*self.c_ped
4193 4203 self.c_ped = self.c_ped +1
4194 4204
4195 ###print("iterador------------->",iterador)
4205 print("iterador------------->",iterador)
4196 4206 if iterador < self.n_Muestras_p:
4197 4207 self.nro_file = self.nro_file
4198 4208 else:
4199 4209 self.nro_file = self.nro_file+1
4200 4210 print("PRUEBA-------------")
4201 4211 utc_ped_setnext=self.gettimeutcfromDirFilename(path=self.path_ped,file=self.list_pedestal[self.nro_file])
4202 4212 utc_adq_setnext=dataOut.utctime
4203 4213 print("utc_pedestal",utc_ped_setnext)
4204 4214 print("utc_adq",utc_adq_setnext)
4205 4215
4216 print("self.c_ped",self.c_ped)
4217 #dif = self.blocksPerfile-(self.nro_key_p+self.f_a_p*(self.c_ped-2))
4218 dif = self.n_Muestras_p-(self.nro_key_p+self.f_a_p*(self.c_ped-2))
4206 4219
4207 dif = self.blocksPerfile-(self.nro_key_p+self.f_a_p*(self.c_ped-2))
4208 4220 self.c_ped = 1
4209 4221 ##tmp = j
4210 4222 ##print("tmp else",tmp)
4211 4223 self.nro_key_p= self.f_a_p-dif
4212 4224 iterador = self.nro_key_p
4213 ###print("iterador else",iterador)
4225 print("iterador else",iterador)
4214 4226 #self.c_ped = self.c_ped +1
4215 4227
4216 ###print("nro_file",self.nro_file)
4228 print("nro_file",self.nro_file)
4217 4229 #print("tmp",tmp)
4218 4230 try:
4219 4231 ff_pedestal = self.list_pedestal[self.nro_file]
4220 4232 print("ff_pedestal",ff_pedestal)
4221 4233 except:
4222 4234 print("############# EXCEPCION ######################")
4223 4235 return numpy.ones(self.blocksPerfile)*numpy.nan
4224 4236
4225 4237 #angulo = self.getDatavaluefromDirFilename(path=self.path_ped,file=ff_pedestal,value="azimuth")
4226 4238 angulo = self.getDatavaluefromDirFilename(path=self.path_ped,file=ff_pedestal,value="azi_pos")
4227 4239
4228 4240 self.angulo_adq[j]= angulo[iterador]
4229 4241
4230 4242 return self.angulo_adq
4231 4243
4232 4244 def setNextFileonline(self):
4233 4245 tmp = 0
4234 4246 self.nTries_p = 3
4235 4247 self.delay = 3
4236 4248 ready = 1
4237 4249 for j in range(self.blocksPerfile):
4238 4250 iterador = self.nro_key_p +self.f_a_p*(j-tmp)
4239 4251 if iterador < self.n_Muestras_p:
4240 4252 self.nro_file = self.nro_file
4241 4253 else:
4242 4254 self.nro_file = self.nro_file+1
4243 4255 dif = self.blocksPerfile-(self.nro_key_p+self.f_a_p*(j-tmp-1))
4244 4256 tmp = j
4245 4257 self.nro_key_p= self.f_a_p-dif
4246 4258 iterador = self.nro_key_p
4247 4259 #print("nro_file---------------- :",self.nro_file)
4248 4260 try:
4249 4261 # update list_pedestal
4250 4262 self.list_pedestal = self.getfirstFilefromPath(path=self.path_ped,meta="PE",ext=".hdf5")
4251 4263 ff_pedestal = self.list_pedestal[self.nro_file]
4252 4264 except:
4253 4265 ff_pedestal = None
4254 4266 ready = 0
4255 4267 for nTries_p in range(self.nTries_p):
4256 4268 try:
4257 4269 # update list_pedestal
4258 4270 self.list_pedestal = self.getfirstFilefromPath(path=self.path_ped,meta="PE",ext=".hdf5")
4259 4271 ff_pedestal = self.list_pedestal[self.nro_file]
4260 4272 except:
4261 4273 ff_pedestal = None
4262 4274 if ff_pedestal is not None:
4263 4275 ready=1
4264 4276 break
4265 4277 log.warning("Waiting %0.2f sec for the next file: \"%s\" , try %02d ..." % (self.delay, self.nro_file, nTries_p + 1))
4266 4278 time.sleep(self.delay)
4267 4279 continue
4268 4280 #return numpy.ones(self.blocksPerfile)*numpy.nan
4269 4281
4270 4282 if ready == 1:
4271 4283 #angulo = self.getDatavaluefromDirFilename(path=self.path_ped,file=ff_pedestal,value="azimuth")
4272 4284 angulo = self.getDatavaluefromDirFilename(path=self.path_ped,file=ff_pedestal,value="azi_pos")
4273 4285
4274 4286 else:
4275 4287 print("there is no pedestal file")
4276 4288 angulo = numpy.ones(self.n_Muestras_p)*numpy.nan
4277 4289 self.angulo_adq[j]= angulo[iterador]
4278 4290 ####print("Angulo",self.angulo_adq)
4279 4291 ####print("Angulo",len(self.angulo_adq))
4280 4292 #self.nro_key_p=iterador + self.f_a_p
4281 4293 #if self.nro_key_p< self.n_Muestras_p:
4282 4294 # self.nro_file = self.nro_file
4283 4295 #else:
4284 4296 # self.nro_file = self.nro_file+1
4285 4297 # self.nro_key_p= self.nro_key_p
4286 4298 return self.angulo_adq
4287 4299
4288 4300
4289 def run(self, dataOut,path_ped,path_adq,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online):
4301 #def run(self, dataOut,path_ped,path_adq,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online):
4302 def run(self, dataOut,path_ped,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online):
4303
4290 4304 if not self.isConfig:
4291 4305 print("######################SETUP#########################################")
4292 self.setup( dataOut, path_ped,path_adq,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online)
4306 #self.setup( dataOut, path_ped,path_adq,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online)
4307 self.setup( dataOut, path_ped,t_Interval_p,n_Muestras_p,blocksPerfile,f_a_p,online)
4293 4308 self.isConfig = True
4294 4309
4295 4310 dataOut.flagNoData = True
4296 #print("profIndex",self.__profIndex)
4311 print("profIndex",self.__profIndex)
4297 4312
4298 4313 if self.__profIndex==0:
4299 4314 angulo_adq = self.setNextFileP(dataOut)
4300 4315 dataOut.azimuth = angulo_adq
4301 4316 print("TIEMPO:",dataOut.utctime)
4302 4317 ##print("####################################################################")
4303 ##print("angulos",dataOut.azimuth,len(dataOut.azimuth))
4318 print("angulos",dataOut.azimuth,len(dataOut.azimuth))
4304 4319 self.__dataReady = True
4305 4320 self.__profIndex += 1
4306 4321 print("TIEMPO_bucle:",dataOut.utctime)
4322 print("profIndex",self.__profIndex)
4307 4323 if self.__profIndex== blocksPerfile:
4308 4324 self.__profIndex = 0
4309 4325 if self.__dataReady:
4310 4326 #print(self.__profIndex,dataOut.azimuth[:10])
4311 4327 dataOut.flagNoData = False
4312 4328 return dataOut
4313 4329
4314 4330
4315 4331 class Block360(Operation):
4316 4332 '''
4317 4333 '''
4318 4334 isConfig = False
4319 4335 __profIndex = 0
4320 4336 __initime = None
4321 4337 __lastdatatime = None
4322 4338 __buffer = None
4323 4339 __dataReady = False
4324 4340 n = None
4325 4341 __nch = 0
4326 4342 __nHeis = 0
4327 4343 index = 0
4344 mode = 0
4328 4345
4329 4346 def __init__(self,**kwargs):
4330 4347 Operation.__init__(self,**kwargs)
4331 4348
4332 def setup(self, dataOut, n = None):
4349 def setup(self, dataOut, n = None, mode = None):
4333 4350 '''
4334 4351 n= Numero de PRF's de entrada
4335 4352 '''
4336 4353 self.__initime = None
4337 4354 self.__lastdatatime = 0
4338 4355 self.__dataReady = False
4339 4356 self.__buffer = 0
4340 4357 self.__buffer_1D = 0
4341 4358 self.__profIndex = 0
4342 4359 self.index = 0
4343 4360 self.__nch = dataOut.nChannels
4344 4361 self.__nHeis = dataOut.nHeights
4345 4362 ##print("ELVALOR DE n es:", n)
4346 4363 if n == None:
4347 4364 raise ValueError("n should be specified.")
4348 4365
4366 if mode == None:
4367 raise ValueError("mode should be specified.")
4368
4349 4369 if n != None:
4350 4370 if n<1:
4351 4371 print("n should be greater than 2")
4352 4372 raise ValueError("n should be greater than 2")
4353 4373
4354 4374 self.n = n
4375 self.mode = mode
4376 print("self.mode",self.mode)
4355 4377 #print("nHeights")
4356 4378 self.__buffer = numpy.zeros(( dataOut.nChannels,n, dataOut.nHeights))
4357 4379 self.__buffer2= numpy.zeros(n)
4358 4380
4359 def putData(self,data):
4381 def putData(self,data,mode):
4360 4382 '''
4361 4383 Add a profile to he __buffer and increase in one the __profiel Index
4362 4384 '''
4363 4385 #print("line 4049",data.dataPP_POW.shape,data.dataPP_POW[:10])
4364 4386 #print("line 4049",data.azimuth.shape,data.azimuth)
4365 self.__buffer[:,self.__profIndex,:]= data.dataPP_POW
4387 if self.mode==0:
4388 self.__buffer[:,self.__profIndex,:]= data.dataPP_POW
4389 if self.mode==1:
4390 self.__buffer[:,self.__profIndex,:]= data.data_pow
4366 4391 #print("me casi",self.index,data.azimuth[self.index])
4367 4392 #print(self.__profIndex, self.index , data.azimuth[self.index] )
4368 4393 #print("magic",data.profileIndex)
4369 4394 #print(data.azimuth[self.index])
4370 4395 #print("index",self.index)
4371 4396
4372 4397 self.__buffer2[self.__profIndex] = data.azimuth[self.index]
4373 4398 #print("q pasa")
4374 4399 self.index+=1
4375 4400 #print("index",self.index,data.azimuth[:10])
4376 4401 self.__profIndex += 1
4377 4402 return #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Remove DCΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
4378 4403
4379 4404 def pushData(self,data):
4380 4405 '''
4381 4406 Return the PULSEPAIR and the profiles used in the operation
4382 4407 Affected : self.__profileIndex
4383 4408 '''
4384 4409 #print("pushData")
4385 4410
4386 4411 data_360 = self.__buffer
4387 4412 data_p = self.__buffer2
4388 4413 n = self.__profIndex
4389 4414
4390 4415 self.__buffer = numpy.zeros((self.__nch, self.n,self.__nHeis))
4391 4416 self.__buffer2 = numpy.zeros(self.n)
4392 4417 self.__profIndex = 0
4393 4418 #print("pushData")
4394 4419 return data_360,n,data_p
4395 4420
4396 4421
4397 4422 def byProfiles(self,dataOut):
4398 4423
4399 4424 self.__dataReady = False
4400 4425 data_360 = None
4401 4426 data_p = None
4402 4427 #print("dataOu",dataOut.dataPP_POW)
4403 self.putData(data=dataOut)
4428 self.putData(data=dataOut,mode = self.mode)
4404 4429 #print("profIndex",self.__profIndex)
4405 4430 if self.__profIndex == self.n:
4406 4431 data_360,n,data_p = self.pushData(data=dataOut)
4407 4432 self.__dataReady = True
4408 4433
4409 4434 return data_360,data_p
4410 4435
4411 4436
4412 4437 def blockOp(self, dataOut, datatime= None):
4413 4438 if self.__initime == None:
4414 4439 self.__initime = datatime
4415 4440 data_360,data_p = self.byProfiles(dataOut)
4416 4441 self.__lastdatatime = datatime
4417 4442
4418 4443 if data_360 is None:
4419 4444 return None, None,None
4420 4445
4421 4446 avgdatatime = self.__initime
4422 4447 deltatime = datatime - self.__lastdatatime
4423 4448 self.__initime = datatime
4424 4449 #print(data_360.shape,avgdatatime,data_p.shape)
4425 4450 return data_360,avgdatatime,data_p
4426 4451
4427 def run(self, dataOut,n = None,**kwargs):
4428
4452 def run(self, dataOut,n = None,mode=None,**kwargs):
4453 print("BLOCK 360 HERE WE GO MOMENTOS")
4429 4454 if not self.isConfig:
4430 self.setup(dataOut = dataOut, n = n , **kwargs)
4455 self.setup(dataOut = dataOut, n = n ,mode= mode ,**kwargs)
4431 4456 self.index = 0
4432 4457 #print("comova",self.isConfig)
4433 4458 self.isConfig = True
4434 4459 if self.index==dataOut.azimuth.shape[0]:
4435 4460 self.index=0
4436 4461 data_360, avgdatatime,data_p = self.blockOp(dataOut, dataOut.utctime)
4437 4462 dataOut.flagNoData = True
4438 4463
4439 4464 if self.__dataReady:
4440 4465 dataOut.data_360 = data_360 # S
4441 4466 ##print("---------------------------------------------------------------------------------")
4442 4467 ##print("---------------------------DATAREADY---------------------------------------------")
4443 4468 ##print("---------------------------------------------------------------------------------")
4444 4469 ##print("data_360",dataOut.data_360.shape)
4445 4470 dataOut.data_azi = data_p
4446 4471 ##print("azi: ",dataOut.data_azi)
4447 4472 #print("jroproc_parameters",data_p[0],data_p[-1])#,data_360.shape,avgdatatime)
4448 4473 dataOut.utctime = avgdatatime
4449 4474 dataOut.flagNoData = False
4450 4475 return dataOut
Show file before | Show file after | Hide comments
@@ -1,75 +1,81
1 1 import os,sys
2 2 import datetime
3 3 import time
4 4 from schainpy.controller import Project
5 5 #path='/DATA_RM/TEST_HDF5/d2021200'
6 6 #path='/DATA_RM/TEST_HDF5/d2021200'
7 7 #path='/DATA_RM/TEST_HDF5/d2021214'
8 8 #path='/DATA_RM/TEST_HDF5/d2021229'
9 9
10 10 #path='/DATA_RM/TEST_HDF5/d2021231'
11 11 #path='/DATA_RM/TEST_HDF5/ADQ_OFFLINE/d2021231'
12 12 path='/DATA_RM/TEST_HDF5/d2021231'
13 13 #path='/DATA_RM/TEST_14_HDF5/d2021257'
14 14 ## TEST ULTIMA PRUEBA 22 DE SEPTIEMBRE
15 15 path = '/DATA_RM/TEST_HDF5_PP_22/d2021265'
16 #path = '/DATA_RM/TEST_HDF5_PP_100/d2021285'
17 path = '/DATA_RM/TEST_HDF5_PP/d2021285'
18
19
16 20 path_adq=path
17 21 #path_ped='/DATA_RM/TEST_PEDESTAL/P2021200'
18 22 #path_ped='/DATA_RM/TEST_PEDESTAL/P2021214'
19 23 #path_ped='/DATA_RM/TEST_PEDESTAL/P2021230'
20 24 #path_ped='/DATA_RM/TEST_PEDESTAL/P20210819'
21 25 #path_ped='/DATA_RM/TEST_PEDESTAL/P20210819-154315'
22 26 #path_ped='/DATA_RM/TEST_PEDESTAL/P20210914-162434'
23 27 #path_ped='/DATA_RM/TEST_PEDESTAL/PEDESTAL_OFFLINE/P20210819-161524'
24 28 #pruebas con perdida de datos
25 29 #path_ped='/DATA_RM/TEST_PEDESTAL/PEDESTAL_OFFLINE/P20210819-161524_TEST'
26 30 ## TEST ULTIMA PRUEBA 22 DE SEPTIEMBRE
27 path_ped='/DATA_RM/TEST_PEDESTAL/P20210922-122731'
31 path_ped='/DATA_RM/TEST_PEDESTAL/P20211012-082745'
28 32
29 33
30 34 figpath = '/home/soporte/Pictures'
31 35 desc = "Simulator Test"
32 36
33 37 controllerObj = Project()
34 38 controllerObj.setup(id='10',name='Test Simulator',description=desc)
35 39 readUnitConfObj = controllerObj.addReadUnit(datatype='HDFReader',
36 40 path=path,
37 41 startDate="2021/01/01", #"2020/01/01",#today,
38 42 endDate= "2021/12/01", #"2020/12/30",#today,
39 43 startTime='00:00:00',
40 44 endTime='23:59:59',
41 45 t_Interval_p=0.01,
42 46 n_Muestras_p=100,
43 47 delay=30,
44 48 #set=0,
45 49 online=0,
46 50 walk=0,
47 51 nTries=6)#1
48 52
49 53 procUnitConfObjA = controllerObj.addProcUnit(datatype='ParametersProc',inputId=readUnitConfObj.getId())
50 V=2
51 blocksPerfile=100
54 V=10
55 blocksPerfile=360
52 56 print("Velocidad del Pedestal",V)
53 f_a_p= int(blocksPerfile/V)
57 tmp_blocksPerfile=100
58 f_a_p= int(tmp_blocksPerfile/V)
54 59
55 60 opObj11 = procUnitConfObjA.addOperation(name='PedestalInformation')
56 61 opObj11.addParameter(name='path_ped', value=path_ped)
57 opObj11.addParameter(name='path_adq', value=path_adq)
62 #opObj11.addParameter(name='path_adq', value=path_adq)
58 63 opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
64 opObj11.addParameter(name='blocksPerfile', value=blocksPerfile, format='int')
59 65 opObj11.addParameter(name='n_Muestras_p', value='100', format='float')
60 opObj11.addParameter(name='blocksPerfile', value='100', format='int')
61 66 opObj11.addParameter(name='f_a_p', value=f_a_p, format='int')
62 67 opObj11.addParameter(name='online', value='0', format='int')# habilitar el enable aqui tambien
63 68
64 69
65 70 opObj11 = procUnitConfObjA.addOperation(name='Block360')
66 71 opObj11.addParameter(name='n', value='10', format='int')
72 opObj11.addParameter(name='mode', value=0, format='int')
67 73 # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
68 74
69 75 opObj11= procUnitConfObjA.addOperation(name='WeatherPlot',optype='other')
70 76 #opObj11.addParameter(name='save', value=figpath)
71 77 #opObj11.addParameter(name='save_period', value=1)
72 78
73 79 controllerObj.start()
74 80 #online 1 utc_adq 1617490240.48
75 81 #online 0 utc_adq 1617489815.4804
General Comments 0
You need to be logged in to leave comments. Login now