diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index a87ea94..e0840e2 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -399,15 +399,10 @@ class WeatherPlot(Plot): meta = {} if hasattr(dataOut, 'dataPP_POWER'): factor = 1 - if hasattr(dataOut, 'nFFTPoints'): factor = dataOut.normFactor - - ####print("factor",factor) data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor)) - ####print("weather",data['weather']) data['azi'] = dataOut.data_azi - data['ele'] = dataOut.data_ele return data, meta @@ -428,7 +423,6 @@ class WeatherPlot(Plot): return ang_ return ang_ - def fixData360HL(self,angulos): vec = numpy.where(angulos>=360) angulos[vec]=angulos[vec]-360 @@ -466,7 +460,6 @@ class WeatherPlot(Plot): c=c+1 return ang_new,ang_new2 - def globalCheckPED(self,angulos): l1,l2 = self.get2List(angulos) if len(l1)>0: @@ -476,7 +469,6 @@ class WeatherPlot(Plot): ang1_,ang2_ = self.fixDataComp(ang_=angulos2,list1_=l1,list2_=l2) ang1_ = self.fixData360HL(ang1_) ang2_ = self.fixData360HL(ang2_) - else: ang1_= angulos ang2_= angulos @@ -512,7 +504,6 @@ class WeatherPlot(Plot): position=list1[i]+1 for j in range(list2[i]): new_data_azi[position+j]=new_data_azi[position+j-1]+1 - return new_data_azi def fixDATA(self,data_azi): @@ -542,17 +533,14 @@ class WeatherPlot(Plot): def const_ploteo(self,data_weather,data_azi,step,res): if self.ini==0: - #------- AZIMUTH + #------- n = (360/res)-len(data_azi) #--------------------- new ------------------------- - ####data_azi_old = data_azi data_azi_new ,data_azi_old= self.globalCheckPED(data_azi) #------------------------ - ####data_azi_new = self.fixDATA(data_azi) - #ata_azi_new = self.fixDATANEW(data_azi) start = data_azi_new[-1] + res end = data_azi_new[0] - res - ##### new + #------ new self.last_data_azi = end if start>end: end = end + 360 @@ -571,30 +559,16 @@ class WeatherPlot(Plot): start_azi = self.res_azi[0] #-----------new------------ data_azi ,data_azi_old= self.globalCheckPED(data_azi) - print("---------------------------------------------------") - print("data_azi",data_azi) - print("data_azi_old",data_azi_old) data_weather = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean) #-------------------------- - ####data_azi_old = data_azi - ### weather ### - ####data_weather = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean) - - ####if numpy.isnan(data_azi[0]): - #### data_azi[0]=self.last_data_azi+1 - ####data_azi = self.fixDATA(data_azi) start = data_azi[0] end = data_azi[-1] self.last_data_azi= end - ####print("start",start) - ####print("end",end) if start< start_azi: start = start +360 if end =0)[0] - r = numpy.arange(len(r_mask))*delta_height - #print("2",r) - self.y = 2*r - ######self.y = self.data.yrange + r_mask = numpy.where(r>=0)[0] + r = numpy.arange(len(r_mask))*delta_height + self.y = 2*r # RADAR #data_weather = data['weather'] # PEDESTAL #data_azi = data['azi'] - res = 1 + res = 1 # STEP - step = (360/(res*data['weather'].shape[0])) - #print("shape wr_data", wr_data.shape) - #print("shape wr_azi",wr_azi.shape) - #print("step",step) - ####print("Time---->",self.data.times[-1],thisDatetime) - #print("alturas", len(self.y))numpy.where(r>=0) + step = (360/(res*data['weather'].shape[0])) + self.res_weather, self.res_azi = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_azi=data['azi'],step=step,res=res) - #numpy.set_printoptions(suppress=True) - #print("resultado",self.res_azi) - self.res_ele =numpy.mean(data['ele']) - ###########################/DATA_RM/10_tmp/ch0############################### + self.res_ele = numpy.mean(data['ele']) ################# PLOTEO ################### - ########################################################## for i,ax in enumerate(self.axes): if ax.firsttime: @@ -685,3 +632,76 @@ class WeatherPlot(Plot): plt.text(1.0, 1.05, 'Azimuth '+str(thisDatetime)+" Step "+str(self.ini)+ " Elev: "+str(round(self.res_ele,2)), transform=caax.transAxes, va='bottom',ha='right') self.ini= self.ini+1 + + +class WeatherRHIPlot(Plot): + CODE = 'weather' + plot_name = 'weather' + plot_type = 'rhistyle' + buffering = False + + def setup(self): + self.ncols = 1 + self.nrows = 1 + self.nplots= 1 + self.ylabel= 'Range [Km]' + self.titles= ['Weather'] + self.colorbar=False + self.width =8 + self.height =8 + self.ini =0 + self.len_azi =0 + self.buffer_ini = None + self.buffer_azi = 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_azi = 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 + data['weather'] = 10*numpy.log10(dataOut.data_360[1]/(factor)) + data['azi'] = dataOut.data_azi + data['ele'] = dataOut.data_ele + return data, meta + + 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] + r = numpy.arange(len(r_mask))*delta_height + self.y = 2*r + ###self.res_weather, self.res_ele = self.const_ploteo(data_weather=data['weather'][:,r_mask],data_azi=data['ele'],step=step,res=res) + ###self.res_azi = numpy.mean(data['azi']) + #------------- + # 90 angulos en el axis 0 + # 1000 step en el axis 1 + self.res_weather = numpy.ones([90,1000]) + r = numpy.linspace(0,1999,1000) + self.res_ele = numpy.arange(0,90) + self.res_azi = 240 + #------------- + for i,ax in enumerate(self.axes): + if ax.firsttime: + plt.clf() + cgax, pm = wrl.vis.plot_rhi(self.res_weather,r=r,th=self.res_ele,fig=self.figures[0], proj='cg') + else: + plt.clf() + cgax, pm = wrl.vis.plot_rhi(self.res_weather,r=r,th=self.res_ele,fig=self.figures[0], proj='cg') + 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') + + self.ini= self.ini+1 diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index fb05f1c..296e2c8 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -3944,12 +3944,13 @@ class WeatherRadar(Operation): def setMoments(self,dataOut,i): - type = dataOut.inputUnit - nCh = dataOut.nChannels - nHeis= dataOut.nHeights + type = dataOut.inputUnit + nCh = dataOut.nChannels + nHeis = dataOut.nHeights data_param = numpy.zeros((nCh,4,nHeis)) if type == "Voltage": - data_param[:,0,:] = dataOut.dataPP_POW/(dataOut.nCohInt**2) + factor = dataOut.normFactor + data_param[:,0,:] = dataOut.dataPP_POW/(factor) data_param[:,1,:] = dataOut.dataPP_DOP data_param[:,2,:] = dataOut.dataPP_WIDTH data_param[:,3,:] = dataOut.dataPP_SNR @@ -3957,7 +3958,6 @@ class WeatherRadar(Operation): data_param[:,0,:] = dataOut.data_POW data_param[:,1,:] = dataOut.data_DOP data_param[:,2,:] = dataOut.data_WIDTH - def setMoments(self,dataOut,i): data_param[:,3,:] = dataOut.data_SNR return data_param[:,i,:] diff --git a/schainpy/scripts/basic_proc00003.py b/schainpy/scripts/basic_proc00003.py index 4ca3942..856561b 100644 --- a/schainpy/scripts/basic_proc00003.py +++ b/schainpy/scripts/basic_proc00003.py @@ -95,7 +95,7 @@ if mode_proc==0: opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other') opObj11.addParameter(name='n', value=int(n), format='int') procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId()) - + # REVISAR EL test_sim00013.py if plot_rti==1: opObj11 = procUnitConfObjB.addOperation(name='GenericRTIPlot',optype='external') opObj11.addParameter(name='attr_data', value='dataPP_POW') diff --git a/schainpy/scripts/basic_proc00004.py b/schainpy/scripts/basic_proc00004.py new file mode 100644 index 0000000..f993d61 --- /dev/null +++ b/schainpy/scripts/basic_proc00004.py @@ -0,0 +1,122 @@ +# Ing. AVP +# 04/01/2022 +# ARCHIVO DE LECTURA +#---- DATA RHI --- 23 DE NOVIEMBRE DEL 2021 --- 23/11/2021--- +#---- PEDESTAL ---------------------------------------------- +#------- HORA 143826 /DATA_RM/TEST_PEDESTAL/P20211123-143826 14:38-15:10 +#---- RADAR ---------------------------------------------- +#------- 14:26-15:00 +#------- /DATA_RM/DRONE/2MHZ_5V_ELEVACION/ +#------- /DATA_RM/DRONE/2MHZ_5V_ELEVACION/ch0/2021-11-23T19-00-00 + +import os, sys +import datetime +import time +import numpy +from ext_met import getfirstFilefromPath,getDatavaluefromDirFilename +from schainpy.controller import Project +#----------------------------------------------------------------------------------------- +print("[SETUP]-RADAR METEOROLOGICO-") +path_ped = "/DATA_RM/TEST_PEDESTAL/P20211123-143826" +print("PATH PEDESTAL :",path_ped) +path_adq = "/DATA_RM/DRONE/2MHZ_5V_ELEVACION/" +print("PATH DATA :",path_adq) +figpath_pp_rti = "/home/soporte/Pictures/TEST_PP_RHI" +print("PATH PP RTI :",figpath_pp_rti) +figpath_pp_rhi = "/home/soporte/Pictures/TEST_PP_RHI" +print("PATH PP RHI :",figpath_pp_rhi) +path_pp_save_int = "/DATA_RM/TEST_SAVE_PP_INT_RHI" +print("PATH SAVE PP INT :",path_pp_save_int) +print(" ") +#------------------------------------------------------------------------------------------- +print("SELECCIONAR MODO: PPI (0) O RHI (1)") +mode_wr = 1 +if mode_wr==0: + print("[ ON ] MODE PPI") + list_ped = getfirstFilefromPath(path=path_ped,meta="PE",ext=".hdf5") + ff_pedestal = list_ped[2] + azi_vel = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="azi_vel") + V = round(azi_vel[0]) + print("VELOCIDAD AZI :", int(numpy.mean(azi_vel)),"°/seg") +else: + print("[ ON ] MODE RHI") + list_ped = getfirstFilefromPath(path=path_ped,meta="PE",ext=".hdf5") + ff_pedestal = list_ped[2] + ele_vel = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="ele_vel") + V = round(ele_vel[0]) + V = 10.0 + print("VELOCIDAD ELE :", int(numpy.mean(ele_vel)),"°/seg") +print(" ") +#--------------------------------------------------------------------------------------- +print("SELECCIONAR MODO: PULSE PAIR (0) O FREQUENCY (1)") +mode_proc = 0 +if mode_proc==0: + print("[ ON ] MODE PULSEPAIR") +else: + print("[ ON ] MODE FREQUENCY") +ipp = 60.0 +print("IPP(Km.) : %1.2f"%ipp) +ipp_sec = (ipp*1.0e3/150.0)*1.0e-6 +print("IPP(useg.) : %1.2f"%(ipp_sec*(1.0e6))) +VEL=V +n= int(1/(VEL*ipp_sec)) +print("N° Profiles : ", n) +#-------------------------------------------- +plot_rti = 0 +plot_rhi = 1 +integration = 1 +save = 0 +#---------------------------RANGO DE PLOTEO---------------------------------- +dBmin = '1' +dBmax = '85' +xmin = '17' +xmax = '17.25' +ymin = '0' +ymax = '600' +#---------------------------------------------------------------------------- +time.sleep(3) +#---------------------SIGNAL CHAIN ------------------------------------ +desc = "USRP_WEATHER_RADAR" +filename = "USRP_processing.xml" +controllerObj = Project() +controllerObj.setup(id = '191', name='Test_USRP', description=desc) +#---------------------UNIDAD DE LECTURA-------------------------------- +readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader', + path=path_adq, + startDate="2021/11/10",#today, + endDate="2021/12/30",#today, + startTime='17:10:25', + endTime='23:59:59', + delay=0, + #set=0, + online=0, + walk=1, + ippKm=ipp) + +procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc',inputId=readUnitConfObj.getId()) + +opObj11 = procUnitConfObjA.addOperation(name='selectHeights') +opObj11.addParameter(name='minIndex', value='1', format='int') +# opObj11.addParameter(name='maxIndex', value='10000', format='int') +opObj11.addParameter(name='maxIndex', value='400', format='int') + +if mode_proc==0: + opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other') + opObj11.addParameter(name='n', value=int(n), format='int') + procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId()) + + if integration==1: + opObj11 = procUnitConfObjB.addOperation(name='PedestalInformation') + opObj11.addParameter(name='path_ped', value=path_ped) + opObj11.addParameter(name='t_Interval_p', value='0.01', format='float') + + if plot_rhi==1: + opObj11 = procUnitConfObjB.addOperation(name='Block360') + opObj11.addParameter(name='n', value='10', format='int') + opObj11.addParameter(name='mode', value=mode_proc, format='int') + # este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180 + opObj11= procUnitConfObjB.addOperation(name='WeatherRHIPlot',optype='other') + opObj11.addParameter(name='save', value=figpath_pp_rhi) + opObj11.addParameter(name='save_period', value=1) + +controllerObj.start() diff --git a/schainpy/scripts/test_wradlib_RHI.py b/schainpy/scripts/test_wradlib_RHI.py new file mode 100644 index 0000000..6739943 --- /dev/null +++ b/schainpy/scripts/test_wradlib_RHI.py @@ -0,0 +1,57 @@ +import numpy as np +import matplotlib.pyplot as plt +import wradlib as wrl +import warnings +# libreia nueva +from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter +warnings.filterwarnings('ignore') +# lectura de gaMIC hdf5 file +filename = wrl.util.get_wradlib_data_file("/home/soporte/Downloads/2014-06-09--185000.rhi.mvol") +data1, metadata = wrl.io.read_gamic_hdf5(filename) +print(data1) +data1 = data1['SCAN0']['ZH']['data'] +print(data1) +print("SHAPE Data",np.array(data1).shape) +r = metadata['SCAN0']['r'] +print("r",r) +print("longitud r",len(r)) +th = metadata['SCAN0']['el'] +print("th",th) +print("longitud th",len(th)) +az = metadata['SCAN0']['az'] +print("az",az) +site = (metadata['VOL']['Longitude'], metadata['VOL']['Latitude'], + metadata['VOL']['Height']) + +print("Longitud,Latitud,Altura",site) +ma1 = np.array(data1) +''' +mask_ind = np.where(data1 <= np.nanmin(data1)) +data1[mask_ind] = np.nan +ma1 = np.ma.array(data1, mask=np.isnan(data1)) +''' +#cgax, pm = wrl.vis.plot_rhi(ma1,r=r,th=th,rf=1e3) +fig = plt.figure(figsize=(10,8)) +cgax, pm = wrl.vis.plot_rhi(ma1,r=r,th=th,rf=1e3,fig=fig, ax=111,proj='cg') +caax = cgax.parasites[0] +paax = cgax.parasites[1] +cgax.set_ylim(0, 14) +#caax = cgax.parasites[0] +#paax = cgax.parasites[1] +#cgax, pm = wrl.vis.plot_rhi(ma1, r=r, th=th, rf=1e3, fig=fig, ax=111, proj='cg') +txt = plt.title('Simple RHI',y=1.05) +#cbar = plt.gcf().colorbar(pm, pad=0.05, ax=paax) +cbar = plt.gcf().colorbar(pm, pad=0.05) +cbar.set_label('reflectivity [dBZ]') +caax.set_xlabel('x_range [km]') +caax.set_ylabel('y_range [km]') +plt.text(1.0, 1.05, 'azimuth', transform=caax.transAxes, va='bottom',ha='right') +gh = cgax.get_grid_helper() + +# set theta to some nice values +locs = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., + 15., 16., 17., 18., 20., 22., 25., 30., 35., 40., 50., 60., 70., 80., 90.] +gh.grid_finder.grid_locator1 = FixedLocator(locs) +gh.grid_finder.tick_formatter1 = DictFormatter(dict([(i, r"${0:.0f}^\circ$".format(i)) for i in locs])) + +plt.show()