schain Commit - r1411:e3b451d0f24e · Jicamarca Repository

update y revision RHI

avaldezp -

r1411:e3b451d0f24e

parent child

Context file:

r1411:e3b451d0f24e

Collapse all files

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

schainpy/scripts/basic_proc00004.py created 644 +122 0

@@ -0,0 +1,122
	1	# Ing. AVP
	2	# 04/01/2022
	3	# ARCHIVO DE LECTURA
	4	#---- DATA RHI --- 23 DE NOVIEMBRE DEL 2021 --- 23/11/2021---
	5	#---- PEDESTAL ----------------------------------------------
	6	#------- HORA 143826 /DATA_RM/TEST_PEDESTAL/P20211123-143826 14:38-15:10
	7	#---- RADAR ----------------------------------------------
	8	#------- 14:26-15:00
	9	#------- /DATA_RM/DRONE/2MHZ_5V_ELEVACION/
	10	#------- /DATA_RM/DRONE/2MHZ_5V_ELEVACION/ch0/2021-11-23T19-00-00
	11
	12	import os, sys
	13	import datetime
	14	import time
	15	import numpy
	16	from ext_met import getfirstFilefromPath,getDatavaluefromDirFilename
	17	from schainpy.controller import Project
	18	#-----------------------------------------------------------------------------------------
	19	print("[SETUP]-RADAR METEOROLOGICO-")
	20	path_ped = "/DATA_RM/TEST_PEDESTAL/P20211123-143826"
	21	print("PATH PEDESTAL :",path_ped)
	22	path_adq = "/DATA_RM/DRONE/2MHZ_5V_ELEVACION/"
	23	print("PATH DATA :",path_adq)
	24	figpath_pp_rti = "/home/soporte/Pictures/TEST_PP_RHI"
	25	print("PATH PP RTI :",figpath_pp_rti)
	26	figpath_pp_rhi = "/home/soporte/Pictures/TEST_PP_RHI"
	27	print("PATH PP RHI :",figpath_pp_rhi)
	28	path_pp_save_int = "/DATA_RM/TEST_SAVE_PP_INT_RHI"
	29	print("PATH SAVE PP INT :",path_pp_save_int)
	30	print(" ")
	31	#-------------------------------------------------------------------------------------------
	32	print("SELECCIONAR MODO: PPI (0) O RHI (1)")
	33	mode_wr = 1
	34	if mode_wr==0:
	35	print("[ ON ] MODE PPI")
	36	list_ped = getfirstFilefromPath(path=path_ped,meta="PE",ext=".hdf5")
	37	ff_pedestal = list_ped[2]
	38	azi_vel = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="azi_vel")
	39	V = round(azi_vel[0])
	40	print("VELOCIDAD AZI :", int(numpy.mean(azi_vel)),"°/seg")
	41	else:
	42	print("[ ON ] MODE RHI")
	43	list_ped = getfirstFilefromPath(path=path_ped,meta="PE",ext=".hdf5")
	44	ff_pedestal = list_ped[2]
	45	ele_vel = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="ele_vel")
	46	V = round(ele_vel[0])
	47	V = 10.0
	48	print("VELOCIDAD ELE :", int(numpy.mean(ele_vel)),"°/seg")
	49	print(" ")
	50	#---------------------------------------------------------------------------------------
	51	print("SELECCIONAR MODO: PULSE PAIR (0) O FREQUENCY (1)")
	52	mode_proc = 0
	53	if mode_proc==0:
	54	print("[ ON ] MODE PULSEPAIR")
	55	else:
	56	print("[ ON ] MODE FREQUENCY")
	57	ipp = 60.0
	58	print("IPP(Km.) : %1.2f"%ipp)
	59	ipp_sec = (ipp1.0e3/150.0)1.0e-6
	60	print("IPP(useg.) : %1.2f"%(ipp_sec*(1.0e6)))
	61	VEL=V
	62	n= int(1/(VEL*ipp_sec))
	63	print("N° Profiles : ", n)
	64	#--------------------------------------------
	65	plot_rti = 0
	66	plot_rhi = 1
	67	integration = 1
	68	save = 0
	69	#---------------------------RANGO DE PLOTEO----------------------------------
	70	dBmin = '1'
	71	dBmax = '85'
	72	xmin = '17'
	73	xmax = '17.25'
	74	ymin = '0'
	75	ymax = '600'
	76	#----------------------------------------------------------------------------
	77	time.sleep(3)
	78	#---------------------SIGNAL CHAIN ------------------------------------
	79	desc = "USRP_WEATHER_RADAR"
	80	filename = "USRP_processing.xml"
	81	controllerObj = Project()
	82	controllerObj.setup(id = '191', name='Test_USRP', description=desc)
	83	#---------------------UNIDAD DE LECTURA--------------------------------
	84	readUnitConfObj = controllerObj.addReadUnit(datatype='DigitalRFReader',
	85	path=path_adq,
	86	startDate="2021/11/10",#today,
	87	endDate="2021/12/30",#today,
	88	startTime='17:10:25',
	89	endTime='23:59:59',
	90	delay=0,
	91	#set=0,
	92	online=0,
	93	walk=1,
	94	ippKm=ipp)
	95
	96	procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc',inputId=readUnitConfObj.getId())
	97
	98	opObj11 = procUnitConfObjA.addOperation(name='selectHeights')
	99	opObj11.addParameter(name='minIndex', value='1', format='int')
	100	# opObj11.addParameter(name='maxIndex', value='10000', format='int')
	101	opObj11.addParameter(name='maxIndex', value='400', format='int')
	102
	103	if mode_proc==0:
	104	opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other')
	105	opObj11.addParameter(name='n', value=int(n), format='int')
	106	procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
	107
	108	if integration==1:
	109	opObj11 = procUnitConfObjB.addOperation(name='PedestalInformation')
	110	opObj11.addParameter(name='path_ped', value=path_ped)
	111	opObj11.addParameter(name='t_Interval_p', value='0.01', format='float')
	112
	113	if plot_rhi==1:
	114	opObj11 = procUnitConfObjB.addOperation(name='Block360')
	115	opObj11.addParameter(name='n', value='10', format='int')
	116	opObj11.addParameter(name='mode', value=mode_proc, format='int')
	117	# este bloque funciona bien con divisores de 360 no olvidar 0 10 20 30 40 60 90 120 180
	118	opObj11= procUnitConfObjB.addOperation(name='WeatherRHIPlot',optype='other')
	119	opObj11.addParameter(name='save', value=figpath_pp_rhi)
	120	opObj11.addParameter(name='save_period', value=1)
	121
	122	controllerObj.start()

schainpy/scripts/test_wradlib_RHI.py created 644 +57 0

@@ -0,0 +1,57
	1	import numpy as np
	2	import matplotlib.pyplot as plt
	3	import wradlib as wrl
	4	import warnings
	5	# libreia nueva
	6	from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
	7	warnings.filterwarnings('ignore')
	8	# lectura de gaMIC hdf5 file
	9	filename = wrl.util.get_wradlib_data_file("/home/soporte/Downloads/2014-06-09--185000.rhi.mvol")
	10	data1, metadata = wrl.io.read_gamic_hdf5(filename)
	11	print(data1)
	12	data1 = data1['SCAN0']['ZH']['data']
	13	print(data1)
	14	print("SHAPE Data",np.array(data1).shape)
	15	r = metadata['SCAN0']['r']
	16	print("r",r)
	17	print("longitud r",len(r))
	18	th = metadata['SCAN0']['el']
	19	print("th",th)
	20	print("longitud th",len(th))
	21	az = metadata['SCAN0']['az']
	22	print("az",az)
	23	site = (metadata['VOL']['Longitude'], metadata['VOL']['Latitude'],
	24	metadata['VOL']['Height'])
	25
	26	print("Longitud,Latitud,Altura",site)
	27	ma1 = np.array(data1)
	28	'''
	29	mask_ind = np.where(data1 <= np.nanmin(data1))
	30	data1[mask_ind] = np.nan
	31	ma1 = np.ma.array(data1, mask=np.isnan(data1))
	32	'''
	33	#cgax, pm = wrl.vis.plot_rhi(ma1,r=r,th=th,rf=1e3)
	34	fig = plt.figure(figsize=(10,8))
	35	cgax, pm = wrl.vis.plot_rhi(ma1,r=r,th=th,rf=1e3,fig=fig, ax=111,proj='cg')
	36	caax = cgax.parasites[0]
	37	paax = cgax.parasites[1]
	38	cgax.set_ylim(0, 14)
	39	#caax = cgax.parasites[0]
	40	#paax = cgax.parasites[1]
	41	#cgax, pm = wrl.vis.plot_rhi(ma1, r=r, th=th, rf=1e3, fig=fig, ax=111, proj='cg')
	42	txt = plt.title('Simple RHI',y=1.05)
	43	#cbar = plt.gcf().colorbar(pm, pad=0.05, ax=paax)
	44	cbar = plt.gcf().colorbar(pm, pad=0.05)
	45	cbar.set_label('reflectivity [dBZ]')
	46	caax.set_xlabel('x_range [km]')
	47	caax.set_ylabel('y_range [km]')
	48	plt.text(1.0, 1.05, 'azimuth', transform=caax.transAxes, va='bottom',ha='right')
	49	gh = cgax.get_grid_helper()
	50
	51	# set theta to some nice values
	52	locs = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14.,
	53	15., 16., 17., 18., 20., 22., 25., 30., 35., 40., 50., 60., 70., 80., 90.]
	54	gh.grid_finder.grid_locator1 = FixedLocator(locs)
	55	gh.grid_finder.tick_formatter1 = DictFormatter(dict([(i, r"${0:.0f}^\circ$".format(i)) for i in locs]))
	56
	57	plt.show()

schainpy/model/graphics/jroplot_parameters.py +77 -57

             import os
             import datetime
             import numpy
             from schainpy.model.graphics.jroplot_base import Plot, plt
             from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
             from schainpy.utils import log
             # libreria wradlib
             import wradlib as wrl
             EARTH_RADIUS = 6.3710e3
             def ll2xy(lat1, lon1, lat2, lon2):
                 p = 0.017453292519943295
                 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
                     numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
                 r = 12742 * numpy.arcsin(numpy.sqrt(a))
                 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
                                       * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
                 theta = -theta + numpy.pi/2
                 return r*numpy.cos(theta), r*numpy.sin(theta)
             def km2deg(km):
                 '''
                 Convert distance in km to degrees
                 '''
                 return numpy.rad2deg(km/EARTH_RADIUS)
             class SpectralMomentsPlot(SpectraPlot):
                 '''
                 Plot for Spectral Moments
                 '''
                 CODE = 'spc_moments'
                 # colormap = 'jet'
                 # plot_type = 'pcolor'
             class DobleGaussianPlot(SpectraPlot):
                 '''
                 Plot for Double Gaussian Plot
                 '''
                 CODE = 'gaussian_fit'
                 # colormap = 'jet'
                 # plot_type = 'pcolor'
             class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
                 '''
                 Plot SpectraCut with Double Gaussian Fit
                 '''
                 CODE = 'cut_gaussian_fit'
             class SnrPlot(RTIPlot):
                 '''
                 Plot for SNR Data
                 '''
                 CODE = 'snr'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'snr': 10*numpy.log10(dataOut.data_snr)
                     }
                     return data, {}
             class DopplerPlot(RTIPlot):
                 '''
                 Plot for DOPPLER Data (1st moment)
                 '''
                 CODE = 'dop'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'dop': 10*numpy.log10(dataOut.data_dop)
                     }
                     return data, {}
             class PowerPlot(RTIPlot):
                 '''
                 Plot for Power Data (0 moment)
                 '''
                 CODE = 'pow'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
                     }
                     return data, {}
             class SpectralWidthPlot(RTIPlot):
                 '''
                 Plot for Spectral Width Data (2nd moment)
                 '''
                 CODE = 'width'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'width': dataOut.data_width
                     }
                     return data, {}
             class SkyMapPlot(Plot):
                 '''
                 Plot for meteors detection data
                 '''
                 CODE = 'param'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.width = 7.2
                     self.height = 7.2
                     self.nplots = 1
                     self.xlabel = 'Zonal Zenith Angle (deg)'
                     self.ylabel = 'Meridional Zenith Angle (deg)'
                     self.polar = True
                     self.ymin = -180
                     self.ymax = 180
                     self.colorbar = False
                 def plot(self):
                     arrayParameters = numpy.concatenate(self.data['param'])
                     error = arrayParameters[:, -1]
                     indValid = numpy.where(error == 0)[0]
                     finalMeteor = arrayParameters[indValid, :]
                     finalAzimuth = finalMeteor[:, 3]
                     finalZenith = finalMeteor[:, 4]
                     x = finalAzimuth * numpy.pi / 180
                     y = finalZenith
                     ax = self.axes[0]
                     if ax.firsttime:
                         ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
                     else:
                         ax.plot.set_data(x, y)
                     dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
                     dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
                     title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
                                                                                                  dt2,
                                                                                                  len(x))
                     self.titles[0] = title
             class GenericRTIPlot(Plot):
                 '''
                 Plot for data_xxxx object
                 '''
                 CODE = 'param'
                 colormap = 'viridis'
                 plot_type = 'pcolorbuffer'
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     self.nrows = self.data.shape('param')[0]
                     self.nplots = self.nrows
                     self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
                     if not self.xlabel:
                         self.xlabel = 'Time'
                     self.ylabel = 'Range [km]'
                     if not self.titles:
                         self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
                 def update(self, dataOut):
                     data = {
                         'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
                     }
                     meta = {}
                     return data, meta
                 def plot(self):
                     # self.data.normalize_heights()
                     self.x = self.data.times
                     self.y = self.data.yrange
                     self.z = self.data['param']
                     self.z = 10*numpy.log10(self.z)
                     self.z = numpy.ma.masked_invalid(self.z)
                     if self.decimation is None:
                         x, y, z = self.fill_gaps(self.x, self.y, self.z)
                     else:
                         x, y, z = self.fill_gaps(*self.decimate())
                     for n, ax in enumerate(self.axes):
                         self.zmax = self.zmax if self.zmax is not None else numpy.max(
                             self.z[n])
                         self.zmin = self.zmin if self.zmin is not None else numpy.min(
                             self.z[n])
                         if ax.firsttime:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
                                                    vmin=self.zmin,
                                                    vmax=self.zmax,
                                                    cmap=self.cmaps[n]
                                                    )
                         else:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.collections.remove(ax.collections[0])
                             ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
                                                    vmin=self.zmin,
                                                    vmax=self.zmax,
                                                    cmap=self.cmaps[n]
                                                    )
             class PolarMapPlot(Plot):
                 '''
                 Plot for weather radar
                 '''
                 CODE = 'param'
                 colormap = 'seismic'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.width = 9
                     self.height = 8
                     self.mode = self.data.meta['mode']
                     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))
                     if self.mode == 'E':
                         self.xlabel = 'Longitude'
                         self.ylabel = 'Latitude'
                     else:
                         self.xlabel = 'Range (km)'
                         self.ylabel = 'Height (km)'
                     self.bgcolor = 'white'
                     self.cb_labels = self.data.meta['units']
                     self.lat = self.data.meta['latitude']
                     self.lon = self.data.meta['longitude']
                     self.xmin, self.xmax = float(
                         km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
                     self.ymin, self.ymax = float(
                         km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
                     #  self.polar = True
                 def plot(self):
                     for n, ax in enumerate(self.axes):
                         data = self.data['param'][self.channels[n]]
                         zeniths = numpy.linspace(
 , self.data.meta['max_range'], data.shape[1])
                         if self.mode == 'E':
                             azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
                             r, theta = numpy.meshgrid(zeniths, azimuths)
                             x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
                                 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
                             x = km2deg(x) + self.lon
                             y = km2deg(y) + self.lat
                         else:
                             azimuths = numpy.radians(self.data.yrange)
                             r, theta = numpy.meshgrid(zeniths, azimuths)
                             x, y = r*numpy.cos(theta), r*numpy.sin(theta)
                         self.y = zeniths
                         if ax.firsttime:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.plt = ax.pcolormesh(  # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 vmin=self.zmin,
                                 vmax=self.zmax,
                                 cmap=self.cmaps[n])
                         else:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.collections.remove(ax.collections[0])
                             ax.plt = ax.pcolormesh(  # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 vmin=self.zmin,
                                 vmax=self.zmax,
                                 cmap=self.cmaps[n])
                         if self.mode == 'A':
                             continue
                         # plot district names
                         f = open('/data/workspace/schain_scripts/distrito.csv')
                         for line in f:
                             label, lon, lat = [s.strip() for s in line.split(',') if s]
                             lat = float(lat)
                             lon = float(lon)
                             # ax.plot(lon, lat, '.b', ms=2)
                             ax.text(lon, lat, label.decode('utf8'), ha='center',
                                     va='bottom', size='8', color='black')
                         # plot limites
                         limites = []
                         tmp = []
                         for line in open('/data/workspace/schain_scripts/lima.csv'):
                             if '#' in line:
                                 if tmp:
                                     limites.append(tmp)
                                 tmp = []
                                 continue
                             values = line.strip().split(',')
                             tmp.append((float(values[0]), float(values[1])))
                         for points in limites:
                             ax.add_patch(
                                 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
                         # plot Cuencas
                         for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
                             f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
                             values = [line.strip().split(',') for line in f]
                             points = [(float(s[0]), float(s[1])) for s in values]
                             ax.add_patch(Polygon(points, ec='b', fc='none'))
                         # plot grid
                         for r in (15, 30, 45, 60):
                             ax.add_artist(plt.Circle((self.lon, self.lat),
                                                      km2deg(r), color='0.6', fill=False, lw=0.2))
                             ax.text(
                                 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
                                 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
                                 '{}km'.format(r),
                                 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
                     if self.mode == 'E':
                         title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
                         label = 'E{:02d}'.format(int(self.data.meta['elevation']))
                     else:
                         title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
                         label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
                     self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
                     self.titles = ['{} {}'.format(
                         self.data.parameters[x], title) for x in self.channels]
             class WeatherPlot(Plot):
                 CODE = 'weather'
                 plot_name = 'weather'
                 plot_type = 'ppistyle'
                 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
-                    ####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
                 def get2List(self,angulos):
                     list1=[]
                     list2=[]
                     for i in reversed(range(len(angulos))):
                         diff_ = angulos[i]-angulos[i-1]
                         if diff_ >1.5:
                             list1.append(i-1)
                             list2.append(diff_)
                     return list(reversed(list1)),list(reversed(list2))
                 def fixData360(self,list_,ang_):
                     if list_[0]==-1:
                         vec = numpy.where(ang_<ang_[0])
                         ang_[vec] = ang_[vec]+360
                         return ang_
                     return ang_
                 def fixData360HL(self,angulos):
                     vec = numpy.where(angulos>=360)
                     angulos[vec]=angulos[vec]-360
                     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_):
                     size  = len(ang_)
                     size2 = 0
                     for i in range(len(list2_)):
                         size2=size2+round(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(list2_[value])-1):
                                 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):
                     l1,l2 = self.get2List(angulos)
                     if len(l1)>0:
                         angulos2 = self.fixData360(list_=l1,ang_=angulos)
                         l1,l2 = self.get2List(angulos2)
                         ang1_,ang2_ = self.fixDataComp(ang_=angulos2,list1_=l1,list2_=l2)
                         ang1_ = self.fixData360HL(ang1_)
                         ang2_ = self.fixData360HL(ang2_)
                     else:
                         ang1_= angulos
                         ang2_= angulos
                     return ang1_,ang2_
                 def analizeDATA(self,data_azi):
                     list1 = []
                     list2 = []
                     dat   = data_azi
                     for i in reversed(range(1,len(dat))):
                         if dat[i]>dat[i-1]:
                             diff = int(dat[i])-int(dat[i-1])
                         else:
                             diff = 360+int(dat[i])-int(dat[i-1])
                         if diff > 1:
                             list1.append(i-1)
                             list2.append(diff-1)
                     return list1,list2
                 def fixDATANEW(self,data_azi,data_weather):
                     list1,list2 = self.analizeDATA(data_azi)
                     if len(list1)== 0:
                         return data_azi,data_weather
                     else:
                         resize = 0
                         for i in range(len(list2)):
                             resize=  resize + list2[i]
                         new_data_azi = numpy.resize(data_azi,resize)
                         new_data_weather= numpy.resize(date_weather,resize)
                         for i in range(len(list2)):
                             j=0
                             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):
                     data=data_azi
                     for i in range(len(data)):
                         if numpy.isnan(data[i]):
                            data[i]=data[i-1]+1
                     return data
                 def replaceNAN(self,data_weather,data_azi,val):
                     data= data_azi
                     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,:]
                                 sc=c+1
                     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,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
                         azi_vacia = numpy.linspace(start,end,int(n))
                         azi_vacia = numpy.where(azi_vacia>360,azi_vacia-360,azi_vacia)
                         data_azi  = numpy.hstack((data_azi_new,azi_vacia))
                         # RADAR
                         val_mean         = numpy.mean(data_weather[:,-1])
                         self.val_mean    = val_mean
                         data_weather_cmp = numpy.ones([(360-data_weather.shape[0]),data_weather.shape[1]])*val_mean
                         data_weather     = self.replaceNAN(data_weather=data_weather,data_azi=data_azi_old,val=self.val_mean)
                         data_weather     = numpy.vstack((data_weather,data_weather_cmp))
                     else:
                         # azimuth
                         flag=0
                         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 <start_azi:
                             end  = end +360
-                        ####print("start",start)
-                        ####print("end",end)
-                        #### AQUI SERA LA MAGIA
                         pos_ini = int((start-start_azi)/res)
                         len_azi = len(data_azi)
                         if (360-pos_ini)<len_azi:
                             if pos_ini+1==360:
                                 pos_ini=0
                             else:
                                 flag=1
                                 dif= 360-pos_ini
                                 comp= len_azi-dif
                         #-----------------
-                        ####print(pos_ini)
-                        ####print(len_azi)
-                        ####print("shape",self.res_azi.shape)
                         if flag==0:
                             # AZIMUTH
                             self.res_azi[pos_ini:pos_ini+len_azi] = data_azi
                             # RADAR
                             self.res_weather[pos_ini:pos_ini+len_azi,:] = data_weather
                         else:
                             # AZIMUTH
                             self.res_azi[pos_ini:pos_ini+dif] = data_azi[0:dif]
                             self.res_azi[0:comp]              = data_azi[dif:]
                             # RADAR
                             self.res_weather[pos_ini:pos_ini+dif,:] = data_weather[0:dif,:]
                             self.res_weather[0:comp,:]              = data_weather[dif:,:]
                             flag=0
                         data_azi                                    = self.res_azi
                         data_weather                                = self.res_weather
                     return data_weather,data_azi
                 def plot(self):
-                    #print("--------------------------------------",self.ini,"-----------------------------------")
-                    #numpy.set_printoptions(suppress=True)
-                    ####print("times: ",self.data.times)
                     thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]).strftime('%Y-%m-%d %H:%M:%S')
-                    #print("times: ",thisDatetime)
                     data         = self.data[-1]
-                    ####ALTURA altura_tmp_h
-                    ###print("Y RANGES",self.data.yrange,len(self.data.yrange))
-                    ###altura_h     = (data['weather'].shape[1])/10.0
-                    ###stoprange = float(altura_h*0.3)#stoprange = float(33*1.5) por ahora 400
-                    ###rangestep = float(0.03)
-                    ###r      = numpy.arange(0, stoprange, rangestep)
-                    ###print("r",r,len(r))
-                    #-----------------------------update----------------------
                     r            = self.data.yrange
                     delta_height = r[1]-r[0]
-                    #print("1",r)
                     r_mask       = numpy.where(r>=0)[0]
                     r            = numpy.arange(len(r_mask))*delta_height
-                    #print("2",r)
                     self.y       = 2*r
-                    ######self.y = self.data.yrange
                     # RADAR
                     #data_weather = data['weather']
                     # PEDESTAL
                     #data_azi  = data['azi']
                     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)
                     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###############################
                     #################    PLOTEO           ###################
-                    ##########################################################
                     for i,ax in enumerate(self.axes):
                         if ax.firsttime:
                             plt.clf()
                             cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=8, vmax=35)
                         else:
                             plt.clf()
                             cgax, pm = wrl.vis.plot_ppi(self.res_weather,r=r,az=self.res_azi,fig=self.figures[0], proj='cg', vmin=8, vmax=35)
                     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, '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

schainpy/model/proc/jroproc_parameters.py 0 0

	1	NO CONTENT: modified file			NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff

schainpy/scripts/basic_proc00003.py +1 -1

             # Ing. AVP
             # 04/01/2022
             # ARCHIVO DE LECTURA
             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/P20211110-171003"
             print("PATH PEDESTAL :",path_ped)
             path_adq  = "/DATA_RM/10"
             print("PATH DATA     :",path_adq)
             figpath_pp_rti  = "/home/soporte/Pictures/TEST_PP_RTI"
             print("PATH PP RTI   :",figpath_pp_rti)
             figpath_pp_ppi  = "/home/soporte/Pictures/TEST_PP_PPI"
             print("PATH PP PPI   :",figpath_pp_ppi)
             path_pp_save_int  = "/DATA_RM/TEST_SAVE_PP_INT"
             print("PATH SAVE PP INT   :",path_pp_save_int)
             print(" ")
             #-------------------------------------------------------------------------------------------
             print("SELECCIONAR MODO: PPI        (0) O RHI       (1)")
             mode_wr     = 0
             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]
                 V            = round(ele_vel[0])
                 ele_vel      = getDatavaluefromDirFilename(path=path_ped,file=ff_pedestal,value="ele_vel")
                 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_ppi    = 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())
+                # REVISAR EL test_sim00013.py
                 if plot_rti==1:
                     opObj11 = procUnitConfObjB.addOperation(name='GenericRTIPlot',optype='external')
                     opObj11.addParameter(name='attr_data', value='dataPP_POW')
                     opObj11.addParameter(name='colormap', value='jet')
                     opObj11.addParameter(name='xmin', value=xmin)
                     opObj11.addParameter(name='xmax', value=xmax)
                     opObj11.addParameter(name='zmin', value=dBmin)
                     opObj11.addParameter(name='zmax', value=dBmax)
                     opObj11.addParameter(name='save', value=figpath_pp_rti)
                     opObj11.addParameter(name='showprofile', value=0)
                     opObj11.addParameter(name='save_period', value=50)
                 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_ppi==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='WeatherPlot',optype='other')
                     opObj11.addParameter(name='save', value=figpath_pp_ppi)
                     opObj11.addParameter(name='save_period', value=1)
                 if save==1:
                     opObj10 = procUnitConfObjB.addOperation(name='HDFWriter')
                     opObj10.addParameter(name='path',value=path_pp_save_int)
                     opObj10.addParameter(name='mode',value="weather")
                     opObj10.addParameter(name='blocksPerFile',value='360',format='int')
                     opObj10.addParameter(name='metadataList',value='utctimeInit,timeZone,paramInterval,profileIndex,channelList,heightList,flagDataAsBlock',format='list')
                     opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,azimuth,elevation,utctime',format='list')#,format='list'
             controllerObj.start()

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages