schain-jc Commit - r513:19f921674eb5

First Spectral Fitting and EW Drifts operative module inside Signal Chain TRUNK

Julio Valdez -

r513:19f921674eb5

parent child

Context file:

r513:19f921674eb5

Collapse all files

schainpy/test/EWDrifts_estimation01.py created 644 +135 0

@@ -0,0 +1,135
	1	# DIAS 19 Y 20 FEB 2014
	2	# Comprobacion de Resultados DBS con SA
	3
	4	import os, sys
	5
	6	path = os.path.split(os.getcwd())[0]
	7	sys.path.append(path)
	8
	9	from controller import *
	10
	11	desc = "DBS Experiment Test"
	12	filename = "DBStest.xml"
	13
	14	controllerObj = Project()
	15
	16	controllerObj.setup(id = '191', name='test01', description=desc)
	17
	18	#Experimentos
	19
	20	path = '/host/Jicamarca/EW_Drifts/d2012248'
	21	pathFigure = '/home/propietario/workspace/Graficos/drifts'
	22
	23
	24	path = "/home/soporte/Data/drifts"
	25	pathFigure = '/home/soporte/workspace/Graficos/drifts/prueba'
	26
	27	xmin = 11.75
	28	xmax = 14.75
	29	#------------------------------------------------------------------------------------------------
	30	readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
	31	path=path,
	32	startDate='2012/01/01',
	33	endDate='2012/12/31',
	34	startTime='00:00:00',
	35	endTime='23:59:59',
	36	online=0,
	37	walk=1)
	38
	39	opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
	40
	41	#--------------------------------------------------------------------------------------------------
	42
	43	procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
	44
	45	opObj11 = procUnitConfObj0.addOperation(name='ProfileSelector', optype='other')
	46	opObj11.addParameter(name='profileRangeList', value='0,127', format='intlist')
	47
	48	opObj11 = procUnitConfObj0.addOperation(name='filterByHeights')
	49	opObj11.addParameter(name='window', value='3', format='int')
	50
	51	opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
	52	# opObj11.addParameter(name='code', value='1,-1', format='floatlist')
	53	# opObj11.addParameter(name='nCode', value='2', format='int')
	54	# opObj11.addParameter(name='nBaud', value='1', format='int')
	55
	56	procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
	57	procUnitConfObj1.addParameter(name='nFFTPoints', value='128', format='int')
	58	procUnitConfObj1.addParameter(name='nProfiles', value='128', format='int')
	59	procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(2,3)', format='pairsList')#,(2,3)
	60
	61	opObj11 = procUnitConfObj1.addOperation(name='selectHeights')
	62	# # opObj11.addParameter(name='minHei', value='320.0', format='float')
	63	# # opObj11.addParameter(name='maxHei', value='350.0', format='float')
	64	opObj11.addParameter(name='minHei', value='200.0', format='float')
	65	opObj11.addParameter(name='maxHei', value='600.0', format='float')
	66
	67	opObj11 = procUnitConfObj1.addOperation(name='selectChannels')
	68	opObj11.addParameter(name='channelList', value='0,1,2,3', format='intlist')
	69
	70	opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
	71	opObj11.addParameter(name='timeInterval', value='300.0', format='float')
	72
	73	opObj13 = procUnitConfObj1.addOperation(name='removeDC')
	74
	75	# opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
	76	# opObj14.addParameter(name='id', value='1', format='int')
	77	# # opObj14.addParameter(name='wintitle', value='Con interf', format='str')
	78	# opObj14.addParameter(name='save', value='1', format='bool')
	79	# opObj14.addParameter(name='figpath', value=pathFigure, format='str')
	80	# # opObj14.addParameter(name='zmin', value='5', format='int')
	81	# opObj14.addParameter(name='zmax', value='30', format='int')
	82	#
	83	# opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
	84	# opObj12.addParameter(name='id', value='2', format='int')
	85	# opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
	86	# opObj12.addParameter(name='save', value='1', format='bool')
	87	# opObj12.addParameter(name='figpath', value = pathFigure, format='str')
	88	# opObj12.addParameter(name='xmin', value=xmin, format='float')
	89	# opObj12.addParameter(name='xmax', value=xmax, format='float')
	90	# # opObj12.addParameter(name='zmin', value='5', format='int')
	91	# opObj12.addParameter(name='zmax', value='30', format='int')
	92
	93	#--------------------------------------------------------------------------------------------------
	94
	95	procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
	96	opObj20 = procUnitConfObj2.addOperation(name='SpectralFitting')
	97	opObj20.addParameter(name='path', value='/home/soporte/workspace/RemoteSystemsTempFiles', format='str')
	98	opObj20.addParameter(name='file', value='modelSpectralFitting', format='str')
	99	opObj20.addParameter(name='groupList', value='(0,1),(2,3)',format='multiList')
	100
	101	opObj11 = procUnitConfObj2.addOperation(name='SpectralFittingPlot', optype='other')
	102	opObj11.addParameter(name='id', value='3', format='int')
	103	opObj11.addParameter(name='wintitle', value='DopplerPlot', format='str')
	104	opObj11.addParameter(name='cutHeight', value='350', format='int')
	105	opObj11.addParameter(name='fit', value='1', format='int')#1--True/include fit
	106	opObj11.addParameter(name='save', value='1', format='bool')
	107	opObj11.addParameter(name='figpath', value = pathFigure, format='str')
	108
	109	opObj11 = procUnitConfObj2.addOperation(name='EWDriftsEstimation', optype='other')
	110	opObj11.addParameter(name='zenith', value='-3.80208,3.10658', format='floatlist')
	111	opObj11.addParameter(name='zenithCorrection', value='0.183201', format='float')
	112
	113	opObj23 = procUnitConfObj2.addOperation(name='EWDriftsPlot', optype='other')
	114	opObj23.addParameter(name='id', value='4', format='int')
	115	opObj23.addParameter(name='wintitle', value='EW Drifts', format='str')
	116	opObj23.addParameter(name='save', value='1', format='bool')
	117	opObj23.addParameter(name='figpath', value = pathFigure, format='str')
	118	opObj23.addParameter(name='zminZonal', value='-150', format='int')
	119	opObj23.addParameter(name='zmaxZonal', value='150', format='int')
	120	opObj23.addParameter(name='zminVertical', value='-30', format='float')
	121	opObj23.addParameter(name='zmaxVertical', value='30', format='float')
	122	opObj23.addParameter(name='SNR_1', value='1', format='bool')
	123	opObj23.addParameter(name='SNRmax', value='5', format='int')
	124	# opObj23.addParameter(name='SNRthresh', value='-50', format='float')
	125	opObj23.addParameter(name='xmin', value=xmin, format='float')
	126	opObj23.addParameter(name='xmax', value=xmax, format='float')
	127	#--------------------------------------------------------------------------------------------------
	128	print "Escribiendo el archivo XML"
	129	controllerObj.writeXml(filename)
	130	print "Leyendo el archivo XML"
	131	controllerObj.readXml(filename)
	132
	133	controllerObj.createObjects()
	134	controllerObj.connectObjects()
	135	controllerObj.run() No newline at end of file

schainpy/model/data/jrodata.py +22 -6

             class Parameters(JROData):
+                #Information from previous data
                 inputUnit = None        #Type of data to be processed
                 operation = None        #Type of operation to parametrize
+                normFactor = None       #Normalization Factor
+                groupList = None        #List of Pairs, Groups, etc
+                #Parameters
                 data_param = None       #Parameters obtained
                 data_pre = None         #Data Pre Parametrization
                 SNR = None              #Signal to Noise Ratio
-                pairsList = None        #List of Pairs for Cross correlations or Cross spectrum
                 initUtcTime = None      #Initial UTC time
                 paramInterval = None    #Time interval to calculate Parameters in seconds
-                windsInterval = None     #Time interval to calculate Winds in seconds
+                #Fitting
+                constants = None
+                error = None
+                library = None
+                #Output signal
+                outputInterval = None     #Time interval to calculate output signal in seconds
+                data_output = None       #Out signal
-                normFactor = None       #Normalization Factor
-                winds = None            #Wind estimations
                 def __init__(self):
                     '''
                     datatime = []
                     datatime.append(self.initUtcTime)
-                    datatime.append(self.initUtcTime + self.windsInterval - 1)
+                    datatime.append(self.initUtcTime + self.outputInterval - 1)
                     datatime = numpy.array(datatime)

schainpy/model/graphics/jroplot_parameters.py +405 -1

@@ -455,7 +455,7 class WindProfilerPlot(Figure):
455	# y = dataOut.heightRange	455	# y = dataOut.heightRange
456	y = dataOut.heightRange	456	y = dataOut.heightRange
457		457
458	z = dataOut.~~winds~~.copy()	458	z = dataOut.data_output.copy()
459	nplots = z.shape[0] #Number of wind dimensions estimated	459	nplots = z.shape[0] #Number of wind dimensions estimated
460	nplotsw = nplots	460	nplotsw = nplots
461		461
@@ -771,4 +771,408 class ParametersPlot(Figure):
771	if x[1] >= self.axesList[0].xmax:	771	if x[1] >= self.axesList[0].xmax:
772	self.counter_imagwr = wr_period	772	self.counter_imagwr = wr_period
773	self.__isConfig = False	773	self.__isConfig = False
		774	self.figfile = None
		775
		776
		777	class SpectralFittingPlot(Figure):
		778
		779	__isConfig = None
		780	__nsubplots = None
		781
		782	WIDTHPROF = None
		783	HEIGHTPROF = None
		784	PREFIX = 'prm'
		785
		786
		787	N = None
		788	ippSeconds = None
		789
		790	def __init__(self):
		791	self.__isConfig = False
		792	self.__nsubplots = 1
		793
		794	self.WIDTH = 450
		795	self.HEIGHT = 250
		796	self.WIDTHPROF = 0
		797	self.HEIGHTPROF = 0
		798
		799	def getSubplots(self):
		800
		801	ncol = int(numpy.sqrt(self.nplots)+0.9)
		802	nrow = int(self.nplots*1./ncol + 0.9)
		803
		804	return nrow, ncol
		805
		806	def setup(self, id, nplots, wintitle, showprofile=False, show=True):
		807
		808	showprofile = False
		809	self.__showprofile = showprofile
		810	self.nplots = nplots
		811
		812	ncolspan = 5
		813	colspan = 4
		814	if showprofile:
		815	ncolspan = 5
		816	colspan = 4
		817	self.__nsubplots = 2
		818
		819	self.createFigure(id = id,
		820	wintitle = wintitle,
		821	widthplot = self.WIDTH + self.WIDTHPROF,
		822	heightplot = self.HEIGHT + self.HEIGHTPROF,
		823	show=show)
		824
		825	nrow, ncol = self.getSubplots()
		826
		827	counter = 0
		828	for y in range(nrow):
		829	for x in range(ncol):
		830
		831	if counter >= self.nplots:
		832	break
		833
		834	self.addAxes(nrow, ncolncolspan, y, xncolspan, colspan, 1)
		835
		836	if showprofile:
		837	self.addAxes(nrow, ncolncolspan, y, xncolspan+colspan, 1, 1)
		838
		839	counter += 1
		840
		841	def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
		842	xmin=None, xmax=None, ymin=None, ymax=None,
		843	save=False, figpath='./', figfile=None, show=True):
		844
		845	"""
		846
		847	Input:
		848	dataOut :
		849	id :
		850	wintitle :
		851	channelList :
		852	showProfile :
		853	xmin : None,
		854	xmax : None,
		855	zmin : None,
		856	zmax : None
		857	"""
		858
		859	if cutHeight==None:
		860	h=270
		861	heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
		862	cutHeight = dataOut.heightList[heightindex]
		863
		864	factor = dataOut.normFactor
		865	x = dataOut.abscissaRange[:-1]
		866	#y = dataOut.getHeiRange()
		867
		868	z = dataOut.data_pre[:,:,heightindex]/factor
		869	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
		870	avg = numpy.average(z, axis=1)
		871	listChannels = z.shape[0]
		872
		873	#Reconstruct Function
		874	if fit==True:
		875	groupArray = dataOut.groupList
		876	listChannels = groupArray.reshape((groupArray.size))
		877	listChannels.sort()
		878	spcFitLine = numpy.zeros(z.shape)
		879	constants = dataOut.constants
		880
		881	nGroups = groupArray.shape[0]
		882	nChannels = groupArray.shape[1]
		883	nProfiles = z.shape[1]
		884
		885	for f in range(nGroups):
		886	groupChann = groupArray[f,:]
		887	p = dataOut.data_param[f,:,heightindex]
		888	# p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
		889	fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
		890	fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
		891	spcFitLine[groupChann,:] = fitLineAux
		892	# spcFitLine = spcFitLine/factor
		893
		894	z = z[listChannels,:]
		895	spcFitLine = spcFitLine[listChannels,:]
		896	spcFitLinedB = 10*numpy.log10(spcFitLine)
		897
		898	zdB = 10*numpy.log10(z)
		899	#thisDatetime = dataOut.datatime
		900	thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
		901	title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
		902	xlabel = "Velocity (m/s)"
		903	ylabel = "Spectrum"
		904
		905	if not self.__isConfig:
		906
		907	nplots = listChannels.size
		908
		909	self.setup(id=id,
		910	nplots=nplots,
		911	wintitle=wintitle,
		912	showprofile=showprofile,
		913	show=show)
		914
		915	if xmin == None: xmin = numpy.nanmin(x)
		916	if xmax == None: xmax = numpy.nanmax(x)
		917	if ymin == None: ymin = numpy.nanmin(zdB)
		918	if ymax == None: ymax = numpy.nanmax(zdB)+2
		919
		920	self.__isConfig = True
		921
		922	self.setWinTitle(title)
		923	for i in range(self.nplots):
		924	# title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
		925	title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i]+1)
		926	axes = self.axesList[i*self.__nsubplots]
		927	if fit == False:
		928	axes.pline(x, zdB[i,:],
		929	xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
		930	xlabel=xlabel, ylabel=ylabel, title=title
		931	)
		932	if fit == True:
		933	fitline=spcFitLinedB[i,:]
		934	y=numpy.vstack([zdB[i,:],fitline] )
		935	legendlabels=['Data','Fitting']
		936	axes.pmultilineyaxis(x, y,
		937	xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
		938	xlabel=xlabel, ylabel=ylabel, title=title,
		939	legendlabels=legendlabels, marker=None,
		940	linestyle='solid', grid='both')
		941
		942	self.draw()
		943
		944	if save:
		945	date = thisDatetime.strftime("%Y%m%d_%H%M%S")
		946	if figfile == None:
		947	figfile = self.getFilename(name = date)
		948
		949	self.saveFigure(figpath, figfile)
		950
		951
		952	class EWDriftsPlot(Figure):
		953
		954	__isConfig = None
		955	__nsubplots = None
		956
		957	WIDTHPROF = None
		958	HEIGHTPROF = None
		959	PREFIX = 'drift'
		960
		961	def __init__(self):
		962
		963	self.timerange = 26060
		964	self.isConfig = False
		965	self.__nsubplots = 1
		966
		967	self.WIDTH = 800
		968	self.HEIGHT = 150
		969	self.WIDTHPROF = 120
		970	self.HEIGHTPROF = 0
		971	self.counter_imagwr = 0
		972
		973	self.PLOT_CODE = 0
		974	self.FTP_WEI = None
		975	self.EXP_CODE = None
		976	self.SUB_EXP_CODE = None
		977	self.PLOT_POS = None
		978	self.tmin = None
		979	self.tmax = None
		980
		981	self.xmin = None
		982	self.xmax = None
		983
		984	self.figfile = None
		985
		986	def getSubplots(self):
		987
		988	ncol = 1
		989	nrow = self.nplots
		990
		991	return nrow, ncol
		992
		993	def setup(self, id, nplots, wintitle, showprofile=True, show=True):
		994
		995	self.__showprofile = showprofile
		996	self.nplots = nplots
		997
		998	ncolspan = 1
		999	colspan = 1
		1000
		1001	self.createFigure(id = id,
		1002	wintitle = wintitle,
		1003	widthplot = self.WIDTH + self.WIDTHPROF,
		1004	heightplot = self.HEIGHT + self.HEIGHTPROF,
		1005	show=show)
		1006
		1007	nrow, ncol = self.getSubplots()
		1008
		1009	counter = 0
		1010	for y in range(nrow):
		1011	if counter >= self.nplots:
		1012	break
		1013
		1014	self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
		1015	counter += 1
		1016
		1017	def run(self, dataOut, id, wintitle="", channelList=None,
		1018	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
		1019	zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
		1020	timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
		1021	save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
		1022	server=None, folder=None, username=None, password=None,
		1023	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
		1024	"""
		1025
		1026	Input:
		1027	dataOut :
		1028	id :
		1029	wintitle :
		1030	channelList :
		1031	showProfile :
		1032	xmin : None,
		1033	xmax : None,
		1034	ymin : None,
		1035	ymax : None,
		1036	zmin : None,
		1037	zmax : None
		1038	"""
		1039
		1040	if channelList == None:
		1041	channelIndexList = dataOut.channelIndexList
		1042	else:
		1043	channelIndexList = []
		1044	for channel in channelList:
		1045	if channel not in dataOut.channelList:
		1046	raise ValueError, "Channel %d is not in dataOut.channelList"
		1047	channelIndexList.append(dataOut.channelList.index(channel))
		1048
		1049	if timerange != None:
		1050	self.timerange = timerange
		1051
		1052	tmin = None
		1053	tmax = None
		1054
		1055	x = dataOut.getTimeRange1()
		1056	# y = dataOut.heightRange
		1057	y = dataOut.heightList
		1058
		1059	z = dataOut.data_output
		1060	nplots = z.shape[0] #Number of wind dimensions estimated
		1061	nplotsw = nplots
		1062
		1063	#If there is a SNR function defined
		1064	if dataOut.SNR != None:
		1065	nplots += 1
		1066	SNR = dataOut.SNR
		1067
		1068	if SNR_1:
		1069	SNR += 1
		1070
		1071	SNRavg = numpy.average(SNR, axis=0)
		1072
		1073	SNRdB = 10*numpy.log10(SNR)
		1074	SNRavgdB = 10*numpy.log10(SNRavg)
		1075
		1076	ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
		1077
		1078	for i in range(nplotsw):
		1079	z[i,ind] = numpy.nan
		1080
		1081
		1082	showprofile = False
		1083	# thisDatetime = dataOut.datatime
		1084	thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
		1085	title = wintitle + " EW Drifts"
		1086	xlabel = ""
		1087	ylabel = "Height (Km)"
		1088
		1089	if not self.__isConfig:
		1090
		1091	self.setup(id=id,
		1092	nplots=nplots,
		1093	wintitle=wintitle,
		1094	showprofile=showprofile,
		1095	show=show)
		1096
		1097	self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
		1098
		1099	if ymin == None: ymin = numpy.nanmin(y)
		1100	if ymax == None: ymax = numpy.nanmax(y)
		1101
		1102	if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
		1103	if zminZonal == None: zminZonal = -zmaxZonal
		1104	if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
		1105	if zminVertical == None: zminVertical = -zmaxVertical
		1106
		1107	if dataOut.SNR != None:
		1108	if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
		1109	if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
		1110
		1111	self.FTP_WEI = ftp_wei
		1112	self.EXP_CODE = exp_code
		1113	self.SUB_EXP_CODE = sub_exp_code
		1114	self.PLOT_POS = plot_pos
		1115
		1116	self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
		1117	self.__isConfig = True
		1118
		1119
		1120	self.setWinTitle(title)
		1121
		1122	if ((self.xmax - x[1]) < (x[1]-x[0])):
		1123	x[1] = self.xmax
		1124
		1125	strWind = ['Zonal','Vertical']
		1126	strCb = 'Velocity (m/s)'
		1127	zmaxVector = [zmaxZonal, zmaxVertical]
		1128	zminVector = [zminZonal, zminVertical]
		1129
		1130	for i in range(nplotsw):
		1131
		1132	title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
		1133	axes = self.axesList[i*self.__nsubplots]
		1134
		1135	z1 = z[i,:].reshape((1,-1))
		1136
		1137	axes.pcolorbuffer(x, y, z1,
		1138	xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
		1139	xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
		1140	ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
		1141
		1142	if dataOut.SNR != None:
		1143	i += 1
		1144	if SNR_1:
		1145	title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
		1146	else:
		1147	title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
		1148	axes = self.axesList[i*self.__nsubplots]
		1149	SNRavgdB = SNRavgdB.reshape((1,-1))
		1150
		1151	axes.pcolorbuffer(x, y, SNRavgdB,
		1152	xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
		1153	xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
		1154	ticksize=9, cblabel='', cbsize="1%", colormap="jet")
		1155
		1156	self.draw()
		1157
		1158	if self.figfile == None:
		1159	str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
		1160	self.figfile = self.getFilename(name = str_datetime)
		1161
		1162	if figpath != '':
		1163
		1164	self.counter_imagwr += 1
		1165	if (self.counter_imagwr>=wr_period):
		1166	# store png plot to local folder
		1167	self.saveFigure(figpath, self.figfile)
		1168	# store png plot to FTP server according to RT-Web format
		1169	name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
		1170	ftp_filename = os.path.join(figpath, name)
		1171	self.saveFigure(figpath, ftp_filename)
		1172
		1173	self.counter_imagwr = 0
		1174
		1175	if x[1] >= self.axesList[0].xmax:
		1176	self.counter_imagwr = wr_period
		1177	self.__isConfig = False
774	self.figfile = None No newline at end of file	1178	self.figfile = None

schainpy/model/graphics/mpldriver.py +1 -1

                 matplotlib.pyplot.ioff()
             #    lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
-                lines = ax.plot(x, y.T, linestyle='None', marker='.', markersize=markersize)
+                lines = ax.plot(x, y.T, linestyle=linestyle, marker=marker, markersize=markersize)
                 leg = ax.legend(lines, legendlabels, loc='upper left', bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
                                 handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)

schainpy/model/proc/jroproc_parameters.py +228 -18

@@ -7,7 +7,9 from scipy import stats
7	import re	7	import re
8	import datetime	8	import datetime
9	import copy	9	import copy
10		10	import sys
		11	import importlib
		12	import itertools
11		13
12	from jroproc_base import ProcessingUnit, Operation	14	from jroproc_base import ProcessingUnit, Operation
13	from model.data.jrodata import Parameters	15	from model.data.jrodata import Parameters
@@ -60,7 +62,7 class ParametersProc(ProcessingUnit):
60		62
61	def run(self, nSeconds = None, nProfiles = None):	63	def run(self, nSeconds = None, nProfiles = None):
62		64
63	self.dataOut.flagNoData = True	65
64		66
65	if self.firstdatatime == None:	67	if self.firstdatatime == None:
66	self.firstdatatime = self.dataIn.utctime	68	self.firstdatatime = self.dataIn.utctime
@@ -68,6 +70,7 class ParametersProc(ProcessingUnit):
68	#---------------------- Voltage Data ---------------------------	70	#---------------------- Voltage Data ---------------------------
69		71
70	if self.dataIn.type == "Voltage":	72	if self.dataIn.type == "Voltage":
		73	self.dataOut.flagNoData = True
71	if nSeconds != None:	74	if nSeconds != None:
72	self.nSeconds = nSeconds	75	self.nSeconds = nSeconds
73	self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))	76	self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
@@ -100,6 +103,7 class ParametersProc(ProcessingUnit):
100	self.dataOut.abscissaRange = self.dataIn.getVelRange(1)	103	self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
101	self.dataOut.noise = self.dataIn.getNoise()	104	self.dataOut.noise = self.dataIn.getNoise()
102	self.dataOut.normFactor = self.dataIn.normFactor	105	self.dataOut.normFactor = self.dataIn.normFactor
		106	self.dataOut.flagNoData = False
103		107
104	#---------------------- Correlation Data ---------------------------	108	#---------------------- Correlation Data ---------------------------
105		109
@@ -112,14 +116,14 class ParametersProc(ProcessingUnit):
112	self.dataOut.noise = self.dataIn.noise	116	self.dataOut.noise = self.dataIn.noise
113	self.dataOut.normFactor = self.dataIn.normFactor	117	self.dataOut.normFactor = self.dataIn.normFactor
114	self.dataOut.SNR = self.dataIn.SNR	118	self.dataOut.SNR = self.dataIn.SNR
115	self.dataOut.~~pairs~~List = self.dataIn.pairsList	119	self.dataOut.groupList = self.dataIn.pairsList
		120	self.dataOut.flagNoData = False
116		121
117		122
118	self.__updateObjFromInput()	123	self.__updateObjFromInput()
119	self.dataOut.flagNoData = False
120	self.firstdatatime = None	124	self.firstdatatime = None
121	self.dataOut.initUtcTime = self.dataIn.ltctime	125	self.dataOut.initUtcTime = self.dataIn.ltctime
122	self.dataOut.~~winds~~Interval = self.dataIn.timeInterval	126	self.dataOut.outputInterval = self.dataIn.timeInterval
123		127
124	#------------------- Get Moments ----------------------------------	128	#------------------- Get Moments ----------------------------------
125	def GetMoments(self, channelList = None):	129	def GetMoments(self, channelList = None):
@@ -238,7 +242,7 class ParametersProc(ProcessingUnit):
238	self.dataOut.noise	242	self.dataOut.noise
239	self.dataOut.normFactor	243	self.dataOut.normFactor
240	self.dataOut.SNR	244	self.dataOut.SNR
241	self.dataOut.~~pairs~~List	245	self.dataOut.groupList
242	self.dataOut.nChannels	246	self.dataOut.nChannels
243		247
244	Affected:	248	Affected:
@@ -251,7 +255,7 class ParametersProc(ProcessingUnit):
251	absc = self.dataOut.abscissaRange[:-1]	255	absc = self.dataOut.abscissaRange[:-1]
252	noise = self.dataOut.noise	256	noise = self.dataOut.noise
253	SNR = self.dataOut.SNR	257	SNR = self.dataOut.SNR
254	pairsList = self.dataOut.~~pairs~~List	258	pairsList = self.dataOut.groupList
255	nChannels = self.dataOut.nChannels	259	nChannels = self.dataOut.nChannels
256	pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)	260	pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
257	self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))	261	self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
@@ -1109,6 +1113,137 class ParametersProc(ProcessingUnit):
1109		1113
1110	return heights, error	1114	return heights, error
1111		1115
		1116	def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
		1117
		1118	'''
		1119	Function GetMoments()
		1120
		1121	Input:
		1122	Output:
		1123	Variables modified:
		1124	'''
		1125	if path != None:
		1126	sys.path.append(path)
		1127	self.dataOut.library = importlib.import_module(file)
		1128
		1129	#To be inserted as a parameter
		1130	groupArray = numpy.array(groupList)
		1131	# groupArray = numpy.array([[0,1],[2,3]])
		1132	self.dataOut.groupList = groupArray
		1133
		1134	nGroups = groupArray.shape[0]
		1135	nChannels = self.dataIn.nChannels
		1136	nHeights=self.dataIn.heightList.size
		1137
		1138	#Parameters Array
		1139	self.dataOut.data_param = None
		1140
		1141	#Set constants
		1142	constants = self.dataOut.library.setConstants(self.dataIn)
		1143	self.dataOut.constants = constants
		1144	M = self.dataIn.normFactor
		1145	N = self.dataIn.nFFTPoints
		1146	ippSeconds = self.dataIn.ippSeconds
		1147	K = self.dataIn.nIncohInt
		1148	pairsArray = numpy.array(self.dataIn.pairsList)
		1149
		1150	#List of possible combinations
		1151	listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
		1152	indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
		1153
		1154	if getSNR:
		1155	listChannels = groupArray.reshape((groupArray.size))
		1156	listChannels.sort()
		1157	noise = self.dataIn.getNoise()
		1158	self.dataOut.SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
		1159
		1160	for i in range(nGroups):
		1161	coord = groupArray[i,:]
		1162
		1163	#Input data array
		1164	data = self.dataIn.data_spc[coord,:,:]/(M*N)
		1165	data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
		1166
		1167	#Cross Spectra data array for Covariance Matrixes
		1168	ind = 0
		1169	for pairs in listComb:
		1170	pairsSel = numpy.array([coord[x],coord[y]])
		1171	indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
		1172	ind += 1
		1173	dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
		1174	dataCross = dataCross**2/K
		1175
		1176	for h in range(nHeights):
		1177	# print self.dataOut.heightList[h]
		1178
		1179	#Input
		1180	d = data[:,h]
		1181
		1182	#Covariance Matrix
		1183	D = numpy.diag(d**2/K)
		1184	ind = 0
		1185	for pairs in listComb:
		1186	#Coordinates in Covariance Matrix
		1187	x = pairs[0]
		1188	y = pairs[1]
		1189	#Channel Index
		1190	S12 = dataCross[ind,:,h]
		1191	D12 = numpy.diag(S12)
		1192	#Completing Covariance Matrix with Cross Spectras
		1193	D[xN:(x+1)N,yN:(y+1)N] = D12
		1194	D[yN:(y+1)N,xN:(x+1)N] = D12
		1195	ind += 1
		1196	Dinv=numpy.linalg.inv(D)
		1197	L=numpy.linalg.cholesky(Dinv)
		1198	LT=L.T
		1199
		1200	dp = numpy.dot(LT,d)
		1201
		1202	#Initial values
		1203	data_spc = self.dataIn.data_spc[coord,:,h]
		1204	p0 = self.dataOut.library.initialValuesFunction(data_spc, constants)
		1205
		1206	#Least Squares
		1207	minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
		1208	# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
		1209	#Chi square error
		1210	error0 = numpy.sum(infodict['fvec']*2)/(2N)
		1211	# error0 = 0
		1212	#Error with Jacobian
		1213	error1 = self.dataOut.library.errorFunction(minp,constants,LT)
		1214	#Save
		1215	if self.dataOut.data_param == None:
		1216	self.dataOut.data_param = numpy.zeros((nGroups, minp.size, nHeights))*numpy.nan
		1217	self.dataOut.error = numpy.zeros((nGroups, error1.size + 1, nHeights))*numpy.nan
		1218
		1219	self.dataOut.error[i,:,h] = numpy.hstack((error0,error1))
		1220	self.dataOut.data_param[i,:,h] = minp
		1221	return
		1222
		1223
		1224	def __residFunction(self, p, dp, LT, constants):
		1225
		1226	fm = self.dataOut.library.modelFunction(p, constants)
		1227	fmp=numpy.dot(LT,fm)
		1228
		1229	return dp-fmp
		1230
		1231	def __getSNR(self, z, noise):
		1232
		1233	avg = numpy.average(z, axis=1)
		1234	SNR = (avg.T-noise)/noise
		1235	SNR = SNR.T
		1236	return SNR
		1237
		1238	def __chisq(p,chindex,hindex):
		1239	#similar to Resid but calculates CHI**2
		1240	[LT,d,fm]=setupLTdfm(p,chindex,hindex)
		1241	dp=numpy.dot(LT,d)
		1242	fmp=numpy.dot(LT,fm)
		1243	chisq=numpy.dot((dp-fmp).T,(dp-fmp))
		1244	return chisq
		1245
		1246
1112		1247
1113	class WindProfiler(Operation):	1248	class WindProfiler(Operation):
1114		1249
@@ -1359,12 +1494,12 class WindProfiler(Operation):
1359	winds = correctFactor*winds	1494	winds = correctFactor*winds
1360	return winds	1495	return winds
1361		1496
1362	def __checkTime(self, currentTime, paramInterval, ~~winds~~Interval):	1497	def __checkTime(self, currentTime, paramInterval, outputInterval):
1363		1498
1364	dataTime = currentTime + paramInterval	1499	dataTime = currentTime + paramInterval
1365	deltaTime = dataTime - self.__initime	1500	deltaTime = dataTime - self.__initime
1366		1501
1367	if deltaTime >= ~~winds~~Interval or deltaTime < 0:	1502	if deltaTime >= outputInterval or deltaTime < 0:
1368	self.__dataReady = True	1503	self.__dataReady = True
1369	return	1504	return
1370		1505
@@ -1462,7 +1597,7 class WindProfiler(Operation):
1462	theta_y = theta_y[arrayChannel]	1597	theta_y = theta_y[arrayChannel]
1463		1598
1464	velRadial0 = param[:,1,:] #Radial velocity	1599	velRadial0 = param[:,1,:] #Radial velocity
1465	dataOut.~~winds~~, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function	1600	dataOut.data_output, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1466		1601
1467	elif technique == 'SA':	1602	elif technique == 'SA':
1468		1603
@@ -1483,12 +1618,12 class WindProfiler(Operation):
1483		1618
1484	tau = dataOut.data_param	1619	tau = dataOut.data_param
1485	_lambda = dataOut.C/dataOut.frequency	1620	_lambda = dataOut.C/dataOut.frequency
1486	pairsList = dataOut.~~pairs~~List	1621	pairsList = dataOut.groupList
1487	nChannels = dataOut.nChannels	1622	nChannels = dataOut.nChannels
1488		1623
1489	dataOut.~~winds~~ = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)	1624	dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1490	dataOut.initUtcTime = dataOut.ltctime	1625	dataOut.initUtcTime = dataOut.ltctime
1491	dataOut.~~winds~~Interval = dataOut.timeInterval	1626	dataOut.outputInterval = dataOut.timeInterval
1492		1627
1493	elif technique == 'Meteors':	1628	elif technique == 'Meteors':
1494	dataOut.flagNoData = True	1629	dataOut.flagNoData = True
@@ -1511,7 +1646,7 class WindProfiler(Operation):
1511	hmax = kwargs['hmax']	1646	hmax = kwargs['hmax']
1512	else: hmax = 110	1647	else: hmax = 110
1513		1648
1514	dataOut.~~winds~~Interval = nHours*3600	1649	dataOut.outputInterval = nHours*3600
1515		1650
1516	if self.__isConfig == False:	1651	if self.__isConfig == False:
1517	# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)	1652	# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
@@ -1526,14 +1661,89 class WindProfiler(Operation):
1526	else:	1661	else:
1527	self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))	1662	self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1528		1663
1529	self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.~~winds~~Interval) #Check if the buffer is ready	1664	self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1530		1665
1531	if self.__dataReady:	1666	if self.__dataReady:
1532	dataOut.initUtcTime = self.__initime	1667	dataOut.initUtcTime = self.__initime
1533	self.__initime = self.__initime + dataOut.~~winds~~Interval #to erase time offset	1668	self.__initime = self.__initime + dataOut.outputInterval #to erase time offset
1534		1669
1535	dataOut.~~winds~~, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)	1670	dataOut.data_output, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1536	dataOut.flagNoData = False	1671	dataOut.flagNoData = False
1537	self.__buffer = None	1672	self.__buffer = None
1538		1673
1539	return No newline at end of file	1674	return
		1675
		1676	class EWDriftsEstimation(Operation):
		1677
		1678
		1679	def __init__(self):
		1680	Operation.__init__(self)
		1681
		1682	def __correctValues(self, heiRang, phi, velRadial, SNR):
		1683	listPhi = phi.tolist()
		1684	maxid = listPhi.index(max(listPhi))
		1685	minid = listPhi.index(min(listPhi))
		1686
		1687	rango = range(len(phi))
		1688	# rango = numpy.delete(rango,maxid)
		1689
		1690	heiRang1 = heiRang*math.cos(phi[maxid])
		1691	heiRangAux = heiRang*math.cos(phi[minid])
		1692	indOut = (heiRang1 < heiRangAux[0]).nonzero()
		1693	heiRang1 = numpy.delete(heiRang1,indOut)
		1694
		1695	velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
		1696	SNR1 = numpy.zeros([len(phi),len(heiRang1)])
		1697
		1698	for i in rango:
		1699	x = heiRang*math.cos(phi[i])
		1700	y1 = velRadial[i,:]
		1701	f1 = interpolate.interp1d(x,y1,kind = 'cubic')
		1702
		1703	x1 = heiRang1
		1704	y11 = f1(x1)
		1705
		1706	y2 = SNR[i,:]
		1707	f2 = interpolate.interp1d(x,y2,kind = 'cubic')
		1708	y21 = f2(x1)
		1709
		1710	velRadial1[i,:] = y11
		1711	SNR1[i,:] = y21
		1712
		1713	return heiRang1, velRadial1, SNR1
		1714
		1715	def run(self, dataOut, zenith, zenithCorrection):
		1716	heiRang = dataOut.heightList
		1717	velRadial = dataOut.data_param[:,3,:]
		1718	SNR = dataOut.SNR
		1719
		1720	zenith = numpy.array(zenith)
		1721	zenith -= zenithCorrection
		1722	zenith *= numpy.pi/180
		1723
		1724	heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
		1725
		1726	alp = zenith[0]
		1727	bet = zenith[1]
		1728
		1729	w_w = velRadial1[0,:]
		1730	w_e = velRadial1[1,:]
		1731
		1732	w = (w_wnumpy.sin(bet) - w_enumpy.sin(alp))/(numpy.cos(alp)numpy.sin(bet) - numpy.cos(bet)numpy.sin(alp))
		1733	u = (w_wnumpy.cos(bet) - w_enumpy.cos(alp))/(numpy.sin(alp)numpy.cos(bet) - numpy.sin(bet)numpy.cos(alp))
		1734
		1735	winds = numpy.vstack((u,w))
		1736
		1737	dataOut.heightList = heiRang1
		1738	dataOut.data_output = winds
		1739	dataOut.SNR = SNR1
		1740
		1741	dataOut.initUtcTime = dataOut.ltctime
		1742	dataOut.outputInterval = dataOut.timeInterval
		1743	return
		1744
		1745
		1746
		1747
		1748
		1749	No newline at end of file

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