diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index 449c772..0502c7e 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -691,9 +691,9 @@ class WeatherRHIPlot(Plot): #data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor)) data['azi'] = dataOut.data_azi data['ele'] = dataOut.data_ele - print("UPDATE") - print("data[weather]",data['weather'].shape) - print("data[azi]",data['azi']) + #print("UPDATE") + #print("data[weather]",data['weather'].shape) + #print("data[azi]",data['azi']) return data, meta def get2List(self,angulos): @@ -797,6 +797,13 @@ class WeatherRHIPlot(Plot): end = data_ele[-1] number = (end-start) len_ang=len(data_ele) + print("start",start) + print("end",end) + print("number",number) + + print("len_ang",len_ang) + + #exit(1) if start=len_ang or (numpy.argmin(data_ele)==0)):#caso subida return 0 @@ -829,6 +836,7 @@ class WeatherRHIPlot(Plot): #---------------------------------------------------------- tipo_case = self.check_case(data_ele,ang_max,ang_min) print("check_case",tipo_case) + #exit(1) #--------------------- new ------------------------- data_ele_new ,data_ele_old= self.globalCheckPED(data_ele,tipo_case) @@ -842,15 +850,21 @@ class WeatherRHIPlot(Plot): if tipo_case==0 or tipo_case==3: # SUBIDA n1= round(self.start_data_ele)- start n2= end - round(self.end_data_ele) + print(self.start_data_ele) + print(self.end_data_ele) if n1>0: ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1) ele1_nan= numpy.ones(n1)*numpy.nan data_ele = numpy.hstack((ele1,data_ele_new)) + print("ele1_nan",ele1_nan.shape) + print("data_ele_old",data_ele_old.shape) data_ele_old = numpy.hstack((ele1_nan,data_ele_old)) if n2>0: ele2= numpy.linspace(self.end_data_ele+1,end,n2) ele2_nan= numpy.ones(n2)*numpy.nan data_ele = numpy.hstack((data_ele,ele2)) + print("ele2_nan",ele2_nan.shape) + print("data_ele_old",data_ele_old.shape) data_ele_old = numpy.hstack((data_ele_old,ele2_nan)) if tipo_case==1 or tipo_case==2: # BAJADA @@ -1065,12 +1079,781 @@ class WeatherRHIPlot(Plot): for i,ax in enumerate(self.axes): self.res_weather[i], self.res_ele = self.const_ploteo(val_ch=i, data_weather=data['weather'][i][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min) self.res_azi = numpy.mean(data['azi']) + if i==0: + print("*****************************************************************************to plot**************************",self.res_weather[i].shape) if ax.firsttime: #plt.clf() cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj='cg',vmin=20, vmax=80) #fig=self.figures[0] else: #plt.clf() + if i==0: + print(self.res_weather[i]) + print(self.res_ele) + cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj='cg',vmin=20, vmax=80) + caax = cgax.parasites[0] + paax = cgax.parasites[1] + cbar = plt.gcf().colorbar(pm, pad=0.075) + caax.set_xlabel('x_range [km]') + caax.set_ylabel('y_range [km]') + plt.text(1.0, 1.05, 'Elevacion '+str(thisDatetime)+" Step "+str(self.ini)+ " Azi: "+str(round(self.res_azi,2)), transform=caax.transAxes, va='bottom',ha='right') + print("***************************self.ini****************************",self.ini) + self.ini= self.ini+1 + +class WeatherRHI_vRF2_Plot(Plot): + CODE = 'weather' + plot_name = 'weather' + plot_type = 'rhistyle' + buffering = False + data_ele_tmp = None + + def setup(self): + print("********************") + print("********************") + print("********************") + print("SETUP WEATHER PLOT") + self.ncols = 1 + self.nrows = 1 + self.nplots= 1 + self.ylabel= 'Range [Km]' + self.titles= ['Weather'] + if self.channels is not None: + self.nplots = len(self.channels) + self.nrows = len(self.channels) + else: + self.nplots = self.data.shape(self.CODE)[0] + self.nrows = self.nplots + self.channels = list(range(self.nplots)) + print("channels",self.channels) + print("que saldra", self.data.shape(self.CODE)[0]) + self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)] + print("self.titles",self.titles) + self.colorbar=False + self.width =8 + self.height =8 + self.ini =0 + self.len_azi =0 + self.buffer_ini = None + self.buffer_ele = None + self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08}) + self.flag =0 + self.indicador= 0 + self.last_data_ele = None + self.val_mean = None + + def update(self, dataOut): + + data = {} + meta = {} + if hasattr(dataOut, 'dataPP_POWER'): + factor = 1 + if hasattr(dataOut, 'nFFTPoints'): + factor = dataOut.normFactor + print("dataOut",dataOut.data_360.shape) + # + data['weather'] = 10*numpy.log10(dataOut.data_360/(factor)) + # + #data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor)) + data['azi'] = dataOut.data_azi + data['ele'] = dataOut.data_ele + data['case_flag'] = dataOut.case_flag + #print("UPDATE") + #print("data[weather]",data['weather'].shape) + #print("data[azi]",data['azi']) + return data, meta + + def get2List(self,angulos): + list1=[] + list2=[] + for i in reversed(range(len(angulos))): + if not i==0:#el caso de i=0 evalula el primero de la lista con el ultimo y no es relevante + diff_ = angulos[i]-angulos[i-1] + if abs(diff_) >1.5: + list1.append(i-1) + list2.append(diff_) + return list(reversed(list1)),list(reversed(list2)) + + def fixData90(self,list_,ang_): + if list_[0]==-1: + vec = numpy.where(ang_=90) + angulos[vec]=angulos[vec]-90 + return angulos + + + def search_pos(self,pos,list_): + for i in range(len(list_)): + if pos == list_[i]: + return True,i + i=None + return False,i + + def fixDataComp(self,ang_,list1_,list2_,tipo_case): + size = len(ang_) + size2 = 0 + for i in range(len(list2_)): + size2=size2+round(abs(list2_[i]))-1 + new_size= size+size2 + ang_new = numpy.zeros(new_size) + ang_new2 = numpy.zeros(new_size) + + tmp = 0 + c = 0 + for i in range(len(ang_)): + ang_new[tmp +c] = ang_[i] + ang_new2[tmp+c] = ang_[i] + condition , value = self.search_pos(i,list1_) + if condition: + pos = tmp + c + 1 + for k in range(round(abs(list2_[value]))-1): + if tipo_case==0 or tipo_case==3:#subida + ang_new[pos+k] = ang_new[pos+k-1]+1 + ang_new2[pos+k] = numpy.nan + elif tipo_case==1 or tipo_case==2:#bajada + ang_new[pos+k] = ang_new[pos+k-1]-1 + ang_new2[pos+k] = numpy.nan + + tmp = pos +k + c = 0 + c=c+1 + return ang_new,ang_new2 + + def globalCheckPED(self,angulos,tipo_case): + l1,l2 = self.get2List(angulos) + ##print("l1",l1) + ##print("l2",l2) + if len(l1)>0: + #angulos2 = self.fixData90(list_=l1,ang_=angulos) + #l1,l2 = self.get2List(angulos2) + ang1_,ang2_ = self.fixDataComp(ang_=angulos,list1_=l1,list2_=l2,tipo_case=tipo_case) + #ang1_ = self.fixData90HL(ang1_) + #ang2_ = self.fixData90HL(ang2_) + else: + ang1_= angulos + ang2_= angulos + return ang1_,ang2_ + + + def replaceNAN(self,data_weather,data_ele,val): + data= data_ele + data_T= data_weather + if data.shape[0]> data_T.shape[0]: + data_N = numpy.ones( [data.shape[0],data_T.shape[1]]) + c = 0 + for i in range(len(data)): + if numpy.isnan(data[i]): + data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan + else: + data_N[i,:]=data_T[c,:] + c=c+1 + return data_N + else: + for i in range(len(data)): + if numpy.isnan(data[i]): + data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan + return data_T + + def check_case(self,data_ele,ang_max,ang_min): + start = data_ele[0] + end = data_ele[-1] + number = (end-start) + len_ang=len(data_ele) + print("start",start) + print("end",end) + print("number",number) + + print("len_ang",len_ang) + + #exit(1) + + if start=len_ang or (numpy.argmin(data_ele)==0)):#caso subida + return 0 + #elif start>end and (round(abs(number)+1)>=len_ang or(numpy.argmax(data_ele)==0)):#caso bajada + # return 1 + elif round(abs(number)+1)>=len_ang and (start>end or(numpy.argmax(data_ele)==0)):#caso bajada + return 1 + elif round(abs(number)+1)data_ele[-1]:# caso BAJADA CAMBIO ANG MAX + return 2 + elif round(abs(number)+1)0: + ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1) + ele1_nan= numpy.ones(n1)*numpy.nan + data_ele = numpy.hstack((ele1,data_ele_new)) + print("ele1_nan",ele1_nan.shape) + print("data_ele_old",data_ele_old.shape) + data_ele_old = numpy.hstack((ele1_nan,data_ele_old)) + if n2>0: + ele2= numpy.linspace(self.end_data_ele+1,end,n2) + ele2_nan= numpy.ones(n2)*numpy.nan + data_ele = numpy.hstack((data_ele,ele2)) + print("ele2_nan",ele2_nan.shape) + print("data_ele_old",data_ele_old.shape) + data_ele_old = numpy.hstack((data_ele_old,ele2_nan)) + + if tipo_case==1 or tipo_case==2: # BAJADA + data_ele_new = data_ele_new[::-1] # reversa + data_ele_old = data_ele_old[::-1]# reversa + data_weather = data_weather[::-1,:]# reversa + vec= numpy.where(data_ele_new0: + ele1= numpy.linspace(ang_min+1,self.start_data_ele-1,n1) + ele1_nan= numpy.ones(n1)*numpy.nan + data_ele = numpy.hstack((ele1,data_ele_new)) + data_ele_old = numpy.hstack((ele1_nan,data_ele_old)) + if n2>0: + ele2= numpy.linspace(self.end_data_ele+1,end,n2) + ele2_nan= numpy.ones(n2)*numpy.nan + data_ele = numpy.hstack((data_ele,ele2)) + data_ele_old = numpy.hstack((data_ele_old,ele2_nan)) + # RADAR + # NOTA data_ele y data_weather es la variable que retorna + val_mean = numpy.mean(data_weather[:,-1]) + self.val_mean = val_mean + data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean) + print("eleold",data_ele_old) + print(self.data_ele_tmp[val_ch]) + print(data_ele_old.shape[0]) + print(self.data_ele_tmp[val_ch].shape[0]) + if (data_ele_old.shape[0]==91 or self.data_ele_tmp[val_ch].shape[0]==91): + import sys + print("EXIT",self.ini) + + sys.exit(1) + self.data_ele_tmp[val_ch]= data_ele_old + else: + #print("**********************************************") + #print("****************VARIABLE**********************") + #-------------------------CAMBIOS RHI--------------------------------- + #--------------------------------------------------------------------- + ##print("INPUT data_ele",data_ele) + flag=0 + start_ele = self.res_ele[0] + #tipo_case = self.check_case(data_ele,ang_max,ang_min) + tipo_case = case_flag[-1] + #print("TIPO DE DATA",tipo_case) + #-----------new------------ + data_ele ,data_ele_old = self.globalCheckPED(data_ele,tipo_case) + data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean) + + #-------------------------------NEW RHI ITERATIVO------------------------- + + if tipo_case==0 : # SUBIDA + vec = numpy.where(data_ele=ang_max-1: + self.data_ele_tmp[val_ch] = numpy.ones(ang_max-ang_min)*numpy.nan + self.res_weather[val_ch] = self.replaceNAN(data_weather=self.res_weather[val_ch],data_ele=self.data_ele_tmp[val_ch],val=self.val_mean) + self.data_ele_tmp[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1]=data_ele_old + self.res_ele[new_i_ele-1:new_i_ele+len(data_ele)-1]= data_ele + self.res_weather[val_ch][new_i_ele-1:new_i_ele+len(data_ele)-1,:]= data_weather + data_ele = self.res_ele + data_weather = self.res_weather[val_ch] + + elif tipo_case==2: #bajada + vec = numpy.where(data_ele0: + ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1) + ele1_nan= numpy.ones(n1)*numpy.nan + data_ele = numpy.hstack((ele1,data_ele_new)) + data_ele_old = numpy.hstack((ele1_nan,data_ele_new)) + if n2>0: + ele2= numpy.linspace(new_f_ele+1,ang_max,n2) + ele2_nan= numpy.ones(n2)*numpy.nan + data_ele = numpy.hstack((data_ele,ele2)) + data_ele_old = numpy.hstack((data_ele_old,ele2_nan)) + + self.data_ele_tmp[val_ch] = data_ele_old + self.res_ele = data_ele + self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean) + data_ele = self.res_ele + data_weather = self.res_weather[val_ch] + + elif tipo_case==3:#subida + vec = numpy.where(00: + len_vec= len(data_ele) + vec3 = numpy.linspace(pos_ini,len_vec-1,len_vec-pos_ini).astype(int) + #print(vec3) + data_ele= data_ele[vec3] + data_ele_new = data_ele + data_ele_old= data_ele_old[vec3] + data_weather= data_weather[vec3] + + new_i_ele = int(data_ele_new[0]) + new_f_ele = int(data_ele_new[-1]) + n1= new_i_ele- ang_min + n2= ang_max - new_f_ele-1 + if n1>0: + ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1) + ele1_nan= numpy.ones(n1)*numpy.nan + data_ele = numpy.hstack((ele1,data_ele_new)) + data_ele_old = numpy.hstack((ele1_nan,data_ele_new)) + if n2>0: + ele2= numpy.linspace(new_f_ele+1,ang_max,n2) + ele2_nan= numpy.ones(n2)*numpy.nan + data_ele = numpy.hstack((data_ele,ele2)) + data_ele_old = numpy.hstack((data_ele_old,ele2_nan)) + + self.data_ele_tmp[val_ch] = data_ele_old + self.res_ele = data_ele + self.res_weather[val_ch] = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean) + data_ele = self.res_ele + data_weather = self.res_weather[val_ch] + #print("self.data_ele_tmp",self.data_ele_tmp) + return data_weather,data_ele + + + def plot(self): + thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S') + data = self.data[-1] + r = self.data.yrange + delta_height = r[1]-r[0] + r_mask = numpy.where(r>=0)[0] + ##print("delta_height",delta_height) + #print("r_mask",r_mask,len(r_mask)) + r = numpy.arange(len(r_mask))*delta_height + self.y = 2*r + res = 1 + ###print("data['weather'].shape[0]",data['weather'].shape[0]) + ang_max = self.ang_max + ang_min = self.ang_min + var_ang =ang_max - ang_min + step = (int(var_ang)/(res*data['weather'].shape[0])) + ###print("step",step) + #-------------------------------------------------------- + ##print('weather',data['weather'].shape) + ##print('ele',data['ele'].shape) + + ###self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min) + ###self.res_azi = numpy.mean(data['azi']) + ###print("self.res_ele",self.res_ele) + plt.clf() + subplots = [121, 122] + try: + if self.data[-2]['ele'].max()1.5: + list1.append(i-1) + list2.append(diff_) + return list(reversed(list1)),list(reversed(list2)) + + def fixData90(self,list_,ang_): + if list_[0]==-1: + vec = numpy.where(ang_=90) + angulos[vec]=angulos[vec]-90 + return angulos + + + def search_pos(self,pos,list_): + for i in range(len(list_)): + if pos == list_[i]: + return True,i + i=None + return False,i + + def fixDataComp(self,ang_,list1_,list2_,tipo_case): + size = len(ang_) + size2 = 0 + for i in range(len(list2_)): + size2=size2+round(abs(list2_[i]))-1 + new_size= size+size2 + ang_new = numpy.zeros(new_size) + ang_new2 = numpy.zeros(new_size) + + tmp = 0 + c = 0 + for i in range(len(ang_)): + ang_new[tmp +c] = ang_[i] + ang_new2[tmp+c] = ang_[i] + condition , value = self.search_pos(i,list1_) + if condition: + pos = tmp + c + 1 + for k in range(round(abs(list2_[value]))-1): + if tipo_case==0 or tipo_case==3:#subida + ang_new[pos+k] = ang_new[pos+k-1]+1 + ang_new2[pos+k] = numpy.nan + elif tipo_case==1 or tipo_case==2:#bajada + ang_new[pos+k] = ang_new[pos+k-1]-1 + ang_new2[pos+k] = numpy.nan + + tmp = pos +k + c = 0 + c=c+1 + return ang_new,ang_new2 + + def globalCheckPED(self,angulos,tipo_case): + l1,l2 = self.get2List(angulos) + print("l1",l1) + print("l2",l2) + if len(l1)>0: + #angulos2 = self.fixData90(list_=l1,ang_=angulos) + #l1,l2 = self.get2List(angulos2) + ang1_,ang2_ = self.fixDataComp(ang_=angulos,list1_=l1,list2_=l2,tipo_case=tipo_case) + #ang1_ = self.fixData90HL(ang1_) + #ang2_ = self.fixData90HL(ang2_) + else: + ang1_= angulos + ang2_= angulos + return ang1_,ang2_ + + + def replaceNAN(self,data_weather,data_ele,val): + data= data_ele + data_T= data_weather + #print(data.shape[0]) + #print(data_T.shape[0]) + #exit(1) + if data.shape[0]> data_T.shape[0]: + data_N = numpy.ones( [data.shape[0],data_T.shape[1]]) + c = 0 + for i in range(len(data)): + if numpy.isnan(data[i]): + data_N[i,:]=numpy.ones(data_T.shape[1])*numpy.nan + else: + data_N[i,:]=data_T[c,:] + c=c+1 + return data_N + else: + for i in range(len(data)): + if numpy.isnan(data[i]): + data_T[i,:]=numpy.ones(data_T.shape[1])*numpy.nan + return data_T + + + def const_ploteo(self,val_ch,data_weather,data_ele,step,res,ang_max,ang_min,case_flag): + ang_max= ang_max + ang_min= ang_min + data_weather=data_weather + val_ch=val_ch + ##print("*********************DATA WEATHER**************************************") + ##print(data_weather) + + ''' + print("**********************************************") + print("**********************************************") + print("***************ini**************") + print("**********************************************") + print("**********************************************") + ''' + #print("data_ele",data_ele) + #---------------------------------------------------------- + + #exit(1) + tipo_case = case_flag[-1] + print("tipo_case",tipo_case) + #--------------------- new ------------------------- + data_ele_new ,data_ele_old= self.globalCheckPED(data_ele,tipo_case) + + #-------------------------CAMBIOS RHI--------------------------------- + + vec = numpy.where(data_ele0: + ele1= numpy.linspace(ang_min+1,new_i_ele-1,n1) + ele1_nan= numpy.ones(n1)*numpy.nan + data_ele = numpy.hstack((ele1,data_ele_new)) + data_ele_old = numpy.hstack((ele1_nan,data_ele_new)) + if n2>0: + ele2= numpy.linspace(new_f_ele+1,ang_max,n2) + ele2_nan= numpy.ones(n2)*numpy.nan + data_ele = numpy.hstack((data_ele,ele2)) + data_ele_old = numpy.hstack((data_ele_old,ele2_nan)) + + + print("ele shape",data_ele.shape) + print(data_ele) + + #print("self.data_ele_tmp",self.data_ele_tmp) + val_mean = numpy.mean(data_weather[:,-1]) + self.val_mean = val_mean + data_weather = self.replaceNAN(data_weather=data_weather,data_ele=data_ele_old,val=self.val_mean) + self.data_ele_tmp[val_ch]= data_ele_old + + + print("data_weather shape",data_weather.shape) + print(data_weather) + #exit(1) + return data_weather,data_ele + + + def plot(self): + thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S') + data = self.data[-1] + r = self.data.yrange + delta_height = r[1]-r[0] + r_mask = numpy.where(r>=0)[0] + ##print("delta_height",delta_height) + #print("r_mask",r_mask,len(r_mask)) + r = numpy.arange(len(r_mask))*delta_height + self.y = 2*r + res = 1 + ###print("data['weather'].shape[0]",data['weather'].shape[0]) + ang_max = self.ang_max + ang_min = self.ang_min + var_ang =ang_max - ang_min + step = (int(var_ang)/(res*data['weather'].shape[0])) + ###print("step",step) + #-------------------------------------------------------- + ##print('weather',data['weather'].shape) + ##print('ele',data['ele'].shape) + + ###self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min) + ###self.res_azi = numpy.mean(data['azi']) + ###print("self.res_ele",self.res_ele) + plt.clf() + subplots = [121, 122] + if self.ini==0: + self.data_ele_tmp = numpy.ones([self.nplots,int(var_ang)])*numpy.nan + self.res_weather= numpy.ones([self.nplots,int(var_ang),len(r_mask)])*numpy.nan + print("SHAPE",self.data_ele_tmp.shape) + + for i,ax in enumerate(self.axes): + self.res_weather[i], self.res_ele = self.const_ploteo(val_ch=i, data_weather=data['weather'][i][:,r_mask],data_ele=data['ele'],step=step,res=res,ang_max=ang_max,ang_min=ang_min,case_flag=self.data['case_flag']) + self.res_azi = numpy.mean(data['azi']) + + print(self.res_ele) + #exit(1) + if ax.firsttime: + #plt.clf() + cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj='cg',vmin=20, vmax=80) + #fig=self.figures[0] + else: + + #plt.clf() cgax, pm = wrl.vis.plot_rhi(self.res_weather[i],r=r,th=self.res_ele,ax=subplots[i], proj='cg',vmin=20, vmax=80) caax = cgax.parasites[0] paax = cgax.parasites[1] diff --git a/schainpy/model/io/jroIO_digitalRF.py b/schainpy/model/io/jroIO_digitalRF.py index bbf17ff..ba5cca8 100644 --- a/schainpy/model/io/jroIO_digitalRF.py +++ b/schainpy/model/io/jroIO_digitalRF.py @@ -565,7 +565,7 @@ class DigitalRFReader(ProcessingUnit): return True def __isBufferEmpty(self): - + return self.__bufferIndex > self.__samples_to_read - self.__nSamples # 40960 - 40 def getData(self, seconds=30, nTries=5): @@ -636,6 +636,7 @@ class DigitalRFReader(ProcessingUnit): self.dataOut.flagNoData = False buffer = self.__data_buffer[:,self.__bufferIndex:self.__bufferIndex + self.__samples_to_read] buffer = buffer.reshape((self.__nChannels, self.nProfileBlocks, int(self.__samples_to_read/self.nProfileBlocks))) + self.dataOut.nProfileBlocks = self.nProfileBlocks self.dataOut.data = buffer self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate self.profileIndex += self.__samples_to_read diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 6713ab2..f44f1aa 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -4071,6 +4071,7 @@ class PedestalInformation(Operation): # utc_ped_list.append(self.gettimeutcfromDirFilename(path=self.path_ped,file=list_pedestal[i])) nro_file,utc_ped,utc_ped_1 =self.getNROFile(utc_adq,utc_ped_list) #print("NROFILE************************************", nro_file,utc_ped) + #print(nro_file) if nro_file < 0: return numpy.NaN,numpy.NaN else: @@ -4090,6 +4091,8 @@ class PedestalInformation(Operation): return numpy.NaN,numpy.NaN #------------------------------------------------------------------------------------------------------ ''' + #print(int(self.samp_rate_ped)) + #print(nro_key_p) if int(self.samp_rate_ped)-1>=nro_key_p>0: #print("angulo_array :",angulo[nro_key_p]) return angulo[nro_key_p],angulo_ele[nro_key_p] @@ -4195,10 +4198,13 @@ class PedestalInformation(Operation): self.samp_rate_ped= samp_rate_ped self.t_Interval_p = t_Interval_p self.list_pedestal = self.getfirstFilefromPath(path=self.path_ped,meta="pos@",ext=".h5") + self.utc_ped_list= [] for i in range(len(self.list_pedestal)): #print(i,self.gettimeutcfromDirFilename(path=self.path_ped,file=self.list_pedestal[i]))# OJO IDENTIFICADOR DE SINCRONISMO self.utc_ped_list.append(self.gettimeutcfromDirFilename(path=self.path_ped,file=self.list_pedestal[i])) + #print(self.utc_ped_list) + #exit(1) #print("que paso") dataOut.wr_exp = wr_exp #print("SETUP READY") @@ -4218,6 +4224,8 @@ class PedestalInformation(Operation): def run(self, dataOut,path_ped,samp_rate_ped,t_Interval_p,wr_exp): #print("INTEGRATION -----") + #print("PEDESTAL") + if not self.isConfig: self.setup(dataOut, path_ped,samp_rate_ped,t_Interval_p,wr_exp) self.__dataReady = True @@ -4231,11 +4239,18 @@ class PedestalInformation(Operation): angulo,angulo_ele = self.getAnguloProfile(utc_adq=utc_adq,utc_ped_list=self.utc_ped_list) #print("angulo**********",angulo) dataOut.flagNoData = False + if numpy.isnan(angulo) or numpy.isnan(angulo_ele) : + #print("PEDESTAL 3") + #exit(1) dataOut.flagNoData = True return dataOut dataOut.azimuth = angulo dataOut.elevation = angulo_ele + #print("PEDESTAL END") + #print(dataOut.azimuth) + #print(dataOut.elevation) + #exit(1) return dataOut class Block360(Operation): @@ -4373,6 +4388,8 @@ class Block360(Operation): def run(self, dataOut,n = None,mode=None,**kwargs): #print("BLOCK 360 HERE WE GO MOMENTOS") + print("Block 360") + #exit(1) if not self.isConfig: self.setup(dataOut = dataOut, n = n ,mode= mode ,**kwargs) ####self.index = 0 @@ -4388,7 +4405,7 @@ class Block360(Operation): #print("DATA 360") #print(dataOut.data_360) #print("---------------------------------------------------------------------------------") - #print("---------------------------DATAREADY---------------------------------------------") + print("---------------------------DATAREADY---------------------------------------------") #print("---------------------------------------------------------------------------------") #print("data_360",dataOut.data_360.shape) dataOut.data_azi = data_p @@ -4399,3 +4416,207 @@ class Block360(Operation): dataOut.utctime = avgdatatime dataOut.flagNoData = False return dataOut + +class Block360_vRF(Operation): + ''' + ''' + isConfig = False + __profIndex = 0 + __initime = None + __lastdatatime = None + __buffer = None + __dataReady = False + n = None + __nch = 0 + __nHeis = 0 + index = 0 + mode = 0 + + def __init__(self,**kwargs): + Operation.__init__(self,**kwargs) + + def setup(self, dataOut, n = None, mode = None): + ''' + n= Numero de PRF's de entrada + ''' + self.__initime = None + self.__lastdatatime = 0 + self.__dataReady = False + self.__buffer = 0 + self.__buffer_1D = 0 + self.__profIndex = 0 + self.index = 0 + self.__nch = dataOut.nChannels + self.__nHeis = dataOut.nHeights + ##print("ELVALOR DE n es:", n) + if n == None: + raise ValueError("n should be specified.") + + if mode == None: + raise ValueError("mode should be specified.") + + if n != None: + if n<1: + print("n should be greater than 2") + raise ValueError("n should be greater than 2") + + self.n = n + self.mode = mode + #print("self.mode",self.mode) + #print("nHeights") + self.__buffer = numpy.zeros(( dataOut.nChannels,n, dataOut.nHeights)) + self.__buffer2 = numpy.zeros(n) + self.__buffer3 = numpy.zeros(n) + + + + + def putData(self,data,mode): + ''' + Add a profile to he __buffer and increase in one the __profiel Index + ''' + #print("line 4049",data.dataPP_POW.shape,data.dataPP_POW[:10]) + #print("line 4049",data.azimuth.shape,data.azimuth) + if self.mode==0: + self.__buffer[:,self.__profIndex,:]= data.dataPP_POWER# PRIMER MOMENTO + if self.mode==1: + self.__buffer[:,self.__profIndex,:]= data.data_pow + #print("me casi",self.index,data.azimuth[self.index]) + #print(self.__profIndex, self.index , data.azimuth[self.index] ) + #print("magic",data.profileIndex) + #print(data.azimuth[self.index]) + #print("index",self.index) + + #####self.__buffer2[self.__profIndex] = data.azimuth[self.index] + self.__buffer2[self.__profIndex] = data.azimuth + self.__buffer3[self.__profIndex] = data.elevation + #print("q pasa") + #####self.index+=1 + #print("index",self.index,data.azimuth[:10]) + self.__profIndex += 1 + return #················· Remove DC··································· + + def pushData(self,data): + ''' + Return the PULSEPAIR and the profiles used in the operation + Affected : self.__profileIndex + ''' + #print("pushData") + + data_360 = self.__buffer + data_p = self.__buffer2 + data_e = self.__buffer3 + n = self.__profIndex + + self.__buffer = numpy.zeros((self.__nch, self.n,self.__nHeis)) + self.__buffer2 = numpy.zeros(self.n) + self.__buffer3 = numpy.zeros(self.n) + self.__profIndex = 0 + #print("pushData") + return data_360,n,data_p,data_e + + + def byProfiles(self,dataOut): + + self.__dataReady = False + data_360 = None + data_p = None + data_e = None + #print("dataOu",dataOut.dataPP_POW) + self.putData(data=dataOut,mode = self.mode) + ##### print("profIndex",self.__profIndex) + if self.__profIndex == self.n: + data_360,n,data_p,data_e = self.pushData(data=dataOut) + self.__dataReady = True + + return data_360,data_p,data_e + + + def blockOp(self, dataOut, datatime= None): + if self.__initime == None: + self.__initime = datatime + data_360,data_p,data_e = self.byProfiles(dataOut) + self.__lastdatatime = datatime + + if data_360 is None: + return None, None,None,None + + + avgdatatime = self.__initime + if self.n==1: + avgdatatime = datatime + deltatime = datatime - self.__lastdatatime + self.__initime = datatime + #print(data_360.shape,avgdatatime,data_p.shape) + return data_360,avgdatatime,data_p,data_e + + def checkcase(self,data_ele): + start = data_ele[0] + end = data_ele[-1] + diff_angle = (end-start) + len_ang=len(data_ele) + print("start",start) + print("end",end) + print("number",diff_angle) + + print("len_ang",len_ang) + + aux = (data_ele<0).any(axis=0) + + #exit(1) + if diff_angle<0 and aux!=1: #Bajada + return 1 + elif diff_angle<0 and aux==1: #Bajada con angulos negativos + return 0 + elif diff_angle == 0: # This case happens when the angle reaches the max_angle if n = 2 + self.flagEraseFirstData = 1 + print("ToDO this case") + exit(1) + elif diff_angle>0: #Subida + return 0 + + def run(self, dataOut,n = None,mode=None,**kwargs): + #print("BLOCK 360 HERE WE GO MOMENTOS") + print("Block 360") + + #exit(1) + if not self.isConfig: + if n == 1: + print("*******************Min Value is 2. Setting n = 2*******************") + n = 2 + #exit(1) + print(n) + self.setup(dataOut = dataOut, n = n ,mode= mode ,**kwargs) + ####self.index = 0 + #print("comova",self.isConfig) + self.isConfig = True + ####if self.index==dataOut.azimuth.shape[0]: + #### self.index=0 + data_360, avgdatatime,data_p,data_e = self.blockOp(dataOut, dataOut.utctime) + dataOut.flagNoData = True + + if self.__dataReady: + dataOut.data_360 = data_360 # S + #print("DATA 360") + #print(dataOut.data_360) + #print("---------------------------------------------------------------------------------") + print("---------------------------DATAREADY---------------------------------------------") + #print("---------------------------------------------------------------------------------") + #print("data_360",dataOut.data_360.shape) + dataOut.data_azi = data_p + dataOut.data_ele = data_e + ###print("azi: ",dataOut.data_azi) + #print("ele: ",dataOut.data_ele) + #print("jroproc_parameters",data_p[0],data_p[-1])#,data_360.shape,avgdatatime) + dataOut.utctime = avgdatatime + + dataOut.case_flag = self.checkcase(dataOut.data_ele) + if dataOut.case_flag: #Si está de bajada empieza a plotear + print("INSIDE CASE FLAG BAJADA") + dataOut.flagNoData = False + else: + print("CASE SUBIDA") + dataOut.flagNoData = True + + #dataOut.flagNoData = False + return dataOut diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 7abc55b..32f17a4 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -216,7 +216,7 @@ class selectHeights(Operation): if (maxIndex >= self.dataOut.nHeights): maxIndex = self.dataOut.nHeights - + #print("shapeeee",self.dataOut.data.shape) #voltage if self.dataOut.flagDataAsBlock: """ @@ -1472,7 +1472,10 @@ class PulsePair(Operation): return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf, avgdatatime def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs): - + #print("hey") + #print(dataOut.data.shape) + #exit(1) + #print(self.__profIndex) if not self.isConfig: self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs) self.isConfig = True @@ -1494,7 +1497,233 @@ class PulsePair(Operation): dataOut.flagNoData = False return dataOut +class PulsePair_vRF(Operation): + ''' + Function PulsePair(Signal Power, Velocity) + The real component of Lag[0] provides Intensity Information + The imag component of Lag[1] Phase provides Velocity Information + + Configuration Parameters: + nPRF = Number of Several PRF + theta = Degree Azimuth angel Boundaries + + Input: + self.dataOut + lag[N] + Affected: + self.dataOut.spc + ''' + isConfig = False + __profIndex = 0 + __initime = None + __lastdatatime = None + __buffer = None + noise = None + __dataReady = False + n = None + __nch = 0 + __nHeis = 0 + removeDC = False + ipp = None + lambda_ = 0 + + def __init__(self,**kwargs): + Operation.__init__(self,**kwargs) + def setup(self, dataOut, n = None, removeDC=False): + ''' + n= Numero de PRF's de entrada + ''' + self.__initime = None + ####print("[INICIO]-setup del METODO PULSE PAIR") + self.__lastdatatime = 0 + self.__dataReady = False + self.__buffer = 0 + self.__profIndex = 0 + self.noise = None + self.__nch = dataOut.nChannels + self.__nHeis = dataOut.nHeights + self.removeDC = removeDC + self.lambda_ = 3.0e8/(9345.0e6) + self.ippSec = dataOut.ippSeconds + self.nCohInt = dataOut.nCohInt + ####print("IPPseconds",dataOut.ippSeconds) + ####print("ELVALOR DE n es:", n) + if n == None: + raise ValueError("n should be specified.") + + if n != None: + if n<2: + raise ValueError("n should be greater than 2") + + self.n = n + self.__nProf = n + + self.__buffer = numpy.zeros((dataOut.nChannels, + n, + dataOut.nHeights), + dtype='complex') + + def putData(self,data): + ''' + Add a profile to he __buffer and increase in one the __profiel Index + ''' + self.__buffer[:,self.__profIndex,:]= data + self.__profIndex += 1 + return + + def putDataByBlock(self,data,n): + ''' + Add a profile to he __buffer and increase in one the __profiel Index + ''' + self.__buffer[:]= data + self.__profIndex = n + return + + def pushData(self,dataOut): + ''' + Return the PULSEPAIR and the profiles used in the operation + Affected : self.__profileIndex + ''' + #----------------- Remove DC----------------------------------- + if self.removeDC==True: + mean = numpy.mean(self.__buffer,1) + tmp = mean.reshape(self.__nch,1,self.__nHeis) + dc= numpy.tile(tmp,[1,self.__nProf,1]) + self.__buffer = self.__buffer - dc + #------------------Calculo de Potencia ------------------------ + pair0 = self.__buffer*numpy.conj(self.__buffer) + pair0 = pair0.real + lag_0 = numpy.sum(pair0,1) + #-----------------Calculo de Cscp------------------------------ New + cspc_pair01 = self.__buffer[0]*self.__buffer[1] + #------------------Calculo de Ruido x canal-------------------- + self.noise = numpy.zeros(self.__nch) + for i in range(self.__nch): + daux = numpy.sort(pair0[i,:,:],axis= None) + self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt) + + self.noise = self.noise.reshape(self.__nch,1) + self.noise = numpy.tile(self.noise,[1,self.__nHeis]) + noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis) + noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1]) + #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N-- + #------------------ P= S+N ,P=lag_0/N --------------------------------- + #-------------------- Power -------------------------------------------------- + data_power = lag_0/(self.n*self.nCohInt) + #--------------------CCF------------------------------------------------------ + data_ccf =numpy.sum(cspc_pair01,axis=0)/(self.n*self.nCohInt) + #------------------ Senal -------------------------------------------------- + data_intensity = pair0 - noise_buffer + data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt) + #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt) + for i in range(self.__nch): + for j in range(self.__nHeis): + if data_intensity[i][j] < 0: + data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j])) + + #----------------- Calculo de Frecuencia y Velocidad doppler-------- + pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:]) + lag_1 = numpy.sum(pair1,1) + data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1) + data_velocity = (self.lambda_/2.0)*data_freq + + #---------------- Potencia promedio estimada de la Senal----------- + lag_0 = lag_0/self.n + S = lag_0-self.noise + + #---------------- Frecuencia Doppler promedio --------------------- + lag_1 = lag_1/(self.n-1) + R1 = numpy.abs(lag_1) + + #---------------- Calculo del SNR---------------------------------- + data_snrPP = S/self.noise + for i in range(self.__nch): + for j in range(self.__nHeis): + if data_snrPP[i][j] < 1.e-20: + data_snrPP[i][j] = 1.e-20 + + #----------------- Calculo del ancho espectral ---------------------- + L = S/R1 + L = numpy.where(L<0,1,L) + L = numpy.log(L) + tmp = numpy.sqrt(numpy.absolute(L)) + data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L) + n = self.__profIndex + + self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex') + self.__profIndex = 0 + return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,data_ccf,n + + + def pulsePairbyProfiles(self,dataOut,n): + + self.__dataReady = False + data_power = None + data_intensity = None + data_velocity = None + data_specwidth = None + data_snrPP = None + data_ccf = None + + if dataOut.flagDataAsBlock: + self.putDataByBlock(data=dataOut.data,n=n) + else: + self.putData(data=dataOut.data) + if self.__profIndex == self.n: + data_power,data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, n = self.pushData(dataOut=dataOut) + self.__dataReady = True + + return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf + + + def pulsePairOp(self, dataOut, n, datatime= None): + + if self.__initime == None: + self.__initime = datatime + data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf = self.pulsePairbyProfiles(dataOut,n) + self.__lastdatatime = datatime + + if data_power is None: + return None, None, None,None,None,None,None + + avgdatatime = self.__initime + deltatime = datatime - self.__lastdatatime + self.__initime = datatime + + return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf, avgdatatime + + def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs): + #print("hey") + #print(dataOut.data.shape) + #exit(1) + if dataOut.flagDataAsBlock: + n = dataOut.nProfileBlocks + #print(self.__profIndex) + if not self.isConfig: + self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs) + self.isConfig = True + + + data_power, data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, avgdatatime = self.pulsePairOp(dataOut, n, dataOut.utctime) + + + dataOut.flagNoData = True + + if self.__dataReady: + ###print("READY ----------------------------------") + dataOut.nCohInt *= self.n + dataOut.dataPP_POW = data_intensity # S + dataOut.dataPP_POWER = data_power # P valor que corresponde a POTENCIA MOMENTO + dataOut.dataPP_DOP = data_velocity + dataOut.dataPP_SNR = data_snrPP + dataOut.dataPP_WIDTH = data_specwidth + dataOut.dataPP_CCF = data_ccf + dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo. + dataOut.nProfiles = int(dataOut.nProfiles/n) + dataOut.utctime = avgdatatime + dataOut.flagNoData = False + return dataOut # import collections # from scipy.stats import mode