##// END OF EJS Templates
Merge branch 'v3.0-devel' of http://jro-dev.igp.gob.pe/rhodecode/schain into v3.0-devel
Juan C. Espinoza -
r1323:133a5771abd8 merge
parent child
Show More

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

@@ -0,0 +1,87
1 import numpy
2 import sys
3 import zmq
4 import time
5 import h5py
6 import os
7
8 path="/home/alex/Downloads/pedestal"
9 ext=".hdf5"
10
11 port ="5556"
12 if len(sys.argv)>1:
13 port = sys.argv[1]
14 int(port)
15
16 if len(sys.argv)>2:
17 port1 = sys.argv[2]
18 int(port1)
19
20 #Socket to talk to server
21 context = zmq.Context()
22 socket = context.socket(zmq.SUB)
23
24 print("Collecting updates from weather server...")
25 socket.connect("tcp://localhost:%s"%port)
26
27 if len(sys.argv)>2:
28 socket.connect("tcp://localhost:%s"%port1)
29
30 #Subscribe to zipcode, default is NYC,10001
31 topicfilter = "10001"
32 socket.setsockopt_string(zmq.SUBSCRIBE,topicfilter)
33 #Process 5 updates
34 total_value=0
35 count= -1
36 azi= []
37 elev=[]
38 time0=[]
39 #for update_nbr in range(250):
40 while(True):
41 string= socket.recv()
42 topic,ang_elev,ang_elev_dec,ang_azi,ang_azi_dec,seconds,seconds_dec= string.split()
43 ang_azi =float(ang_azi)+1e-3*float(ang_azi_dec)
44 ang_elev =float(ang_elev)+1e-3*float(ang_elev_dec)
45 seconds =float(seconds) +1e-6*float(seconds_dec)
46 azi.append(ang_azi)
47 elev.append(ang_elev)
48 time0.append(seconds)
49 count +=1
50 if count == 100:
51 timetuple=time.localtime()
52 epoc = time.mktime(timetuple)
53 #print(epoc)
54 fullpath = path + ("/" if path[-1]!="/" else "")
55
56 if not os.path.exists(fullpath):
57 os.mkdir(fullpath)
58
59 azi_array = numpy.array(azi)
60 elev_array = numpy.array(elev)
61 time0_array= numpy.array(time0)
62 pedestal_array=numpy.array([azi,elev,time0])
63 count=0
64 azi= []
65 elev=[]
66 time0=[]
67 #print(pedestal_array[0])
68 #print(pedestal_array[1])
69
70 meta='PE'
71 filex="%s%4.4d%3.3d%10.4d%s"%(meta,timetuple.tm_year,timetuple.tm_yday,epoc,ext)
72 filename = os.path.join(fullpath,filex)
73 fp = h5py.File(filename,'w')
74 #print("Escribiendo HDF5...",epoc)
75 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· DataΒ·....Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
76 grp = fp.create_group("Data")
77 dset = grp.create_dataset("azimuth" , data=pedestal_array[0])
78 dset = grp.create_dataset("elevacion", data=pedestal_array[1])
79 dset = grp.create_dataset("utc" , data=pedestal_array[2])
80 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· MetadataΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
81 grp = fp.create_group("Metadata")
82 dset = grp.create_dataset("utctimeInit", data=pedestal_array[2][0])
83 timeInterval = pedestal_array[2][1]-pedestal_array[2][0]
84 dset = grp.create_dataset("timeInterval", data=timeInterval)
85 fp.close()
86
87 #print ("Average messagedata value for topic '%s' was %dF" % ( topicfilter,total_value / update_nbr))
@@ -0,0 +1,48
1 ###########################################################################
2 ############################### SERVIDOR###################################
3 ######################### SIMULADOR DE PEDESTAL############################
4 ###########################################################################
5 import time
6 import math
7 import numpy
8 import struct
9 from time import sleep
10 import zmq
11 import pickle
12 port="5556"
13 context = zmq.Context()
14 socket = context.socket(zmq.PUB)
15 socket.bind("tcp://*:%s"%port)
16 ###### PARAMETROS DE ENTRADA################################
17 print("PEDESTAL RESOLUCION 0.01")
18 print("MAXIMA VELOCIDAD DEL PEDESTAL")
19 ang_elev = 4.12
20 ang_azi = 30
21 velocidad= input ("Ingresa velocidad:")
22 velocidad= float(velocidad)
23 print (velocidad)
24 ############################################################
25 sleep(3)
26 print("Start program")
27 t1 = time.time()
28 count=0
29 while(True):
30 tmp_vuelta = int(360/velocidad)
31 t1=t1+tmp_vuelta*count
32 count= count+1
33 muestras_seg = 100
34 t2 = time.time()
35 for i in range(tmp_vuelta):
36 for j in range(muestras_seg):
37 tmp_variable = (i+j/100.0)
38 ang_azi = (tmp_variable)*float(velocidad)
39 seconds = t1+ tmp_variable
40 topic=10001
41 print ("AzimΒ°: ","%.4f"%ang_azi,"Time:" ,"%.5f"%seconds)
42 seconds_dec=(seconds-int(seconds))*1e6
43 ang_azi_dec= (ang_azi-int(ang_azi))*1e3
44 ang_elev_dec=(ang_elev-int(ang_elev))*1e3
45 sleep(0.0088)
46 socket.send_string("%d %d %d %d %d %d %d"%(topic,ang_elev,ang_elev_dec,ang_azi,ang_azi_dec,seconds,seconds_dec))
47 t3 = time.time()
48 print ("Total time for 1 vuelta in Seconds",t3-t2)
@@ -0,0 +1,82
1 import os, sys
2 import datetime
3 import time
4 from schainpy.controller import Project
5
6 desc = "USRP_test"
7 filename = "USRP_processing.xml"
8 controllerObj = Project()
9 controllerObj.setup(id = '191', name='Test_USRP', description=desc)
10
11 ############## USED TO PLOT IQ VOLTAGE, POWER AND SPECTRA #############
12 ######PATH DE LECTURA, ESCRITURA, GRAFICOS Y ENVIO WEB#################
13 path = '/home/alex/Downloads/test_rawdata'
14 figpath = '/home/alex/Downloads/hdf5_test'
15 ######################## UNIDAD DE LECTURA#############################
16 '''
17 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
18 path=path,
19 startDate="2020/01/01", #"2020/01/01",#today,
20 endDate= "2020/12/01", #"2020/12/30",#today,
21 startTime='00:00:00',
22 endTime='23:59:59',
23 delay=0,
24 #set=0,
25 online=0,
26 walk=1)
27
28 '''
29 readUnitConfObj = controllerObj.addReadUnit(datatype='SimulatorReader',
30 frequency=9.345e9,
31 FixRCP_IPP= 60,
32 Tau_0 = 30,
33 AcqH0_0=0,
34 samples=330,
35 AcqDH_0=0.15,
36 FixRCP_TXA=0.15,
37 FixRCP_TXB=0.15,
38 Fdoppler=600.0,
39 Hdoppler=36,
40 Adoppler=300,#300
41 delay=0,
42 online=0,
43 walk=0,
44 profilesPerBlock=625,
45 dataBlocksPerFile=100)
46 #nTotalReadFiles=2)
47
48
49 #opObj11 = readUnitConfObj.addOperation(name='printInfo')
50
51 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
52
53 procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
54 procUnitConfObjB.addParameter(name='nFFTPoints', value=625, format='int')
55 procUnitConfObjB.addParameter(name='nProfiles', value=625, format='int')
56
57 opObj11 = procUnitConfObjB.addOperation(name='removeDC')
58 opObj11.addParameter(name='mode', value=2)
59 #opObj11 = procUnitConfObjB.addOperation(name='SpectraPlot')
60 #opObj11 = procUnitConfObjB.addOperation(name='PowerProfilePlot')
61
62 procUnitConfObjC= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjB.getId())
63 procUnitConfObjC.addOperation(name='SpectralMoments')
64 #opObj11 = procUnitConfObjC.addOperation(name='PowerPlot')
65
66 '''
67 opObj11 = procUnitConfObjC.addOperation(name='SpectralMomentsPlot')
68 #opObj11.addParameter(name='xmin', value=14)
69 #opObj11.addParameter(name='xmax', value=15)
70 #opObj11.addParameter(name='save', value=figpath)
71 opObj11.addParameter(name='showprofile', value=1)
72 #opObj11.addParameter(name='save_period', value=10)
73 '''
74
75 opObj10 = procUnitConfObjC.addOperation(name='ParameterWriter')
76 opObj10.addParameter(name='path',value=figpath)
77 #opObj10.addParameter(name='mode',value=0)
78 opObj10.addParameter(name='blocksPerFile',value='100',format='int')
79 opObj10.addParameter(name='metadataList',value='utctimeInit,timeInterval',format='list')
80 opObj10.addParameter(name='dataList',value='data_POW,data_DOP,data_WIDTH,data_SNR')#,format='list'
81
82 controllerObj.start()
@@ -0,0 +1,162
1 import os,numpy,h5py
2 from shutil import copyfile
3
4 def isNumber(str):
5 try:
6 float(str)
7 return True
8 except:
9 return False
10
11 def getfirstFilefromPath(path,meta,ext):
12 validFilelist = []
13 fileList = os.listdir(path)
14 if len(fileList)<1:
15 return None
16 # meta 1234 567 8-18 BCDE
17 # H,D,PE YYYY DDD EPOC .ext
18
19 for thisFile in fileList:
20 if meta =="PE":
21 try:
22 number= int(thisFile[len(meta)+7:len(meta)+17])
23 except:
24 print("There is a file or folder with different format")
25 if meta == "D":
26 try:
27 number= int(thisFile[8:11])
28 except:
29 print("There is a file or folder with different format")
30
31 if not isNumber(str=number):
32 continue
33 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
34 continue
35 validFilelist.sort()
36 validFilelist.append(thisFile)
37 if len(validFilelist)>0:
38 validFilelist = sorted(validFilelist,key=str.lower)
39 return validFilelist
40 return None
41
42 def gettimeutcfromDirFilename(path,file):
43 dir_file= path+"/"+file
44 fp = h5py.File(dir_file,'r')
45 epoc = fp['Metadata'].get('utctimeInit')[()]
46 fp.close()
47 return epoc
48
49 def getDatavaluefromDirFilename(path,file,value):
50 dir_file= path+"/"+file
51 fp = h5py.File(dir_file,'r')
52 array = fp['Data'].get(value)[()]
53 fp.close()
54 return array
55
56
57 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Velocidad de PedestalΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
58 w = input ("Ingresa velocidad de Pedestal: ")
59 w = 4
60 w = float(w)
61 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Resolucion minimo en gradosΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
62 alfa = input ("Ingresa resolucion minima en grados: ")
63 alfa = 1
64 alfa = float(alfa)
65 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· IPP del Experimento Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
66 IPP = input ("Ingresa el IPP del experimento: ")
67 IPP = 0.0004
68 IPP = float(IPP)
69 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· MODE Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
70 mode = input ("Ingresa el MODO del experimento T or F: ")
71 mode = "T"
72 mode = str(mode)
73
74 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Tiempo en generar la resolucion minΒ·Β·Β·
75 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· MCU Β·Β· var_ang = w * (var_tiempo)Β·Β·Β·
76 var_tiempo = alfa/w
77 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Tiempo Equivalente en perfilesΒ·Β·Β·Β·Β·Β·Β·Β·
78 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· var_tiempo = IPP * ( num_perfiles )Β·
79 num_perfiles = int(var_tiempo/IPP)
80
81 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·DATA PEDESTALΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
82 dir_pedestal = "/home/alex/Downloads/pedestal"
83 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· DATA ADQΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
84 if mode=="T":
85 dir_adq = "/home/alex/Downloads/hdf5_testPP/d2020194" # Time domain
86 else:
87 dir_adq = "/home/alex/Downloads/hdf5_test/d2020194" # Frequency domain
88
89 print( "Velocidad angular :", w)
90 print( "Resolucion minima en grados :", alfa)
91 print( "Numero de perfiles equivalente:", num_perfiles)
92 print( "Mode :", mode)
93
94 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· First FileΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
95 list_pedestal = getfirstFilefromPath(path=dir_pedestal,meta="PE",ext=".hdf5")
96 list_adq = getfirstFilefromPath(path=dir_adq ,meta="D",ext=".hdf5")
97
98 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· utc time Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
99 utc_pedestal= gettimeutcfromDirFilename(path=dir_pedestal,file=list_pedestal[0])
100 utc_adq = gettimeutcfromDirFilename(path=dir_adq ,file=list_adq[0])
101
102 print("utc_pedestal :",utc_pedestal)
103 print("utc_adq :",utc_adq)
104 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Relacion: utc_adq (+/-) var_tiempo*nro_file= utc_pedestal
105 time_Interval_p = 0.01
106 n_perfiles_p = 100
107 if utc_adq>utc_pedestal:
108 nro_file = int((int(utc_adq) - int(utc_pedestal))/(time_Interval_p*n_perfiles_p))
109 ff_pedestal = list_pedestal[nro_file]
110 utc_pedestal = gettimeutcfromDirFilename(path=dir_pedestal,file=ff_pedestal)
111 nro_key_p = int((utc_adq-utc_pedestal)/time_Interval_p)
112 if utc_adq >utc_pedestal:
113 ff_pedestal = ff_pedestal
114 else:
115 nro_file = nro_file-1
116 ff_pedestal = list_pedestal[nro_file]
117 angulo = getDatavaluefromDirFilename(path=dir_pedestal,file=ff_pedestal,value="azimuth")
118 nro_key_p = int((utc_adq-utc_pedestal)/time_Interval_p)
119 print("nro_file :",nro_file)
120 print("name_file :",ff_pedestal)
121 print("utc_pedestal_file :",utc_pedestal)
122 print("nro_key_p :",nro_key_p)
123 print("utc_pedestal_init :",utc_pedestal+nro_key_p*time_Interval_p)
124 print("angulo_array :",angulo[nro_key_p])
125 #4+25+25+25+21
126 #while True:
127 list_pedestal = getfirstFilefromPath(path=dir_pedestal,meta="PE",ext=".hdf5")
128 list_adq = getfirstFilefromPath(path=dir_adq ,meta="D",ext=".hdf5")
129
130 nro_file = nro_file #10
131 nro_key_perfil = nro_key_p
132 blocksPerFile = 100
133 wr_path = "/home/alex/Downloads/hdf5_wr/"
134 # Lectura de archivos de adquisicion para adicion de azimuth
135 for thisFile in range(len(list_adq)):
136 print("thisFileAdq",thisFile)
137 angulo_adq = numpy.zeros(blocksPerFile)
138 tmp = 0
139 for j in range(blocksPerFile):
140 iterador = nro_key_perfil + 25*(j-tmp)
141 if iterador < n_perfiles_p:
142 nro_file = nro_file
143 else:
144 nro_file = nro_file+1
145 tmp = j
146 iterador = nro_key_perfil
147 ff_pedestal = list_pedestal[nro_file]
148 angulo = getDatavaluefromDirFilename(path=dir_pedestal,file=ff_pedestal,value="azimuth")
149 angulo_adq[j]= angulo[iterador]
150 copyfile(dir_adq+"/"+list_adq[thisFile],wr_path+list_adq[thisFile])
151 fp = h5py.File(wr_path+list_adq[thisFile],'a')
152 grp = fp.create_group("Pedestal")
153 dset = grp.create_dataset("azimuth" , data=angulo_adq)
154 fp.close()
155 print("Angulo",angulo_adq)
156 print("Angulo",len(angulo_adq))
157 nro_key_perfil=iterador + 25
158 if nro_key_perfil< n_perfiles_p:
159 nro_file = nro_file
160 else:
161 nro_file = nro_file+1
162 nro_key_perfil= nro_key_p
@@ -1,1397 +1,1403
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10 import json
11 11
12 12 import schainpy.admin
13 13 from schainpy.utils import log
14 14 from .jroheaderIO import SystemHeader, RadarControllerHeader
15 15 from schainpy.model.data import _noise
16 16
17 17
18 18 def getNumpyDtype(dataTypeCode):
19 19
20 20 if dataTypeCode == 0:
21 21 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
22 22 elif dataTypeCode == 1:
23 23 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
24 24 elif dataTypeCode == 2:
25 25 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
26 26 elif dataTypeCode == 3:
27 27 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
28 28 elif dataTypeCode == 4:
29 29 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
30 30 elif dataTypeCode == 5:
31 31 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
32 32 else:
33 33 raise ValueError('dataTypeCode was not defined')
34 34
35 35 return numpyDtype
36 36
37 37
38 38 def getDataTypeCode(numpyDtype):
39 39
40 40 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
41 41 datatype = 0
42 42 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
43 43 datatype = 1
44 44 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
45 45 datatype = 2
46 46 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
47 47 datatype = 3
48 48 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
49 49 datatype = 4
50 50 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
51 51 datatype = 5
52 52 else:
53 53 datatype = None
54 54
55 55 return datatype
56 56
57 57
58 58 def hildebrand_sekhon(data, navg):
59 59 """
60 60 This method is for the objective determination of the noise level in Doppler spectra. This
61 61 implementation technique is based on the fact that the standard deviation of the spectral
62 62 densities is equal to the mean spectral density for white Gaussian noise
63 63
64 64 Inputs:
65 65 Data : heights
66 66 navg : numbers of averages
67 67
68 68 Return:
69 69 mean : noise's level
70 70 """
71 71
72 72 sortdata = numpy.sort(data, axis=None)
73 73 '''
74 74 lenOfData = len(sortdata)
75 75 nums_min = lenOfData*0.2
76 76
77 77 if nums_min <= 5:
78 78
79 79 nums_min = 5
80 80
81 81 sump = 0.
82 82 sumq = 0.
83 83
84 84 j = 0
85 85 cont = 1
86 86
87 87 while((cont == 1)and(j < lenOfData)):
88 88
89 89 sump += sortdata[j]
90 90 sumq += sortdata[j]**2
91 91
92 92 if j > nums_min:
93 93 rtest = float(j)/(j-1) + 1.0/navg
94 94 if ((sumq*j) > (rtest*sump**2)):
95 95 j = j - 1
96 96 sump = sump - sortdata[j]
97 97 sumq = sumq - sortdata[j]**2
98 98 cont = 0
99 99
100 100 j += 1
101 101
102 102 lnoise = sump / j
103 103 '''
104 104 return _noise.hildebrand_sekhon(sortdata, navg)
105 105
106 106
107 107 class Beam:
108 108
109 109 def __init__(self):
110 110 self.codeList = []
111 111 self.azimuthList = []
112 112 self.zenithList = []
113 113
114 114
115 115 class GenericData(object):
116 116
117 117 flagNoData = True
118 118
119 119 def copy(self, inputObj=None):
120 120
121 121 if inputObj == None:
122 122 return copy.deepcopy(self)
123 123
124 124 for key in list(inputObj.__dict__.keys()):
125 125
126 126 attribute = inputObj.__dict__[key]
127 127
128 128 # If this attribute is a tuple or list
129 129 if type(inputObj.__dict__[key]) in (tuple, list):
130 130 self.__dict__[key] = attribute[:]
131 131 continue
132 132
133 133 # If this attribute is another object or instance
134 134 if hasattr(attribute, '__dict__'):
135 135 self.__dict__[key] = attribute.copy()
136 136 continue
137 137
138 138 self.__dict__[key] = inputObj.__dict__[key]
139 139
140 140 def deepcopy(self):
141 141
142 142 return copy.deepcopy(self)
143 143
144 144 def isEmpty(self):
145 145
146 146 return self.flagNoData
147 147
148 148 def isReady(self):
149 149
150 150 return not self.flagNoData
151 151
152 152
153 153 class JROData(GenericData):
154 154
155 155 # m_BasicHeader = BasicHeader()
156 156 # m_ProcessingHeader = ProcessingHeader()
157 157
158 158 systemHeaderObj = SystemHeader()
159 159 radarControllerHeaderObj = RadarControllerHeader()
160 160 # data = None
161 161 type = None
162 162 datatype = None # dtype but in string
163 163 # dtype = None
164 164 # nChannels = None
165 165 # nHeights = None
166 166 nProfiles = None
167 167 heightList = None
168 168 channelList = None
169 169 flagDiscontinuousBlock = False
170 170 useLocalTime = False
171 171 utctime = None
172 172 timeZone = None
173 173 dstFlag = None
174 174 errorCount = None
175 175 blocksize = None
176 176 # nCode = None
177 177 # nBaud = None
178 178 # code = None
179 179 flagDecodeData = False # asumo q la data no esta decodificada
180 180 flagDeflipData = False # asumo q la data no esta sin flip
181 181 flagShiftFFT = False
182 182 # ippSeconds = None
183 183 # timeInterval = None
184 184 nCohInt = None
185 185 # noise = None
186 186 windowOfFilter = 1
187 187 # Speed of ligth
188 188 C = 3e8
189 189 frequency = 49.92e6
190 190 realtime = False
191 191 beacon_heiIndexList = None
192 192 last_block = None
193 193 blocknow = None
194 194 azimuth = None
195 195 zenith = None
196 196 beam = Beam()
197 197 profileIndex = None
198 198 error = None
199 199 data = None
200 200 nmodes = None
201 201
202 202 def __str__(self):
203 203
204 204 return '{} - {}'.format(self.type, self.getDatatime())
205 205
206 206 def getNoise(self):
207 207
208 208 raise NotImplementedError
209 209
210 210 def getNChannels(self):
211 211
212 212 return len(self.channelList)
213 213
214 214 def getChannelIndexList(self):
215 215
216 216 return list(range(self.nChannels))
217 217
218 218 def getNHeights(self):
219 219
220 220 return len(self.heightList)
221 221
222 222 def getHeiRange(self, extrapoints=0):
223 223
224 224 heis = self.heightList
225 225 # deltah = self.heightList[1] - self.heightList[0]
226 226 #
227 227 # heis.append(self.heightList[-1])
228 228
229 229 return heis
230 230
231 231 def getDeltaH(self):
232 232
233 233 delta = self.heightList[1] - self.heightList[0]
234 234
235 235 return delta
236 236
237 237 def getltctime(self):
238 238
239 239 if self.useLocalTime:
240 240 return self.utctime - self.timeZone * 60
241 241
242 242 return self.utctime
243 243
244 244 def getDatatime(self):
245 245
246 246 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
247 247 return datatimeValue
248 248
249 249 def getTimeRange(self):
250 250
251 251 datatime = []
252 252
253 253 datatime.append(self.ltctime)
254 254 datatime.append(self.ltctime + self.timeInterval + 1)
255 255
256 256 datatime = numpy.array(datatime)
257 257
258 258 return datatime
259 259
260 260 def getFmaxTimeResponse(self):
261 261
262 262 period = (10**-6) * self.getDeltaH() / (0.15)
263 263
264 264 PRF = 1. / (period * self.nCohInt)
265 265
266 266 fmax = PRF
267 267
268 268 return fmax
269 269
270 270 def getFmax(self):
271 271 PRF = 1. / (self.ippSeconds * self.nCohInt)
272 272
273 273 fmax = PRF
274 274 return fmax
275 275
276 276 def getVmax(self):
277 277
278 278 _lambda = self.C / self.frequency
279 279
280 280 vmax = self.getFmax() * _lambda / 2
281 281
282 282 return vmax
283 283
284 284 def get_ippSeconds(self):
285 285 '''
286 286 '''
287 287 return self.radarControllerHeaderObj.ippSeconds
288 288
289 289 def set_ippSeconds(self, ippSeconds):
290 290 '''
291 291 '''
292 292
293 293 self.radarControllerHeaderObj.ippSeconds = ippSeconds
294 294
295 295 return
296 296
297 297 def get_dtype(self):
298 298 '''
299 299 '''
300 300 return getNumpyDtype(self.datatype)
301 301
302 302 def set_dtype(self, numpyDtype):
303 303 '''
304 304 '''
305 305
306 306 self.datatype = getDataTypeCode(numpyDtype)
307 307
308 308 def get_code(self):
309 309 '''
310 310 '''
311 311 return self.radarControllerHeaderObj.code
312 312
313 313 def set_code(self, code):
314 314 '''
315 315 '''
316 316 self.radarControllerHeaderObj.code = code
317 317
318 318 return
319 319
320 320 def get_ncode(self):
321 321 '''
322 322 '''
323 323 return self.radarControllerHeaderObj.nCode
324 324
325 325 def set_ncode(self, nCode):
326 326 '''
327 327 '''
328 328 self.radarControllerHeaderObj.nCode = nCode
329 329
330 330 return
331 331
332 332 def get_nbaud(self):
333 333 '''
334 334 '''
335 335 return self.radarControllerHeaderObj.nBaud
336 336
337 337 def set_nbaud(self, nBaud):
338 338 '''
339 339 '''
340 340 self.radarControllerHeaderObj.nBaud = nBaud
341 341
342 342 return
343 343
344 344 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
345 345 channelIndexList = property(
346 346 getChannelIndexList, "I'm the 'channelIndexList' property.")
347 347 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
348 348 #noise = property(getNoise, "I'm the 'nHeights' property.")
349 349 datatime = property(getDatatime, "I'm the 'datatime' property")
350 350 ltctime = property(getltctime, "I'm the 'ltctime' property")
351 351 ippSeconds = property(get_ippSeconds, set_ippSeconds)
352 352 dtype = property(get_dtype, set_dtype)
353 353 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
354 354 code = property(get_code, set_code)
355 355 nCode = property(get_ncode, set_ncode)
356 356 nBaud = property(get_nbaud, set_nbaud)
357 357
358 358
359 359 class Voltage(JROData):
360 360
361 361 # data es un numpy array de 2 dmensiones (canales, alturas)
362 362 data = None
363 data_intensity = None
364 data_velocity = None
365 data_specwidth = None
363 dataPP_POW = None
364 dataPP_DOP = None
365 dataPP_WIDTH = None
366 dataPP_SNR = None
367
366 368 def __init__(self):
367 369 '''
368 370 Constructor
369 371 '''
370 372
371 373 self.useLocalTime = True
372 374 self.radarControllerHeaderObj = RadarControllerHeader()
373 375 self.systemHeaderObj = SystemHeader()
374 376 self.type = "Voltage"
375 377 self.data = None
376 378 # self.dtype = None
377 379 # self.nChannels = 0
378 380 # self.nHeights = 0
379 381 self.nProfiles = None
380 382 self.heightList = None
381 383 self.channelList = None
382 384 # self.channelIndexList = None
383 385 self.flagNoData = True
384 386 self.flagDiscontinuousBlock = False
385 387 self.utctime = None
386 388 self.timeZone = None
387 389 self.dstFlag = None
388 390 self.errorCount = None
389 391 self.nCohInt = None
390 392 self.blocksize = None
391 393 self.flagCohInt = False
392 394 self.flagDecodeData = False # asumo q la data no esta decodificada
393 395 self.flagDeflipData = False # asumo q la data no esta sin flip
394 396 self.flagShiftFFT = False
395 397 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
396 398 self.profileIndex = 0
397 399
398 400 def getNoisebyHildebrand(self, channel=None):
399 401 """
400 402 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
401 403
402 404 Return:
403 405 noiselevel
404 406 """
405 407
406 408 if channel != None:
407 409 data = self.data[channel]
408 410 nChannels = 1
409 411 else:
410 412 data = self.data
411 413 nChannels = self.nChannels
412 414
413 415 noise = numpy.zeros(nChannels)
414 416 power = data * numpy.conjugate(data)
415 417
416 418 for thisChannel in range(nChannels):
417 419 if nChannels == 1:
418 420 daux = power[:].real
419 421 else:
420 422 daux = power[thisChannel, :].real
421 423 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
422 424
423 425 return noise
424 426
425 427 def getNoise(self, type=1, channel=None):
426 428
427 429 if type == 1:
428 430 noise = self.getNoisebyHildebrand(channel)
429 431
430 432 return noise
431 433
432 434 def getPower(self, channel=None):
433 435
434 436 if channel != None:
435 437 data = self.data[channel]
436 438 else:
437 439 data = self.data
438 440
439 441 power = data * numpy.conjugate(data)
440 442 powerdB = 10 * numpy.log10(power.real)
441 443 powerdB = numpy.squeeze(powerdB)
442 444
443 445 return powerdB
444 446
445 447 def getTimeInterval(self):
446 448
447 449 timeInterval = self.ippSeconds * self.nCohInt
448 450
449 451 return timeInterval
450 452
451 453 noise = property(getNoise, "I'm the 'nHeights' property.")
452 454 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
453 455
454 456
455 457 class Spectra(JROData):
456 458
457 459 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
458 460 data_spc = None
459 461 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
460 462 data_cspc = None
461 463 # data dc es un numpy array de 2 dmensiones (canales, alturas)
462 464 data_dc = None
463 465 # data power
464 466 data_pwr = None
465 467 nFFTPoints = None
466 468 # nPairs = None
467 469 pairsList = None
468 470 nIncohInt = None
469 471 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
470 472 nCohInt = None # se requiere para determinar el valor de timeInterval
471 473 ippFactor = None
472 474 profileIndex = 0
473 475 plotting = "spectra"
474 476
475 477 def __init__(self):
476 478 '''
477 479 Constructor
478 480 '''
479 481
480 482 self.useLocalTime = True
481 483 self.radarControllerHeaderObj = RadarControllerHeader()
482 484 self.systemHeaderObj = SystemHeader()
483 485 self.type = "Spectra"
484 486 # self.data = None
485 487 # self.dtype = None
486 488 # self.nChannels = 0
487 489 # self.nHeights = 0
488 490 self.nProfiles = None
489 491 self.heightList = None
490 492 self.channelList = None
491 493 # self.channelIndexList = None
492 494 self.pairsList = None
493 495 self.flagNoData = True
494 496 self.flagDiscontinuousBlock = False
495 497 self.utctime = None
496 498 self.nCohInt = None
497 499 self.nIncohInt = None
498 500 self.blocksize = None
499 501 self.nFFTPoints = None
500 502 self.wavelength = None
501 503 self.flagDecodeData = False # asumo q la data no esta decodificada
502 504 self.flagDeflipData = False # asumo q la data no esta sin flip
503 505 self.flagShiftFFT = False
504 506 self.ippFactor = 1
505 507 #self.noise = None
506 508 self.beacon_heiIndexList = []
507 509 self.noise_estimation = None
508 510
509 511 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
510 512 """
511 513 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
512 514
513 515 Return:
514 516 noiselevel
515 517 """
516 518
517 519 noise = numpy.zeros(self.nChannels)
518 520
519 521 for channel in range(self.nChannels):
520 522 daux = self.data_spc[channel,
521 523 xmin_index:xmax_index, ymin_index:ymax_index]
522 524 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
523 525
524 526 return noise
525 527
526 528 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
527 529
528 530 if self.noise_estimation is not None:
529 531 # this was estimated by getNoise Operation defined in jroproc_spectra.py
530 532 return self.noise_estimation
531 533 else:
532 534 noise = self.getNoisebyHildebrand(
533 535 xmin_index, xmax_index, ymin_index, ymax_index)
534 536 return noise
535 537
536 538 def getFreqRangeTimeResponse(self, extrapoints=0):
537 539
538 540 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
539 541 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
540 542
541 543 return freqrange
542 544
543 545 def getAcfRange(self, extrapoints=0):
544 546
545 547 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
546 548 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
547 549
548 550 return freqrange
549 551
550 552 def getFreqRange(self, extrapoints=0):
551 553
552 554 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
553 555 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
554 556
555 557 return freqrange
556 558
557 559 def getVelRange(self, extrapoints=0):
558 560
559 561 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
560 562 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
561 563
562 564 if self.nmodes:
563 565 return velrange/self.nmodes
564 566 else:
565 567 return velrange
566 568
567 569 def getNPairs(self):
568 570
569 571 return len(self.pairsList)
570 572
571 573 def getPairsIndexList(self):
572 574
573 575 return list(range(self.nPairs))
574 576
575 577 def getNormFactor(self):
576 578
577 579 pwcode = 1
578 580
579 581 if self.flagDecodeData:
580 582 pwcode = numpy.sum(self.code[0]**2)
581 583 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
582 584 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
583 585
584 586 return normFactor
585 587
586 588 def getFlagCspc(self):
587 589
588 590 if self.data_cspc is None:
589 591 return True
590 592
591 593 return False
592 594
593 595 def getFlagDc(self):
594 596
595 597 if self.data_dc is None:
596 598 return True
597 599
598 600 return False
599 601
600 602 def getTimeInterval(self):
601 603
602 604 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
603 605 if self.nmodes:
604 606 return self.nmodes*timeInterval
605 607 else:
606 608 return timeInterval
607 609
608 610 def getPower(self):
609 611
610 612 factor = self.normFactor
611 613 z = self.data_spc / factor
612 614 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
613 615 avg = numpy.average(z, axis=1)
614 616
615 617 return 10 * numpy.log10(avg)
616 618
617 619 def getCoherence(self, pairsList=None, phase=False):
618 620
619 621 z = []
620 622 if pairsList is None:
621 623 pairsIndexList = self.pairsIndexList
622 624 else:
623 625 pairsIndexList = []
624 626 for pair in pairsList:
625 627 if pair not in self.pairsList:
626 628 raise ValueError("Pair %s is not in dataOut.pairsList" % (
627 629 pair))
628 630 pairsIndexList.append(self.pairsList.index(pair))
629 631 for i in range(len(pairsIndexList)):
630 632 pair = self.pairsList[pairsIndexList[i]]
631 633 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
632 634 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
633 635 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
634 636 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
635 637 if phase:
636 638 data = numpy.arctan2(avgcoherenceComplex.imag,
637 639 avgcoherenceComplex.real) * 180 / numpy.pi
638 640 else:
639 641 data = numpy.abs(avgcoherenceComplex)
640 642
641 643 z.append(data)
642 644
643 645 return numpy.array(z)
644 646
645 647 def setValue(self, value):
646 648
647 649 print("This property should not be initialized")
648 650
649 651 return
650 652
651 653 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
652 654 pairsIndexList = property(
653 655 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
654 656 normFactor = property(getNormFactor, setValue,
655 657 "I'm the 'getNormFactor' property.")
656 658 flag_cspc = property(getFlagCspc, setValue)
657 659 flag_dc = property(getFlagDc, setValue)
658 660 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
659 661 timeInterval = property(getTimeInterval, setValue,
660 662 "I'm the 'timeInterval' property")
661 663
662 664
663 665 class SpectraHeis(Spectra):
664 666
665 667 data_spc = None
666 668 data_cspc = None
667 669 data_dc = None
668 670 nFFTPoints = None
669 671 # nPairs = None
670 672 pairsList = None
671 673 nCohInt = None
672 674 nIncohInt = None
673 675
674 676 def __init__(self):
675 677
676 678 self.radarControllerHeaderObj = RadarControllerHeader()
677 679
678 680 self.systemHeaderObj = SystemHeader()
679 681
680 682 self.type = "SpectraHeis"
681 683
682 684 # self.dtype = None
683 685
684 686 # self.nChannels = 0
685 687
686 688 # self.nHeights = 0
687 689
688 690 self.nProfiles = None
689 691
690 692 self.heightList = None
691 693
692 694 self.channelList = None
693 695
694 696 # self.channelIndexList = None
695 697
696 698 self.flagNoData = True
697 699
698 700 self.flagDiscontinuousBlock = False
699 701
700 702 # self.nPairs = 0
701 703
702 704 self.utctime = None
703 705
704 706 self.blocksize = None
705 707
706 708 self.profileIndex = 0
707 709
708 710 self.nCohInt = 1
709 711
710 712 self.nIncohInt = 1
711 713
712 714 def getNormFactor(self):
713 715 pwcode = 1
714 716 if self.flagDecodeData:
715 717 pwcode = numpy.sum(self.code[0]**2)
716 718
717 719 normFactor = self.nIncohInt * self.nCohInt * pwcode
718 720
719 721 return normFactor
720 722
721 723 def getTimeInterval(self):
722 724
723 725 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
724 726
725 727 return timeInterval
726 728
727 729 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
728 730 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
729 731
730 732
731 733 class Fits(JROData):
732 734
733 735 heightList = None
734 736 channelList = None
735 737 flagNoData = True
736 738 flagDiscontinuousBlock = False
737 739 useLocalTime = False
738 740 utctime = None
739 741 timeZone = None
740 742 # ippSeconds = None
741 743 # timeInterval = None
742 744 nCohInt = None
743 745 nIncohInt = None
744 746 noise = None
745 747 windowOfFilter = 1
746 748 # Speed of ligth
747 749 C = 3e8
748 750 frequency = 49.92e6
749 751 realtime = False
750 752
751 753 def __init__(self):
752 754
753 755 self.type = "Fits"
754 756
755 757 self.nProfiles = None
756 758
757 759 self.heightList = None
758 760
759 761 self.channelList = None
760 762
761 763 # self.channelIndexList = None
762 764
763 765 self.flagNoData = True
764 766
765 767 self.utctime = None
766 768
767 769 self.nCohInt = 1
768 770
769 771 self.nIncohInt = 1
770 772
771 773 self.useLocalTime = True
772 774
773 775 self.profileIndex = 0
774 776
775 777 # self.utctime = None
776 778 # self.timeZone = None
777 779 # self.ltctime = None
778 780 # self.timeInterval = None
779 781 # self.header = None
780 782 # self.data_header = None
781 783 # self.data = None
782 784 # self.datatime = None
783 785 # self.flagNoData = False
784 786 # self.expName = ''
785 787 # self.nChannels = None
786 788 # self.nSamples = None
787 789 # self.dataBlocksPerFile = None
788 790 # self.comments = ''
789 791 #
790 792
791 793 def getltctime(self):
792 794
793 795 if self.useLocalTime:
794 796 return self.utctime - self.timeZone * 60
795 797
796 798 return self.utctime
797 799
798 800 def getDatatime(self):
799 801
800 802 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
801 803 return datatime
802 804
803 805 def getTimeRange(self):
804 806
805 807 datatime = []
806 808
807 809 datatime.append(self.ltctime)
808 810 datatime.append(self.ltctime + self.timeInterval)
809 811
810 812 datatime = numpy.array(datatime)
811 813
812 814 return datatime
813 815
814 816 def getHeiRange(self):
815 817
816 818 heis = self.heightList
817 819
818 820 return heis
819 821
820 822 def getNHeights(self):
821 823
822 824 return len(self.heightList)
823 825
824 826 def getNChannels(self):
825 827
826 828 return len(self.channelList)
827 829
828 830 def getChannelIndexList(self):
829 831
830 832 return list(range(self.nChannels))
831 833
832 834 def getNoise(self, type=1):
833 835
834 836 #noise = numpy.zeros(self.nChannels)
835 837
836 838 if type == 1:
837 839 noise = self.getNoisebyHildebrand()
838 840
839 841 if type == 2:
840 842 noise = self.getNoisebySort()
841 843
842 844 if type == 3:
843 845 noise = self.getNoisebyWindow()
844 846
845 847 return noise
846 848
847 849 def getTimeInterval(self):
848 850
849 851 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
850 852
851 853 return timeInterval
852 854
853 855 def get_ippSeconds(self):
854 856 '''
855 857 '''
856 858 return self.ipp_sec
857 859
858 860
859 861 datatime = property(getDatatime, "I'm the 'datatime' property")
860 862 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
861 863 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
862 864 channelIndexList = property(
863 865 getChannelIndexList, "I'm the 'channelIndexList' property.")
864 866 noise = property(getNoise, "I'm the 'nHeights' property.")
865 867
866 868 ltctime = property(getltctime, "I'm the 'ltctime' property")
867 869 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
868 870 ippSeconds = property(get_ippSeconds, '')
869 871
870 872 class Correlation(JROData):
871 873
872 874 noise = None
873 875 SNR = None
874 876 #--------------------------------------------------
875 877 mode = None
876 878 split = False
877 879 data_cf = None
878 880 lags = None
879 881 lagRange = None
880 882 pairsList = None
881 883 normFactor = None
882 884 #--------------------------------------------------
883 885 # calculateVelocity = None
884 886 nLags = None
885 887 nPairs = None
886 888 nAvg = None
887 889
888 890 def __init__(self):
889 891 '''
890 892 Constructor
891 893 '''
892 894 self.radarControllerHeaderObj = RadarControllerHeader()
893 895
894 896 self.systemHeaderObj = SystemHeader()
895 897
896 898 self.type = "Correlation"
897 899
898 900 self.data = None
899 901
900 902 self.dtype = None
901 903
902 904 self.nProfiles = None
903 905
904 906 self.heightList = None
905 907
906 908 self.channelList = None
907 909
908 910 self.flagNoData = True
909 911
910 912 self.flagDiscontinuousBlock = False
911 913
912 914 self.utctime = None
913 915
914 916 self.timeZone = None
915 917
916 918 self.dstFlag = None
917 919
918 920 self.errorCount = None
919 921
920 922 self.blocksize = None
921 923
922 924 self.flagDecodeData = False # asumo q la data no esta decodificada
923 925
924 926 self.flagDeflipData = False # asumo q la data no esta sin flip
925 927
926 928 self.pairsList = None
927 929
928 930 self.nPoints = None
929 931
930 932 def getPairsList(self):
931 933
932 934 return self.pairsList
933 935
934 936 def getNoise(self, mode=2):
935 937
936 938 indR = numpy.where(self.lagR == 0)[0][0]
937 939 indT = numpy.where(self.lagT == 0)[0][0]
938 940
939 941 jspectra0 = self.data_corr[:, :, indR, :]
940 942 jspectra = copy.copy(jspectra0)
941 943
942 944 num_chan = jspectra.shape[0]
943 945 num_hei = jspectra.shape[2]
944 946
945 947 freq_dc = jspectra.shape[1] / 2
946 948 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
947 949
948 950 if ind_vel[0] < 0:
949 951 ind_vel[list(range(0, 1))] = ind_vel[list(
950 952 range(0, 1))] + self.num_prof
951 953
952 954 if mode == 1:
953 955 jspectra[:, freq_dc, :] = (
954 956 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
955 957
956 958 if mode == 2:
957 959
958 960 vel = numpy.array([-2, -1, 1, 2])
959 961 xx = numpy.zeros([4, 4])
960 962
961 963 for fil in range(4):
962 964 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
963 965
964 966 xx_inv = numpy.linalg.inv(xx)
965 967 xx_aux = xx_inv[0, :]
966 968
967 969 for ich in range(num_chan):
968 970 yy = jspectra[ich, ind_vel, :]
969 971 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
970 972
971 973 junkid = jspectra[ich, freq_dc, :] <= 0
972 974 cjunkid = sum(junkid)
973 975
974 976 if cjunkid.any():
975 977 jspectra[ich, freq_dc, junkid.nonzero()] = (
976 978 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
977 979
978 980 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
979 981
980 982 return noise
981 983
982 984 def getTimeInterval(self):
983 985
984 986 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
985 987
986 988 return timeInterval
987 989
988 990 def splitFunctions(self):
989 991
990 992 pairsList = self.pairsList
991 993 ccf_pairs = []
992 994 acf_pairs = []
993 995 ccf_ind = []
994 996 acf_ind = []
995 997 for l in range(len(pairsList)):
996 998 chan0 = pairsList[l][0]
997 999 chan1 = pairsList[l][1]
998 1000
999 1001 # Obteniendo pares de Autocorrelacion
1000 1002 if chan0 == chan1:
1001 1003 acf_pairs.append(chan0)
1002 1004 acf_ind.append(l)
1003 1005 else:
1004 1006 ccf_pairs.append(pairsList[l])
1005 1007 ccf_ind.append(l)
1006 1008
1007 1009 data_acf = self.data_cf[acf_ind]
1008 1010 data_ccf = self.data_cf[ccf_ind]
1009 1011
1010 1012 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1011 1013
1012 1014 def getNormFactor(self):
1013 1015 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1014 1016 acf_pairs = numpy.array(acf_pairs)
1015 1017 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1016 1018
1017 1019 for p in range(self.nPairs):
1018 1020 pair = self.pairsList[p]
1019 1021
1020 1022 ch0 = pair[0]
1021 1023 ch1 = pair[1]
1022 1024
1023 1025 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1024 1026 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1025 1027 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1026 1028
1027 1029 return normFactor
1028 1030
1029 1031 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1030 1032 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1031 1033
1032 1034
1033 1035 class Parameters(Spectra):
1034 1036
1035 1037 experimentInfo = None # Information about the experiment
1036 1038 # Information from previous data
1037 1039 inputUnit = None # Type of data to be processed
1038 1040 operation = None # Type of operation to parametrize
1039 1041 # normFactor = None #Normalization Factor
1040 1042 groupList = None # List of Pairs, Groups, etc
1041 1043 # Parameters
1042 1044 data_param = None # Parameters obtained
1043 1045 data_pre = None # Data Pre Parametrization
1044 1046 data_SNR = None # Signal to Noise Ratio
1045 1047 # heightRange = None #Heights
1046 1048 abscissaList = None # Abscissa, can be velocities, lags or time
1047 1049 # noise = None #Noise Potency
1048 1050 utctimeInit = None # Initial UTC time
1049 1051 paramInterval = None # Time interval to calculate Parameters in seconds
1050 1052 useLocalTime = True
1051 1053 # Fitting
1052 1054 data_error = None # Error of the estimation
1053 1055 constants = None
1054 1056 library = None
1055 1057 # Output signal
1056 1058 outputInterval = None # Time interval to calculate output signal in seconds
1057 1059 data_output = None # Out signal
1058 1060 nAvg = None
1059 1061 noise_estimation = None
1060 1062 GauSPC = None # Fit gaussian SPC
1061 1063
1062 1064 def __init__(self):
1063 1065 '''
1064 1066 Constructor
1065 1067 '''
1066 1068 self.radarControllerHeaderObj = RadarControllerHeader()
1067 1069
1068 1070 self.systemHeaderObj = SystemHeader()
1069 1071
1070 1072 self.type = "Parameters"
1071 1073
1072 1074 def getTimeRange1(self, interval):
1073 1075
1074 1076 datatime = []
1075 1077
1076 1078 if self.useLocalTime:
1077 1079 time1 = self.utctimeInit - self.timeZone * 60
1078 1080 else:
1079 1081 time1 = self.utctimeInit
1080 1082
1081 1083 datatime.append(time1)
1082 1084 datatime.append(time1 + interval)
1083 1085 datatime = numpy.array(datatime)
1084 1086
1085 1087 return datatime
1086 1088
1087 1089 def getTimeInterval(self):
1088 1090
1089 1091 if hasattr(self, 'timeInterval1'):
1090 1092 return self.timeInterval1
1091 1093 else:
1092 1094 return self.paramInterval
1093 1095
1094 1096 def setValue(self, value):
1095 1097
1096 1098 print("This property should not be initialized")
1097 1099
1098 1100 return
1099 1101
1100 1102 def getNoise(self):
1101 1103
1102 1104 return self.spc_noise
1103 1105
1104 1106 timeInterval = property(getTimeInterval)
1105 1107 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1106 1108
1107 1109
1108 1110 class PlotterData(object):
1109 1111 '''
1110 1112 Object to hold data to be plotted
1111 1113 '''
1112 1114
1113 1115 MAXNUMX = 200
1114 1116 MAXNUMY = 200
1115 1117
1116 1118 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1117 1119
1118 1120 self.key = code
1119 1121 self.throttle = throttle_value
1120 1122 self.exp_code = exp_code
1121 1123 self.buffering = buffering
1122 1124 self.ready = False
1123 1125 self.flagNoData = False
1124 1126 self.localtime = False
1125 1127 self.data = {}
1126 1128 self.meta = {}
1127 1129 self.__heights = []
1128 1130
1129 1131 if 'snr' in code:
1130 1132 self.plottypes = ['snr']
1131 1133 elif code == 'spc':
1132 1134 self.plottypes = ['spc', 'noise', 'rti']
1133 1135 elif code == 'cspc':
1134 1136 self.plottypes = ['cspc', 'spc', 'noise', 'rti']
1135 1137 elif code == 'rti':
1136 1138 self.plottypes = ['noise', 'rti']
1137 1139 else:
1138 1140 self.plottypes = [code]
1139 1141
1140 1142 if 'snr' not in self.plottypes and snr:
1141 1143 self.plottypes.append('snr')
1142 1144
1143 1145 for plot in self.plottypes:
1144 1146 self.data[plot] = {}
1145 1147
1146 1148
1147 1149 def __str__(self):
1148 1150 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1149 1151 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1150 1152
1151 1153 def __len__(self):
1152 1154 return len(self.data[self.key])
1153 1155
1154 1156 def __getitem__(self, key):
1155 1157
1156 1158 if key not in self.data:
1157 1159 raise KeyError(log.error('Missing key: {}'.format(key)))
1158 1160 if 'spc' in key or not self.buffering:
1159 1161 ret = self.data[key][self.tm]
1160 1162 elif 'scope' in key:
1161 1163 ret = numpy.array(self.data[key][float(self.tm)])
1162 1164 else:
1163 1165 ret = numpy.array([self.data[key][x] for x in self.times])
1164 1166 if ret.ndim > 1:
1165 1167 ret = numpy.swapaxes(ret, 0, 1)
1166 1168 return ret
1167 1169
1168 1170 def __contains__(self, key):
1169 1171 return key in self.data
1170 1172
1171 1173 def setup(self):
1172 1174 '''
1173 1175 Configure object
1174 1176 '''
1175 1177 self.type = ''
1176 1178 self.ready = False
1177 1179 del self.data
1178 1180 self.data = {}
1179 1181 self.__heights = []
1180 1182 self.__all_heights = set()
1181 1183 for plot in self.plottypes:
1182 1184 if 'snr' in plot:
1183 1185 plot = 'snr'
1184 1186 elif 'spc_moments' == plot:
1185 1187 plot = 'moments'
1186 1188 self.data[plot] = {}
1187 1189
1188 1190 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1189 1191 self.data['noise'] = {}
1190 1192 self.data['rti'] = {}
1191 1193 if 'noise' not in self.plottypes:
1192 1194 self.plottypes.append('noise')
1193 1195 if 'rti' not in self.plottypes:
1194 1196 self.plottypes.append('rti')
1195 1197
1196 1198 def shape(self, key):
1197 1199 '''
1198 1200 Get the shape of the one-element data for the given key
1199 1201 '''
1200 1202
1201 1203 if len(self.data[key]):
1202 1204 if 'spc' in key or not self.buffering:
1203 1205 return self.data[key].shape
1204 1206 return self.data[key][self.times[0]].shape
1205 1207 return (0,)
1206 1208
1207 1209 def update(self, dataOut, tm):
1208 1210 '''
1209 1211 Update data object with new dataOut
1210 1212 '''
1211
1213
1212 1214 self.profileIndex = dataOut.profileIndex
1213 1215 self.tm = tm
1214 1216 self.type = dataOut.type
1215 1217 self.parameters = getattr(dataOut, 'parameters', [])
1216 1218
1217 1219 if hasattr(dataOut, 'meta'):
1218 1220 self.meta.update(dataOut.meta)
1219 1221
1220 1222 if hasattr(dataOut, 'pairsList'):
1221 1223 self.pairs = dataOut.pairsList
1222 1224
1223 1225 self.interval = dataOut.getTimeInterval()
1224 1226 self.localtime = dataOut.useLocalTime
1225 1227 if True in ['spc' in ptype for ptype in self.plottypes]:
1226 1228 self.xrange = (dataOut.getFreqRange(1)/1000.,
1227 1229 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1228 1230 self.__heights.append(dataOut.heightList)
1229 1231 self.__all_heights.update(dataOut.heightList)
1230 1232
1231 1233 for plot in self.plottypes:
1232 1234 if plot in ('spc', 'spc_moments', 'spc_cut'):
1233 1235 z = dataOut.data_spc/dataOut.normFactor
1234 1236 buffer = 10*numpy.log10(z)
1235 1237 if plot == 'cspc':
1236 1238 buffer = (dataOut.data_spc, dataOut.data_cspc)
1237 1239 if plot == 'noise':
1238 1240 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1239 1241 if plot in ('rti', 'spcprofile'):
1240 1242 buffer = dataOut.getPower()
1241 1243 if plot == 'snr_db':
1242 1244 buffer = dataOut.data_SNR
1243 1245 if plot == 'snr':
1244 1246 buffer = 10*numpy.log10(dataOut.data_SNR)
1245 1247 if plot == 'dop':
1246 1248 buffer = dataOut.data_DOP
1247 1249 if plot == 'pow':
1248 1250 buffer = 10*numpy.log10(dataOut.data_POW)
1249 1251 if plot == 'width':
1250 1252 buffer = dataOut.data_WIDTH
1251 1253 if plot == 'coh':
1252 1254 buffer = dataOut.getCoherence()
1253 1255 if plot == 'phase':
1254 1256 buffer = dataOut.getCoherence(phase=True)
1255 1257 if plot == 'output':
1256 1258 buffer = dataOut.data_output
1257 1259 if plot == 'param':
1258 1260 buffer = dataOut.data_param
1259 1261 if plot == 'scope':
1260 1262 buffer = dataOut.data
1261 1263 self.flagDataAsBlock = dataOut.flagDataAsBlock
1262 1264 self.nProfiles = dataOut.nProfiles
1263 1265 if plot == 'pp_power':
1264 buffer = dataOut.data_intensity
1266 buffer = dataOut.dataPP_POWER
1267 self.flagDataAsBlock = dataOut.flagDataAsBlock
1268 self.nProfiles = dataOut.nProfiles
1269 if plot == 'pp_signal':
1270 buffer = dataOut.dataPP_POW
1265 1271 self.flagDataAsBlock = dataOut.flagDataAsBlock
1266 1272 self.nProfiles = dataOut.nProfiles
1267 1273 if plot == 'pp_velocity':
1268 buffer = dataOut.data_velocity
1274 buffer = dataOut.dataPP_DOP
1269 1275 self.flagDataAsBlock = dataOut.flagDataAsBlock
1270 1276 self.nProfiles = dataOut.nProfiles
1271 1277 if plot == 'pp_specwidth':
1272 buffer = dataOut.data_specwidth
1278 buffer = dataOut.dataPP_WIDTH
1273 1279 self.flagDataAsBlock = dataOut.flagDataAsBlock
1274 1280 self.nProfiles = dataOut.nProfiles
1275 1281
1276 1282 if plot == 'spc':
1277 1283 self.data['spc'][tm] = buffer
1278 1284 elif plot == 'cspc':
1279 1285 self.data['cspc'][tm] = buffer
1280 1286 elif plot == 'spc_moments':
1281 1287 self.data['spc'][tm] = buffer
1282 1288 self.data['moments'][tm] = dataOut.moments
1283 1289 else:
1284 1290 if self.buffering:
1285 1291 self.data[plot][tm] = buffer
1286 1292 else:
1287 1293 self.data[plot][tm] = buffer
1288 1294
1289 1295 if dataOut.channelList is None:
1290 1296 self.channels = range(buffer.shape[0])
1291 1297 else:
1292 1298 self.channels = dataOut.channelList
1293 1299
1294 1300 if buffer is None:
1295 1301 self.flagNoData = True
1296 1302 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1297 1303
1298 1304 def normalize_heights(self):
1299 1305 '''
1300 1306 Ensure same-dimension of the data for different heighList
1301 1307 '''
1302 1308
1303 1309 H = numpy.array(list(self.__all_heights))
1304 1310 H.sort()
1305 1311 for key in self.data:
1306 1312 shape = self.shape(key)[:-1] + H.shape
1307 1313 for tm, obj in list(self.data[key].items()):
1308 1314 h = self.__heights[self.times.tolist().index(tm)]
1309 1315 if H.size == h.size:
1310 1316 continue
1311 1317 index = numpy.where(numpy.in1d(H, h))[0]
1312 1318 dummy = numpy.zeros(shape) + numpy.nan
1313 1319 if len(shape) == 2:
1314 1320 dummy[:, index] = obj
1315 1321 else:
1316 1322 dummy[index] = obj
1317 1323 self.data[key][tm] = dummy
1318 1324
1319 1325 self.__heights = [H for tm in self.times]
1320 1326
1321 1327 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1322 1328 '''
1323 1329 Convert data to json
1324 1330 '''
1325 1331
1326 1332 dy = int(self.heights.size/self.MAXNUMY) + 1
1327 1333 if self.key in ('spc', 'cspc'):
1328 1334 dx = int(self.data[self.key][tm].shape[1]/self.MAXNUMX) + 1
1329 1335 data = self.roundFloats(
1330 1336 self.data[self.key][tm][::, ::dx, ::dy].tolist())
1331 1337 else:
1332 1338 if self.key is 'noise':
1333 1339 data = [[x] for x in self.roundFloats(self.data[self.key][tm].tolist())]
1334 1340 else:
1335 1341 data = self.roundFloats(self.data[self.key][tm][::, ::dy].tolist())
1336 1342
1337 1343 meta = {}
1338 1344 ret = {
1339 1345 'plot': plot_name,
1340 1346 'code': self.exp_code,
1341 1347 'time': float(tm),
1342 1348 'data': data,
1343 1349 }
1344 1350 meta['type'] = plot_type
1345 1351 meta['interval'] = float(self.interval)
1346 1352 meta['localtime'] = self.localtime
1347 1353 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1348 1354 if 'spc' in self.data or 'cspc' in self.data:
1349 1355 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1350 1356 else:
1351 1357 meta['xrange'] = []
1352 1358
1353 1359 meta.update(self.meta)
1354 1360 ret['metadata'] = meta
1355 1361 return json.dumps(ret)
1356 1362
1357 1363 @property
1358 1364 def times(self):
1359 1365 '''
1360 1366 Return the list of times of the current data
1361 1367 '''
1362 1368
1363 1369 ret = numpy.array([*self.data[self.key]])
1364 1370 if self:
1365 1371 ret.sort()
1366 1372 return ret
1367 1373
1368 1374 @property
1369 1375 def min_time(self):
1370 1376 '''
1371 1377 Return the minimun time value
1372 1378 '''
1373 1379
1374 1380 return self.times[0]
1375 1381
1376 1382 @property
1377 1383 def max_time(self):
1378 1384 '''
1379 1385 Return the maximun time value
1380 1386 '''
1381 1387
1382 1388 return self.times[-1]
1383 1389
1384 1390 @property
1385 1391 def heights(self):
1386 1392 '''
1387 1393 Return the list of heights of the current data
1388 1394 '''
1389 1395
1390 1396 return numpy.array(self.__heights[-1])
1391 1397
1392 1398 @staticmethod
1393 1399 def roundFloats(obj):
1394 1400 if isinstance(obj, list):
1395 1401 return list(map(PlotterData.roundFloats, obj))
1396 1402 elif isinstance(obj, float):
1397 1403 return round(obj, 2)
@@ -1,276 +1,302
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9
10 10 from schainpy.model.graphics.jroplot_base import Plot, plt
11 11
12 12
13 13 class ScopePlot(Plot):
14 14
15 15 '''
16 16 Plot for Scope
17 17 '''
18 18
19 19 CODE = 'scope'
20 20 plot_name = 'Scope'
21 21 plot_type = 'scatter'
22 22
23 23 def setup(self):
24 24
25 25 self.xaxis = 'Range (Km)'
26 26 self.ncols = 1
27 27 self.nrows = 1
28 28 self.nplots = 1
29 29 self.ylabel = 'Intensity [dB]'
30 30 self.titles = ['Scope']
31 31 self.colorbar = False
32 32 self.width = 6
33 33 self.height = 4
34 34
35 35 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
36 36
37 37 yreal = y[channelIndexList,:].real
38 38 yimag = y[channelIndexList,:].imag
39 39 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
40 40 self.xlabel = "Range (Km)"
41 41 self.ylabel = "Intensity - IQ"
42 42
43 43 self.y = yreal
44 44 self.x = x
45 45 self.xmin = min(x)
46 46 self.xmax = max(x)
47 47
48 48
49 49 self.titles[0] = title
50 50
51 51 for i,ax in enumerate(self.axes):
52 52 title = "Channel %d" %(i)
53 53 if ax.firsttime:
54 54 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
55 55 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
56 56 else:
57 57 ax.plt_r.set_data(x, yreal[i,:])
58 58 ax.plt_i.set_data(x, yimag[i,:])
59 59
60 60 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
61 61 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
62 62 yreal = y.real
63 63 yreal = 10*numpy.log10(yreal)
64 64 self.y = yreal
65 65 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
66 66 self.xlabel = "Range (Km)"
67 67 self.ylabel = "Intensity"
68 68 self.xmin = min(x)
69 69 self.xmax = max(x)
70 70
71 71
72 72 self.titles[0] = title
73 73
74 74 for i,ax in enumerate(self.axes):
75 75 title = "Channel %d" %(i)
76 76
77 77 ychannel = yreal[i,:]
78 78
79 79 if ax.firsttime:
80 80 ax.plt_r = ax.plot(x, ychannel)[0]
81 81 else:
82 82 #pass
83 83 ax.plt_r.set_data(x, ychannel)
84 84
85 85 def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle):
86 86
87 87
88 88 y = y[channelIndexList,:]
89 89 yreal = y.real
90 90 yreal = 10*numpy.log10(yreal)
91 91 self.y = yreal
92 92 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
93 93 self.xlabel = "Range (Km)"
94 94 self.ylabel = "Intensity"
95 95 self.xmin = min(x)
96 96 self.xmax = max(x)
97 97
98 98 self.titles[0] =title
99 99 for i,ax in enumerate(self.axes):
100 100 title = "Channel %d" %(i)
101 101
102 102 ychannel = yreal[i,:]
103 103
104 104 if ax.firsttime:
105 105 ax.plt_r = ax.plot(x, ychannel)[0]
106 106 else:
107 107 #pass
108 108 ax.plt_r.set_data(x, ychannel)
109 109
110 110 def plot_weathervelocity(self, x, y, channelIndexList, thisDatetime, wintitle):
111 111
112 112 x = x[channelIndexList,:]
113 113 yreal = y
114 114 self.y = yreal
115 115 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
116 116 self.xlabel = "Velocity (m/s)"
117 117 self.ylabel = "Range (Km)"
118 118 self.xmin = numpy.min(x)
119 119 self.xmax = numpy.max(x)
120 120 self.titles[0] =title
121 121 for i,ax in enumerate(self.axes):
122 122 title = "Channel %d" %(i)
123 123 xchannel = x[i,:]
124 124 if ax.firsttime:
125 125 ax.plt_r = ax.plot(xchannel, yreal)[0]
126 126 else:
127 127 #pass
128 128 ax.plt_r.set_data(xchannel, yreal)
129 129
130 130 def plot_weatherspecwidth(self, x, y, channelIndexList, thisDatetime, wintitle):
131 131
132 132 x = x[channelIndexList,:]
133 133 yreal = y
134 134 self.y = yreal
135 135 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
136 136 self.xlabel = "width "
137 137 self.ylabel = "Range (Km)"
138 138 self.xmin = numpy.min(x)
139 139 self.xmax = numpy.max(x)
140 140 self.titles[0] =title
141 141 for i,ax in enumerate(self.axes):
142 142 title = "Channel %d" %(i)
143 143 xchannel = x[i,:]
144 144 if ax.firsttime:
145 145 ax.plt_r = ax.plot(xchannel, yreal)[0]
146 146 else:
147 147 #pass
148 148 ax.plt_r.set_data(xchannel, yreal)
149 149
150 150 def plot(self):
151 151 if self.channels:
152 152 channels = self.channels
153 153 else:
154 154 channels = self.data.channels
155 155
156 156 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
157 157 if self.CODE == "pp_power":
158 158 scope = self.data['pp_power']
159 elif self.CODE == "pp_signal":
160 scope = self.data["pp_signal"]
159 161 elif self.CODE == "pp_velocity":
160 162 scope = self.data["pp_velocity"]
161 163 elif self.CODE == "pp_specwidth":
162 164 scope = self.data["pp_specwidth"]
163 165 else:
164 166 scope =self.data["scope"]
165 167
166 168 if self.data.flagDataAsBlock:
167 169
168 170 for i in range(self.data.nProfiles):
169 171
170 172 wintitle1 = " [Profile = %d] " %i
171 173 if self.CODE =="scope":
172 174 if self.type == "power":
173 175 self.plot_power(self.data.heights,
174 176 scope[:,i,:],
175 177 channels,
176 178 thisDatetime,
177 179 wintitle1
178 180 )
179 181
180 182 if self.type == "iq":
181 183 self.plot_iq(self.data.heights,
182 184 scope[:,i,:],
183 185 channels,
184 186 thisDatetime,
185 187 wintitle1
186 188 )
187 189 if self.CODE=="pp_power":
188 190 self.plot_weatherpower(self.data.heights,
189 191 scope[:,i,:],
190 192 channels,
191 193 thisDatetime,
192 194 wintitle
193 195 )
196 if self.CODE=="pp_signal":
197 self.plot_weatherpower(self.data.heights,
198 scope[:,i,:],
199 channels,
200 thisDatetime,
201 wintitle
202 )
194 203 if self.CODE=="pp_velocity":
195 204 self.plot_weathervelocity(scope[:,i,:],
196 205 self.data.heights,
197 206 channels,
198 207 thisDatetime,
199 208 wintitle
200 209 )
201 210 if self.CODE=="pp_spcwidth":
202 211 self.plot_weatherspecwidth(scope[:,i,:],
203 212 self.data.heights,
204 213 channels,
205 214 thisDatetime,
206 215 wintitle
207 216 )
208 217 else:
209 218 wintitle = " [Profile = %d] " %self.data.profileIndex
210 219 if self.CODE== "scope":
211 220 if self.type == "power":
212 221 self.plot_power(self.data.heights,
213 222 scope,
214 223 channels,
215 224 thisDatetime,
216 225 wintitle
217 226 )
218 227
219 228 if self.type == "iq":
220 229 self.plot_iq(self.data.heights,
221 230 scope,
222 231 channels,
223 232 thisDatetime,
224 233 wintitle
225 234 )
226 235 if self.CODE=="pp_power":
227 236 self.plot_weatherpower(self.data.heights,
228 237 scope,
229 238 channels,
230 239 thisDatetime,
231 240 wintitle
232 241 )
242 if self.CODE=="pp_signal":
243 self.plot_weatherpower(self.data.heights,
244 scope,
245 channels,
246 thisDatetime,
247 wintitle
248 )
233 249 if self.CODE=="pp_velocity":
234 250 self.plot_weathervelocity(scope,
235 251 self.data.heights,
236 252 channels,
237 253 thisDatetime,
238 254 wintitle
239 255 )
240 256 if self.CODE=="pp_specwidth":
241 257 self.plot_weatherspecwidth(scope,
242 258 self.data.heights,
243 259 channels,
244 260 thisDatetime,
245 261 wintitle
246 262 )
247 263
248 264
249 265
250 266 class PulsepairPowerPlot(ScopePlot):
251 267 '''
252 Plot for
268 Plot for P= S+N
253 269 '''
254 270
255 271 CODE = 'pp_power'
256 272 plot_name = 'PulsepairPower'
257 273 plot_type = 'scatter'
258 274 buffering = False
259 275
260 276 class PulsepairVelocityPlot(ScopePlot):
261 277 '''
262 Plot for
278 Plot for VELOCITY
263 279 '''
264 280 CODE = 'pp_velocity'
265 281 plot_name = 'PulsepairVelocity'
266 282 plot_type = 'scatter'
267 283 buffering = False
268 284
269 285 class PulsepairSpecwidthPlot(ScopePlot):
270 286 '''
271 Plot for
287 Plot for WIDTH
272 288 '''
273 289 CODE = 'pp_specwidth'
274 290 plot_name = 'PulsepairSpecwidth'
275 291 plot_type = 'scatter'
276 292 buffering = False
293
294 class PulsepairSignalPlot(ScopePlot):
295 '''
296 Plot for S
297 '''
298
299 CODE = 'pp_signal'
300 plot_name = 'PulsepairSignal'
301 plot_type = 'scatter'
302 buffering = False
@@ -1,512 +1,519
1 1 import numpy,math,random,time
2 2 #---------------1 Heredamos JRODatareader
3 3 from schainpy.model.io.jroIO_base import *
4 4 #---------------2 Heredamos las propiedades de ProcessingUnit
5 5 from schainpy.model.proc.jroproc_base import ProcessingUnit,Operation,MPDecorator
6 6 #---------------3 Importaremos las clases BascicHeader, SystemHeader, RadarControlHeader, ProcessingHeader
7 7 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader,SystemHeader,RadarControllerHeader, ProcessingHeader
8 8 #---------------4 Importaremos el objeto Voltge
9 9 from schainpy.model.data.jrodata import Voltage
10 10
11 11 class SimulatorReader(JRODataReader, ProcessingUnit):
12 12 incIntFactor = 1
13 13 nFFTPoints = 0
14 14 FixPP_IncInt = 1
15 15 FixRCP_IPP = 1000
16 16 FixPP_CohInt = 1
17 17 Tau_0 = 250
18 18 AcqH0_0 = 70
19 19 H0 = AcqH0_0
20 20 AcqDH_0 = 1.25
21 21 DH0 = AcqDH_0
22 22 Bauds = 32
23 23 BaudWidth = None
24 24 FixRCP_TXA = 40
25 25 FixRCP_TXB = 70
26 26 fAngle = 2.0*math.pi*(1/16)
27 27 DC_level = 500
28 28 stdev = 8
29 29 Num_Codes = 2
30 30 #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1])
31 31 #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0])
32 32 #Dyn_snCode = numpy.array([Num_Codes,Bauds])
33 33 Dyn_snCode = None
34 34 Samples = 200
35 35 channels = 2
36 36 pulses = None
37 37 Reference = None
38 38 pulse_size = None
39 39 prof_gen = None
40 40 Fdoppler = 100
41 41 Hdoppler = 36
42 42 Adoppler = 300
43 43 frequency = 9345
44 44 nTotalReadFiles = 1000
45 45
46 46 def __init__(self):
47 47 """
48 48 Inicializador de la clases SimulatorReader para
49 49 generar datos de voltage simulados.
50 50 Input:
51 51 dataOut: Objeto de la clase Voltage.
52 52 Este Objeto sera utilizado apra almacenar
53 53 un perfil de datos cada vez qe se haga
54 54 un requerimiento (getData)
55 55 """
56 56 ProcessingUnit.__init__(self)
57 57 print(" [ START ] init - Metodo Simulator Reader")
58 58
59 59 self.isConfig = False
60 60 self.basicHeaderObj = BasicHeader(LOCALTIME)
61 61 self.systemHeaderObj = SystemHeader()
62 62 self.radarControllerHeaderObj = RadarControllerHeader()
63 63 self.processingHeaderObj = ProcessingHeader()
64 64 self.profileIndex = 2**32-1
65 65 self.dataOut = Voltage()
66 66 #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1])
67 67 code0 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,1,1,1,-1,1,1,-1,1,-1,-1,-1,1,1,1,-1,1])
68 68 #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0])
69 69 code1 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,-1,-1,1,-1,-1,1,-1,1,1,1,-1,-1,-1,1,-1])
70 70 #self.Dyn_snCode = numpy.array([code0,code1])
71 71 self.Dyn_snCode = None
72 72
73 73 def set_kwargs(self, **kwargs):
74 74 for key, value in kwargs.items():
75 75 setattr(self, key, value)
76 76
77 77 def __hasNotDataInBuffer(self):
78 78
79 79 if self.profileIndex >= self.processingHeaderObj.profilesPerBlock* self.nTxs:
80 80 if self.nReadBlocks>0:
81 81 tmp = self.dataOut.utctime
82 82 tmp_utc = int(self.dataOut.utctime)
83 83 tmp_milisecond = int((tmp-tmp_utc)*1000)
84 84 self.basicHeaderObj.utc = tmp_utc
85 85 self.basicHeaderObj.miliSecond= tmp_milisecond
86 86 return 1
87 87 return 0
88 88
89 89 def setNextFile(self):
90 90 """Set the next file to be readed open it and parse de file header"""
91 91
92 92 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
93 93 self.nReadFiles=self.nReadFiles+1
94 94 if self.nReadFiles > self.nTotalReadFiles:
95 95 self.flagNoMoreFiles=1
96 96 raise schainpy.admin.SchainWarning('No more files to read')
97 97
98 98 print('------------------- [Opening file] ------------------------------',self.nReadFiles)
99 99 self.nReadBlocks = 0
100 100 #if self.nReadBlocks==0:
101 101 # self.readFirstHeader()
102 102
103 103 def __setNewBlock(self):
104 104 self.setNextFile()
105 105 if self.flagIsNewFile:
106 106 return 1
107 107
108 108 def readNextBlock(self):
109 109 while True:
110 110 self.__setNewBlock()
111 111 if not(self.readBlock()):
112 112 return 0
113 113 self.getBasicHeader()
114 114 break
115 115 if self.verbose:
116 116 print("[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
117 117 self.processingHeaderObj.dataBlocksPerFile,
118 118 self.dataOut.datatime.ctime()) )
119 119 return 1
120 120
121 121 def getFirstHeader(self):
122 122 self.getBasicHeader()
123 123 self.dataOut.processingHeaderObj = self.processingHeaderObj.copy()
124 124 self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
125 125 self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
126 126 self.dataOut.dtype = self.dtype
127 127
128 128 self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
129 129 self.dataOut.heightList = numpy.arange(self.processingHeaderObj.nHeights) * self.processingHeaderObj.deltaHeight + self.processingHeaderObj.firstHeight
130 130 self.dataOut.channelList = list(range(self.systemHeaderObj.nChannels))
131 131 self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
132 132 # asumo q la data no esta decodificada
133 133 self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode
134 134 # asumo q la data no esta sin flip
135 135 self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip
136 136 self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft
137 137 self.dataOut.frequency = self.frequency
138 138
139 139 def getBasicHeader(self):
140 140 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
141 141 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
142 142
143 143 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
144 144 self.dataOut.timeZone = self.basicHeaderObj.timeZone
145 145 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
146 146 self.dataOut.errorCount = self.basicHeaderObj.errorCount
147 147 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
148 148 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
149 149
150 150 def readFirstHeader(self):
151 151
152 152 datatype = int(numpy.log2((self.processingHeaderObj.processFlags &
153 153 PROCFLAG.DATATYPE_MASK)) - numpy.log2(PROCFLAG.DATATYPE_CHAR))
154 154 if datatype == 0:
155 155 datatype_str = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
156 156 elif datatype == 1:
157 157 datatype_str = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
158 158 elif datatype == 2:
159 159 datatype_str = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
160 160 elif datatype == 3:
161 161 datatype_str = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
162 162 elif datatype == 4:
163 163 datatype_str = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
164 164 elif datatype == 5:
165 165 datatype_str = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
166 166 else:
167 167 raise ValueError('Data type was not defined')
168 168
169 169 self.dtype = datatype_str
170 170
171 171
172 172 def set_RCH(self, expType=2, nTx=1,ipp=None, txA=0, txB=0,
173 173 nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None,
174 174 numTaus=0, line6Function=0, line5Function=0, fClock=None,
175 175 prePulseBefore=0, prePulseAfter=0,
176 176 codeType=0, nCode=0, nBaud=0, code=None,
177 177 flip1=0, flip2=0,Taus=0):
178 178 self.radarControllerHeaderObj.expType = expType
179 179 self.radarControllerHeaderObj.nTx = nTx
180 180 self.radarControllerHeaderObj.ipp = float(ipp)
181 181 self.radarControllerHeaderObj.txA = float(txA)
182 182 self.radarControllerHeaderObj.txB = float(txB)
183 183 self.radarControllerHeaderObj.rangeIpp = b'A\n'#ipp
184 184 self.radarControllerHeaderObj.rangeTxA = b''
185 185 self.radarControllerHeaderObj.rangeTxB = b''
186 186
187 187 self.radarControllerHeaderObj.nHeights = int(nHeights)
188 188 self.radarControllerHeaderObj.firstHeight = numpy.array([firstHeight])
189 189 self.radarControllerHeaderObj.deltaHeight = numpy.array([deltaHeight])
190 190 self.radarControllerHeaderObj.samplesWin = numpy.array([nHeights])
191 191
192 192
193 193 self.radarControllerHeaderObj.nWindows = nWindows
194 194 self.radarControllerHeaderObj.numTaus = numTaus
195 195 self.radarControllerHeaderObj.codeType = codeType
196 196 self.radarControllerHeaderObj.line6Function = line6Function
197 197 self.radarControllerHeaderObj.line5Function = line5Function
198 198 #self.radarControllerHeaderObj.fClock = fClock
199 199 self.radarControllerHeaderObj.prePulseBefore= prePulseBefore
200 200 self.radarControllerHeaderObj.prePulseAfter = prePulseAfter
201 201
202 202 self.radarControllerHeaderObj.flip1 = flip1
203 203 self.radarControllerHeaderObj.flip2 = flip2
204 204
205 205 self.radarControllerHeaderObj.code_size = 0
206 206 if self.radarControllerHeaderObj.codeType != 0:
207 207 self.radarControllerHeaderObj.nCode = nCode
208 208 self.radarControllerHeaderObj.nBaud = nBaud
209 209 self.radarControllerHeaderObj.code = code
210 210 self.radarControllerHeaderObj.code_size = int(numpy.ceil(nBaud / 32.)) * nCode * 4
211 211
212 212 if fClock is None and deltaHeight is not None:
213 213 self.fClock = 0.15 / (deltaHeight * 1e-6)
214 214 self.radarControllerHeaderObj.fClock = self.fClock
215 215 if numTaus==0:
216 216 self.radarControllerHeaderObj.Taus = numpy.array(0,'<f4')
217 217 else:
218 218 self.radarControllerHeaderObj.Taus = numpy.array(Taus,'<f4')
219 219
220 220 def set_PH(self, dtype=0, blockSize=0, profilesPerBlock=0,
221 221 dataBlocksPerFile=0, nWindows=0, processFlags=0, nCohInt=0,
222 222 nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0,
223 223 deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0,
224 224 code=0, nBaud=None, shif_fft=False, flag_dc=False,
225 225 flag_cspc=False, flag_decode=False, flag_deflip=False):
226 226
227 227 self.processingHeaderObj.dtype = dtype
228 228 self.processingHeaderObj.profilesPerBlock = profilesPerBlock
229 229 self.processingHeaderObj.dataBlocksPerFile = dataBlocksPerFile
230 230 self.processingHeaderObj.nWindows = nWindows
231 231 self.processingHeaderObj.processFlags = processFlags
232 232 self.processingHeaderObj.nCohInt = nCohInt
233 233 self.processingHeaderObj.nIncohInt = nIncohInt
234 234 self.processingHeaderObj.totalSpectra = totalSpectra
235 235
236 236 self.processingHeaderObj.nHeights = int(nHeights)
237 237 self.processingHeaderObj.firstHeight = firstHeight#numpy.array([firstHeight])#firstHeight
238 238 self.processingHeaderObj.deltaHeight = deltaHeight#numpy.array([deltaHeight])#deltaHeight
239 239 self.processingHeaderObj.samplesWin = nHeights#numpy.array([nHeights])#nHeights
240 240
241 241 def set_BH(self, utc = 0, miliSecond = 0, timeZone = 0):
242 242 self.basicHeaderObj.utc = utc
243 243 self.basicHeaderObj.miliSecond = miliSecond
244 244 self.basicHeaderObj.timeZone = timeZone
245 245
246 246 def set_SH(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWidth=32):
247 247 #self.systemHeaderObj.size = size
248 248 self.systemHeaderObj.nSamples = nSamples
249 249 self.systemHeaderObj.nProfiles = nProfiles
250 250 self.systemHeaderObj.nChannels = nChannels
251 251 self.systemHeaderObj.adcResolution = adcResolution
252 252 self.systemHeaderObj.pciDioBusWidth = pciDioBusWidth
253 253
254 254 def init_acquisition(self):
255 255
256 256 if self.nFFTPoints != 0:
257 257 self.incIntFactor = m_nProfilesperBlock/self.nFFTPoints
258 258 if (self.FixPP_IncInt > self.incIntFactor):
259 259 self.incIntFactor = self.FixPP_IncInt/ self.incIntFactor
260 260 elif(self.FixPP_IncInt< self.incIntFactor):
261 261 print("False alert...")
262 262
263 263 ProfilesperBlock = self.processingHeaderObj.profilesPerBlock
264 264
265 265 self.timeperblock =int(((self.FixRCP_IPP
266 266 *ProfilesperBlock
267 267 *self.FixPP_CohInt
268 268 *self.incIntFactor)
269 269 /150.0)
270 270 *0.9
271 271 +0.5)
272 272 # para cada canal
273 273 self.profiles = ProfilesperBlock*self.FixPP_CohInt
274 274 self.profiles = ProfilesperBlock
275 275 self.Reference = int((self.Tau_0-self.AcqH0_0)/(self.AcqDH_0)+0.5)
276 276 self.BaudWidth = int((self.FixRCP_TXA/self.AcqDH_0)/self.Bauds + 0.5 )
277 277
278 278 if (self.BaudWidth==0):
279 279 self.BaudWidth=1
280 280
281 281 def init_pulse(self,Num_Codes=Num_Codes,Bauds=Bauds,BaudWidth=BaudWidth,Dyn_snCode=Dyn_snCode):
282 282
283 283 Num_Codes = Num_Codes
284 284 Bauds = Bauds
285 285 BaudWidth = BaudWidth
286 286 Dyn_snCode = Dyn_snCode
287 287
288 288 if Dyn_snCode:
289 289 print("EXISTE")
290 290 else:
291 291 print("No existe")
292 292
293 293 if Dyn_snCode: # if Bauds:
294 294 pulses = list(range(0,Num_Codes))
295 295 num_codes = Num_Codes
296 296 for i in range(num_codes):
297 297 pulse_size = Bauds*BaudWidth
298 298 pulses[i] = numpy.zeros(pulse_size)
299 299 for j in range(Bauds):
300 300 for k in range(BaudWidth):
301 301 pulses[i][j*BaudWidth+k] = int(Dyn_snCode[i][j]*600)
302 302 else:
303 303 print("sin code")
304 304 pulses = list(range(1))
305 305 if self.AcqDH_0>0.149:
306 306 pulse_size = int(self.FixRCP_TXB/0.15+0.5)
307 307 else:
308 308 pulse_size = int((self.FixRCP_TXB/self.AcqDH_0)+0.5) #0.0375
309 309 pulses[0] = numpy.ones(pulse_size)
310 310 pulses = 600*pulses[0]
311 311
312 312 return pulses,pulse_size
313 313
314 314 def jro_GenerateBlockOfData(self,Samples=Samples,DC_level= DC_level,stdev=stdev,
315 315 Reference= Reference,pulses= pulses,
316 316 Num_Codes= Num_Codes,pulse_size=pulse_size,
317 317 prof_gen= prof_gen,H0 = H0,DH0=DH0,
318 318 Adoppler=Adoppler,Fdoppler= Fdoppler,Hdoppler=Hdoppler):
319 319 Samples = Samples
320 320 DC_level = DC_level
321 321 stdev = stdev
322 322 m_nR = Reference
323 323 pulses = pulses
324 324 num_codes = Num_Codes
325 325 ps = pulse_size
326 326 prof_gen = prof_gen
327 327 channels = self.channels
328 328 H0 = H0
329 329 DH0 = DH0
330 330 ippSec = self.radarControllerHeaderObj.ippSeconds
331 331 Fdoppler = self.Fdoppler
332 332 Hdoppler = self.Hdoppler
333 333 Adoppler = self.Adoppler
334 334
335 335 self.datablock = numpy.zeros([channels,prof_gen,Samples],dtype= numpy.complex64)
336 336 for i in range(channels):
337 337 for k in range(prof_gen):
338 338 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·NOISEΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
339 339 Noise_r = numpy.random.normal(DC_level,stdev,Samples)
340 340 Noise_i = numpy.random.normal(DC_level,stdev,Samples)
341 341 Noise = numpy.zeros(Samples,dtype=complex)
342 342 Noise.real = Noise_r
343 343 Noise.imag = Noise_i
344 344 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·PULSOSΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
345 345 Pulso = numpy.zeros(pulse_size,dtype=complex)
346 346 Pulso.real = pulses[k%num_codes]
347 347 Pulso.imag = pulses[k%num_codes]
348 348 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· PULSES+NOISEΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·
349 349 InBuffer = numpy.zeros(Samples,dtype=complex)
350 350 InBuffer[m_nR:m_nR+ps] = Pulso
351 351 InBuffer = InBuffer+Noise
352 352 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· ANGLE Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
353 353 InBuffer.real[m_nR:m_nR+ps] = InBuffer.real[m_nR:m_nR+ps]*(math.cos( self.fAngle)*5)
354 354 InBuffer.imag[m_nR:m_nR+ps] = InBuffer.imag[m_nR:m_nR+ps]*(math.sin( self.fAngle)*5)
355 355 InBuffer=InBuffer
356 356 self.datablock[i][k]= InBuffer
357 357
358 358 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·DOPPLER SIGNAL...............................................
359 359 time_vec = numpy.linspace(0,(prof_gen-1)*ippSec,int(prof_gen))+self.nReadBlocks*ippSec*prof_gen+(self.nReadFiles-1)*ippSec*prof_gen
360 360 fd = Fdoppler #+(600.0/120)*self.nReadBlocks
361 361 d_signal = Adoppler*numpy.array(numpy.exp(1.0j*2.0*math.pi*fd*time_vec),dtype=numpy.complex64)
362 362 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·SeΓ±al con ancho espectralΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
363 #specw_sig = numpy.linspace(-149,150,300)
364 #w = 8
365 #A = 20
366 #specw_sig = specw_sig/w
367 #specw_sig = numpy.sinc(specw_sig)
368 #specw_sig = A*numpy.array(specw_sig,dtype=numpy.complex64)
363 if prof_gen%2==0:
364 min = int(prof_gen/2.0-1.0)
365 max = int(prof_gen/2.0)
366 else:
367 min = int(prof_gen/2.0)
368 max = int(prof_gen/2.0)
369 specw_sig = numpy.linspace(-min,max,prof_gen)
370 w = 4
371 A = 20
372 specw_sig = specw_sig/w
373 specw_sig = numpy.sinc(specw_sig)
374 specw_sig = A*numpy.array(specw_sig,dtype=numpy.complex64)
369 375 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· DATABLOCK + DOPPLERΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
370 376 HD=int(Hdoppler/self.AcqDH_0)
371 377 for i in range(12):
372 378 self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ d_signal# RESULT
373 379 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· DATABLOCK + DOPPLER*Sinc(x)Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
374 #HD=int(Hdoppler/self.AcqDH_0)
375 #HD=int(HD/2)
376 #for i in range(12):
377 # self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ specw_sig*d_signal# RESULT
380 HD=int(Hdoppler/self.AcqDH_0)
381 HD=int(HD/2)
382 for i in range(12):
383 self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ specw_sig*d_signal# RESULT
378 384
379 385 def readBlock(self):
380 386
381 387 self.jro_GenerateBlockOfData(Samples= self.samples,DC_level=self.DC_level,
382 388 stdev=self.stdev,Reference= self.Reference,
383 389 pulses = self.pulses,Num_Codes=self.Num_Codes,
384 390 pulse_size=self.pulse_size,prof_gen=self.profiles,
385 391 H0=self.H0,DH0=self.DH0)
386 392
387 393 self.profileIndex = 0
388 394 self.flagIsNewFile = 0
389 395 self.flagIsNewBlock = 1
390 396 self.nTotalBlocks += 1
391 397 self.nReadBlocks += 1
392 398
393 399 return 1
394 400
395 401
396 402 def getData(self):
397 403 if self.flagNoMoreFiles:
398 404 self.dataOut.flagNodata = True
399 405 return 0
400 406 self.flagDiscontinuousBlock = 0
401 407 self.flagIsNewBlock = 0
402 408 if self.__hasNotDataInBuffer(): # aqui es verdad
403 409 if not(self.readNextBlock()): # return 1 y por eso el if not salta a getBasic Header
404 410 return 0
405 411 self.getFirstHeader() # atributo
406 412
407 413 if not self.getByBlock:
408 414 self.dataOut.flagDataAsBlock = False
409 415 self.dataOut.data = self.datablock[:, self.profileIndex, :]
410 416 self.dataOut.profileIndex = self.profileIndex
411 417 self.profileIndex += 1
412 418 else:
413 419 pass
414 420 self.dataOut.flagNoData = False
415 421 self.getBasicHeader()
416 422 self.dataOut.realtime = self.online
417 423 return self.dataOut.data
418 424
419 425
420 426 def setup(self,frequency=49.92e6,incIntFactor= 1, nFFTPoints = 0, FixPP_IncInt=1,FixRCP_IPP=1000,
421 427 FixPP_CohInt= 1,Tau_0= 250,AcqH0_0 = 70 ,AcqDH_0=1.25, Bauds= 32,
422 428 FixRCP_TXA = 40, FixRCP_TXB = 50, fAngle = 2.0*math.pi*(1/16),DC_level= 50,
423 429 stdev= 8,Num_Codes = 1 , Dyn_snCode = None, samples=200,
424 channels=2,Fdoppler=20,Hdoppler=36,Adoppler=500,nTotalReadFiles=10000,
430 channels=2,Fdoppler=20,Hdoppler=36,Adoppler=500,
431 profilesPerBlock=300,dataBlocksPerFile=120,nTotalReadFiles=10000,
425 432 **kwargs):
426 433
427 434 self.set_kwargs(**kwargs)
428 435 self.nReadBlocks = 0
429 436 self.nReadFiles = 1
430 437 print('------------------- [Opening file: ] ------------------------------',self.nReadFiles)
431 438
432 439 tmp = time.time()
433 440 tmp_utc = int(tmp)
434 441 tmp_milisecond = int((tmp-tmp_utc)*1000)
435 442 print(" SETUP -basicHeaderObj.utc",datetime.datetime.utcfromtimestamp(tmp))
436 443 if Dyn_snCode is None:
437 444 Num_Codes=1
438 445 Bauds =1
439 446
440 447
441 448
442 449 self.set_BH(utc= tmp_utc,miliSecond= tmp_milisecond,timeZone=300 )
443 450 self.set_RCH( expType=0, nTx=150,ipp=FixRCP_IPP, txA=FixRCP_TXA, txB= FixRCP_TXB,
444 451 nWindows=1 , nHeights=samples, firstHeight=AcqH0_0, deltaHeight=AcqDH_0,
445 452 numTaus=1, line6Function=0, line5Function=0, fClock=None,
446 453 prePulseBefore=0, prePulseAfter=0,
447 454 codeType=0, nCode=Num_Codes, nBaud=32, code=Dyn_snCode,
448 455 flip1=0, flip2=0,Taus=Tau_0)
449 456
450 self.set_PH(dtype=0, blockSize=0, profilesPerBlock=300,
451 dataBlocksPerFile=120, nWindows=1, processFlags=numpy.array([1024]), nCohInt=1,
457 self.set_PH(dtype=0, blockSize=0, profilesPerBlock=profilesPerBlock,
458 dataBlocksPerFile=dataBlocksPerFile, nWindows=1, processFlags=numpy.array([1024]), nCohInt=1,
452 459 nIncohInt=1, totalSpectra=0, nHeights=samples, firstHeight=AcqH0_0,
453 460 deltaHeight=AcqDH_0, samplesWin=samples, spectraComb=0, nCode=0,
454 461 code=0, nBaud=None, shif_fft=False, flag_dc=False,
455 462 flag_cspc=False, flag_decode=False, flag_deflip=False)
456 463
457 self.set_SH(nSamples=samples, nProfiles=300, nChannels=channels)
464 self.set_SH(nSamples=samples, nProfiles=profilesPerBlock, nChannels=channels)
458 465
459 466 self.readFirstHeader()
460 467
461 468 self.frequency = frequency
462 469 self.incIntFactor = incIntFactor
463 470 self.nFFTPoints = nFFTPoints
464 471 self.FixPP_IncInt = FixPP_IncInt
465 472 self.FixRCP_IPP = FixRCP_IPP
466 473 self.FixPP_CohInt = FixPP_CohInt
467 474 self.Tau_0 = Tau_0
468 475 self.AcqH0_0 = AcqH0_0
469 476 self.H0 = AcqH0_0
470 477 self.AcqDH_0 = AcqDH_0
471 478 self.DH0 = AcqDH_0
472 479 self.Bauds = Bauds
473 480 self.FixRCP_TXA = FixRCP_TXA
474 481 self.FixRCP_TXB = FixRCP_TXB
475 482 self.fAngle = fAngle
476 483 self.DC_level = DC_level
477 484 self.stdev = stdev
478 485 self.Num_Codes = Num_Codes
479 486 self.Dyn_snCode = Dyn_snCode
480 487 self.samples = samples
481 488 self.channels = channels
482 489 self.profiles = None
483 490 self.m_nReference = None
484 491 self.Baudwidth = None
485 492 self.Fdoppler = Fdoppler
486 493 self.Hdoppler = Hdoppler
487 494 self.Adoppler = Adoppler
488 495 self.nTotalReadFiles = int(nTotalReadFiles)
489 496
490 497 print("IPP ", self.FixRCP_IPP)
491 498 print("Tau_0 ",self.Tau_0)
492 499 print("AcqH0_0",self.AcqH0_0)
493 500 print("samples,window ",self.samples)
494 501 print("AcqDH_0",AcqDH_0)
495 502 print("FixRCP_TXA",self.FixRCP_TXA)
496 503 print("FixRCP_TXB",self.FixRCP_TXB)
497 504 print("Dyn_snCode",Dyn_snCode)
498 505 print("Fdoppler", Fdoppler)
499 506 print("Hdoppler",Hdoppler)
500 507 print("Vdopplermax",Fdoppler*(3.0e8/self.frequency)/2.0)
501 508 print("nTotalReadFiles", nTotalReadFiles)
502 509
503 510 self.init_acquisition()
504 511 self.pulses,self.pulse_size=self.init_pulse(Num_Codes=self.Num_Codes,Bauds=self.Bauds,BaudWidth=self.BaudWidth,Dyn_snCode=Dyn_snCode)
505 512 print(" [ END ] - SETUP metodo")
506 513 return
507 514
508 515 def run(self,**kwargs): # metodo propio
509 516 if not(self.isConfig):
510 517 self.setup(**kwargs)
511 518 self.isConfig = True
512 519 self.getData()
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
@@ -1,1604 +1,1627
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 from schainpy.model.data.jrodata import Voltage
5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62 62 for channel in channelList:
63 63 if channel not in self.dataOut.channelList:
64 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
65 65
66 66 index = self.dataOut.channelList.index(channel)
67 67 channelIndexList.append(index)
68 68 self.selectChannelsByIndex(channelIndexList)
69 69 return self.dataOut
70 70
71 71 def selectChannelsByIndex(self, channelIndexList):
72 72 """
73 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
74 74
75 75 Input:
76 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
77 77
78 78 Affected:
79 79 self.dataOut.data
80 80 self.dataOut.channelIndexList
81 81 self.dataOut.nChannels
82 82 self.dataOut.m_ProcessingHeader.totalSpectra
83 83 self.dataOut.systemHeaderObj.numChannels
84 84 self.dataOut.m_ProcessingHeader.blockSize
85 85
86 86 Return:
87 87 None
88 88 """
89 89
90 90 for channelIndex in channelIndexList:
91 91 if channelIndex not in self.dataOut.channelIndexList:
92 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
93 93
94 94 if self.dataOut.type == 'Voltage':
95 95 if self.dataOut.flagDataAsBlock:
96 96 """
97 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
98 98 """
99 99 data = self.dataOut.data[channelIndexList,:,:]
100 100 else:
101 101 data = self.dataOut.data[channelIndexList,:]
102 102
103 103 self.dataOut.data = data
104 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
105 105 self.dataOut.channelList = range(len(channelIndexList))
106 106
107 107 elif self.dataOut.type == 'Spectra':
108 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
109 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
110 110
111 111 self.dataOut.data_spc = data_spc
112 112 self.dataOut.data_dc = data_dc
113 113
114 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 self.dataOut.channelList = range(len(channelIndexList))
116 116 self.__selectPairsByChannel(channelIndexList)
117 117
118 118 return 1
119 119
120 120 def __selectPairsByChannel(self, channelList=None):
121 121
122 122 if channelList == None:
123 123 return
124 124
125 125 pairsIndexListSelected = []
126 126 for pairIndex in self.dataOut.pairsIndexList:
127 127 # First pair
128 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
129 129 continue
130 130 # Second pair
131 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
132 132 continue
133 133
134 134 pairsIndexListSelected.append(pairIndex)
135 135
136 136 if not pairsIndexListSelected:
137 137 self.dataOut.data_cspc = None
138 138 self.dataOut.pairsList = []
139 139 return
140 140
141 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
142 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
143 143 for i in pairsIndexListSelected]
144 144
145 145 return
146 146
147 147 class selectHeights(Operation):
148 148
149 149 def run(self, dataOut, minHei=None, maxHei=None):
150 150 """
151 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
152 152 minHei <= height <= maxHei
153 153
154 154 Input:
155 155 minHei : valor minimo de altura a considerar
156 156 maxHei : valor maximo de altura a considerar
157 157
158 158 Affected:
159 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
160 160
161 161 Return:
162 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
163 163 """
164 164
165 165 self.dataOut = dataOut
166 166
167 167 if minHei == None:
168 168 minHei = self.dataOut.heightList[0]
169 169
170 170 if maxHei == None:
171 171 maxHei = self.dataOut.heightList[-1]
172 172
173 173 if (minHei < self.dataOut.heightList[0]):
174 174 minHei = self.dataOut.heightList[0]
175 175
176 176 if (maxHei > self.dataOut.heightList[-1]):
177 177 maxHei = self.dataOut.heightList[-1]
178 178
179 179 minIndex = 0
180 180 maxIndex = 0
181 181 heights = self.dataOut.heightList
182 182
183 183 inda = numpy.where(heights >= minHei)
184 184 indb = numpy.where(heights <= maxHei)
185 185
186 186 try:
187 187 minIndex = inda[0][0]
188 188 except:
189 189 minIndex = 0
190 190
191 191 try:
192 192 maxIndex = indb[0][-1]
193 193 except:
194 194 maxIndex = len(heights)
195 195
196 196 self.selectHeightsByIndex(minIndex, maxIndex)
197 197
198 198 return self.dataOut
199 199
200 200 def selectHeightsByIndex(self, minIndex, maxIndex):
201 201 """
202 202 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
203 203 minIndex <= index <= maxIndex
204 204
205 205 Input:
206 206 minIndex : valor de indice minimo de altura a considerar
207 207 maxIndex : valor de indice maximo de altura a considerar
208 208
209 209 Affected:
210 210 self.dataOut.data
211 211 self.dataOut.heightList
212 212
213 213 Return:
214 214 1 si el metodo se ejecuto con exito caso contrario devuelve 0
215 215 """
216 216
217 217 if self.dataOut.type == 'Voltage':
218 218 if (minIndex < 0) or (minIndex > maxIndex):
219 219 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
220 220
221 221 if (maxIndex >= self.dataOut.nHeights):
222 222 maxIndex = self.dataOut.nHeights
223 223
224 224 #voltage
225 225 if self.dataOut.flagDataAsBlock:
226 226 """
227 227 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
228 228 """
229 229 data = self.dataOut.data[:,:, minIndex:maxIndex]
230 230 else:
231 231 data = self.dataOut.data[:, minIndex:maxIndex]
232 232
233 233 # firstHeight = self.dataOut.heightList[minIndex]
234 234
235 235 self.dataOut.data = data
236 236 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
237 237
238 238 if self.dataOut.nHeights <= 1:
239 239 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
240 240 elif self.dataOut.type == 'Spectra':
241 241 if (minIndex < 0) or (minIndex > maxIndex):
242 242 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
243 243 minIndex, maxIndex))
244 244
245 245 if (maxIndex >= self.dataOut.nHeights):
246 246 maxIndex = self.dataOut.nHeights - 1
247 247
248 248 # Spectra
249 249 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
250 250
251 251 data_cspc = None
252 252 if self.dataOut.data_cspc is not None:
253 253 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
254 254
255 255 data_dc = None
256 256 if self.dataOut.data_dc is not None:
257 257 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
258 258
259 259 self.dataOut.data_spc = data_spc
260 260 self.dataOut.data_cspc = data_cspc
261 261 self.dataOut.data_dc = data_dc
262 262
263 263 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
264 264
265 265 return 1
266 266
267 267
268 268 class filterByHeights(Operation):
269 269
270 270 def run(self, dataOut, window):
271 271
272 272 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
273 273
274 274 if window == None:
275 275 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
276 276
277 277 newdelta = deltaHeight * window
278 278 r = dataOut.nHeights % window
279 279 newheights = (dataOut.nHeights-r)/window
280 280
281 281 if newheights <= 1:
282 282 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
283 283
284 284 if dataOut.flagDataAsBlock:
285 285 """
286 286 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
287 287 """
288 288 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
289 289 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
290 290 buffer = numpy.sum(buffer,3)
291 291
292 292 else:
293 293 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
294 294 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
295 295 buffer = numpy.sum(buffer,2)
296 296
297 297 dataOut.data = buffer
298 298 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
299 299 dataOut.windowOfFilter = window
300 300
301 301 return dataOut
302 302
303 303
304 304 class setH0(Operation):
305 305
306 306 def run(self, dataOut, h0, deltaHeight = None):
307 307
308 308 if not deltaHeight:
309 309 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
310 310
311 311 nHeights = dataOut.nHeights
312 312
313 313 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
314 314
315 315 dataOut.heightList = newHeiRange
316 316
317 317 return dataOut
318 318
319 319
320 320 class deFlip(Operation):
321 321
322 322 def run(self, dataOut, channelList = []):
323 323
324 324 data = dataOut.data.copy()
325 325
326 326 if dataOut.flagDataAsBlock:
327 327 flip = self.flip
328 328 profileList = list(range(dataOut.nProfiles))
329 329
330 330 if not channelList:
331 331 for thisProfile in profileList:
332 332 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
333 333 flip *= -1.0
334 334 else:
335 335 for thisChannel in channelList:
336 336 if thisChannel not in dataOut.channelList:
337 337 continue
338 338
339 339 for thisProfile in profileList:
340 340 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
341 341 flip *= -1.0
342 342
343 343 self.flip = flip
344 344
345 345 else:
346 346 if not channelList:
347 347 data[:,:] = data[:,:]*self.flip
348 348 else:
349 349 for thisChannel in channelList:
350 350 if thisChannel not in dataOut.channelList:
351 351 continue
352 352
353 353 data[thisChannel,:] = data[thisChannel,:]*self.flip
354 354
355 355 self.flip *= -1.
356 356
357 357 dataOut.data = data
358 358
359 359 return dataOut
360 360
361 361
362 362 class setAttribute(Operation):
363 363 '''
364 364 Set an arbitrary attribute(s) to dataOut
365 365 '''
366 366
367 367 def __init__(self):
368 368
369 369 Operation.__init__(self)
370 370 self._ready = False
371 371
372 372 def run(self, dataOut, **kwargs):
373 373
374 374 for key, value in kwargs.items():
375 375 setattr(dataOut, key, value)
376 376
377 377 return dataOut
378 378
379 379
380 380 @MPDecorator
381 381 class printAttribute(Operation):
382 382 '''
383 383 Print an arbitrary attribute of dataOut
384 384 '''
385 385
386 386 def __init__(self):
387 387
388 388 Operation.__init__(self)
389 389
390 390 def run(self, dataOut, attributes):
391 391
392 392 for attr in attributes:
393 393 if hasattr(dataOut, attr):
394 394 log.log(getattr(dataOut, attr), attr)
395 395
396 396
397 397 class interpolateHeights(Operation):
398 398
399 399 def run(self, dataOut, topLim, botLim):
400 400 #69 al 72 para julia
401 401 #82-84 para meteoros
402 402 if len(numpy.shape(dataOut.data))==2:
403 403 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
404 404 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
405 405 #dataOut.data[:,botLim:limSup+1] = sampInterp
406 406 dataOut.data[:,botLim:topLim+1] = sampInterp
407 407 else:
408 408 nHeights = dataOut.data.shape[2]
409 409 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
410 410 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
411 411 f = interpolate.interp1d(x, y, axis = 2)
412 412 xnew = numpy.arange(botLim,topLim+1)
413 413 ynew = f(xnew)
414 414 dataOut.data[:,:,botLim:topLim+1] = ynew
415 415
416 416 return dataOut
417 417
418 418
419 419 class CohInt(Operation):
420 420
421 421 isConfig = False
422 422 __profIndex = 0
423 423 __byTime = False
424 424 __initime = None
425 425 __lastdatatime = None
426 426 __integrationtime = None
427 427 __buffer = None
428 428 __bufferStride = []
429 429 __dataReady = False
430 430 __profIndexStride = 0
431 431 __dataToPutStride = False
432 432 n = None
433 433
434 434 def __init__(self, **kwargs):
435 435
436 436 Operation.__init__(self, **kwargs)
437 437
438 438 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
439 439 """
440 440 Set the parameters of the integration class.
441 441
442 442 Inputs:
443 443
444 444 n : Number of coherent integrations
445 445 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
446 446 overlapping :
447 447 """
448 448
449 449 self.__initime = None
450 450 self.__lastdatatime = 0
451 451 self.__buffer = None
452 452 self.__dataReady = False
453 453 self.byblock = byblock
454 454 self.stride = stride
455 455
456 456 if n == None and timeInterval == None:
457 457 raise ValueError("n or timeInterval should be specified ...")
458 458
459 459 if n != None:
460 460 self.n = n
461 461 self.__byTime = False
462 462 else:
463 463 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
464 464 self.n = 9999
465 465 self.__byTime = True
466 466
467 467 if overlapping:
468 468 self.__withOverlapping = True
469 469 self.__buffer = None
470 470 else:
471 471 self.__withOverlapping = False
472 472 self.__buffer = 0
473 473
474 474 self.__profIndex = 0
475 475
476 476 def putData(self, data):
477 477
478 478 """
479 479 Add a profile to the __buffer and increase in one the __profileIndex
480 480
481 481 """
482 482
483 483 if not self.__withOverlapping:
484 484 self.__buffer += data.copy()
485 485 self.__profIndex += 1
486 486 return
487 487
488 488 #Overlapping data
489 489 nChannels, nHeis = data.shape
490 490 data = numpy.reshape(data, (1, nChannels, nHeis))
491 491
492 492 #If the buffer is empty then it takes the data value
493 493 if self.__buffer is None:
494 494 self.__buffer = data
495 495 self.__profIndex += 1
496 496 return
497 497
498 498 #If the buffer length is lower than n then stakcing the data value
499 499 if self.__profIndex < self.n:
500 500 self.__buffer = numpy.vstack((self.__buffer, data))
501 501 self.__profIndex += 1
502 502 return
503 503
504 504 #If the buffer length is equal to n then replacing the last buffer value with the data value
505 505 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
506 506 self.__buffer[self.n-1] = data
507 507 self.__profIndex = self.n
508 508 return
509 509
510 510
511 511 def pushData(self):
512 512 """
513 513 Return the sum of the last profiles and the profiles used in the sum.
514 514
515 515 Affected:
516 516
517 517 self.__profileIndex
518 518
519 519 """
520 520
521 521 if not self.__withOverlapping:
522 522 data = self.__buffer
523 523 n = self.__profIndex
524 524
525 525 self.__buffer = 0
526 526 self.__profIndex = 0
527 527
528 528 return data, n
529 529
530 530 #Integration with Overlapping
531 531 data = numpy.sum(self.__buffer, axis=0)
532 532 # print data
533 533 # raise
534 534 n = self.__profIndex
535 535
536 536 return data, n
537 537
538 538 def byProfiles(self, data):
539 539
540 540 self.__dataReady = False
541 541 avgdata = None
542 542 # n = None
543 543 # print data
544 544 # raise
545 545 self.putData(data)
546 546
547 547 if self.__profIndex == self.n:
548 548 avgdata, n = self.pushData()
549 549 self.__dataReady = True
550 550
551 551 return avgdata
552 552
553 553 def byTime(self, data, datatime):
554 554
555 555 self.__dataReady = False
556 556 avgdata = None
557 557 n = None
558 558
559 559 self.putData(data)
560 560
561 561 if (datatime - self.__initime) >= self.__integrationtime:
562 562 avgdata, n = self.pushData()
563 563 self.n = n
564 564 self.__dataReady = True
565 565
566 566 return avgdata
567 567
568 568 def integrateByStride(self, data, datatime):
569 569 # print data
570 570 if self.__profIndex == 0:
571 571 self.__buffer = [[data.copy(), datatime]]
572 572 else:
573 573 self.__buffer.append([data.copy(),datatime])
574 574 self.__profIndex += 1
575 575 self.__dataReady = False
576 576
577 577 if self.__profIndex == self.n * self.stride :
578 578 self.__dataToPutStride = True
579 579 self.__profIndexStride = 0
580 580 self.__profIndex = 0
581 581 self.__bufferStride = []
582 582 for i in range(self.stride):
583 583 current = self.__buffer[i::self.stride]
584 584 data = numpy.sum([t[0] for t in current], axis=0)
585 585 avgdatatime = numpy.average([t[1] for t in current])
586 586 # print data
587 587 self.__bufferStride.append((data, avgdatatime))
588 588
589 589 if self.__dataToPutStride:
590 590 self.__dataReady = True
591 591 self.__profIndexStride += 1
592 592 if self.__profIndexStride == self.stride:
593 593 self.__dataToPutStride = False
594 594 # print self.__bufferStride[self.__profIndexStride - 1]
595 595 # raise
596 596 return self.__bufferStride[self.__profIndexStride - 1]
597 597
598 598
599 599 return None, None
600 600
601 601 def integrate(self, data, datatime=None):
602 602
603 603 if self.__initime == None:
604 604 self.__initime = datatime
605 605
606 606 if self.__byTime:
607 607 avgdata = self.byTime(data, datatime)
608 608 else:
609 609 avgdata = self.byProfiles(data)
610 610
611 611
612 612 self.__lastdatatime = datatime
613 613
614 614 if avgdata is None:
615 615 return None, None
616 616
617 617 avgdatatime = self.__initime
618 618
619 619 deltatime = datatime - self.__lastdatatime
620 620
621 621 if not self.__withOverlapping:
622 622 self.__initime = datatime
623 623 else:
624 624 self.__initime += deltatime
625 625
626 626 return avgdata, avgdatatime
627 627
628 628 def integrateByBlock(self, dataOut):
629 629
630 630 times = int(dataOut.data.shape[1]/self.n)
631 631 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
632 632
633 633 id_min = 0
634 634 id_max = self.n
635 635
636 636 for i in range(times):
637 637 junk = dataOut.data[:,id_min:id_max,:]
638 638 avgdata[:,i,:] = junk.sum(axis=1)
639 639 id_min += self.n
640 640 id_max += self.n
641 641
642 642 timeInterval = dataOut.ippSeconds*self.n
643 643 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
644 644 self.__dataReady = True
645 645 return avgdata, avgdatatime
646 646
647 647 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
648 648
649 649 if not self.isConfig:
650 650 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
651 651 self.isConfig = True
652 652
653 653 if dataOut.flagDataAsBlock:
654 654 """
655 655 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
656 656 """
657 657 avgdata, avgdatatime = self.integrateByBlock(dataOut)
658 658 dataOut.nProfiles /= self.n
659 659 else:
660 660 if stride is None:
661 661 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
662 662 else:
663 663 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
664 664
665 665
666 666 # dataOut.timeInterval *= n
667 667 dataOut.flagNoData = True
668 668
669 669 if self.__dataReady:
670 670 dataOut.data = avgdata
671 671 if not dataOut.flagCohInt:
672 672 dataOut.nCohInt *= self.n
673 673 dataOut.flagCohInt = True
674 674 dataOut.utctime = avgdatatime
675 675 # print avgdata, avgdatatime
676 676 # raise
677 677 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
678 678 dataOut.flagNoData = False
679 679 return dataOut
680 680
681 681 class Decoder(Operation):
682 682
683 683 isConfig = False
684 684 __profIndex = 0
685 685
686 686 code = None
687 687
688 688 nCode = None
689 689 nBaud = None
690 690
691 691 def __init__(self, **kwargs):
692 692
693 693 Operation.__init__(self, **kwargs)
694 694
695 695 self.times = None
696 696 self.osamp = None
697 697 # self.__setValues = False
698 698 self.isConfig = False
699 699 self.setupReq = False
700 700 def setup(self, code, osamp, dataOut):
701 701
702 702 self.__profIndex = 0
703 703
704 704 self.code = code
705 705
706 706 self.nCode = len(code)
707 707 self.nBaud = len(code[0])
708 708
709 709 if (osamp != None) and (osamp >1):
710 710 self.osamp = osamp
711 711 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
712 712 self.nBaud = self.nBaud*self.osamp
713 713
714 714 self.__nChannels = dataOut.nChannels
715 715 self.__nProfiles = dataOut.nProfiles
716 716 self.__nHeis = dataOut.nHeights
717 717
718 718 if self.__nHeis < self.nBaud:
719 719 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
720 720
721 721 #Frequency
722 722 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
723 723
724 724 __codeBuffer[:,0:self.nBaud] = self.code
725 725
726 726 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
727 727
728 728 if dataOut.flagDataAsBlock:
729 729
730 730 self.ndatadec = self.__nHeis #- self.nBaud + 1
731 731
732 732 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
733 733
734 734 else:
735 735
736 736 #Time
737 737 self.ndatadec = self.__nHeis #- self.nBaud + 1
738 738
739 739 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
740 740
741 741 def __convolutionInFreq(self, data):
742 742
743 743 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
744 744
745 745 fft_data = numpy.fft.fft(data, axis=1)
746 746
747 747 conv = fft_data*fft_code
748 748
749 749 data = numpy.fft.ifft(conv,axis=1)
750 750
751 751 return data
752 752
753 753 def __convolutionInFreqOpt(self, data):
754 754
755 755 raise NotImplementedError
756 756
757 757 def __convolutionInTime(self, data):
758 758
759 759 code = self.code[self.__profIndex]
760 760 for i in range(self.__nChannels):
761 761 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
762 762
763 763 return self.datadecTime
764 764
765 765 def __convolutionByBlockInTime(self, data):
766 766
767 767 repetitions = int(self.__nProfiles / self.nCode)
768 768 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
769 769 junk = junk.flatten()
770 770 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
771 771 profilesList = range(self.__nProfiles)
772 772
773 773 for i in range(self.__nChannels):
774 774 for j in profilesList:
775 775 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
776 776 return self.datadecTime
777 777
778 778 def __convolutionByBlockInFreq(self, data):
779 779
780 780 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
781 781
782 782
783 783 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
784 784
785 785 fft_data = numpy.fft.fft(data, axis=2)
786 786
787 787 conv = fft_data*fft_code
788 788
789 789 data = numpy.fft.ifft(conv,axis=2)
790 790
791 791 return data
792 792
793 793
794 794 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
795 795
796 796 if dataOut.flagDecodeData:
797 797 print("This data is already decoded, recoding again ...")
798 798
799 799 if not self.isConfig:
800 800
801 801 if code is None:
802 802 if dataOut.code is None:
803 803 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
804 804
805 805 code = dataOut.code
806 806 else:
807 807 code = numpy.array(code).reshape(nCode,nBaud)
808 808 self.setup(code, osamp, dataOut)
809 809
810 810 self.isConfig = True
811 811
812 812 if mode == 3:
813 813 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
814 814
815 815 if times != None:
816 816 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
817 817
818 818 if self.code is None:
819 819 print("Fail decoding: Code is not defined.")
820 820 return
821 821
822 822 self.__nProfiles = dataOut.nProfiles
823 823 datadec = None
824 824
825 825 if mode == 3:
826 826 mode = 0
827 827
828 828 if dataOut.flagDataAsBlock:
829 829 """
830 830 Decoding when data have been read as block,
831 831 """
832 832
833 833 if mode == 0:
834 834 datadec = self.__convolutionByBlockInTime(dataOut.data)
835 835 if mode == 1:
836 836 datadec = self.__convolutionByBlockInFreq(dataOut.data)
837 837 else:
838 838 """
839 839 Decoding when data have been read profile by profile
840 840 """
841 841 if mode == 0:
842 842 datadec = self.__convolutionInTime(dataOut.data)
843 843
844 844 if mode == 1:
845 845 datadec = self.__convolutionInFreq(dataOut.data)
846 846
847 847 if mode == 2:
848 848 datadec = self.__convolutionInFreqOpt(dataOut.data)
849 849
850 850 if datadec is None:
851 851 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
852 852
853 853 dataOut.code = self.code
854 854 dataOut.nCode = self.nCode
855 855 dataOut.nBaud = self.nBaud
856 856
857 857 dataOut.data = datadec
858 858
859 859 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
860 860
861 861 dataOut.flagDecodeData = True #asumo q la data esta decodificada
862 862
863 863 if self.__profIndex == self.nCode-1:
864 864 self.__profIndex = 0
865 865 return dataOut
866 866
867 867 self.__profIndex += 1
868 868
869 869 return dataOut
870 870 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
871 871
872 872
873 873 class ProfileConcat(Operation):
874 874
875 875 isConfig = False
876 876 buffer = None
877 877
878 878 def __init__(self, **kwargs):
879 879
880 880 Operation.__init__(self, **kwargs)
881 881 self.profileIndex = 0
882 882
883 883 def reset(self):
884 884 self.buffer = numpy.zeros_like(self.buffer)
885 885 self.start_index = 0
886 886 self.times = 1
887 887
888 888 def setup(self, data, m, n=1):
889 889 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
890 890 self.nHeights = data.shape[1]#.nHeights
891 891 self.start_index = 0
892 892 self.times = 1
893 893
894 894 def concat(self, data):
895 895
896 896 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
897 897 self.start_index = self.start_index + self.nHeights
898 898
899 899 def run(self, dataOut, m):
900 900 dataOut.flagNoData = True
901 901
902 902 if not self.isConfig:
903 903 self.setup(dataOut.data, m, 1)
904 904 self.isConfig = True
905 905
906 906 if dataOut.flagDataAsBlock:
907 907 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
908 908
909 909 else:
910 910 self.concat(dataOut.data)
911 911 self.times += 1
912 912 if self.times > m:
913 913 dataOut.data = self.buffer
914 914 self.reset()
915 915 dataOut.flagNoData = False
916 916 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
917 917 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
918 918 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
919 919 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
920 920 dataOut.ippSeconds *= m
921 921 return dataOut
922 922
923 923 class ProfileSelector(Operation):
924 924
925 925 profileIndex = None
926 926 # Tamanho total de los perfiles
927 927 nProfiles = None
928 928
929 929 def __init__(self, **kwargs):
930 930
931 931 Operation.__init__(self, **kwargs)
932 932 self.profileIndex = 0
933 933
934 934 def incProfileIndex(self):
935 935
936 936 self.profileIndex += 1
937 937
938 938 if self.profileIndex >= self.nProfiles:
939 939 self.profileIndex = 0
940 940
941 941 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
942 942
943 943 if profileIndex < minIndex:
944 944 return False
945 945
946 946 if profileIndex > maxIndex:
947 947 return False
948 948
949 949 return True
950 950
951 951 def isThisProfileInList(self, profileIndex, profileList):
952 952
953 953 if profileIndex not in profileList:
954 954 return False
955 955
956 956 return True
957 957
958 958 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
959 959
960 960 """
961 961 ProfileSelector:
962 962
963 963 Inputs:
964 964 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
965 965
966 966 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
967 967
968 968 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
969 969
970 970 """
971 971
972 972 if rangeList is not None:
973 973 if type(rangeList[0]) not in (tuple, list):
974 974 rangeList = [rangeList]
975 975
976 976 dataOut.flagNoData = True
977 977
978 978 if dataOut.flagDataAsBlock:
979 979 """
980 980 data dimension = [nChannels, nProfiles, nHeis]
981 981 """
982 982 if profileList != None:
983 983 dataOut.data = dataOut.data[:,profileList,:]
984 984
985 985 if profileRangeList != None:
986 986 minIndex = profileRangeList[0]
987 987 maxIndex = profileRangeList[1]
988 988 profileList = list(range(minIndex, maxIndex+1))
989 989
990 990 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
991 991
992 992 if rangeList != None:
993 993
994 994 profileList = []
995 995
996 996 for thisRange in rangeList:
997 997 minIndex = thisRange[0]
998 998 maxIndex = thisRange[1]
999 999
1000 1000 profileList.extend(list(range(minIndex, maxIndex+1)))
1001 1001
1002 1002 dataOut.data = dataOut.data[:,profileList,:]
1003 1003
1004 1004 dataOut.nProfiles = len(profileList)
1005 1005 dataOut.profileIndex = dataOut.nProfiles - 1
1006 1006 dataOut.flagNoData = False
1007 1007
1008 1008 return dataOut
1009 1009
1010 1010 """
1011 1011 data dimension = [nChannels, nHeis]
1012 1012 """
1013 1013
1014 1014 if profileList != None:
1015 1015
1016 1016 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1017 1017
1018 1018 self.nProfiles = len(profileList)
1019 1019 dataOut.nProfiles = self.nProfiles
1020 1020 dataOut.profileIndex = self.profileIndex
1021 1021 dataOut.flagNoData = False
1022 1022
1023 1023 self.incProfileIndex()
1024 1024 return dataOut
1025 1025
1026 1026 if profileRangeList != None:
1027 1027
1028 1028 minIndex = profileRangeList[0]
1029 1029 maxIndex = profileRangeList[1]
1030 1030
1031 1031 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1032 1032
1033 1033 self.nProfiles = maxIndex - minIndex + 1
1034 1034 dataOut.nProfiles = self.nProfiles
1035 1035 dataOut.profileIndex = self.profileIndex
1036 1036 dataOut.flagNoData = False
1037 1037
1038 1038 self.incProfileIndex()
1039 1039 return dataOut
1040 1040
1041 1041 if rangeList != None:
1042 1042
1043 1043 nProfiles = 0
1044 1044
1045 1045 for thisRange in rangeList:
1046 1046 minIndex = thisRange[0]
1047 1047 maxIndex = thisRange[1]
1048 1048
1049 1049 nProfiles += maxIndex - minIndex + 1
1050 1050
1051 1051 for thisRange in rangeList:
1052 1052
1053 1053 minIndex = thisRange[0]
1054 1054 maxIndex = thisRange[1]
1055 1055
1056 1056 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1057 1057
1058 1058 self.nProfiles = nProfiles
1059 1059 dataOut.nProfiles = self.nProfiles
1060 1060 dataOut.profileIndex = self.profileIndex
1061 1061 dataOut.flagNoData = False
1062 1062
1063 1063 self.incProfileIndex()
1064 1064
1065 1065 break
1066 1066
1067 1067 return dataOut
1068 1068
1069 1069
1070 1070 if beam != None: #beam is only for AMISR data
1071 1071 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1072 1072 dataOut.flagNoData = False
1073 1073 dataOut.profileIndex = self.profileIndex
1074 1074
1075 1075 self.incProfileIndex()
1076 1076
1077 1077 return dataOut
1078 1078
1079 1079 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1080 1080
1081 1081
1082 1082 class Reshaper(Operation):
1083 1083
1084 1084 def __init__(self, **kwargs):
1085 1085
1086 1086 Operation.__init__(self, **kwargs)
1087 1087
1088 1088 self.__buffer = None
1089 1089 self.__nitems = 0
1090 1090
1091 1091 def __appendProfile(self, dataOut, nTxs):
1092 1092
1093 1093 if self.__buffer is None:
1094 1094 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1095 1095 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1096 1096
1097 1097 ini = dataOut.nHeights * self.__nitems
1098 1098 end = ini + dataOut.nHeights
1099 1099
1100 1100 self.__buffer[:, ini:end] = dataOut.data
1101 1101
1102 1102 self.__nitems += 1
1103 1103
1104 1104 return int(self.__nitems*nTxs)
1105 1105
1106 1106 def __getBuffer(self):
1107 1107
1108 1108 if self.__nitems == int(1./self.__nTxs):
1109 1109
1110 1110 self.__nitems = 0
1111 1111
1112 1112 return self.__buffer.copy()
1113 1113
1114 1114 return None
1115 1115
1116 1116 def __checkInputs(self, dataOut, shape, nTxs):
1117 1117
1118 1118 if shape is None and nTxs is None:
1119 1119 raise ValueError("Reshaper: shape of factor should be defined")
1120 1120
1121 1121 if nTxs:
1122 1122 if nTxs < 0:
1123 1123 raise ValueError("nTxs should be greater than 0")
1124 1124
1125 1125 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1126 1126 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1127 1127
1128 1128 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1129 1129
1130 1130 return shape, nTxs
1131 1131
1132 1132 if len(shape) != 2 and len(shape) != 3:
1133 1133 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1134 1134
1135 1135 if len(shape) == 2:
1136 1136 shape_tuple = [dataOut.nChannels]
1137 1137 shape_tuple.extend(shape)
1138 1138 else:
1139 1139 shape_tuple = list(shape)
1140 1140
1141 1141 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1142 1142
1143 1143 return shape_tuple, nTxs
1144 1144
1145 1145 def run(self, dataOut, shape=None, nTxs=None):
1146 1146
1147 1147 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1148 1148
1149 1149 dataOut.flagNoData = True
1150 1150 profileIndex = None
1151 1151
1152 1152 if dataOut.flagDataAsBlock:
1153 1153
1154 1154 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1155 1155 dataOut.flagNoData = False
1156 1156
1157 1157 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1158 1158
1159 1159 else:
1160 1160
1161 1161 if self.__nTxs < 1:
1162 1162
1163 1163 self.__appendProfile(dataOut, self.__nTxs)
1164 1164 new_data = self.__getBuffer()
1165 1165
1166 1166 if new_data is not None:
1167 1167 dataOut.data = new_data
1168 1168 dataOut.flagNoData = False
1169 1169
1170 1170 profileIndex = dataOut.profileIndex*nTxs
1171 1171
1172 1172 else:
1173 1173 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1174 1174
1175 1175 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1176 1176
1177 1177 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1178 1178
1179 1179 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1180 1180
1181 1181 dataOut.profileIndex = profileIndex
1182 1182
1183 1183 dataOut.ippSeconds /= self.__nTxs
1184 1184
1185 1185 return dataOut
1186 1186
1187 1187 class SplitProfiles(Operation):
1188 1188
1189 1189 def __init__(self, **kwargs):
1190 1190
1191 1191 Operation.__init__(self, **kwargs)
1192 1192
1193 1193 def run(self, dataOut, n):
1194 1194
1195 1195 dataOut.flagNoData = True
1196 1196 profileIndex = None
1197 1197
1198 1198 if dataOut.flagDataAsBlock:
1199 1199
1200 1200 #nchannels, nprofiles, nsamples
1201 1201 shape = dataOut.data.shape
1202 1202
1203 1203 if shape[2] % n != 0:
1204 1204 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1205 1205
1206 1206 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1207 1207
1208 1208 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1209 1209 dataOut.flagNoData = False
1210 1210
1211 1211 profileIndex = int(dataOut.nProfiles/n) - 1
1212 1212
1213 1213 else:
1214 1214
1215 1215 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1216 1216
1217 1217 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1218 1218
1219 1219 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1220 1220
1221 1221 dataOut.nProfiles = int(dataOut.nProfiles*n)
1222 1222
1223 1223 dataOut.profileIndex = profileIndex
1224 1224
1225 1225 dataOut.ippSeconds /= n
1226 1226
1227 1227 return dataOut
1228 1228
1229 1229 class CombineProfiles(Operation):
1230 1230 def __init__(self, **kwargs):
1231 1231
1232 1232 Operation.__init__(self, **kwargs)
1233 1233
1234 1234 self.__remData = None
1235 1235 self.__profileIndex = 0
1236 1236
1237 1237 def run(self, dataOut, n):
1238 1238
1239 1239 dataOut.flagNoData = True
1240 1240 profileIndex = None
1241 1241
1242 1242 if dataOut.flagDataAsBlock:
1243 1243
1244 1244 #nchannels, nprofiles, nsamples
1245 1245 shape = dataOut.data.shape
1246 1246 new_shape = shape[0], shape[1]/n, shape[2]*n
1247 1247
1248 1248 if shape[1] % n != 0:
1249 1249 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1250 1250
1251 1251 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1252 1252 dataOut.flagNoData = False
1253 1253
1254 1254 profileIndex = int(dataOut.nProfiles*n) - 1
1255 1255
1256 1256 else:
1257 1257
1258 1258 #nchannels, nsamples
1259 1259 if self.__remData is None:
1260 1260 newData = dataOut.data
1261 1261 else:
1262 1262 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1263 1263
1264 1264 self.__profileIndex += 1
1265 1265
1266 1266 if self.__profileIndex < n:
1267 1267 self.__remData = newData
1268 1268 #continue
1269 1269 return
1270 1270
1271 1271 self.__profileIndex = 0
1272 1272 self.__remData = None
1273 1273
1274 1274 dataOut.data = newData
1275 1275 dataOut.flagNoData = False
1276 1276
1277 1277 profileIndex = dataOut.profileIndex/n
1278 1278
1279 1279
1280 1280 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1281 1281
1282 1282 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1283 1283
1284 1284 dataOut.nProfiles = int(dataOut.nProfiles/n)
1285 1285
1286 1286 dataOut.profileIndex = profileIndex
1287 1287
1288 1288 dataOut.ippSeconds *= n
1289 1289
1290 1290 return dataOut
1291 1291
1292 1292 class PulsePairVoltage(Operation):
1293 1293 '''
1294 1294 Function PulsePair(Signal Power, Velocity)
1295 1295 The real component of Lag[0] provides Intensity Information
1296 1296 The imag component of Lag[1] Phase provides Velocity Information
1297 1297
1298 1298 Configuration Parameters:
1299 1299 nPRF = Number of Several PRF
1300 1300 theta = Degree Azimuth angel Boundaries
1301 1301
1302 1302 Input:
1303 1303 self.dataOut
1304 1304 lag[N]
1305 1305 Affected:
1306 1306 self.dataOut.spc
1307 1307 '''
1308 1308 isConfig = False
1309 1309 __profIndex = 0
1310 1310 __initime = None
1311 1311 __lastdatatime = None
1312 1312 __buffer = None
1313 1313 noise = None
1314 1314 __dataReady = False
1315 1315 n = None
1316 1316 __nch = 0
1317 1317 __nHeis = 0
1318 1318 removeDC = False
1319 1319 ipp = None
1320 1320 lambda_ = 0
1321 1321
1322 1322 def __init__(self,**kwargs):
1323 1323 Operation.__init__(self,**kwargs)
1324 1324
1325 1325 def setup(self, dataOut, n = None, removeDC=False):
1326 1326 '''
1327 1327 n= Numero de PRF's de entrada
1328 1328 '''
1329 1329 self.__initime = None
1330 1330 self.__lastdatatime = 0
1331 1331 self.__dataReady = False
1332 1332 self.__buffer = 0
1333 1333 self.__profIndex = 0
1334 1334 self.noise = None
1335 1335 self.__nch = dataOut.nChannels
1336 1336 self.__nHeis = dataOut.nHeights
1337 1337 self.removeDC = removeDC
1338 1338 self.lambda_ = 3.0e8/(9345.0e6)
1339 1339 self.ippSec = dataOut.ippSeconds
1340 1340 self.nCohInt = dataOut.nCohInt
1341 1341 print("IPPseconds",dataOut.ippSeconds)
1342 1342
1343 1343 print("ELVALOR DE n es:", n)
1344 1344 if n == None:
1345 1345 raise ValueError("n should be specified.")
1346 1346
1347 1347 if n != None:
1348 1348 if n<2:
1349 1349 raise ValueError("n should be greater than 2")
1350 1350
1351 1351 self.n = n
1352 1352 self.__nProf = n
1353 1353
1354 1354 self.__buffer = numpy.zeros((dataOut.nChannels,
1355 1355 n,
1356 1356 dataOut.nHeights),
1357 1357 dtype='complex')
1358 #self.noise = numpy.zeros([self.__nch,self.__nHeis])
1359 #for i in range(self.__nch):
1360 # self.noise[i]=dataOut.getNoise(channel=i)
1361 1358
1362 1359 def putData(self,data):
1363 1360 '''
1364 1361 Add a profile to he __buffer and increase in one the __profiel Index
1365 1362 '''
1366 1363 self.__buffer[:,self.__profIndex,:]= data
1367 1364 self.__profIndex += 1
1368 1365 return
1369 1366
1370 1367 def pushData(self,dataOut):
1371 1368 '''
1372 1369 Return the PULSEPAIR and the profiles used in the operation
1373 1370 Affected : self.__profileIndex
1374 1371 '''
1372 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Remove DCΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1375 1373 if self.removeDC==True:
1376 1374 mean = numpy.mean(self.__buffer,1)
1377 1375 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1378 1376 dc= numpy.tile(tmp,[1,self.__nProf,1])
1379 1377 self.__buffer = self.__buffer - dc
1378 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Calculo de Potencia Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1379 pair0 = self.__buffer*numpy.conj(self.__buffer)
1380 pair0 = pair0.real
1381 lag_0 = numpy.sum(pair0,1)
1382 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Calculo de Ruido x canalΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1383 self.noise = numpy.zeros(self.__nch)
1384 for i in range(self.__nch):
1385 daux = numpy.sort(pair0[i,:,:],axis= None)
1386 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1387
1388 self.noise = self.noise.reshape(self.__nch,1)
1389 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1390 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1391 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1392 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Potencia recibida= P , Potencia senal = S , Ruido= NΒ·Β·
1393 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· P= S+N ,P=lag_0/N Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1394 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Power Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1395 data_power = lag_0/(self.n*self.nCohInt)
1396 #------------------ Senal Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1397 data_intensity = pair0 - noise_buffer
1398 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1399 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1400 for i in range(self.__nch):
1401 for j in range(self.__nHeis):
1402 if data_intensity[i][j] < 0:
1403 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1380 1404
1381 lag_0 = numpy.sum(self.__buffer*numpy.conj(self.__buffer),1)
1382 data_intensity = lag_0/(self.n*self.nCohInt)#*self.nCohInt)
1383
1405 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Calculo de Frecuencia y Velocidad dopplerΒ·Β·Β·Β·Β·Β·Β·Β·
1384 1406 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1385 1407 lag_1 = numpy.sum(pair1,1)
1386 #angle = numpy.angle(numpy.sum(pair1,1))*180/(math.pi)
1387 data_velocity = (-1.0*self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(lag_1)#self.ippSec*self.nCohInt
1408 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1409 data_velocity = (self.lambda_/2.0)*data_freq
1388 1410
1389 self.noise = numpy.zeros([self.__nch,self.__nHeis])
1390 for i in range(self.__nch):
1391 self.noise[i]=dataOut.getNoise(channel=i)
1411 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Potencia promedio estimada de la SenalΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1412 lag_0 = lag_0/self.n
1413 S = lag_0-self.noise
1392 1414
1393 lag_0 = lag_0.real/(self.n)
1415 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Frecuencia Doppler promedio Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1394 1416 lag_1 = lag_1/(self.n-1)
1395 1417 R1 = numpy.abs(lag_1)
1396 S = (lag_0-self.noise)
1397 1418
1419 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Calculo del SNRΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1398 1420 data_snrPP = S/self.noise
1399 data_snrPP = numpy.where(data_snrPP<0,1,data_snrPP)
1421 for i in range(self.__nch):
1422 for j in range(self.__nHeis):
1423 if data_snrPP[i][j] < 1.e-20:
1424 data_snrPP[i][j] = 1.e-20
1400 1425
1426 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· Calculo del ancho espectral Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
1401 1427 L = S/R1
1402 1428 L = numpy.where(L<0,1,L)
1403 1429 L = numpy.log(L)
1404
1405 1430 tmp = numpy.sqrt(numpy.absolute(L))
1406
1407 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec))*tmp*numpy.sign(L)
1408 #data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec))*k
1431 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1409 1432 n = self.__profIndex
1410 1433
1411 1434 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1412 1435 self.__profIndex = 0
1413 return data_intensity,data_velocity,data_snrPP,data_specwidth,n
1436 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1437
1414 1438
1415 1439 def pulsePairbyProfiles(self,dataOut):
1416 1440
1417 1441 self.__dataReady = False
1442 data_power = None
1418 1443 data_intensity = None
1419 1444 data_velocity = None
1420 1445 data_specwidth = None
1421 1446 data_snrPP = None
1422 1447 self.putData(data=dataOut.data)
1423 1448 if self.__profIndex == self.n:
1424 #self.noise = numpy.zeros([self.__nch,self.__nHeis])
1425 #for i in range(self.__nch):
1426 # self.noise[i]=data.getNoise(channel=i)
1427 #print(self.noise.shape)
1428 data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1449 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1429 1450 self.__dataReady = True
1430 1451
1431 return data_intensity, data_velocity,data_snrPP,data_specwidth
1452 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1453
1432 1454
1433 1455 def pulsePairOp(self, dataOut, datatime= None):
1434 1456
1435 1457 if self.__initime == None:
1436 1458 self.__initime = datatime
1437 #print("hola")
1438 data_intensity, data_velocity,data_snrPP,data_specwidth = self.pulsePairbyProfiles(dataOut)
1459 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1439 1460 self.__lastdatatime = datatime
1440 1461
1441 if data_intensity is None:
1442 return None, None,None,None,None
1462 if data_power is None:
1463 return None, None, None,None,None,None
1443 1464
1444 1465 avgdatatime = self.__initime
1445 1466 deltatime = datatime - self.__lastdatatime
1446 1467 self.__initime = datatime
1447 1468
1448 return data_intensity, data_velocity,data_snrPP,data_specwidth,avgdatatime
1469 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1449 1470
1450 1471 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1451 1472
1452 1473 if not self.isConfig:
1453 1474 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1454 1475 self.isConfig = True
1455 data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1476 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1456 1477 dataOut.flagNoData = True
1457 1478
1458 1479 if self.__dataReady:
1459 1480 dataOut.nCohInt *= self.n
1460 dataOut.data_intensity = data_intensity #valor para intensidad
1461 dataOut.data_velocity = data_velocity #valor para velocidad
1462 dataOut.data_snrPP = data_snrPP # valor para snr
1463 dataOut.data_specwidth = data_specwidth
1481 dataOut.dataPP_POW = data_intensity # S
1482 dataOut.dataPP_POWER = data_power # P
1483 dataOut.dataPP_DOP = data_velocity
1484 dataOut.dataPP_SNR = data_snrPP
1485 dataOut.dataPP_WIDTH = data_specwidth
1464 1486 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1465 1487 dataOut.utctime = avgdatatime
1466 1488 dataOut.flagNoData = False
1467 1489 return dataOut
1468 1490
1469 1491
1492
1470 1493 # import collections
1471 1494 # from scipy.stats import mode
1472 1495 #
1473 1496 # class Synchronize(Operation):
1474 1497 #
1475 1498 # isConfig = False
1476 1499 # __profIndex = 0
1477 1500 #
1478 1501 # def __init__(self, **kwargs):
1479 1502 #
1480 1503 # Operation.__init__(self, **kwargs)
1481 1504 # # self.isConfig = False
1482 1505 # self.__powBuffer = None
1483 1506 # self.__startIndex = 0
1484 1507 # self.__pulseFound = False
1485 1508 #
1486 1509 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1487 1510 #
1488 1511 # #Read data
1489 1512 #
1490 1513 # powerdB = dataOut.getPower(channel = channel)
1491 1514 # noisedB = dataOut.getNoise(channel = channel)[0]
1492 1515 #
1493 1516 # self.__powBuffer.extend(powerdB.flatten())
1494 1517 #
1495 1518 # dataArray = numpy.array(self.__powBuffer)
1496 1519 #
1497 1520 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1498 1521 #
1499 1522 # maxValue = numpy.nanmax(filteredPower)
1500 1523 #
1501 1524 # if maxValue < noisedB + 10:
1502 1525 # #No se encuentra ningun pulso de transmision
1503 1526 # return None
1504 1527 #
1505 1528 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1506 1529 #
1507 1530 # if len(maxValuesIndex) < 2:
1508 1531 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1509 1532 # return None
1510 1533 #
1511 1534 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1512 1535 #
1513 1536 # #Seleccionar solo valores con un espaciamiento de nSamples
1514 1537 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1515 1538 #
1516 1539 # if len(pulseIndex) < 2:
1517 1540 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1518 1541 # return None
1519 1542 #
1520 1543 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1521 1544 #
1522 1545 # #remover senales que se distancien menos de 10 unidades o muestras
1523 1546 # #(No deberian existir IPP menor a 10 unidades)
1524 1547 #
1525 1548 # realIndex = numpy.where(spacing > 10 )[0]
1526 1549 #
1527 1550 # if len(realIndex) < 2:
1528 1551 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1529 1552 # return None
1530 1553 #
1531 1554 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1532 1555 # realPulseIndex = pulseIndex[realIndex]
1533 1556 #
1534 1557 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1535 1558 #
1536 1559 # print "IPP = %d samples" %period
1537 1560 #
1538 1561 # self.__newNSamples = dataOut.nHeights #int(period)
1539 1562 # self.__startIndex = int(realPulseIndex[0])
1540 1563 #
1541 1564 # return 1
1542 1565 #
1543 1566 #
1544 1567 # def setup(self, nSamples, nChannels, buffer_size = 4):
1545 1568 #
1546 1569 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1547 1570 # maxlen = buffer_size*nSamples)
1548 1571 #
1549 1572 # bufferList = []
1550 1573 #
1551 1574 # for i in range(nChannels):
1552 1575 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1553 1576 # maxlen = buffer_size*nSamples)
1554 1577 #
1555 1578 # bufferList.append(bufferByChannel)
1556 1579 #
1557 1580 # self.__nSamples = nSamples
1558 1581 # self.__nChannels = nChannels
1559 1582 # self.__bufferList = bufferList
1560 1583 #
1561 1584 # def run(self, dataOut, channel = 0):
1562 1585 #
1563 1586 # if not self.isConfig:
1564 1587 # nSamples = dataOut.nHeights
1565 1588 # nChannels = dataOut.nChannels
1566 1589 # self.setup(nSamples, nChannels)
1567 1590 # self.isConfig = True
1568 1591 #
1569 1592 # #Append new data to internal buffer
1570 1593 # for thisChannel in range(self.__nChannels):
1571 1594 # bufferByChannel = self.__bufferList[thisChannel]
1572 1595 # bufferByChannel.extend(dataOut.data[thisChannel])
1573 1596 #
1574 1597 # if self.__pulseFound:
1575 1598 # self.__startIndex -= self.__nSamples
1576 1599 #
1577 1600 # #Finding Tx Pulse
1578 1601 # if not self.__pulseFound:
1579 1602 # indexFound = self.__findTxPulse(dataOut, channel)
1580 1603 #
1581 1604 # if indexFound == None:
1582 1605 # dataOut.flagNoData = True
1583 1606 # return
1584 1607 #
1585 1608 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1586 1609 # self.__pulseFound = True
1587 1610 # self.__startIndex = indexFound
1588 1611 #
1589 1612 # #If pulse was found ...
1590 1613 # for thisChannel in range(self.__nChannels):
1591 1614 # bufferByChannel = self.__bufferList[thisChannel]
1592 1615 # #print self.__startIndex
1593 1616 # x = numpy.array(bufferByChannel)
1594 1617 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1595 1618 #
1596 1619 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1597 1620 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1598 1621 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1599 1622 #
1600 1623 # dataOut.data = self.__arrayBuffer
1601 1624 #
1602 1625 # self.__startIndex += self.__newNSamples
1603 1626 #
1604 1627 # return
@@ -1,49 +1,73
1 1 import os,sys
2 2 import datetime
3 3 import time
4 4 from schainpy.controller import Project
5 5 path = '/home/alex/Downloads/NEW_WR2/spc16removeDC'
6 6 figpath = path
7 7 desc = "Simulator Test"
8 8
9 9 controllerObj = Project()
10 10
11 11 controllerObj.setup(id='10',name='Test Simulator',description=desc)
12 12
13 13 readUnitConfObj = controllerObj.addReadUnit(datatype='SimulatorReader',
14 14 frequency=9.345e9,
15 15 FixRCP_IPP= 60,
16 16 Tau_0 = 30,
17 17 AcqH0_0=0,
18 18 samples=330,
19 19 AcqDH_0=0.15,
20 20 FixRCP_TXA=0.15,
21 21 FixRCP_TXB=0.15,
22 22 Fdoppler=600.0,
23 23 Hdoppler=36,
24 24 Adoppler=300,#300
25 25 delay=0,
26 26 online=0,
27 27 walk=0,
28 nTotalReadFiles=3)
29
28 profilesPerBlock=625,
29 dataBlocksPerFile=100)#,#nTotalReadFiles=2)
30 '''
31 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
32 path=path,
33 startDate="2020/01/01", #"2020/01/01",#today,
34 endDate= "2020/12/01", #"2020/12/30",#today,
35 startTime='00:00:00',
36 endTime='23:59:59',
37 delay=0,
38 #set=0,
39 online=0,
40 walk=1)
41 '''
30 42 opObj11 = readUnitConfObj.addOperation(name='printInfo')
31 43
32 44 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
33 45 #opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other')
34 46 #opObj11.addParameter(name='n', value='10', format='int')
35 47
36 48 #opObj10 = procUnitConfObjA.addOperation(name='selectChannels')
37 49 #opObj10.addParameter(name='channelList', value=[0])
38 50 opObj11 = procUnitConfObjA.addOperation(name='PulsePairVoltage', optype='other')
39 opObj11.addParameter(name='n', value='300', format='int')#10
51 opObj11.addParameter(name='n', value='625', format='int')#10
40 52 opObj11.addParameter(name='removeDC', value=1, format='int')
41 53
42 54 #opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')
55 #opObj11 = procUnitConfObjA.addOperation(name='PulsepairSignalPlot', optype='other')
56
43 57
44 opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')
58 #opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')
45 59 #opObj11.addParameter(name='xmax', value=8)
46 60
47 opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')
61 #opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')
62
63 procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
64
65
66 opObj10 = procUnitConfObjB.addOperation(name='ParameterWriter')
67 opObj10.addParameter(name='path',value=figpath)
68 #opObj10.addParameter(name='mode',value=0)
69 opObj10.addParameter(name='blocksPerFile',value='100',format='int')
70 opObj10.addParameter(name='metadataList',value='utctimeInit,timeInterval',format='list')
71 opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,dataPP_SNR,dataPP_WIDTH')#,format='list'
48 72
49 73 controllerObj.start()
General Comments 0
You need to be logged in to leave comments. Login now