schain-jc Commit - r1032:c711a14430f7

merge with BLTR and Madrigal modules

Juan C. Espinoza -

r1032:c711a14430f7

parent child

Context file:

r1032:c711a14430f7

Expand all files

schainpy/model/data/BLTRheaderIO.py created 644 +404 0

			@@ -0,0 +1,404
		1	'''
		2
		3	$Author: murco $
		4	$Id: JROHeaderIO.py 151 2012-10-31 19:00:51Z murco $
		5	'''
		6	import sys
		7	import numpy
		8	import copy
		9	import datetime
		10	from __builtin__ import None
		11
		12	SPEED_OF_LIGHT = 299792458
		13	SPEED_OF_LIGHT = 3e8
		14
		15	FILE_STRUCTURE = numpy.dtype([ #HEADER 48bytes
		16	('FileMgcNumber','<u4'), #0x23020100
		17	('nFDTdataRecors','<u4'), #No Of FDT data records in this file (0 or more)
		18	('RadarUnitId','<u4'),
		19	('SiteName','<s32'), #Null terminated
		20	])
		21
		22	RECORD_STRUCTURE = numpy.dtype([ #RECORD HEADER 180+20N bytes
		23	('RecMgcNumber','<u4'), #0x23030001
		24	('RecCounter','<u4'), #Record counter(0,1, ...)
		25	('Off2StartNxtRec','<u4'), #Offset to start of next record form start of this record
		26	('Off2StartData','<u4'), #Offset to start of data from start of this record
		27	('EpTimeStamp','<i4'), #Epoch time stamp of start of acquisition (seconds)
		28	('msCompTimeStamp','<u4'), #Millisecond component of time stamp (0,...,999)
		29	('ExpTagName','<s32'), #Experiment tag name (null terminated)
		30	('ExpComment','<s32'), #Experiment comment (null terminated)
		31	('SiteLatDegrees','<f4'), #Site latitude (from GPS) in degrees (positive implies North)
		32	('SiteLongDegrees','<f4'), #Site longitude (from GPS) in degrees (positive implies East)
		33	('RTCgpsStatus','<u4'), #RTC GPS engine status (0=SEEK, 1=LOCK, 2=NOT FITTED, 3=UNAVAILABLE)
		34	('TransmitFrec','<u4'), #Transmit frequency (Hz)
		35	('ReceiveFrec','<u4'), #Receive frequency
		36	('FirstOsciFrec','<u4'), #First local oscillator frequency (Hz)
		37	('Polarisation','<u4'), #(0="O", 1="E", 2="linear 1", 3="linear2")
		38	('ReceiverFiltSett','<u4'), #Receiver filter settings (0,1,2,3)
		39	('nModesInUse','<u4'), #Number of modes in use (1 or 2)
		40	('DualModeIndex','<u4'), #Dual Mode index number for these data (0 or 1)
		41	('DualModeRange','<u4'), #Dual Mode range correction for these data (m)
		42	('nDigChannels','<u4'), #Number of digital channels acquired (2*N)
		43	('SampResolution','<u4'), #Sampling resolution (meters)
		44	('nRangeGatesSamp','<u4'), #Number of range gates sampled
		45	('StartRangeSamp','<u4'), #Start range of sampling (meters)
		46	('PRFhz','<u4'), #PRF (Hz)
		47	('Integrations','<u4'), #Integrations
		48	('nDataPointsTrsf','<u4'), #Number of data points transformed
		49	('nReceiveBeams','<u4'), #Number of receive beams stored in file (1 or N)
		50	('nSpectAverages','<u4'), #Number of spectral averages
		51	('FFTwindowingInd','<u4'), #FFT windowing index (0 = no window)
		52	('BeamAngleAzim','<f4'), #Beam steer angle (azimuth) in degrees (clockwise from true North)
		53	('BeamAngleZen','<f4'), #Beam steer angle (zenith) in degrees (0=> vertical)
		54	('AntennaCoord','<f24'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
		55	('RecPhaseCalibr','<f12'), #Receiver phase calibration (degrees) - N values
		56	('RecAmpCalibr','<f12'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
		57	('ReceiverGaindB','<u12'), #Receiver gains in dB - N values
		58	])
		59
		60
		61	class Header(object):
		62
		63	def __init__(self):
		64	raise NotImplementedError
		65
		66
		67	def read(self):
		68
		69	raise NotImplementedError
		70
		71	def write(self):
		72
		73	raise NotImplementedError
		74
		75	def printInfo(self):
		76
		77	message = "#"*50 + "\n"
		78	message += self.__class__.__name__.upper() + "\n"
		79	message += "#"*50 + "\n"
		80
		81	keyList = self.__dict__.keys()
		82	keyList.sort()
		83
		84	for key in keyList:
		85	message += "%s = %s" %(key, self.__dict__[key]) + "\n"
		86
		87	if "size" not in keyList:
		88	attr = getattr(self, "size")
		89
		90	if attr:
		91	message += "%s = %s" %("size", attr) + "\n"
		92
		93	print message
		94
		95	class FileHeader(Header):
		96
		97	FileMgcNumber= None
		98	nFDTdataRecors=None #No Of FDT data records in this file (0 or more)
		99	RadarUnitId= None
		100	SiteName= None
		101
		102	#__LOCALTIME = None
		103
		104	def __init__(self, useLocalTime=True):
		105
		106	self.FileMgcNumber= 0 #0x23020100
		107	self.nFDTdataRecors=0 #No Of FDT data records in this file (0 or more)
		108	self.RadarUnitId= 0
		109	self.SiteName= ""
		110	self.size = 48
		111
		112	#self.useLocalTime = useLocalTime
		113
		114	def read(self, fp):
		115
		116	try:
		117	header = numpy.fromfile(fp, FILE_STRUCTURE,1)
		118	''' numpy.fromfile(file, dtype, count, sep='')
		119	file : file or str
		120	Open file object or filename.
		121
		122	dtype : data-type
		123	Data type of the returned array. For binary files, it is used to determine
		124	the size and byte-order of the items in the file.
		125
		126	count : int
		127	Number of items to read. -1 means all items (i.e., the complete file).
		128
		129	sep : str
		130	Separator between items if file is a text file. Empty (“”) separator means
		131	the file should be treated as binary. Spaces (” ”) in the separator match zero
		132	or more whitespace characters. A separator consisting only of spaces must match
		133	at least one whitespace.
		134
		135	'''
		136
		137	except Exception, e:
		138	print "FileHeader: "
		139	print eBasicHeader
		140	return 0
		141
		142	self.FileMgcNumber= byte(header['FileMgcNumber'][0])
		143	self.nFDTdataRecors=int(header['nFDTdataRecors'][0]) #No Of FDT data records in this file (0 or more)
		144	self.RadarUnitId= int(header['RadarUnitId'][0])
		145	self.SiteName= char(header['SiteName'][0])
		146
		147
		148	if self.size <48:
		149	return 0
		150
		151	return 1
		152
		153	def write(self, fp):
		154
		155	headerTuple = (self.FileMgcNumber,
		156	self.nFDTdataRecors,
		157	self.RadarUnitId,
		158	self.SiteName,
		159	self.size)
		160
		161
		162	header = numpy.array(headerTuple, FILE_STRUCTURE)
		163	# numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0)
		164	header.tofile(fp)
		165	''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
		166
		167	fid : file or str
		168	An open file object, or a string containing a filename.
		169
		170	sep : str
		171	Separator between array items for text output. If “” (empty), a binary file is written,
		172	equivalent to file.write(a.tobytes()).
		173
		174	format : str
		175	Format string for text file output. Each entry in the array is formatted to text by
		176	first converting it to the closest Python type, and then using “format” % item.
		177
		178	'''
		179
		180	return 1
		181
		182
		183	class RecordHeader(Header):
		184
		185	RecMgcNumber=None #0x23030001
		186	RecCounter= None
		187	Off2StartNxtRec= None
		188	EpTimeStamp= None
		189	msCompTimeStamp= None
		190	ExpTagName= None
		191	ExpComment=None
		192	SiteLatDegrees=None
		193	SiteLongDegrees= None
		194	RTCgpsStatus= None
		195	TransmitFrec= None
		196	ReceiveFrec= None
		197	FirstOsciFrec= None
		198	Polarisation= None
		199	ReceiverFiltSett= None
		200	nModesInUse= None
		201	DualModeIndex= None
		202	DualModeRange= None
		203	nDigChannels= None
		204	SampResolution= None
		205	nRangeGatesSamp= None
		206	StartRangeSamp= None
		207	PRFhz= None
		208	Integrations= None
		209	nDataPointsTrsf= None
		210	nReceiveBeams= None
		211	nSpectAverages= None
		212	FFTwindowingInd= None
		213	BeamAngleAzim= None
		214	BeamAngleZen= None
		215	AntennaCoord= None
		216	RecPhaseCalibr= None
		217	RecAmpCalibr= None
		218	ReceiverGaindB= None
		219
		220	'''size = None
		221	nSamples = None
		222	nProfiles = None
		223	nChannels = None
		224	adcResolution = None
		225	pciDioBusWidth = None'''
		226
		227	def __init__(self, RecMgcNumber=None, RecCounter= 0, Off2StartNxtRec= 0,
		228	EpTimeStamp= 0, msCompTimeStamp= 0, ExpTagName= None,
		229	ExpComment=None, SiteLatDegrees=0, SiteLongDegrees= 0,
		230	RTCgpsStatus= 0, TransmitFrec= 0, ReceiveFrec= 0,
		231	FirstOsciFrec= 0, Polarisation= 0, ReceiverFiltSett= 0,
		232	nModesInUse= 0, DualModeIndex= 0, DualModeRange= 0,
		233	nDigChannels= 0, SampResolution= 0, nRangeGatesSamp= 0,
		234	StartRangeSamp= 0, PRFhz= 0, Integrations= 0,
		235	nDataPointsTrsf= 0, nReceiveBeams= 0, nSpectAverages= 0,
		236	FFTwindowingInd= 0, BeamAngleAzim= 0, BeamAngleZen= 0,
		237	AntennaCoord= 0, RecPhaseCalibr= 0, RecAmpCalibr= 0,
		238	ReceiverGaindB= 0):
		239
		240	self.RecMgcNumber = RecMgcNumber #0x23030001
		241	self.RecCounter = RecCounter
		242	self.Off2StartNxtRec = Off2StartNxtRec
		243	self.EpTimeStamp = EpTimeStamp
		244	self.msCompTimeStamp = msCompTimeStamp
		245	self.ExpTagName = ExpTagName
		246	self.ExpComment = ExpComment
		247	self.SiteLatDegrees = SiteLatDegrees
		248	self.SiteLongDegrees = SiteLongDegrees
		249	self.RTCgpsStatus = RTCgpsStatus
		250	self.TransmitFrec = TransmitFrec
		251	self.ReceiveFrec = ReceiveFrec
		252	self.FirstOsciFrec = FirstOsciFrec
		253	self.Polarisation = Polarisation
		254	self.ReceiverFiltSett = ReceiverFiltSett
		255	self.nModesInUse = nModesInUse
		256	self.DualModeIndex = DualModeIndex
		257	self.DualModeRange = DualModeRange
		258	self.nDigChannels = nDigChannels
		259	self.SampResolution = SampResolution
		260	self.nRangeGatesSamp = nRangeGatesSamp
		261	self.StartRangeSamp = StartRangeSamp
		262	self.PRFhz = PRFhz
		263	self.Integrations = Integrations
		264	self.nDataPointsTrsf = nDataPointsTrsf
		265	self.nReceiveBeams = nReceiveBeams
		266	self.nSpectAverages = nSpectAverages
		267	self.FFTwindowingInd = FFTwindowingInd
		268	self.BeamAngleAzim = BeamAngleAzim
		269	self.BeamAngleZen = BeamAngleZen
		270	self.AntennaCoord = AntennaCoord
		271	self.RecPhaseCalibr = RecPhaseCalibr
		272	self.RecAmpCalibr = RecAmpCalibr
		273	self.ReceiverGaindB = ReceiverGaindB
		274
		275
		276	def read(self, fp):
		277
		278	startFp = fp.tell() #The method tell() returns the current position of the file read/write pointer within the file.
		279
		280	try:
		281	header = numpy.fromfile(fp,RECORD_STRUCTURE,1)
		282	except Exception, e:
		283	print "System Header: " + e
		284	return 0
		285
		286	self.RecMgcNumber = header['RecMgcNumber'][0] #0x23030001
		287	self.RecCounter = header['RecCounter'][0]
		288	self.Off2StartNxtRec = header['Off2StartNxtRec'][0]
		289	self.EpTimeStamp = header['EpTimeStamp'][0]
		290	self.msCompTimeStamp = header['msCompTimeStamp'][0]
		291	self.ExpTagName = header['ExpTagName'][0]
		292	self.ExpComment = header['ExpComment'][0]
		293	self.SiteLatDegrees = header['SiteLatDegrees'][0]
		294	self.SiteLongDegrees = header['SiteLongDegrees'][0]
		295	self.RTCgpsStatus = header['RTCgpsStatus'][0]
		296	self.TransmitFrec = header['TransmitFrec'][0]
		297	self.ReceiveFrec = header['ReceiveFrec'][0]
		298	self.FirstOsciFrec = header['FirstOsciFrec'][0]
		299	self.Polarisation = header['Polarisation'][0]
		300	self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
		301	self.nModesInUse = header['nModesInUse'][0]
		302	self.DualModeIndex = header['DualModeIndex'][0]
		303	self.DualModeRange = header['DualModeRange'][0]
		304	self.nDigChannels = header['nDigChannels'][0]
		305	self.SampResolution = header['SampResolution'][0]
		306	self.nRangeGatesSamp = header['nRangeGatesSamp'][0]
		307	self.StartRangeSamp = header['StartRangeSamp'][0]
		308	self.PRFhz = header['PRFhz'][0]
		309	self.Integrations = header['Integrations'][0]
		310	self.nDataPointsTrsf = header['nDataPointsTrsf'][0]
		311	self.nReceiveBeams = header['nReceiveBeams'][0]
		312	self.nSpectAverages = header['nSpectAverages'][0]
		313	self.FFTwindowingInd = header['FFTwindowingInd'][0]
		314	self.BeamAngleAzim = header['BeamAngleAzim'][0]
		315	self.BeamAngleZen = header['BeamAngleZen'][0]
		316	self.AntennaCoord = header['AntennaCoord'][0]
		317	self.RecPhaseCalibr = header['RecPhaseCalibr'][0]
		318	self.RecAmpCalibr = header['RecAmpCalibr'][0]
		319	self.ReceiverGaindB = header['ReceiverGaindB'][0]
		320
		321	Self.size = 180+20*3
		322
		323	endFp = self.size + startFp
		324
		325	if fp.tell() > endFp:
		326	sys.stderr.write("Warning %s: Size value read from System Header is lower than it has to be\n" %fp.name)
		327	return 0
		328
		329	if fp.tell() < endFp:
		330	sys.stderr.write("Warning %s: Size value read from System Header size is greater than it has to be\n" %fp.name)
		331	return 0
		332
		333	return 1
		334
		335	def write(self, fp):
		336
		337	headerTuple = (self.RecMgcNumber,
		338	self.RecCounter,
		339	self.Off2StartNxtRec,
		340	self.EpTimeStamp,
		341	self.msCompTimeStamp,
		342	self.ExpTagName,
		343	self.ExpComment,
		344	self.SiteLatDegrees,
		345	self.SiteLongDegrees,
		346	self.RTCgpsStatus,
		347	self.TransmitFrec,
		348	self.ReceiveFrec,
		349	self.FirstOsciFrec,
		350	self.Polarisation,
		351	self.ReceiverFiltSett,
		352	self.nModesInUse,
		353	self.DualModeIndex,
		354	self.DualModeRange,
		355	self.nDigChannels,
		356	self.SampResolution,
		357	self.nRangeGatesSamp,
		358	self.StartRangeSamp,
		359	self.PRFhz,
		360	self.Integrations,
		361	self.nDataPointsTrsf,
		362	self.nReceiveBeams,
		363	self.nSpectAverages,
		364	self.FFTwindowingInd,
		365	self.BeamAngleAzim,
		366	self.BeamAngleZen,
		367	self.AntennaCoord,
		368	self.RecPhaseCalibr,
		369	self.RecAmpCalibr,
		370	self.ReceiverGaindB)
		371
		372	# self.size,self.nSamples,
		373	# self.nProfiles,
		374	# self.nChannels,
		375	# self.adcResolution,
		376	# self.pciDioBusWidth
		377
		378	header = numpy.array(headerTuple,RECORD_STRUCTURE)
		379	header.tofile(fp)
		380
		381	return 1
		382
		383
		384	def get_dtype_index(numpy_dtype):
		385
		386	index = None
		387
		388	for i in range(len(NUMPY_DTYPE_LIST)):
		389	if numpy_dtype == NUMPY_DTYPE_LIST[i]:
		390	index = i
		391	break
		392
		393	return index
		394
		395	def get_numpy_dtype(index):
		396
		397	#dtype4 = numpy.dtype([('real','<f4'),('imag','<f4')])
		398
		399	return NUMPY_DTYPE_LIST[index]
		400
		401
		402	def get_dtype_width(index):
		403
		404	return DTYPE_WIDTH[index] No newline at end of file

schainpy/model/graphics/mpldriver2.py created 644 +469 0

			@@ -0,0 +1,469
		1	import numpy
		2	import datetime
		3	import sys
		4	import matplotlib
		5
		6	if 'linux' in sys.platform:
		7	matplotlib.use("TKAgg")
		8
		9	if 'darwin' in sys.platform:
		10	matplotlib.use('TKAgg')
		11	#Qt4Agg', 'GTK', 'GTKAgg', 'ps', 'agg', 'cairo', 'MacOSX', 'GTKCairo', 'WXAgg', 'template', 'TkAgg', 'GTK3Cairo', 'GTK3Agg', 'svg', 'WebAgg', 'CocoaAgg', 'emf', 'gdk', 'WX'
		12	import matplotlib.pyplot
		13
		14	from mpl_toolkits.axes_grid1 import make_axes_locatable
		15	from matplotlib.ticker import FuncFormatter, LinearLocator
		16
		17	###########################################
		18	#Actualizacion de las funciones del driver
		19	###########################################
		20
		21	jet_values = matplotlib.pyplot.get_cmap("jet", 100)(numpy.arange(100))[10:90]
		22	blu_values = matplotlib.pyplot.get_cmap("seismic_r", 20)(numpy.arange(20))[10:15]
		23	ncmap = matplotlib.colors.LinearSegmentedColormap.from_list("jro", numpy.vstack((blu_values, jet_values)))
		24	matplotlib.pyplot.register_cmap(cmap=ncmap)
		25
		26	def createFigure(id, wintitle, width, height, facecolor="w", show=True, dpi = 80):
		27
		28	matplotlib.pyplot.ioff()
		29
		30	fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor, figsize=(1.0width/dpi, 1.0height/dpi))
		31	fig.canvas.manager.set_window_title(wintitle)
		32	# fig.canvas.manager.resize(width, height)
		33	matplotlib.pyplot.ion()
		34
		35
		36	if show:
		37	matplotlib.pyplot.show()
		38
		39	return fig
		40
		41	def closeFigure(show=False, fig=None):
		42
		43	# matplotlib.pyplot.ioff()
		44	# matplotlib.pyplot.pause(0)
		45
		46	if show:
		47	matplotlib.pyplot.show()
		48
		49	if fig != None:
		50	matplotlib.pyplot.close(fig)
		51	# matplotlib.pyplot.pause(0)
		52	# matplotlib.pyplot.ion()
		53
		54	return
		55
		56	matplotlib.pyplot.close("all")
		57	# matplotlib.pyplot.pause(0)
		58	# matplotlib.pyplot.ion()
		59
		60	return
		61
		62	def saveFigure(fig, filename):
		63
		64	# matplotlib.pyplot.ioff()
		65	fig.savefig(filename, dpi=matplotlib.pyplot.gcf().dpi)
		66	# matplotlib.pyplot.ion()
		67
		68	def clearFigure(fig):
		69
		70	fig.clf()
		71
		72	def setWinTitle(fig, title):
		73
		74	fig.canvas.manager.set_window_title(title)
		75
		76	def setTitle(fig, title):
		77
		78	fig.suptitle(title)
		79
		80	def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
		81
		82	matplotlib.pyplot.ioff()
		83	matplotlib.pyplot.figure(fig.number)
		84	axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
		85	(xpos, ypos),
		86	colspan=colspan,
		87	rowspan=rowspan,
		88	polar=polar)
		89
		90	axes.grid(True)
		91	matplotlib.pyplot.ion()
		92	return axes
		93
		94	def setAxesText(ax, text):
		95
		96	ax.annotate(text,
		97	xy = (.1, .99),
		98	xycoords = 'figure fraction',
		99	horizontalalignment = 'left',
		100	verticalalignment = 'top',
		101	fontsize = 10)
		102
		103	def printLabels(ax, xlabel, ylabel, title):
		104
		105	ax.set_xlabel(xlabel, size=11)
		106	ax.set_ylabel(ylabel, size=11)
		107	ax.set_title(title, size=8)
		108
		109	def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='',
		110	ticksize=9, xtick_visible=True, ytick_visible=True,
		111	nxticks=4, nyticks=10,
		112	grid=None,color='blue'):
		113
		114	"""
		115
		116	Input:
		117	grid : None, 'both', 'x', 'y'
		118	"""
		119
		120	matplotlib.pyplot.ioff()
		121
		122	ax.set_xlim([xmin,xmax])
		123	ax.set_ylim([ymin,ymax])
		124
		125	printLabels(ax, xlabel, ylabel, title)
		126
		127	######################################################
		128	if (xmax-xmin)<=1:
		129	xtickspos = numpy.linspace(xmin,xmax,nxticks)
		130	xtickspos = numpy.array([float("%.1f"%i) for i in xtickspos])
		131	ax.set_xticks(xtickspos)
		132	else:
		133	xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
		134	# xtickspos = numpy.arange(nxticks)*float(xmax-xmin)/float(nxticks) + int(xmin)
		135	ax.set_xticks(xtickspos)
		136
		137	for tick in ax.get_xticklabels():
		138	tick.set_visible(xtick_visible)
		139
		140	for tick in ax.xaxis.get_major_ticks():
		141	tick.label.set_fontsize(ticksize)
		142
		143	######################################################
		144	for tick in ax.get_yticklabels():
		145	tick.set_visible(ytick_visible)
		146
		147	for tick in ax.yaxis.get_major_ticks():
		148	tick.label.set_fontsize(ticksize)
		149
		150	ax.plot(x, y, color=color)
		151	iplot = ax.lines[-1]
		152
		153	######################################################
		154	if '0.' in matplotlib.__version__[0:2]:
		155	print "The matplotlib version has to be updated to 1.1 or newer"
		156	return iplot
		157
		158	if '1.0.' in matplotlib.__version__[0:4]:
		159	print "The matplotlib version has to be updated to 1.1 or newer"
		160	return iplot
		161
		162	if grid != None:
		163	ax.grid(b=True, which='major', axis=grid)
		164
		165	matplotlib.pyplot.tight_layout()
		166
		167	matplotlib.pyplot.ion()
		168
		169	return iplot
		170
		171	def set_linedata(ax, x, y, idline):
		172
		173	ax.lines[idline].set_data(x,y)
		174
		175	def pline(iplot, x, y, xlabel='', ylabel='', title=''):
		176
		177	ax = iplot.get_axes()
		178
		179	printLabels(ax, xlabel, ylabel, title)
		180
		181	set_linedata(ax, x, y, idline=0)
		182
		183	def addpline(ax, x, y, color, linestyle, lw):
		184
		185	ax.plot(x,y,color=color,linestyle=linestyle,lw=lw)
		186
		187
		188	def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax,
		189	xlabel='', ylabel='', title='', ticksize = 9,
		190	colormap='jet',cblabel='', cbsize="5%",
		191	XAxisAsTime=False):
		192
		193	matplotlib.pyplot.ioff()
		194
		195	divider = make_axes_locatable(ax)
		196	ax_cb = divider.new_horizontal(size=cbsize, pad=0.05)
		197	fig = ax.get_figure()
		198	fig.add_axes(ax_cb)
		199
		200	ax.set_xlim([xmin,xmax])
		201	ax.set_ylim([ymin,ymax])
		202
		203	printLabels(ax, xlabel, ylabel, title)
		204
		205	z = numpy.ma.masked_invalid(z)
		206	cmap=matplotlib.pyplot.get_cmap(colormap)
		207	cmap.set_bad('white',1.)
		208	imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=cmap)
		209	cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
		210	cb.set_label(cblabel)
		211
		212	# for tl in ax_cb.get_yticklabels():
		213	# tl.set_visible(True)
		214
		215	for tick in ax.yaxis.get_major_ticks():
		216	tick.label.set_fontsize(ticksize)
		217
		218	for tick in ax.xaxis.get_major_ticks():
		219	tick.label.set_fontsize(ticksize)
		220
		221	for tick in cb.ax.get_yticklabels():
		222	tick.set_fontsize(ticksize)
		223
		224	ax_cb.yaxis.tick_right()
		225
		226	if '0.' in matplotlib.__version__[0:2]:
		227	print "The matplotlib version has to be updated to 1.1 or newer"
		228	return imesh
		229
		230	if '1.0.' in matplotlib.__version__[0:4]:
		231	print "The matplotlib version has to be updated to 1.1 or newer"
		232	return imesh
		233
		234	matplotlib.pyplot.tight_layout()
		235
		236	if XAxisAsTime:
		237
		238	func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
		239	ax.xaxis.set_major_formatter(FuncFormatter(func))
		240	ax.xaxis.set_major_locator(LinearLocator(7))
		241	ax.grid(True)
		242	matplotlib.pyplot.ion()
		243	return imesh
		244
		245	def pcolor(imesh, z, xlabel='', ylabel='', title=''):
		246
		247	z = z.T
		248	ax = imesh.get_axes()
		249	printLabels(ax, xlabel, ylabel, title)
		250	imesh.set_array(z.ravel())
		251	ax.grid(True)
		252
		253	def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
		254
		255	printLabels(ax, xlabel, ylabel, title)
		256	ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
		257	ax.grid(True)
		258
		259	def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
		260
		261	printLabels(ax, xlabel, ylabel, title)
		262
		263	ax.collections.remove(ax.collections[0])
		264
		265	z = numpy.ma.masked_invalid(z)
		266
		267	cmap=matplotlib.pyplot.get_cmap(colormap)
		268	cmap.set_bad('white',1.)
		269
		270	ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=cmap)
		271	ax.grid(True)
		272
		273	def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
		274	ticksize=9, xtick_visible=True, ytick_visible=True,
		275	nxticks=4, nyticks=10,
		276	grid=None):
		277
		278	"""
		279
		280	Input:
		281	grid : None, 'both', 'x', 'y'
		282	"""
		283
		284	matplotlib.pyplot.ioff()
		285
		286	lines = ax.plot(x.T, y)
		287	leg = ax.legend(lines, legendlabels, loc='upper right')
		288	leg.get_frame().set_alpha(0.5)
		289	ax.set_xlim([xmin,xmax])
		290	ax.set_ylim([ymin,ymax])
		291	printLabels(ax, xlabel, ylabel, title)
		292
		293	xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
		294	ax.set_xticks(xtickspos)
		295
		296	for tick in ax.get_xticklabels():
		297	tick.set_visible(xtick_visible)
		298
		299	for tick in ax.xaxis.get_major_ticks():
		300	tick.label.set_fontsize(ticksize)
		301
		302	for tick in ax.get_yticklabels():
		303	tick.set_visible(ytick_visible)
		304
		305	for tick in ax.yaxis.get_major_ticks():
		306	tick.label.set_fontsize(ticksize)
		307
		308	iplot = ax.lines[-1]
		309
		310	if '0.' in matplotlib.__version__[0:2]:
		311	print "The matplotlib version has to be updated to 1.1 or newer"
		312	return iplot
		313
		314	if '1.0.' in matplotlib.__version__[0:4]:
		315	print "The matplotlib version has to be updated to 1.1 or newer"
		316	return iplot
		317
		318	if grid != None:
		319	ax.grid(b=True, which='major', axis=grid)
		320
		321	matplotlib.pyplot.tight_layout()
		322
		323	matplotlib.pyplot.ion()
		324
		325	return iplot
		326
		327
		328	def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''):
		329
		330	ax = iplot.get_axes()
		331
		332	printLabels(ax, xlabel, ylabel, title)
		333
		334	for i in range(len(ax.lines)):
		335	line = ax.lines[i]
		336	line.set_data(x[i,:],y)
		337
		338	def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
		339	ticksize=9, xtick_visible=True, ytick_visible=True,
		340	nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
		341	grid=None, XAxisAsTime=False):
		342
		343	"""
		344
		345	Input:
		346	grid : None, 'both', 'x', 'y'
		347	"""
		348
		349	matplotlib.pyplot.ioff()
		350
		351	# lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
		352	lines = ax.plot(x, y.T)
		353	# leg = ax.legend(lines, legendlabels, loc=2, bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
		354	# handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)
		355
		356	leg = ax.legend(lines, legendlabels,
		357	loc='upper right', bbox_to_anchor=(1.16, 1), borderaxespad=0)
		358
		359	for label in leg.get_texts(): label.set_fontsize(9)
		360
		361	ax.set_xlim([xmin,xmax])
		362	ax.set_ylim([ymin,ymax])
		363	printLabels(ax, xlabel, ylabel, title)
		364
		365	# xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
		366	# ax.set_xticks(xtickspos)
		367
		368	for tick in ax.get_xticklabels():
		369	tick.set_visible(xtick_visible)
		370
		371	for tick in ax.xaxis.get_major_ticks():
		372	tick.label.set_fontsize(ticksize)
		373
		374	for tick in ax.get_yticklabels():
		375	tick.set_visible(ytick_visible)
		376
		377	for tick in ax.yaxis.get_major_ticks():
		378	tick.label.set_fontsize(ticksize)
		379
		380	iplot = ax.lines[-1]
		381
		382	if '0.' in matplotlib.__version__[0:2]:
		383	print "The matplotlib version has to be updated to 1.1 or newer"
		384	return iplot
		385
		386	if '1.0.' in matplotlib.__version__[0:4]:
		387	print "The matplotlib version has to be updated to 1.1 or newer"
		388	return iplot
		389
		390	if grid != None:
		391	ax.grid(b=True, which='major', axis=grid)
		392
		393	matplotlib.pyplot.tight_layout()
		394
		395	if XAxisAsTime:
		396
		397	func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
		398	ax.xaxis.set_major_formatter(FuncFormatter(func))
		399	ax.xaxis.set_major_locator(LinearLocator(7))
		400
		401	matplotlib.pyplot.ion()
		402
		403	return iplot
		404
		405	def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
		406
		407	ax = iplot.get_axes()
		408
		409	printLabels(ax, xlabel, ylabel, title)
		410
		411	for i in range(len(ax.lines)):
		412	line = ax.lines[i]
		413	line.set_data(x,y[i,:])
		414
		415	def createPolar(ax, x, y,
		416	xlabel='', ylabel='', title='', ticksize = 9,
		417	colormap='jet',cblabel='', cbsize="5%",
		418	XAxisAsTime=False):
		419
		420	matplotlib.pyplot.ioff()
		421
		422	ax.plot(x,y,'bo', markersize=5)
		423	# ax.set_rmax(90)
		424	ax.set_ylim(0,90)
		425	ax.set_yticks(numpy.arange(0,90,20))
		426	# ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
		427	# ax.text(0, 50, ylabel, rotation='vertical', va ='center', ha = 'left' ,size='11')
		428	# ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
		429	ax.yaxis.labelpad = 230
		430	printLabels(ax, xlabel, ylabel, title)
		431	iplot = ax.lines[-1]
		432
		433	if '0.' in matplotlib.__version__[0:2]:
		434	print "The matplotlib version has to be updated to 1.1 or newer"
		435	return iplot
		436
		437	if '1.0.' in matplotlib.__version__[0:4]:
		438	print "The matplotlib version has to be updated to 1.1 or newer"
		439	return iplot
		440
		441	# if grid != None:
		442	# ax.grid(b=True, which='major', axis=grid)
		443
		444	matplotlib.pyplot.tight_layout()
		445
		446	matplotlib.pyplot.ion()
		447
		448
		449	return iplot
		450
		451	def polar(iplot, x, y, xlabel='', ylabel='', title=''):
		452
		453	ax = iplot.get_axes()
		454
		455	# ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
		456	printLabels(ax, xlabel, ylabel, title)
		457
		458	set_linedata(ax, x, y, idline=0)
		459
		460	def draw(fig):
		461
		462	if type(fig) == 'int':
		463	raise ValueError, "Error drawing: Fig parameter should be a matplotlib figure object figure"
		464
		465	fig.canvas.draw()
		466
		467	def pause(interval=0.000001):
		468
		469	matplotlib.pyplot.pause(interval)

schainpy/model/io/MIRAtest.py created 644 +321 0

			@@ -0,0 +1,321
		1	import os, sys
		2	import glob
		3	import fnmatch
		4	import datetime
		5	import time
		6	import re
		7	import h5py
		8	import numpy
		9	import matplotlib.pyplot as plt
		10
		11	import pylab as plb
		12	from scipy.optimize import curve_fit
		13	from scipy import asarray as ar,exp
		14	from scipy import stats
		15
		16	from duplicity.path import Path
		17	from numpy.ma.core import getdata
		18
		19	SPEED_OF_LIGHT = 299792458
		20	SPEED_OF_LIGHT = 3e8
		21
		22	try:
		23	from gevent import sleep
		24	except:
		25	from time import sleep
		26
		27	from schainpy.model.data.jrodata import Spectra
		28	#from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader
		29	from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
		30	#from schainpy.model.io.jroIO_bltr import BLTRReader
		31	from numpy import imag, shape, NaN
		32
		33
		34	startFp = open('/home/erick/Documents/MIRA35C/20160117/20160117_0000.zspc',"rb")
		35
		36
		37	FILE_HEADER = numpy.dtype([ #HEADER 1024bytes
		38	('Hname',numpy.str_,32), #Original file name
		39	('Htime',numpy.str_,32), #Date and time when the file was created
		40	('Hoper',numpy.str_,64), #Name of operator who created the file
		41	('Hplace',numpy.str_,128), #Place where the measurements was carried out
		42	('Hdescr',numpy.str_,256), #Description of measurements
		43	('Hdummy',numpy.str_,512), #Reserved space
		44	#Main chunk
		45	('Msign','<i4'), #Main chunk signature FZKF or NUIG
		46	('MsizeData','<i4'), #Size of data block main chunk
		47	#Processing DSP parameters
		48	('PPARsign','<i4'), #PPAR signature
		49	('PPARsize','<i4'), #PPAR size of block
		50	('PPARprf','<i4'), #Pulse repetition frequency
		51	('PPARpdr','<i4'), #Pulse duration
		52	('PPARsft','<i4'), #FFT length
		53	('PPARavc','<i4'), #Number of spectral (in-coherent) averages
		54	('PPARihp','<i4'), #Number of lowest range gate for moment estimation
		55	('PPARchg','<i4'), #Count for gates for moment estimation
		56	('PPARpol','<i4'), #switch on/off polarimetric measurements. Should be 1.
		57	#Service DSP parameters
		58	('SPARatt','<i4'), #STC attenuation on the lowest ranges on/off
		59	('SPARtx','<i4'), #OBSOLETE
		60	('SPARaddGain0','<f4'), #OBSOLETE
		61	('SPARaddGain1','<f4'), #OBSOLETE
		62	('SPARwnd','<i4'), #Debug only. It normal mode it is 0.
		63	('SPARpos','<i4'), #Delay between sync pulse and tx pulse for phase corr, ns
		64	('SPARadd','<i4'), #"add to pulse" to compensate for delay between the leading edge of driver pulse and envelope of the RF signal.
		65	('SPARlen','<i4'), #Time for measuring txn pulse phase. OBSOLETE
		66	('SPARcal','<i4'), #OBSOLETE
		67	('SPARnos','<i4'), #OBSOLETE
		68	('SPARof0','<i4'), #detection threshold
		69	('SPARof1','<i4'), #OBSOLETE
		70	('SPARswt','<i4'), #2nd moment estimation threshold
		71	('SPARsum','<i4'), #OBSOLETE
		72	('SPARosc','<i4'), #flag Oscillosgram mode
		73	('SPARtst','<i4'), #OBSOLETE
		74	('SPARcor','<i4'), #OBSOLETE
		75	('SPARofs','<i4'), #OBSOLETE
		76	('SPARhsn','<i4'), #Hildebrand div noise detection on noise gate
		77	('SPARhsa','<f4'), #Hildebrand div noise detection on all gates
		78	('SPARcalibPow_M','<f4'), #OBSOLETE
		79	('SPARcalibSNR_M','<f4'), #OBSOLETE
		80	('SPARcalibPow_S','<f4'), #OBSOLETE
		81	('SPARcalibSNR_S','<f4'), #OBSOLETE
		82	('SPARrawGate1','<i4'), #Lowest range gate for spectra saving Raw_Gate1 >=5
		83	('SPARrawGate2','<i4'), #Number of range gates with atmospheric signal
		84	('SPARraw','<i4'), #flag - IQ or spectra saving on/off
		85	('SPARprc','<i4'),]) #flag - Moment estimation switched on/off
		86
		87
		88
		89	self.Hname= None
		90	self.Htime= None
		91	self.Hoper= None
		92	self.Hplace= None
		93	self.Hdescr= None
		94	self.Hdummy= None
		95
		96	self.Msign=None
		97	self.MsizeData=None
		98
		99	self.PPARsign=None
		100	self.PPARsize=None
		101	self.PPARprf=None
		102	self.PPARpdr=None
		103	self.PPARsft=None
		104	self.PPARavc=None
		105	self.PPARihp=None
		106	self.PPARchg=None
		107	self.PPARpol=None
		108	#Service DSP parameters
		109	self.SPARatt=None
		110	self.SPARtx=None
		111	self.SPARaddGain0=None
		112	self.SPARaddGain1=None
		113	self.SPARwnd=None
		114	self.SPARpos=None
		115	self.SPARadd=None
		116	self.SPARlen=None
		117	self.SPARcal=None
		118	self.SPARnos=None
		119	self.SPARof0=None
		120	self.SPARof1=None
		121	self.SPARswt=None
		122	self.SPARsum=None
		123	self.SPARosc=None
		124	self.SPARtst=None
		125	self.SPARcor=None
		126	self.SPARofs=None
		127	self.SPARhsn=None
		128	self.SPARhsa=None
		129	self.SPARcalibPow_M=None
		130	self.SPARcalibSNR_M=None
		131	self.SPARcalibPow_S=None
		132	self.SPARcalibSNR_S=None
		133	self.SPARrawGate1=None
		134	self.SPARrawGate2=None
		135	self.SPARraw=None
		136	self.SPARprc=None
		137
		138
		139
		140	header = numpy.fromfile(fp, FILE_HEADER,1)
		141	''' numpy.fromfile(file, dtype, count, sep='')
		142	file : file or str
		143	Open file object or filename.
		144
		145	dtype : data-type
		146	Data type of the returned array. For binary files, it is used to determine
		147	the size and byte-order of the items in the file.
		148
		149	count : int
		150	Number of items to read. -1 means all items (i.e., the complete file).
		151
		152	sep : str
		153	Separator between items if file is a text file. Empty ("") separator means
		154	the file should be treated as binary. Spaces (" ") in the separator match zero
		155	or more whitespace characters. A separator consisting only of spaces must match
		156	at least one whitespace.
		157
		158	'''
		159
		160	Hname= str(header['Hname'][0])
		161	Htime= str(header['Htime'][0])
		162	Hoper= str(header['Hoper'][0])
		163	Hplace= str(header['Hplace'][0])
		164	Hdescr= str(header['Hdescr'][0])
		165	Hdummy= str(header['Hdummy'][0])
		166
		167	Msign=header['Msign'][0]
		168	MsizeData=header['MsizeData'][0]
		169
		170	PPARsign=header['PPARsign'][0]
		171	PPARsize=header['PPARsize'][0]
		172	PPARprf=header['PPARprf'][0]
		173	PPARpdr=header['PPARpdr'][0]
		174	PPARsft=header['PPARsft'][0]
		175	PPARavc=header['PPARavc'][0]
		176	PPARihp=header['PPARihp'][0]
		177	PPARchg=header['PPARchg'][0]
		178	PPARpol=header['PPARpol'][0]
		179	#Service DSP parameters
		180	SPARatt=header['SPARatt'][0]
		181	SPARtx=header['SPARtx'][0]
		182	SPARaddGain0=header['SPARaddGain0'][0]
		183	SPARaddGain1=header['SPARaddGain1'][0]
		184	SPARwnd=header['SPARwnd'][0]
		185	SPARpos=header['SPARpos'][0]
		186	SPARadd=header['SPARadd'][0]
		187	SPARlen=header['SPARlen'][0]
		188	SPARcal=header['SPARcal'][0]
		189	SPARnos=header['SPARnos'][0]
		190	SPARof0=header['SPARof0'][0]
		191	SPARof1=header['SPARof1'][0]
		192	SPARswt=header['SPARswt'][0]
		193	SPARsum=header['SPARsum'][0]
		194	SPARosc=header['SPARosc'][0]
		195	SPARtst=header['SPARtst'][0]
		196	SPARcor=header['SPARcor'][0]
		197	SPARofs=header['SPARofs'][0]
		198	SPARhsn=header['SPARhsn'][0]
		199	SPARhsa=header['SPARhsa'][0]
		200	SPARcalibPow_M=header['SPARcalibPow_M'][0]
		201	SPARcalibSNR_M=header['SPARcalibSNR_M'][0]
		202	SPARcalibPow_S=header['SPARcalibPow_S'][0]
		203	SPARcalibSNR_S=header['SPARcalibSNR_S'][0]
		204	SPARrawGate1=header['SPARrawGate1'][0]
		205	SPARrawGate2=header['SPARrawGate2'][0]
		206	SPARraw=header['SPARraw'][0]
		207	SPARprc=header['SPARprc'][0]
		208
		209
		210
		211	SRVI_STRUCTURE = numpy.dtype([
		212	('frame_cnt','<u4'),#
		213	('time_t','<u4'), #
		214	('tpow','<f4'), #
		215	('npw1','<f4'), #
		216	('npw2','<f4'), #
		217	('cpw1','<f4'), #
		218	('pcw2','<f4'), #
		219	('ps_err','<u4'), #
		220	('te_err','<u4'), #
		221	('rc_err','<u4'), #
		222	('grs1','<u4'), #
		223	('grs2','<u4'), #
		224	('azipos','<f4'), #
		225	('azivel','<f4'), #
		226	('elvpos','<f4'), #
		227	('elvvel','<f4'), #
		228	('northAngle','<f4'), #
		229	('microsec','<u4'), #
		230	('azisetvel','<f4'), #
		231	('elvsetpos','<f4'), #
		232	('RadarConst','<f4'),]) #
		233
		234	JUMP_STRUCTURE = numpy.dtype([
		235	('jump','<u140'),#
		236	('SizeOfDataBlock1',numpy.str_,32),#
		237	('jump','<i4'),#
		238	('DataBlockTitleSRVI1',numpy.str_,32),#
		239	('SizeOfSRVI1','<i4'),])#
		240
		241
		242
		243	#frame_cnt=0, time_t= 0, tpow=0, npw1=0, npw2=0,
		244	#cpw1=0, pcw2=0, ps_err=0, te_err=0, rc_err=0, grs1=0,
		245	#grs2=0, azipos=0, azivel=0, elvpos=0, elvvel=0, northangle=0,
		246	#microsec=0, azisetvel=0, elvsetpos=0, RadarConst=0
		247
		248
		249	frame_cnt = frame_cnt
		250	dwell = time_t
		251	tpow = tpow
		252	npw1 = npw1
		253	npw2 = npw2
		254	cpw1 = cpw1
		255	pcw2 = pcw2
		256	ps_err = ps_err
		257	te_err = te_err
		258	rc_err = rc_err
		259	grs1 = grs1
		260	grs2 = grs2
		261	azipos = azipos
		262	azivel = azivel
		263	elvpos = elvpos
		264	elvvel = elvvel
		265	northAngle = northAngle
		266	microsec = microsec
		267	azisetvel = azisetvel
		268	elvsetpos = elvsetpos
		269	RadarConst5 = RadarConst
		270
		271
		272
		273	#print fp
		274	#startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
		275	#startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
		276	#RecCounter=0
		277	#Off2StartNxtRec=811248
		278	#print 'OffsetStartHeader ',self.OffsetStartHeader,'RecCounter ', self.RecCounter, 'Off2StartNxtRec ' , self.Off2StartNxtRec
		279	#OffRHeader= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
		280	#startFp.seek(OffRHeader, os.SEEK_SET)
		281	print 'debe ser 48, RecCounter*811248', self.OffsetStartHeader,self.RecCounter,self.Off2StartNxtRec
		282	print 'Posicion del bloque: ',OffRHeader
		283
		284	header = numpy.fromfile(startFp,SRVI_STRUCTURE,1)
		285
		286	self.frame_cnt = header['frame_cnt'][0]#
		287	self.time_t = header['frame_cnt'][0] #
		288	self.tpow = header['frame_cnt'][0] #
		289	self.npw1 = header['frame_cnt'][0] #
		290	self.npw2 = header['frame_cnt'][0] #
		291	self.cpw1 = header['frame_cnt'][0] #
		292	self.pcw2 = header['frame_cnt'][0] #
		293	self.ps_err = header['frame_cnt'][0] #
		294	self.te_err = header['frame_cnt'][0] #
		295	self.rc_err = header['frame_cnt'][0] #
		296	self.grs1 = header['frame_cnt'][0] #
		297	self.grs2 = header['frame_cnt'][0] #
		298	self.azipos = header['frame_cnt'][0] #
		299	self.azivel = header['frame_cnt'][0] #
		300	self.elvpos = header['frame_cnt'][0] #
		301	self.elvvel = header['frame_cnt'][0] #
		302	self.northAngle = header['frame_cnt'][0] #
		303	self.microsec = header['frame_cnt'][0] #
		304	self.azisetvel = header['frame_cnt'][0] #
		305	self.elvsetpos = header['frame_cnt'][0] #
		306	self.RadarConst = header['frame_cnt'][0] #
		307
		308
		309	self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
		310
		311	self.RHsize = 180+20*self.nChannels
		312	self.Datasize= self.nProfilesself.nChannelsself.nHeights24
		313	#print 'Datasize',self.Datasize
		314	endFp = self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
		315
		316	print '=============================================='
		317
		318	print '=============================================='
		319
		320
		321	No newline at end of file

schainpy/model/io/bltrIO_param.py created 644 +362 0

			@@ -0,0 +1,362
		1	'''
		2	Created on Nov 9, 2016
		3
		4	@author: roj- LouVD
		5	'''
		6
		7
		8	import os
		9	import sys
		10	import time
		11	import glob
		12	import datetime
		13
		14	import numpy
		15
		16	from schainpy.model.proc.jroproc_base import ProcessingUnit
		17	from schainpy.model.data.jrodata import Parameters
		18	from schainpy.model.io.jroIO_base import JRODataReader, isNumber
		19
		20	FILE_HEADER_STRUCTURE = numpy.dtype([
		21	('FMN', '<u4'),
		22	('nrec', '<u4'),
		23	('fr_offset', '<u4'),
		24	('id', '<u4'),
		25	('site', 'u1', (32,))
		26	])
		27
		28	REC_HEADER_STRUCTURE = numpy.dtype([
		29	('rmn', '<u4'),
		30	('rcounter', '<u4'),
		31	('nr_offset', '<u4'),
		32	('tr_offset', '<u4'),
		33	('time', '<u4'),
		34	('time_msec', '<u4'),
		35	('tag', 'u1', (32,)),
		36	('comments', 'u1', (32,)),
		37	('lat', '<f4'),
		38	('lon', '<f4'),
		39	('gps_status', '<u4'),
		40	('freq', '<u4'),
		41	('freq0', '<u4'),
		42	('nchan', '<u4'),
		43	('delta_r', '<u4'),
		44	('nranges', '<u4'),
		45	('r0', '<u4'),
		46	('prf', '<u4'),
		47	('ncoh', '<u4'),
		48	('npoints', '<u4'),
		49	('polarization', '<i4'),
		50	('rx_filter', '<u4'),
		51	('nmodes', '<u4'),
		52	('dmode_index', '<u4'),
		53	('dmode_rngcorr', '<u4'),
		54	('nrxs', '<u4'),
		55	('acf_length', '<u4'),
		56	('acf_lags', '<u4'),
		57	('sea_to_atmos', '<f4'),
		58	('sea_notch', '<u4'),
		59	('lh_sea', '<u4'),
		60	('hh_sea', '<u4'),
		61	('nbins_sea', '<u4'),
		62	('min_snr', '<f4'),
		63	('min_cc', '<f4'),
		64	('max_time_diff', '<f4')
		65	])
		66
		67	DATA_STRUCTURE = numpy.dtype([
		68	('range', '<u4'),
		69	('status', '<u4'),
		70	('zonal', '<f4'),
		71	('meridional', '<f4'),
		72	('vertical', '<f4'),
		73	('zonal_a', '<f4'),
		74	('meridional_a', '<f4'),
		75	('corrected_fading', '<f4'), # seconds
		76	('uncorrected_fading', '<f4'), # seconds
		77	('time_diff', '<f4'),
		78	('major_axis', '<f4'),
		79	('axial_ratio', '<f4'),
		80	('orientation', '<f4'),
		81	('sea_power', '<u4'),
		82	('sea_algorithm', '<u4')
		83	])
		84
		85	class BLTRParamReader(JRODataReader, ProcessingUnit):
		86	'''
		87	Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR from *.sswma files
		88	'''
		89
		90	ext = '.sswma'
		91
		92	def __init__(self, **kwargs):
		93
		94	ProcessingUnit.__init__(self , **kwargs)
		95
		96	self.dataOut = Parameters()
		97	self.counter_records = 0
		98	self.flagNoMoreFiles = 0
		99	self.isConfig = False
		100	self.filename = None
		101
		102	def setup(self,
		103	path=None,
		104	startDate=None,
		105	endDate=None,
		106	ext=None,
		107	startTime=datetime.time(0, 0, 0),
		108	endTime=datetime.time(23, 59, 59),
		109	timezone=0,
		110	status_value=0,
		111	**kwargs):
		112
		113	self.path = path
		114	self.startTime = startTime
		115	self.endTime = endTime
		116	self.status_value = status_value
		117
		118	if self.path is None:
		119	raise ValueError, "The path is not valid"
		120
		121	if ext is None:
		122	ext = self.ext
		123
		124	self.search_files(self.path, startDate, endDate, ext)
		125	self.timezone = timezone
		126	self.fileIndex = 0
		127
		128	if not self.fileList:
		129	raise Warning, "There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' "%(path)
		130
		131	self.setNextFile()
		132
		133	def search_files(self, path, startDate, endDate, ext):
		134	'''
		135	Searching for BLTR rawdata file in path
		136	Creating a list of file to proces included in [startDate,endDate]
		137
		138	Input:
		139	path - Path to find BLTR rawdata files
		140	startDate - Select file from this date
		141	enDate - Select file until this date
		142	ext - Extension of the file to read
		143
		144	'''
		145
		146	print 'Searching file in %s ' % (path)
		147	foldercounter = 0
		148	fileList0 = glob.glob1(path, "*%s" % ext)
		149	fileList0.sort()
		150
		151	self.fileList = []
		152	self.dateFileList = []
		153
		154	for thisFile in fileList0:
		155	year = thisFile[-14:-10]
		156	if not isNumber(year):
		157	continue
		158
		159	month = thisFile[-10:-8]
		160	if not isNumber(month):
		161	continue
		162
		163	day = thisFile[-8:-6]
		164	if not isNumber(day):
		165	continue
		166
		167	year, month, day = int(year), int(month), int(day)
		168	dateFile = datetime.date(year, month, day)
		169
		170	if (startDate > dateFile) or (endDate < dateFile):
		171	continue
		172
		173	self.fileList.append(thisFile)
		174	self.dateFileList.append(dateFile)
		175
		176	return
		177
		178	def setNextFile(self):
		179
		180	file_id = self.fileIndex
		181
		182	if file_id == len(self.fileList):
		183	print '\nNo more files in the folder'
		184	print 'Total number of file(s) read : {}'.format(self.fileIndex + 1)
		185	self.flagNoMoreFiles = 1
		186	return 0
		187
		188	print '\n[Setting file] (%s) ...' % self.fileList[file_id]
		189	filename = os.path.join(self.path, self.fileList[file_id])
		190
		191	dirname, name = os.path.split(filename)
		192	self.siteFile = name.split('.')[0] # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya
		193	if self.filename is not None:
		194	self.fp.close()
		195	self.filename = filename
		196	self.fp = open(self.filename, 'rb')
		197	self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
		198	self.nrecords = self.header_file['nrec'][0]
		199	self.sizeOfFile = os.path.getsize(self.filename)
		200	self.counter_records = 0
		201	self.flagIsNewFile = 0
		202	self.fileIndex += 1
		203
		204	return 1
		205
		206	def readNextBlock(self):
		207
		208	while True:
		209	if self.counter_records == self.nrecords:
		210	self.flagIsNewFile = 1
		211	if not self.setNextFile():
		212	return 0
		213
		214	self.readBlock()
		215
		216	if (self.datatime.time() < self.startTime) or (self.datatime.time() > self.endTime):
		217	print "[Reading] Record No. %d/%d -> %s [Skipping]" %(
		218	self.counter_records,
		219	self.nrecords,
		220	self.datatime.ctime())
		221	continue
		222	break
		223
		224	print "[Reading] Record No. %d/%d -> %s" %(
		225	self.counter_records,
		226	self.nrecords,
		227	self.datatime.ctime())
		228
		229	return 1
		230
		231	def readBlock(self):
		232
		233	pointer = self.fp.tell()
		234	header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
		235	self.nchannels = header_rec['nchan'][0]/2
		236	self.kchan = header_rec['nrxs'][0]
		237	self.nmodes = header_rec['nmodes'][0]
		238	self.nranges = header_rec['nranges'][0]
		239	self.fp.seek(pointer)
		240	self.height = numpy.empty((self.nmodes, self.nranges))
		241	self.snr = numpy.empty((self.nmodes, self.nchannels, self.nranges))
		242	self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
		243
		244	for mode in range(self.nmodes):
		245	self.readHeader()
		246	data = self.readData()
		247	self.height[mode] = (data[0] - self.correction) / 1000.
		248	self.buffer[mode] = data[1]
		249	self.snr[mode] = data[2]
		250
		251	self.counter_records = self.counter_records + self.nmodes
		252
		253	return
		254
		255	def readHeader(self):
		256	'''
		257	RecordHeader of BLTR rawdata file
		258	'''
		259
		260	header_structure = numpy.dtype(
		261	REC_HEADER_STRUCTURE.descr + [
		262	('antenna_coord', 'f4', (2, self.nchannels)),
		263	('rx_gains', 'u4', (self.nchannels,)),
		264	('rx_analysis', 'u4', (self.nchannels,))
		265	]
		266	)
		267
		268	self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
		269	self.lat = self.header_rec['lat'][0]
		270	self.lon = self.header_rec['lon'][0]
		271	self.delta = self.header_rec['delta_r'][0]
		272	self.correction = self.header_rec['dmode_rngcorr'][0]
		273	self.imode = self.header_rec['dmode_index'][0]
		274	self.antenna = self.header_rec['antenna_coord']
		275	self.rx_gains = self.header_rec['rx_gains']
		276	self.time = self.header_rec['time'][0]
		277	tseconds = self.header_rec['time'][0]
		278	local_t1 = time.localtime(tseconds)
		279	self.year = local_t1.tm_year
		280	self.month = local_t1.tm_mon
		281	self.day = local_t1.tm_mday
		282	self.t = datetime.datetime(self.year, self.month, self.day)
		283	self.datatime = datetime.datetime.utcfromtimestamp(self.time)
		284
		285	def readData(self):
		286	'''
		287	Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
		288
		289	Input:
		290	status_value - Array data is set to NAN for values that are not equal to status_value
		291
		292	'''
		293
		294	data_structure = numpy.dtype(
		295	DATA_STRUCTURE.descr + [
		296	('rx_saturation', 'u4', (self.nchannels,)),
		297	('chan_offset', 'u4', (2 * self.nchannels,)),
		298	('rx_amp', 'u4', (self.nchannels,)),
		299	('rx_snr', 'f4', (self.nchannels,)),
		300	('cross_snr', 'f4', (self.kchan,)),
		301	('sea_power_relative', 'f4', (self.kchan,))]
		302	)
		303
		304	data = numpy.fromfile(self.fp, data_structure, self.nranges)
		305
		306	height = data['range']
		307	winds = numpy.array((data['zonal'], data['meridional'], data['vertical']))
		308	snr = data['rx_snr'].T
		309
		310	winds[numpy.where(winds == -9999.)] = numpy.nan
		311	winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
		312	snr[numpy.where(snr == -9999.)] = numpy.nan
		313	snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
		314	snr = numpy.power(10, snr / 10)
		315
		316	return height, winds, snr
		317
		318	def set_output(self):
		319	'''
		320	Storing data from databuffer to dataOut object
		321	'''
		322
		323	self.dataOut.data_SNR = self.snr
		324	self.dataOut.height = self.height
		325	self.dataOut.data_output = self.buffer
		326	self.dataOut.utctimeInit = self.time
		327	self.dataOut.utctime = self.dataOut.utctimeInit
		328	self.dataOut.useLocalTime = False
		329	self.dataOut.paramInterval = 157
		330	self.dataOut.timezone = self.timezone
		331	self.dataOut.site = self.siteFile
		332	self.dataOut.nrecords = self.nrecords/self.nmodes
		333	self.dataOut.sizeOfFile = self.sizeOfFile
		334	self.dataOut.lat = self.lat
		335	self.dataOut.lon = self.lon
		336	self.dataOut.channelList = range(self.nchannels)
		337	self.dataOut.kchan = self.kchan
		338	# self.dataOut.nHeights = self.nranges
		339	self.dataOut.delta = self.delta
		340	self.dataOut.correction = self.correction
		341	self.dataOut.nmodes = self.nmodes
		342	self.dataOut.imode = self.imode
		343	self.dataOut.antenna = self.antenna
		344	self.dataOut.rx_gains = self.rx_gains
		345	self.dataOut.flagNoData = False
		346
		347	def getData(self):
		348	'''
		349	Storing data from databuffer to dataOut object
		350	'''
		351	if self.flagNoMoreFiles:
		352	self.dataOut.flagNoData = True
		353	print 'No file left to process'
		354	return 0
		355
		356	if not self.readNextBlock():
		357	self.dataOut.flagNoData = True
		358	return 0
		359
		360	self.set_output()
		361
		362	return 1

schainpy/model/io/jroIO_bltr.py created 644 +1154 0

This diff has been collapsed as it changes many lines, (1154 lines changed) Show them Hide them
		@@ -0,0 +1,1154
	1	import os, sys
	2	import glob
	3	import fnmatch
	4	import datetime
	5	import time
	6	import re
	7	import h5py
	8	import numpy
	9	import matplotlib.pyplot as plt
	10
	11	import pylab as plb
	12	from scipy.optimize import curve_fit
	13	from scipy import asarray as ar, exp
	14	from scipy import stats
	15
	16	from numpy.ma.core import getdata
	17
	18	SPEED_OF_LIGHT = 299792458
	19	SPEED_OF_LIGHT = 3e8
	20
	21	try:
	22	from gevent import sleep
	23	except:
	24	from time import sleep
	25
	26	from schainpy.model.data.jrodata import Spectra
	27	#from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader
	28	from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
	29	#from schainpy.model.io.jroIO_bltr import BLTRReader
	30	from numpy import imag, shape, NaN
	31
	32	from jroIO_base import JRODataReader
	33
	34
	35	class Header(object):
	36
	37	def __init__(self):
	38	raise NotImplementedError
	39
	40
	41	def read(self):
	42
	43	raise NotImplementedError
	44
	45	def write(self):
	46
	47	raise NotImplementedError
	48
	49	def printInfo(self):
	50
	51	message = "#"*50 + "\n"
	52	message += self.__class__.__name__.upper() + "\n"
	53	message += "#"*50 + "\n"
	54
	55	keyList = self.__dict__.keys()
	56	keyList.sort()
	57
	58	for key in keyList:
	59	message += "%s = %s" %(key, self.__dict__[key]) + "\n"
	60
	61	if "size" not in keyList:
	62	attr = getattr(self, "size")
	63
	64	if attr:
	65	message += "%s = %s" %("size", attr) + "\n"
	66
	67	#print message
	68
	69
	70
	71
	72
	73	FILE_STRUCTURE = numpy.dtype([ #HEADER 48bytes
	74	('FileMgcNumber','<u4'), #0x23020100
	75	('nFDTdataRecors','<u4'), #No Of FDT data records in this file (0 or more)
	76	('OffsetStartHeader','<u4'),
	77	('RadarUnitId','<u4'),
	78	('SiteName',numpy.str_,32), #Null terminated
	79	])
	80
	81	class FileHeaderBLTR(Header):
	82
	83	def __init__(self):
	84
	85	self.FileMgcNumber= 0 #0x23020100
	86	self.nFDTdataRecors=0 #No Of FDT data records in this file (0 or more)
	87	self.RadarUnitId= 0
	88	self.OffsetStartHeader=0
	89	self.SiteName= ""
	90	self.size = 48
	91
	92	def FHread(self, fp):
	93	#try:
	94	startFp = open(fp,"rb")
	95
	96	header = numpy.fromfile(startFp, FILE_STRUCTURE,1)
	97
	98	print ' '
	99	print 'puntero file header', startFp.tell()
	100	print ' '
	101
	102
	103	''' numpy.fromfile(file, dtype, count, sep='')
	104	file : file or str
	105	Open file object or filename.
	106
	107	dtype : data-type
	108	Data type of the returned array. For binary files, it is used to determine
	109	the size and byte-order of the items in the file.
	110
	111	count : int
	112	Number of items to read. -1 means all items (i.e., the complete file).
	113
	114	sep : str
	115	Separator between items if file is a text file. Empty ("") separator means
	116	the file should be treated as binary. Spaces (" ") in the separator match zero
	117	or more whitespace characters. A separator consisting only of spaces must match
	118	at least one whitespace.
	119
	120	'''
	121
	122
	123
	124	self.FileMgcNumber= hex(header['FileMgcNumber'][0])
	125	self.nFDTdataRecors=int(header['nFDTdataRecors'][0]) #No Of FDT data records in this file (0 or more)
	126	self.RadarUnitId= int(header['RadarUnitId'][0])
	127	self.OffsetStartHeader= int(header['OffsetStartHeader'][0])
	128	self.SiteName= str(header['SiteName'][0])
	129
	130	#print 'Numero de bloques', self.nFDTdataRecors
	131
	132
	133	if self.size <48:
	134	return 0
	135
	136	return 1
	137
	138
	139	def write(self, fp):
	140
	141	headerTuple = (self.FileMgcNumber,
	142	self.nFDTdataRecors,
	143	self.RadarUnitId,
	144	self.SiteName,
	145	self.size)
	146
	147
	148	header = numpy.array(headerTuple, FILE_STRUCTURE)
	149	# numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0)
	150	header.tofile(fp)
	151	''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
	152
	153	fid : file or str
	154	An open file object, or a string containing a filename.
	155
	156	sep : str
	157	Separator between array items for text output. If "" (empty), a binary file is written,
	158	equivalent to file.write(a.tobytes()).
	159
	160	format : str
	161	Format string for text file output. Each entry in the array is formatted to text by
	162	first converting it to the closest Python type, and then using "format" % item.
	163
	164	'''
	165
	166	return 1
	167
	168
	169
	170
	171
	172	RECORD_STRUCTURE = numpy.dtype([ #RECORD HEADER 180+20N bytes
	173	('RecMgcNumber','<u4'), #0x23030001
	174	('RecCounter','<u4'), #Record counter(0,1, ...)
	175	('Off2StartNxtRec','<u4'), #Offset to start of next record form start of this record
	176	('Off2StartData','<u4'), #Offset to start of data from start of this record
	177	('nUtime','<i4'), #Epoch time stamp of start of acquisition (seconds)
	178	('nMilisec','<u4'), #Millisecond component of time stamp (0,...,999)
	179	('ExpTagName',numpy.str_,32), #Experiment tag name (null terminated)
	180	('ExpComment',numpy.str_,32), #Experiment comment (null terminated)
	181	('SiteLatDegrees','<f4'), #Site latitude (from GPS) in degrees (positive implies North)
	182	('SiteLongDegrees','<f4'), #Site longitude (from GPS) in degrees (positive implies East)
	183	('RTCgpsStatus','<u4'), #RTC GPS engine status (0=SEEK, 1=LOCK, 2=NOT FITTED, 3=UNAVAILABLE)
	184	('TransmitFrec','<u4'), #Transmit frequency (Hz)
	185	('ReceiveFrec','<u4'), #Receive frequency
	186	('FirstOsciFrec','<u4'), #First local oscillator frequency (Hz)
	187	('Polarisation','<u4'), #(0="O", 1="E", 2="linear 1", 3="linear2")
	188	('ReceiverFiltSett','<u4'), #Receiver filter settings (0,1,2,3)
	189	('nModesInUse','<u4'), #Number of modes in use (1 or 2)
	190	('DualModeIndex','<u4'), #Dual Mode index number for these data (0 or 1)
	191	('DualModeRange','<u4'), #Dual Mode range correction for these data (m)
	192	('nDigChannels','<u4'), #Number of digital channels acquired (2*N)
	193	('SampResolution','<u4'), #Sampling resolution (meters)
	194	('nHeights','<u4'), #Number of range gates sampled
	195	('StartRangeSamp','<u4'), #Start range of sampling (meters)
	196	('PRFhz','<u4'), #PRF (Hz)
	197	('nCohInt','<u4'), #Integrations
	198	('nProfiles','<u4'), #Number of data points transformed
	199	('nChannels','<u4'), #Number of receive beams stored in file (1 or N)
	200	('nIncohInt','<u4'), #Number of spectral averages
	201	('FFTwindowingInd','<u4'), #FFT windowing index (0 = no window)
	202	('BeamAngleAzim','<f4'), #Beam steer angle (azimuth) in degrees (clockwise from true North)
	203	('BeamAngleZen','<f4'), #Beam steer angle (zenith) in degrees (0=> vertical)
	204	('AntennaCoord0','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
	205	('AntennaAngl0','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
	206	('AntennaCoord1','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
	207	('AntennaAngl1','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
	208	('AntennaCoord2','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
	209	('AntennaAngl2','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
	210	('RecPhaseCalibr0','<f4'), #Receiver phase calibration (degrees) - N values
	211	('RecPhaseCalibr1','<f4'), #Receiver phase calibration (degrees) - N values
	212	('RecPhaseCalibr2','<f4'), #Receiver phase calibration (degrees) - N values
	213	('RecAmpCalibr0','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
	214	('RecAmpCalibr1','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
	215	('RecAmpCalibr2','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
	216	('ReceiverGaindB0','<i4'), #Receiver gains in dB - N values
	217	('ReceiverGaindB1','<i4'), #Receiver gains in dB - N values
	218	('ReceiverGaindB2','<i4'), #Receiver gains in dB - N values
	219	])
	220
	221
	222	class RecordHeaderBLTR(Header):
	223
	224	def __init__(self, RecMgcNumber=None, RecCounter= 0, Off2StartNxtRec= 811248,
	225	nUtime= 0, nMilisec= 0, ExpTagName= None,
	226	ExpComment=None, SiteLatDegrees=0, SiteLongDegrees= 0,
	227	RTCgpsStatus= 0, TransmitFrec= 0, ReceiveFrec= 0,
	228	FirstOsciFrec= 0, Polarisation= 0, ReceiverFiltSett= 0,
	229	nModesInUse= 0, DualModeIndex= 0, DualModeRange= 0,
	230	nDigChannels= 0, SampResolution= 0, nHeights= 0,
	231	StartRangeSamp= 0, PRFhz= 0, nCohInt= 0,
	232	nProfiles= 0, nChannels= 0, nIncohInt= 0,
	233	FFTwindowingInd= 0, BeamAngleAzim= 0, BeamAngleZen= 0,
	234	AntennaCoord0= 0, AntennaCoord1= 0, AntennaCoord2= 0,
	235	RecPhaseCalibr0= 0, RecPhaseCalibr1= 0, RecPhaseCalibr2= 0,
	236	RecAmpCalibr0= 0, RecAmpCalibr1= 0, RecAmpCalibr2= 0,
	237	AntennaAngl0=0, AntennaAngl1=0, AntennaAngl2=0,
	238	ReceiverGaindB0= 0, ReceiverGaindB1= 0, ReceiverGaindB2= 0, Off2StartData=0, OffsetStartHeader=0):
	239
	240	self.RecMgcNumber = RecMgcNumber #0x23030001
	241	self.RecCounter = RecCounter
	242	self.Off2StartNxtRec = Off2StartNxtRec
	243	self.Off2StartData = Off2StartData
	244	self.nUtime = nUtime
	245	self.nMilisec = nMilisec
	246	self.ExpTagName = ExpTagName
	247	self.ExpComment = ExpComment
	248	self.SiteLatDegrees = SiteLatDegrees
	249	self.SiteLongDegrees = SiteLongDegrees
	250	self.RTCgpsStatus = RTCgpsStatus
	251	self.TransmitFrec = TransmitFrec
	252	self.ReceiveFrec = ReceiveFrec
	253	self.FirstOsciFrec = FirstOsciFrec
	254	self.Polarisation = Polarisation
	255	self.ReceiverFiltSett = ReceiverFiltSett
	256	self.nModesInUse = nModesInUse
	257	self.DualModeIndex = DualModeIndex
	258	self.DualModeRange = DualModeRange
	259	self.nDigChannels = nDigChannels
	260	self.SampResolution = SampResolution
	261	self.nHeights = nHeights
	262	self.StartRangeSamp = StartRangeSamp
	263	self.PRFhz = PRFhz
	264	self.nCohInt = nCohInt
	265	self.nProfiles = nProfiles
	266	self.nChannels = nChannels
	267	self.nIncohInt = nIncohInt
	268	self.FFTwindowingInd = FFTwindowingInd
	269	self.BeamAngleAzim = BeamAngleAzim
	270	self.BeamAngleZen = BeamAngleZen
	271	self.AntennaCoord0 = AntennaCoord0
	272	self.AntennaAngl0 = AntennaAngl0
	273	self.AntennaAngl1 = AntennaAngl1
	274	self.AntennaAngl2 = AntennaAngl2
	275	self.AntennaCoord1 = AntennaCoord1
	276	self.AntennaCoord2 = AntennaCoord2
	277	self.RecPhaseCalibr0 = RecPhaseCalibr0
	278	self.RecPhaseCalibr1 = RecPhaseCalibr1
	279	self.RecPhaseCalibr2 = RecPhaseCalibr2
	280	self.RecAmpCalibr0 = RecAmpCalibr0
	281	self.RecAmpCalibr1 = RecAmpCalibr1
	282	self.RecAmpCalibr2 = RecAmpCalibr2
	283	self.ReceiverGaindB0 = ReceiverGaindB0
	284	self.ReceiverGaindB1 = ReceiverGaindB1
	285	self.ReceiverGaindB2 = ReceiverGaindB2
	286	self.OffsetStartHeader = 48
	287
	288
	289
	290	def RHread(self, fp):
	291	#print fp
	292	#startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
	293	startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
	294	#RecCounter=0
	295	#Off2StartNxtRec=811248
	296	OffRHeader= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
	297	print ' '
	298	print 'puntero Record Header', startFp.tell()
	299	print ' '
	300
	301
	302	startFp.seek(OffRHeader, os.SEEK_SET)
	303
	304	print ' '
	305	print 'puntero Record Header con seek', startFp.tell()
	306	print ' '
	307
	308	#print 'Posicion del bloque: ',OffRHeader
	309
	310	header = numpy.fromfile(startFp,RECORD_STRUCTURE,1)
	311
	312	print ' '
	313	print 'puntero Record Header con seek', startFp.tell()
	314	print ' '
	315
	316	print ' '
	317	#
	318	#print 'puntero Record Header despues de seek', header.tell()
	319	print ' '
	320
	321	self.RecMgcNumber = hex(header['RecMgcNumber'][0]) #0x23030001
	322	self.RecCounter = int(header['RecCounter'][0])
	323	self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
	324	self.Off2StartData = int(header['Off2StartData'][0])
	325	self.nUtime = header['nUtime'][0]
	326	self.nMilisec = header['nMilisec'][0]
	327	self.ExpTagName = str(header['ExpTagName'][0])
	328	self.ExpComment = str(header['ExpComment'][0])
	329	self.SiteLatDegrees = header['SiteLatDegrees'][0]
	330	self.SiteLongDegrees = header['SiteLongDegrees'][0]
	331	self.RTCgpsStatus = header['RTCgpsStatus'][0]
	332	self.TransmitFrec = header['TransmitFrec'][0]
	333	self.ReceiveFrec = header['ReceiveFrec'][0]
	334	self.FirstOsciFrec = header['FirstOsciFrec'][0]
	335	self.Polarisation = header['Polarisation'][0]
	336	self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
	337	self.nModesInUse = header['nModesInUse'][0]
	338	self.DualModeIndex = header['DualModeIndex'][0]
	339	self.DualModeRange = header['DualModeRange'][0]
	340	self.nDigChannels = header['nDigChannels'][0]
	341	self.SampResolution = header['SampResolution'][0]
	342	self.nHeights = header['nHeights'][0]
	343	self.StartRangeSamp = header['StartRangeSamp'][0]
	344	self.PRFhz = header['PRFhz'][0]
	345	self.nCohInt = header['nCohInt'][0]
	346	self.nProfiles = header['nProfiles'][0]
	347	self.nChannels = header['nChannels'][0]
	348	self.nIncohInt = header['nIncohInt'][0]
	349	self.FFTwindowingInd = header['FFTwindowingInd'][0]
	350	self.BeamAngleAzim = header['BeamAngleAzim'][0]
	351	self.BeamAngleZen = header['BeamAngleZen'][0]
	352	self.AntennaCoord0 = header['AntennaCoord0'][0]
	353	self.AntennaAngl0 = header['AntennaAngl0'][0]
	354	self.AntennaCoord1 = header['AntennaCoord1'][0]
	355	self.AntennaAngl1 = header['AntennaAngl1'][0]
	356	self.AntennaCoord2 = header['AntennaCoord2'][0]
	357	self.AntennaAngl2 = header['AntennaAngl2'][0]
	358	self.RecPhaseCalibr0 = header['RecPhaseCalibr0'][0]
	359	self.RecPhaseCalibr1 = header['RecPhaseCalibr1'][0]
	360	self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0]
	361	self.RecAmpCalibr0 = header['RecAmpCalibr0'][0]
	362	self.RecAmpCalibr1 = header['RecAmpCalibr1'][0]
	363	self.RecAmpCalibr2 = header['RecAmpCalibr2'][0]
	364	self.ReceiverGaindB0 = header['ReceiverGaindB0'][0]
	365	self.ReceiverGaindB1 = header['ReceiverGaindB1'][0]
	366	self.ReceiverGaindB2 = header['ReceiverGaindB2'][0]
	367
	368	self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
	369
	370	self.RHsize = 180+20*self.nChannels
	371	self.Datasize= self.nProfilesself.nChannelsself.nHeights24
	372	#print 'Datasize',self.Datasize
	373	endFp = self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
	374
	375	print '=============================================='
	376	print 'RecMgcNumber ',self.RecMgcNumber
	377	print 'RecCounter ',self.RecCounter
	378	print 'Off2StartNxtRec ',self.Off2StartNxtRec
	379	print 'Off2StartData ',self.Off2StartData
	380	print 'Range Resolution ',self.SampResolution
	381	print 'First Height ',self.StartRangeSamp
	382	print 'PRF (Hz) ',self.PRFhz
	383	print 'Heights (K) ',self.nHeights
	384	print 'Channels (N) ',self.nChannels
	385	print 'Profiles (J) ',self.nProfiles
	386	print 'iCoh ',self.nCohInt
	387	print 'iInCoh ',self.nIncohInt
	388	print 'BeamAngleAzim ',self.BeamAngleAzim
	389	print 'BeamAngleZen ',self.BeamAngleZen
	390
	391	#print 'ModoEnUso ',self.DualModeIndex
	392	#print 'UtcTime ',self.nUtime
	393	#print 'MiliSec ',self.nMilisec
	394	#print 'Exp TagName ',self.ExpTagName
	395	#print 'Exp Comment ',self.ExpComment
	396	#print 'FFT Window Index ',self.FFTwindowingInd
	397	#print 'N Dig. Channels ',self.nDigChannels
	398	print 'Size de bloque ',self.RHsize
	399	print 'DataSize ',self.Datasize
	400	print 'BeamAngleAzim ',self.BeamAngleAzim
	401	#print 'AntennaCoord0 ',self.AntennaCoord0
	402	#print 'AntennaAngl0 ',self.AntennaAngl0
	403	#print 'AntennaCoord1 ',self.AntennaCoord1
	404	#print 'AntennaAngl1 ',self.AntennaAngl1
	405	#print 'AntennaCoord2 ',self.AntennaCoord2
	406	#print 'AntennaAngl2 ',self.AntennaAngl2
	407	print 'RecPhaseCalibr0 ',self.RecPhaseCalibr0
	408	print 'RecPhaseCalibr1 ',self.RecPhaseCalibr1
	409	print 'RecPhaseCalibr2 ',self.RecPhaseCalibr2
	410	print 'RecAmpCalibr0 ',self.RecAmpCalibr0
	411	print 'RecAmpCalibr1 ',self.RecAmpCalibr1
	412	print 'RecAmpCalibr2 ',self.RecAmpCalibr2
	413	print 'ReceiverGaindB0 ',self.ReceiverGaindB0
	414	print 'ReceiverGaindB1 ',self.ReceiverGaindB1
	415	print 'ReceiverGaindB2 ',self.ReceiverGaindB2
	416	print '=============================================='
	417
	418	if OffRHeader > endFp:
	419	sys.stderr.write("Warning %s: Size value read from System Header is lower than it has to be\n" %fp)
	420	return 0
	421
	422	if OffRHeader < endFp:
	423	sys.stderr.write("Warning %s: Size value read from System Header size is greater than it has to be\n" %fp)
	424	return 0
	425
	426	return 1
	427
	428
	429	class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODataReader):
	430
	431	path = None
	432	startDate = None
	433	endDate = None
	434	startTime = None
	435	endTime = None
	436	walk = None
	437	isConfig = False
	438
	439
	440	fileList= None
	441
	442	#metadata
	443	TimeZone= None
	444	Interval= None
	445	heightList= None
	446
	447	#data
	448	data= None
	449	utctime= None
	450
	451
	452
	453	def __init__(self, **kwargs):
	454
	455	#Eliminar de la base la herencia
	456	ProcessingUnit.__init__(self, **kwargs)
	457
	458	#self.isConfig = False
	459
	460	#self.pts2read_SelfSpectra = 0
	461	#self.pts2read_CrossSpectra = 0
	462	#self.pts2read_DCchannels = 0
	463	#self.datablock = None
	464	self.utc = None
	465	self.ext = ".fdt"
	466	self.optchar = "P"
	467	self.fpFile=None
	468	self.fp = None
	469	self.BlockCounter=0
	470	self.dtype = None
	471	self.fileSizeByHeader = None
	472	self.filenameList = []
	473	self.fileSelector = 0
	474	self.Off2StartNxtRec=0
	475	self.RecCounter=0
	476	self.flagNoMoreFiles = 0
	477	self.data_spc=None
	478	self.data_cspc=None
	479	self.data_output=None
	480	self.path = None
	481	self.OffsetStartHeader=0
	482	self.Off2StartData=0
	483	self.ipp = 0
	484	self.nFDTdataRecors=0
	485	self.blocksize = 0
	486	self.dataOut = Spectra()
	487	self.profileIndex = 1 #Always
	488	self.dataOut.flagNoData=False
	489	self.dataOut.nRdPairs = 0
	490	self.dataOut.pairsList = []
	491	self.dataOut.data_spc=None
	492	self.dataOut.noise=[]
	493	self.dataOut.velocityX=[]
	494	self.dataOut.velocityY=[]
	495	self.dataOut.velocityV=[]
	496
	497
	498
	499	def Files2Read(self, fp):
	500	'''
	501	Function that indicates the number of .fdt files that exist in the folder to be read.
	502	It also creates an organized list with the names of the files to read.
	503	'''
	504	#self.__checkPath()
	505
	506	ListaData=os.listdir(fp) #Gets the list of files within the fp address
	507	ListaData=sorted(ListaData) #Sort the list of files from least to largest by names
	508	nFiles=0 #File Counter
	509	FileList=[] #A list is created that will contain the .fdt files
	510	for IndexFile in ListaData :
	511	if '.fdt' in IndexFile:
	512	FileList.append(IndexFile)
	513	nFiles+=1
	514
	515	#print 'Files2Read'
	516	#print 'Existen '+str(nFiles)+' archivos .fdt'
	517
	518	self.filenameList=FileList #List of files from least to largest by names
	519
	520
	521	def run(self, **kwargs):
	522	'''
	523	This method will be the one that will initiate the data entry, will be called constantly.
	524	You should first verify that your Setup () is set up and then continue to acquire
	525	the data to be processed with getData ().
	526	'''
	527	if not self.isConfig:
	528	self.setup(**kwargs)
	529	self.isConfig = True
	530
	531	self.getData()
	532	#print 'running'
	533
	534
	535	def setup(self, path=None,
	536	startDate=None,
	537	endDate=None,
	538	startTime=None,
	539	endTime=None,
	540	walk=True,
	541	timezone='utc',
	542	code = None,
	543	online=False,
	544	ReadMode=None,
	545	**kwargs):
	546
	547	self.isConfig = True
	548
	549	self.path=path
	550	self.startDate=startDate
	551	self.endDate=endDate
	552	self.startTime=startTime
	553	self.endTime=endTime
	554	self.walk=walk
	555	self.ReadMode=int(ReadMode)
	556
	557	pass
	558
	559
	560	def getData(self):
	561	'''
	562	Before starting this function, you should check that there is still an unread file,
	563	If there are still blocks to read or if the data block is empty.
	564
	565	You should call the file "read".
	566
	567	'''
	568
	569	if self.flagNoMoreFiles:
	570	self.dataOut.flagNoData = True
	571	print 'NoData se vuelve true'
	572	return 0
	573
	574	self.fp=self.path
	575	self.Files2Read(self.fp)
	576	self.readFile(self.fp)
	577	self.dataOut.data_spc = self.data_spc
	578	self.dataOut.data_cspc =self.data_cspc
	579	self.dataOut.data_output=self.data_output
	580
	581	print 'self.dataOut.data_output', shape(self.dataOut.data_output)
	582
	583	#self.removeDC()
	584	return self.dataOut.data_spc
	585
	586
	587	def readFile(self,fp):
	588	'''
	589	You must indicate if you are reading in Online or Offline mode and load the
	590	The parameters for this file reading mode.
	591
	592	Then you must do 2 actions:
	593
	594	1. Get the BLTR FileHeader.
	595	2. Start reading the first block.
	596	'''
	597
	598	#The address of the folder is generated the name of the .fdt file that will be read
	599	print "File: ",self.fileSelector+1
	600
	601	if self.fileSelector < len(self.filenameList):
	602
	603	self.fpFile=str(fp)+'/'+str(self.filenameList[self.fileSelector])
	604	#print self.fpFile
	605	fheader = FileHeaderBLTR()
	606	fheader.FHread(self.fpFile) #Bltr FileHeader Reading
	607	self.nFDTdataRecors=fheader.nFDTdataRecors
	608
	609	self.readBlock() #Block reading
	610	else:
	611	print 'readFile FlagNoData becomes true'
	612	self.flagNoMoreFiles=True
	613	self.dataOut.flagNoData = True
	614	return 0
	615
	616	def getVelRange(self, extrapoints=0):
	617	Lambda= SPEED_OF_LIGHT/50000000
	618	PRF = self.dataOut.PRF#1./(self.dataOut.ippSeconds * self.dataOut.nCohInt)
	619	Vmax=-Lambda/(4.(1./PRF)self.dataOut.nCohInt*2.)
	620	deltafreq = PRF / (self.nProfiles)
	621	deltavel = (Vmax*2) / (self.nProfiles)
	622	freqrange = deltafreq*(numpy.arange(self.nProfiles)-self.nProfiles/2.) - deltafreq/2
	623	velrange = deltavel*(numpy.arange(self.nProfiles)-self.nProfiles/2.)
	624	return velrange
	625
	626	def readBlock(self):
	627	'''
	628	It should be checked if the block has data, if it is not passed to the next file.
	629
	630	Then the following is done:
	631
	632	1. Read the RecordHeader
	633	2. Fill the buffer with the current block number.
	634
	635	'''
	636
	637	if self.BlockCounter < self.nFDTdataRecors-2:
	638	print self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
	639	if self.ReadMode==1:
	640	rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1)
	641	elif self.ReadMode==0:
	642	rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter)
	643
	644	rheader.RHread(self.fpFile) #Bltr FileHeader Reading
	645
	646	self.OffsetStartHeader=rheader.OffsetStartHeader
	647	self.RecCounter=rheader.RecCounter
	648	self.Off2StartNxtRec=rheader.Off2StartNxtRec
	649	self.Off2StartData=rheader.Off2StartData
	650	self.nProfiles=rheader.nProfiles
	651	self.nChannels=rheader.nChannels
	652	self.nHeights=rheader.nHeights
	653	self.frequency=rheader.TransmitFrec
	654	self.DualModeIndex=rheader.DualModeIndex
	655
	656	self.pairsList =[(0,1),(0,2),(1,2)]
	657	self.dataOut.pairsList = self.pairsList
	658
	659	self.nRdPairs=len(self.dataOut.pairsList)
	660	self.dataOut.nRdPairs = self.nRdPairs
	661
	662	self.__firstHeigth=rheader.StartRangeSamp
	663	self.__deltaHeigth=rheader.SampResolution
	664	self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth
	665	self.dataOut.channelList = range(self.nChannels)
	666	self.dataOut.nProfiles=rheader.nProfiles
	667	self.dataOut.nIncohInt=rheader.nIncohInt
	668	self.dataOut.nCohInt=rheader.nCohInt
	669	self.dataOut.ippSeconds= 1/float(rheader.PRFhz)
	670	self.dataOut.PRF=rheader.PRFhz
	671	self.dataOut.nFFTPoints=rheader.nProfiles
	672	self.dataOut.utctime=rheader.nUtime
	673	self.dataOut.timeZone=0
	674	self.dataOut.normFactor= self.dataOut.nProfilesself.dataOut.nIncohIntself.dataOut.nCohInt
	675	self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
	676
	677	self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN
	678	print 'self.data_output', shape(self.data_output)
	679	self.dataOut.velocityX=[]
	680	self.dataOut.velocityY=[]
	681	self.dataOut.velocityV=[]
	682
	683	'''Block Reading, the Block Data is received and Reshape is used to give it
	684	shape.
	685	'''
	686
	687	#Procedure to take the pointer to where the date block starts
	688	startDATA = open(self.fpFile,"rb")
	689	OffDATA= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec+self.Off2StartData
	690	startDATA.seek(OffDATA, os.SEEK_SET)
	691
	692	def moving_average(x, N=2):
	693	return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
	694
	695	def gaus(xSamples,a,x0,sigma):
	696	return aexp(-(xSamples-x0)2/(2sigma**2))
	697
	698	def Find(x,value):
	699	for index in range(len(x)):
	700	if x[index]==value:
	701	return index
	702
	703	def pol2cart(rho, phi):
	704	x = rho * numpy.cos(phi)
	705	y = rho * numpy.sin(phi)
	706	return(x, y)
	707
	708
	709
	710
	711	if self.DualModeIndex==self.ReadMode:
	712
	713	self.data_fft = numpy.fromfile( startDATA, [('complex','<c8')],self.nProfilesself.nChannelsself.nHeights )
	714
	715	self.data_fft=self.data_fft.astype(numpy.dtype('complex'))
	716
	717	self.data_block=numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles ))
	718
	719	self.data_block = numpy.transpose(self.data_block, (1,2,0))
	720
	721	copy = self.data_block.copy()
	722	spc = copy * numpy.conjugate(copy)
	723
	724	self.data_spc = numpy.absolute(spc) # valor absoluto o magnitud
	725
	726	factor = self.dataOut.normFactor
	727
	728
	729	z = self.data_spc.copy()#/factor
	730	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
	731	#zdB = 10*numpy.log10(z)
	732	print ' '
	733	print 'Z: '
	734	print shape(z)
	735	print ' '
	736	print ' '
	737
	738	self.dataOut.data_spc=self.data_spc
	739
	740	self.noise = self.dataOut.getNoise(ymin_index=80, ymax_index=132)#/factor
	741	#noisedB = 10*numpy.log10(self.noise)
	742
	743
	744	ySamples=numpy.ones([3,self.nProfiles])
	745	phase=numpy.ones([3,self.nProfiles])
	746	CSPCSamples=numpy.ones([3,self.nProfiles],dtype=numpy.complex_)
	747	coherence=numpy.ones([3,self.nProfiles])
	748	PhaseSlope=numpy.ones(3)
	749	PhaseInter=numpy.ones(3)
	750
	751	'''**** Getting CrossSpectra ****'''
	752	cspc=self.data_block.copy()
	753	self.data_cspc=self.data_block.copy()
	754
	755	xFrec=self.getVelRange(1)
	756	VelRange=self.getVelRange(1)
	757	self.dataOut.VelRange=VelRange
	758	#print ' '
	759	#print ' '
	760	#print 'xFrec',xFrec
	761	#print ' '
	762	#print ' '
	763	#Height=35
	764	for i in range(self.nRdPairs):
	765
	766	chan_index0 = self.dataOut.pairsList[i][0]
	767	chan_index1 = self.dataOut.pairsList[i][1]
	768
	769	self.data_cspc[i,:,:]=cspc[chan_index0,:,:] * numpy.conjugate(cspc[chan_index1,:,:])
	770
	771
	772	'''Getting Eij and Nij'''
	773	(AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180)
	774	(AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180)
	775	(AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180)
	776
	777	E01=AntennaX0-AntennaX1
	778	N01=AntennaY0-AntennaY1
	779
	780	E02=AntennaX0-AntennaX2
	781	N02=AntennaY0-AntennaY2
	782
	783	E12=AntennaX1-AntennaX2
	784	N12=AntennaY1-AntennaY2
	785
	786	self.ChanDist= numpy.array([[E01, N01],[E02,N02],[E12,N12]])
	787
	788	self.dataOut.ChanDist = self.ChanDist
	789
	790
	791	# for Height in range(self.nHeights):
	792	#
	793	# for i in range(self.nRdPairs):
	794	#
	795	# '''**** Line of Data SPC ****'''
	796	# zline=z[i,:,Height]
	797	#
	798	# '''**** DC is removed ****'''
	799	# DC=Find(zline,numpy.amax(zline))
	800	# zline[DC]=(zline[DC-1]+zline[DC+1])/2
	801	#
	802	#
	803	# '''**** SPC is normalized ****'''
	804	# FactNorm= zline.copy() / numpy.sum(zline.copy())
	805	# FactNorm= FactNorm/numpy.sum(FactNorm)
	806	#
	807	# SmoothSPC=moving_average(FactNorm,N=3)
	808	#
	809	# xSamples = ar(range(len(SmoothSPC)))
	810	# ySamples[i] = SmoothSPC-self.noise[i]
	811	#
	812	# for i in range(self.nRdPairs):
	813	#
	814	# '''**** Line of Data CSPC ****'''
	815	# cspcLine=self.data_cspc[i,:,Height].copy()
	816	#
	817	#
	818	#
	819	# '''**** CSPC is normalized ****'''
	820	# chan_index0 = self.dataOut.pairsList[i][0]
	821	# chan_index1 = self.dataOut.pairsList[i][1]
	822	# CSPCFactor= numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])
	823	#
	824	#
	825	# CSPCNorm= cspcLine.copy() / numpy.sqrt(CSPCFactor)
	826	#
	827	#
	828	# CSPCSamples[i] = CSPCNorm-self.noise[i]
	829	# coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
	830	#
	831	# '''**** DC is removed ****'''
	832	# DC=Find(coherence[i],numpy.amax(coherence[i]))
	833	# coherence[i][DC]=(coherence[i][DC-1]+coherence[i][DC+1])/2
	834	# coherence[i]= moving_average(coherence[i],N=2)
	835	#
	836	# phase[i] = moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
	837	#
	838	#
	839	# '''**** Getting fij width ****'''
	840	#
	841	# yMean=[]
	842	# yMean2=[]
	843	#
	844	# for j in range(len(ySamples[1])):
	845	# yMean=numpy.append(yMean,numpy.average([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
	846	#
	847	# '''***** Getting fitting Gaussian ****'''
	848	# meanGauss=sum(xSamples*yMean) / len(xSamples)
	849	# sigma=sum(yMean(xSamples-meanGauss)*2) / len(xSamples)
	850	# #print 'Height',Height,'SNR', meanGauss/sigma**2
	851	#
	852	# if (abs(meanGauss/sigma**2) > 0.0001) :
	853	#
	854	# try:
	855	# popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
	856	#
	857	# if numpy.amax(popt)>numpy.amax(yMean)*0.3:
	858	# FitGauss=gaus(xSamples,*popt)
	859	#
	860	# else:
	861	# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
	862	# print 'Verificador: Dentro', Height
	863	# except RuntimeError:
	864	#
	865	# try:
	866	# for j in range(len(ySamples[1])):
	867	# yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]]))
	868	# popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma])
	869	# FitGauss=gaus(xSamples,*popt)
	870	# print 'Verificador: Exepcion1', Height
	871	# except RuntimeError:
	872	#
	873	# try:
	874	# popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma])
	875	# FitGauss=gaus(xSamples,*popt)
	876	# print 'Verificador: Exepcion2', Height
	877	# except RuntimeError:
	878	# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
	879	# print 'Verificador: Exepcion3', Height
	880	# else:
	881	# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
	882	# #print 'Verificador: Fuera', Height
	883	#
	884	#
	885	#
	886	# Maximun=numpy.amax(yMean)
	887	# eMinus1=Maximun*numpy.exp(-1)
	888	#
	889	# HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
	890	# HalfWidth= xFrec[HWpos]
	891	# GCpos=Find(FitGauss, numpy.amax(FitGauss))
	892	# Vpos=Find(FactNorm, numpy.amax(FactNorm))
	893	# #Vpos=numpy.sum(FactNorm)/len(FactNorm)
	894	# #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) )))
	895	# #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos
	896	# '''**** Getting Fij ****'''
	897	#
	898	# GaussCenter=xFrec[GCpos]
	899	# if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
	900	# Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
	901	# else:
	902	# Fij=abs(GaussCenter-HalfWidth)+0.0000001
	903	#
	904	# '''**** Getting Frecuency range of significant data ****'''
	905	#
	906	# Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
	907	#
	908	# if Rangpos<GCpos:
	909	# Range=numpy.array([Rangpos,2*GCpos-Rangpos])
	910	# else:
	911	# Range=numpy.array([2*GCpos-Rangpos,Rangpos])
	912	#
	913	# FrecRange=xFrec[Range[0]:Range[1]]
	914	#
	915	# #print 'FrecRange', FrecRange
	916	# '''**** Getting SCPC Slope ****'''
	917	#
	918	# for i in range(self.nRdPairs):
	919	#
	920	# if len(FrecRange)>5 and len(FrecRange)<self.nProfiles*0.5:
	921	# PhaseRange=moving_average(phase[i,Range[0]:Range[1]],N=3)
	922	#
	923	# slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
	924	# PhaseSlope[i]=slope
	925	# PhaseInter[i]=intercept
	926	# else:
	927	# PhaseSlope[i]=0
	928	# PhaseInter[i]=0
	929	#
	930	# # plt.figure(i+15)
	931	# # plt.title('FASE ( CH%s*CH%s )' %(self.dataOut.pairsList[i][0],self.dataOut.pairsList[i][1]))
	932	# # plt.xlabel('Frecuencia (KHz)')
	933	# # plt.ylabel('Magnitud')
	934	# # #plt.subplot(311+i)
	935	# # plt.plot(FrecRange,PhaseRange,'b')
	936	# # plt.plot(FrecRange,FrecRange*PhaseSlope[i]+PhaseInter[i],'r')
	937	#
	938	# #plt.axis([-0.6, 0.2, -3.2, 3.2])
	939	#
	940	#
	941	# '''Getting constant C'''
	942	# cC=(Fijnumpy.pi)*2
	943	#
	944	# # '''Getting Eij and Nij'''
	945	# # (AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180)
	946	# # (AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180)
	947	# # (AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180)
	948	# #
	949	# # E01=AntennaX0-AntennaX1
	950	# # N01=AntennaY0-AntennaY1
	951	# #
	952	# # E02=AntennaX0-AntennaX2
	953	# # N02=AntennaY0-AntennaY2
	954	# #
	955	# # E12=AntennaX1-AntennaX2
	956	# # N12=AntennaY1-AntennaY2
	957	#
	958	# '''**** Getting constants F and G ****'''
	959	# MijEijNij=numpy.array([[E02,N02], [E12,N12]])
	960	# MijResult0=(-PhaseSlope[1]cC) / (2numpy.pi)
	961	# MijResult1=(-PhaseSlope[2]cC) / (2numpy.pi)
	962	# MijResults=numpy.array([MijResult0,MijResult1])
	963	# (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
	964	#
	965	# '''**** Getting constants A, B and H ****'''
	966	# W01=numpy.amax(coherence[0])
	967	# W02=numpy.amax(coherence[1])
	968	# W12=numpy.amax(coherence[2])
	969	#
	970	# WijResult0=((cFE01+cGN01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
	971	# WijResult1=((cFE02+cGN02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
	972	# WijResult2=((cFE12+cGN12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
	973	#
	974	# WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
	975	#
	976	# WijEijNij=numpy.array([ [E012, N012, 2E01N01] , [E022, N022, 2E02N02] , [E122, N122, 2E12N12] ])
	977	# (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
	978	#
	979	# VxVy=numpy.array([[cA,cH],[cH,cB]])
	980	#
	981	# VxVyResults=numpy.array([-cF,-cG])
	982	# (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
	983	# Vzon = Vy
	984	# Vmer = Vx
	985	# Vmag=numpy.sqrt(Vzon2+Vmer2)
	986	# Vang=numpy.arctan2(Vmer,Vzon)
	987	#
	988	# if abs(Vy)<100 and abs(Vy)> 0.:
	989	# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag
	990	# #print 'Vmag',Vmag
	991	# else:
	992	# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN)
	993	#
	994	# if abs(Vx)<100 and abs(Vx) > 0.:
	995	# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang
	996	# #print 'Vang',Vang
	997	# else:
	998	# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN)
	999	#
	1000	# if abs(GaussCenter)<2:
	1001	# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos])
	1002	#
	1003	# else:
	1004	# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN)
	1005	#
	1006	#
	1007	# # print '********************************************'
	1008	# # print 'HalfWidth ', HalfWidth
	1009	# # print 'Maximun ', Maximun
	1010	# # print 'eMinus1 ', eMinus1
	1011	# # print 'Rangpos ', Rangpos
	1012	# # print 'GaussCenter ',GaussCenter
	1013	# # print 'E01 ',E01
	1014	# # print 'N01 ',N01
	1015	# # print 'E02 ',E02
	1016	# # print 'N02 ',N02
	1017	# # print 'E12 ',E12
	1018	# # print 'N12 ',N12
	1019	# #print 'self.dataOut.velocityX ', self.dataOut.velocityX
	1020	# # print 'Fij ', Fij
	1021	# # print 'cC ', cC
	1022	# # print 'cF ', cF
	1023	# # print 'cG ', cG
	1024	# # print 'cA ', cA
	1025	# # print 'cB ', cB
	1026	# # print 'cH ', cH
	1027	# # print 'Vx ', Vx
	1028	# # print 'Vy ', Vy
	1029	# # print 'Vmag ', Vmag
	1030	# # print 'Vang ', Vang*180/numpy.pi
	1031	# # print 'PhaseSlope ',PhaseSlope[0]
	1032	# # print 'PhaseSlope ',PhaseSlope[1]
	1033	# # print 'PhaseSlope ',PhaseSlope[2]
	1034	# # print '********************************************'
	1035	# #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY)
	1036	#
	1037	# #print 'self.dataOut.velocityX', len(self.dataOut.velocityX)
	1038	# #print 'self.dataOut.velocityY', len(self.dataOut.velocityY)
	1039	# #print 'self.dataOut.velocityV', self.dataOut.velocityV
	1040	#
	1041	# self.data_output[0]=numpy.array(self.dataOut.velocityX)
	1042	# self.data_output[1]=numpy.array(self.dataOut.velocityY)
	1043	# self.data_output[2]=numpy.array(self.dataOut.velocityV)
	1044	#
	1045	# prin= self.data_output[0][~numpy.isnan(self.data_output[0])]
	1046	# print ' '
	1047	# print 'VmagAverage',numpy.mean(prin)
	1048	# print ' '
	1049	# # plt.figure(5)
	1050	# # plt.subplot(211)
	1051	# # plt.plot(self.dataOut.velocityX,'yo:')
	1052	# # plt.subplot(212)
	1053	# # plt.plot(self.dataOut.velocityY,'yo:')
	1054	#
	1055	# # plt.figure(1)
	1056	# # # plt.subplot(121)
	1057	# # # plt.plot(xFrec,ySamples[0],'k',label='Ch0')
	1058	# # # plt.plot(xFrec,ySamples[1],'g',label='Ch1')
	1059	# # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
	1060	# # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
	1061	# # # plt.legend()
	1062	# # plt.title('DATOS A ALTURA DE 2850 METROS')
	1063	# #
	1064	# # plt.xlabel('Frecuencia (KHz)')
	1065	# # plt.ylabel('Magnitud')
	1066	# # # plt.subplot(122)
	1067	# # # plt.title('Fit for Time Constant')
	1068	# # #plt.plot(xFrec,zline)
	1069	# # #plt.plot(xFrec,SmoothSPC,'g')
	1070	# # plt.plot(xFrec,FactNorm)
	1071	# # plt.axis([-4, 4, 0, 0.15])
	1072	# # # plt.xlabel('SelfSpectra KHz')
	1073	# #
	1074	# # plt.figure(10)
	1075	# # # plt.subplot(121)
	1076	# # plt.plot(xFrec,ySamples[0],'b',label='Ch0')
	1077	# # plt.plot(xFrec,ySamples[1],'y',label='Ch1')
	1078	# # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
	1079	# # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
	1080	# # plt.legend()
	1081	# # plt.title('SELFSPECTRA EN CANALES')
	1082	# #
	1083	# # plt.xlabel('Frecuencia (KHz)')
	1084	# # plt.ylabel('Magnitud')
	1085	# # # plt.subplot(122)
	1086	# # # plt.title('Fit for Time Constant')
	1087	# # #plt.plot(xFrec,zline)
	1088	# # #plt.plot(xFrec,SmoothSPC,'g')
	1089	# # # plt.plot(xFrec,FactNorm)
	1090	# # # plt.axis([-4, 4, 0, 0.15])
	1091	# # # plt.xlabel('SelfSpectra KHz')
	1092	# #
	1093	# # plt.figure(9)
	1094	# #
	1095	# #
	1096	# # plt.title('DATOS SUAVIZADOS')
	1097	# # plt.xlabel('Frecuencia (KHz)')
	1098	# # plt.ylabel('Magnitud')
	1099	# # plt.plot(xFrec,SmoothSPC,'g')
	1100	# #
	1101	# # #plt.plot(xFrec,FactNorm)
	1102	# # plt.axis([-4, 4, 0, 0.15])
	1103	# # # plt.xlabel('SelfSpectra KHz')
	1104	# # #
	1105	# # plt.figure(2)
	1106	# # # #plt.subplot(121)
	1107	# # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra')
	1108	# # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano')
	1109	# # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo')
	1110	# # # #plt.plot(xFrec,phase)
	1111	# # # plt.xlabel('Suavizado, promediado KHz')
	1112	# # plt.title('SELFSPECTRA PROMEDIADO')
	1113	# # # #plt.subplot(122)
	1114	# # # #plt.plot(xSamples,zline)
	1115	# # plt.xlabel('Frecuencia (KHz)')
	1116	# # plt.ylabel('Magnitud')
	1117	# # plt.legend()
	1118	# # #
	1119	# # # plt.figure(3)
	1120	# # # plt.subplot(311)
	1121	# # # #plt.plot(xFrec,phase[0])
	1122	# # # plt.plot(xFrec,phase[0],'g')
	1123	# # # plt.subplot(312)
	1124	# # # plt.plot(xFrec,phase[1],'g')
	1125	# # # plt.subplot(313)
	1126	# # # plt.plot(xFrec,phase[2],'g')
	1127	# # # #plt.plot(xFrec,phase[2])
	1128	# # #
	1129	# # # plt.figure(4)
	1130	# # #
	1131	# # # plt.plot(xSamples,coherence[0],'b')
	1132	# # # plt.plot(xSamples,coherence[1],'r')
	1133	# # # plt.plot(xSamples,coherence[2],'g')
	1134	# # plt.show()
	1135	# # #
	1136	# # # plt.clf()
	1137	# # # plt.cla()
	1138	# # # plt.close()
	1139	#
	1140	# print ' '
	1141
	1142
	1143
	1144	self.BlockCounter+=2
	1145
	1146	else:
	1147	self.fileSelector+=1
	1148	self.BlockCounter=0
	1149	print "Next File"
	1150
	1151
	1152
	1153
	1154

schainpy/model/io/jroIO_madrigal.py created 644 +243 0

			@@ -0,0 +1,243
		1	'''
		2	Created on Aug 1, 2017
		3
		4	@author: Juan C. Espinoza
		5	'''
		6
		7	import os
		8	import sys
		9	import time
		10	import json
		11	import datetime
		12
		13	import numpy
		14
		15	try:
		16	import madrigal
		17	import madrigal.cedar
		18	except:
		19	print 'You should install "madrigal library" module if you want to read/write Madrigal data'
		20
		21	from schainpy.model.proc.jroproc_base import Operation
		22	from schainpy.model.data.jrodata import Parameters
		23
		24	MISSING = -32767
		25	DEF_CATALOG = {
		26	'principleInvestigator': 'Marco Milla',
		27	'expPurpose': None,
		28	'expMode': None,
		29	'cycleTime': None,
		30	'correlativeExp': None,
		31	'sciRemarks': None,
		32	'instRemarks': None
		33	}
		34	DEF_HEADER = {
		35	'kindatDesc': None,
		36	'analyst': 'Jicamarca User',
		37	'comments': None,
		38	'history': None
		39	}
		40	MNEMONICS = {
		41	10: 'jro',
		42	11: 'jbr',
		43	840: 'jul',
		44	13: 'jas',
		45	1000: 'pbr',
		46	1001: 'hbr',
		47	1002: 'obr',
		48	}
		49
		50	def load_json(obj):
		51	'''
		52	Parse json as string instead of unicode
		53	'''
		54
		55	if isinstance(obj, str):
		56	obj = json.loads(obj)
		57
		58	return {str(k): load_json(v) if isinstance(v, dict) else str(v) if isinstance(v, unicode) else v
		59	for k, v in obj.items()}
		60
		61
		62	class MAD2Writer(Operation):
		63
		64	def __init__(self, **kwargs):
		65
		66	Operation.__init__(self, **kwargs)
		67	self.dataOut = Parameters()
		68	self.path = None
		69	self.dataOut = None
		70	self.ext = '.dat'
		71
		72	return
		73
		74	def run(self, dataOut, path, oneDList, twoDParam='', twoDList='{}', metadata='{}', **kwargs):
		75	'''
		76	Inputs:
		77	path - path where files will be created
		78	oneDList - json of one-dimensional parameters in record where keys
		79	are Madrigal codes (integers or mnemonics) and values the corresponding
		80	dataOut attribute e.g: {
		81	'gdlatr': 'lat',
		82	'gdlonr': 'lon',
		83	'gdlat2':'lat',
		84	'glon2':'lon'}
		85	twoDParam - independent parameter to get the number of rows e.g:
		86	heighList
		87	twoDList - json of two-dimensional parameters in record where keys
		88	are Madrigal codes (integers or mnemonics) and values the corresponding
		89	dataOut attribute if multidimensional array specify as tupple
		90	('attr', pos) e.g: {
		91	'gdalt': 'heightList',
		92	'vn1p2': ('data_output', 0),
		93	'vn2p2': ('data_output', 1),
		94	'vn3': ('data_output', 2),
		95	'snl': ('data_SNR', 'db')
		96	}
		97	metadata - json of madrigal metadata (kinst, kindat, catalog and header)
		98	'''
		99	if not self.isConfig:
		100	self.setup(dataOut, path, oneDList, twoDParam, twoDList, metadata, **kwargs)
		101	self.isConfig = True
		102
		103	self.putData()
		104	return
		105
		106	def setup(self, dataOut, path, oneDList, twoDParam, twoDList, metadata, **kwargs):
		107	'''
		108	Configure Operation
		109	'''
		110
		111	self.dataOut = dataOut
		112	self.nmodes = self.dataOut.nmodes
		113	self.path = path
		114	self.blocks = kwargs.get('blocks', None)
		115	self.counter = 0
		116	self.oneDList = load_json(oneDList)
		117	self.twoDList = load_json(twoDList)
		118	self.twoDParam = twoDParam
		119	meta = load_json(metadata)
		120	self.kinst = meta.get('kinst')
		121	self.kindat = meta.get('kindat')
		122	self.catalog = meta.get('catalog', DEF_CATALOG)
		123	self.header = meta.get('header', DEF_HEADER)
		124
		125	return
		126
		127	def setFile(self):
		128	'''
		129	Create new cedar file object
		130	'''
		131
		132	self.mnemonic = MNEMONICS[self.kinst] #TODO get mnemonic from madrigal
		133	date = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
		134
		135	filename = '%s%s_%s%s' % (self.mnemonic,
		136	date.strftime('%Y%m%d_%H%M%S'),
		137	self.dataOut.mode,
		138	self.ext)
		139
		140	self.fullname = os.path.join(self.path, filename)
		141
		142	if os.path.isfile(self.fullname) :
		143	print "Destination path '%s' already exists. Previous file deleted. " %self.fullname
		144	os.remove(self.fullname)
		145
		146	try:
		147	print '[Writing] creating file : %s' % (self.fullname)
		148	self.cedarObj = madrigal.cedar.MadrigalCedarFile(self.fullname, True)
		149	except ValueError, e:
		150	print '[Error]: Impossible to create a cedar object with "madrigal.cedar.MadrigalCedarFile" '
		151	return
		152
		153	return 1
		154
		155	def writeBlock(self):
		156	'''
		157	Add data records to cedar file taking data from oneDList and twoDList
		158	attributes.
		159	Allowed parameters in: parcodes.tab
		160	'''
		161
		162	startTime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
		163	endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
		164	nrows = len(getattr(self.dataOut, self.twoDParam))
		165
		166	rec = madrigal.cedar.MadrigalDataRecord(
		167	self.kinst,
		168	self.kindat,
		169	startTime.year,
		170	startTime.month,
		171	startTime.day,
		172	startTime.hour,
		173	startTime.minute,
		174	startTime.second,
		175	startTime.microsecond/10000,
		176	endTime.year,
		177	endTime.month,
		178	endTime.day,
		179	endTime.hour,
		180	endTime.minute,
		181	endTime.second,
		182	endTime.microsecond/10000,
		183	self.oneDList.keys(),
		184	self.twoDList.keys(),
		185	nrows
		186	)
		187
		188	# Setting 1d values
		189	for key in self.oneDList:
		190	rec.set1D(key, getattr(self.dataOut, self.oneDList[key]))
		191
		192	# Setting 2d values
		193	invalid = numpy.isnan(self.dataOut.data_output)
		194	self.dataOut.data_output[invalid] = MISSING
		195	out = {}
		196	for key, value in self.twoDList.items():
		197	if isinstance(value, str):
		198	out[key] = getattr(self.dataOut, value)
		199	elif isinstance(value, tuple):
		200	attr, x = value
		201	if isinstance(x, (int, float)):
		202	out[key] = getattr(self.dataOut, attr)[int(x)]
		203	elif x.lower()=='db':
		204	tmp = getattr(self.dataOut, attr)
		205	SNRavg = numpy.average(tmp, axis=0)
		206	out[key] = 10*numpy.log10(SNRavg)
		207
		208	for n in range(nrows):
		209	for key in out:
		210	rec.set2D(key, n, out[key][n])
		211
		212	self.cedarObj.append(rec)
		213	self.cedarObj.dump()
		214	print '[Writing] Record No. {} (mode {}).'.format(
		215	self.counter,
		216	self.dataOut.mode
		217	)
		218
		219	def setHeader(self):
		220	'''
		221	Create an add catalog and header to cedar file
		222	'''
		223
		224	header = madrigal.cedar.CatalogHeaderCreator(self.fullname)
		225	header.createCatalog(**self.catalog)
		226	header.createHeader(**self.header)
		227	header.write()
		228
		229	def putData(self):
		230
		231	if self.dataOut.flagNoData:
		232	return 0
		233
		234	if self.counter == 0:
		235	self.setFile()
		236
		237	if self.counter <= self.dataOut.nrecords:
		238	self.writeBlock()
		239	self.counter += 1
		240
		241	if self.counter == self.dataOut.nrecords or self.counter == self.blocks:
		242	self.setHeader()
		243	self.counter = 0

schainpy/model/io/jroIO_mira35c.py created 644 +803 0

This diff has been collapsed as it changes many lines, (803 lines changed) Show them Hide them
		@@ -0,0 +1,803
	1	import os, sys
	2	import glob
	3	import fnmatch
	4	import datetime
	5	import time
	6	import re
	7	import h5py
	8	import numpy
	9	import matplotlib.pyplot as plt
	10
	11	import pylab as plb
	12	from scipy.optimize import curve_fit
	13	from scipy import asarray as ar,exp
	14	from scipy import stats
	15
	16	from numpy.ma.core import getdata
	17
	18	SPEED_OF_LIGHT = 299792458
	19	SPEED_OF_LIGHT = 3e8
	20
	21	try:
	22	from gevent import sleep
	23	except:
	24	from time import sleep
	25
	26	from schainpy.model.data.jrodata import Spectra
	27	#from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader
	28	from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
	29	#from schainpy.model.io.jroIO_bltr import BLTRReader
	30	from numpy import imag, shape, NaN, empty
	31
	32
	33
	34	class Header(object):
	35
	36	def __init__(self):
	37	raise NotImplementedError
	38
	39
	40	def read(self):
	41
	42	raise NotImplementedError
	43
	44	def write(self):
	45
	46	raise NotImplementedError
	47
	48	def printInfo(self):
	49
	50	message = "#"*50 + "\n"
	51	message += self.__class__.__name__.upper() + "\n"
	52	message += "#"*50 + "\n"
	53
	54	keyList = self.__dict__.keys()
	55	keyList.sort()
	56
	57	for key in keyList:
	58	message += "%s = %s" %(key, self.__dict__[key]) + "\n"
	59
	60	if "size" not in keyList:
	61	attr = getattr(self, "size")
	62
	63	if attr:
	64	message += "%s = %s" %("size", attr) + "\n"
	65
	66	#print message
	67
	68
	69	FILE_HEADER = numpy.dtype([ #HEADER 1024bytes
	70	('Hname','a32'), #Original file name
	71	('Htime',numpy.str_,32), #Date and time when the file was created
	72	('Hoper',numpy.str_,64), #Name of operator who created the file
	73	('Hplace',numpy.str_,128), #Place where the measurements was carried out
	74	('Hdescr',numpy.str_,256), #Description of measurements
	75	('Hdummy',numpy.str_,512), #Reserved space
	76	#Main chunk 8bytes
	77	('Msign',numpy.str_,4), #Main chunk signature FZKF or NUIG
	78	('MsizeData','<i4'), #Size of data block main chunk
	79	#Processing DSP parameters 36bytes
	80	('PPARsign',numpy.str_,4), #PPAR signature
	81	('PPARsize','<i4'), #PPAR size of block
	82	('PPARprf','<i4'), #Pulse repetition frequency
	83	('PPARpdr','<i4'), #Pulse duration
	84	('PPARsft','<i4'), #FFT length
	85	('PPARavc','<i4'), #Number of spectral (in-coherent) averages
	86	('PPARihp','<i4'), #Number of lowest range gate for moment estimation
	87	('PPARchg','<i4'), #Count for gates for moment estimation
	88	('PPARpol','<i4'), #switch on/off polarimetric measurements. Should be 1.
	89	#Service DSP parameters 112bytes
	90	('SPARatt','<i4'), #STC attenuation on the lowest ranges on/off
	91	('SPARtx','<i4'), #OBSOLETE
	92	('SPARaddGain0','<f4'), #OBSOLETE
	93	('SPARaddGain1','<f4'), #OBSOLETE
	94	('SPARwnd','<i4'), #Debug only. It normal mode it is 0.
	95	('SPARpos','<i4'), #Delay between sync pulse and tx pulse for phase corr, ns
	96	('SPARadd','<i4'), #"add to pulse" to compensate for delay between the leading edge of driver pulse and envelope of the RF signal.
	97	('SPARlen','<i4'), #Time for measuring txn pulse phase. OBSOLETE
	98	('SPARcal','<i4'), #OBSOLETE
	99	('SPARnos','<i4'), #OBSOLETE
	100	('SPARof0','<i4'), #detection threshold
	101	('SPARof1','<i4'), #OBSOLETE
	102	('SPARswt','<i4'), #2nd moment estimation threshold
	103	('SPARsum','<i4'), #OBSOLETE
	104	('SPARosc','<i4'), #flag Oscillosgram mode
	105	('SPARtst','<i4'), #OBSOLETE
	106	('SPARcor','<i4'), #OBSOLETE
	107	('SPARofs','<i4'), #OBSOLETE
	108	('SPARhsn','<i4'), #Hildebrand div noise detection on noise gate
	109	('SPARhsa','<f4'), #Hildebrand div noise detection on all gates
	110	('SPARcalibPow_M','<f4'), #OBSOLETE
	111	('SPARcalibSNR_M','<f4'), #OBSOLETE
	112	('SPARcalibPow_S','<f4'), #OBSOLETE
	113	('SPARcalibSNR_S','<f4'), #OBSOLETE
	114	('SPARrawGate1','<i4'), #Lowest range gate for spectra saving Raw_Gate1 >=5
	115	('SPARrawGate2','<i4'), #Number of range gates with atmospheric signal
	116	('SPARraw','<i4'), #flag - IQ or spectra saving on/off
	117	('SPARprc','<i4'),]) #flag - Moment estimation switched on/off
	118
	119
	120
	121	class FileHeaderMIRA35c(Header):
	122
	123	def __init__(self):
	124
	125	self.Hname= None
	126	self.Htime= None
	127	self.Hoper= None
	128	self.Hplace= None
	129	self.Hdescr= None
	130	self.Hdummy= None
	131
	132	self.Msign=None
	133	self.MsizeData=None
	134
	135	self.PPARsign=None
	136	self.PPARsize=None
	137	self.PPARprf=None
	138	self.PPARpdr=None
	139	self.PPARsft=None
	140	self.PPARavc=None
	141	self.PPARihp=None
	142	self.PPARchg=None
	143	self.PPARpol=None
	144	#Service DSP parameters
	145	self.SPARatt=None
	146	self.SPARtx=None
	147	self.SPARaddGain0=None
	148	self.SPARaddGain1=None
	149	self.SPARwnd=None
	150	self.SPARpos=None
	151	self.SPARadd=None
	152	self.SPARlen=None
	153	self.SPARcal=None
	154	self.SPARnos=None
	155	self.SPARof0=None
	156	self.SPARof1=None
	157	self.SPARswt=None
	158	self.SPARsum=None
	159	self.SPARosc=None
	160	self.SPARtst=None
	161	self.SPARcor=None
	162	self.SPARofs=None
	163	self.SPARhsn=None
	164	self.SPARhsa=None
	165	self.SPARcalibPow_M=None
	166	self.SPARcalibSNR_M=None
	167	self.SPARcalibPow_S=None
	168	self.SPARcalibSNR_S=None
	169	self.SPARrawGate1=None
	170	self.SPARrawGate2=None
	171	self.SPARraw=None
	172	self.SPARprc=None
	173
	174	self.FHsize=1180
	175
	176	def FHread(self, fp):
	177
	178	header = numpy.fromfile(fp, FILE_HEADER,1)
	179	''' numpy.fromfile(file, dtype, count, sep='')
	180	file : file or str
	181	Open file object or filename.
	182
	183	dtype : data-type
	184	Data type of the returned array. For binary files, it is used to determine
	185	the size and byte-order of the items in the file.
	186
	187	count : int
	188	Number of items to read. -1 means all items (i.e., the complete file).
	189
	190	sep : str
	191	Separator between items if file is a text file. Empty ("") separator means
	192	the file should be treated as binary. Spaces (" ") in the separator match zero
	193	or more whitespace characters. A separator consisting only of spaces must match
	194	at least one whitespace.
	195
	196	'''
	197
	198
	199	self.Hname= str(header['Hname'][0])
	200	self.Htime= str(header['Htime'][0])
	201	self.Hoper= str(header['Hoper'][0])
	202	self.Hplace= str(header['Hplace'][0])
	203	self.Hdescr= str(header['Hdescr'][0])
	204	self.Hdummy= str(header['Hdummy'][0])
	205	#1024
	206
	207	self.Msign=str(header['Msign'][0])
	208	self.MsizeData=header['MsizeData'][0]
	209	#8
	210
	211	self.PPARsign=str(header['PPARsign'][0])
	212	self.PPARsize=header['PPARsize'][0]
	213	self.PPARprf=header['PPARprf'][0]
	214	self.PPARpdr=header['PPARpdr'][0]
	215	self.PPARsft=header['PPARsft'][0]
	216	self.PPARavc=header['PPARavc'][0]
	217	self.PPARihp=header['PPARihp'][0]
	218	self.PPARchg=header['PPARchg'][0]
	219	self.PPARpol=header['PPARpol'][0]
	220	#Service DSP parameters
	221	#36
	222
	223	self.SPARatt=header['SPARatt'][0]
	224	self.SPARtx=header['SPARtx'][0]
	225	self.SPARaddGain0=header['SPARaddGain0'][0]
	226	self.SPARaddGain1=header['SPARaddGain1'][0]
	227	self.SPARwnd=header['SPARwnd'][0]
	228	self.SPARpos=header['SPARpos'][0]
	229	self.SPARadd=header['SPARadd'][0]
	230	self.SPARlen=header['SPARlen'][0]
	231	self.SPARcal=header['SPARcal'][0]
	232	self.SPARnos=header['SPARnos'][0]
	233	self.SPARof0=header['SPARof0'][0]
	234	self.SPARof1=header['SPARof1'][0]
	235	self.SPARswt=header['SPARswt'][0]
	236	self.SPARsum=header['SPARsum'][0]
	237	self.SPARosc=header['SPARosc'][0]
	238	self.SPARtst=header['SPARtst'][0]
	239	self.SPARcor=header['SPARcor'][0]
	240	self.SPARofs=header['SPARofs'][0]
	241	self.SPARhsn=header['SPARhsn'][0]
	242	self.SPARhsa=header['SPARhsa'][0]
	243	self.SPARcalibPow_M=header['SPARcalibPow_M'][0]
	244	self.SPARcalibSNR_M=header['SPARcalibSNR_M'][0]
	245	self.SPARcalibPow_S=header['SPARcalibPow_S'][0]
	246	self.SPARcalibSNR_S=header['SPARcalibSNR_S'][0]
	247	self.SPARrawGate1=header['SPARrawGate1'][0]
	248	self.SPARrawGate2=header['SPARrawGate2'][0]
	249	self.SPARraw=header['SPARraw'][0]
	250	self.SPARprc=header['SPARprc'][0]
	251	#112
	252	#1180
	253	#print 'Pointer fp header', fp.tell()
	254	#print ' '
	255	#print 'SPARrawGate'
	256	#print self.SPARrawGate2 - self.SPARrawGate1
	257
	258	#print ' '
	259	#print 'Hname'
	260	#print self.Hname
	261
	262	#print ' '
	263	#print 'Msign'
	264	#print self.Msign
	265
	266	def write(self, fp):
	267
	268	headerTuple = (self.Hname,
	269	self.Htime,
	270	self.Hoper,
	271	self.Hplace,
	272	self.Hdescr,
	273	self.Hdummy)
	274
	275
	276	header = numpy.array(headerTuple, FILE_HEADER)
	277	# numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0)
	278	header.tofile(fp)
	279	''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
	280
	281	fid : file or str
	282	An open file object, or a string containing a filename.
	283
	284	sep : str
	285	Separator between array items for text output. If "" (empty), a binary file is written,
	286	equivalent to file.write(a.tobytes()).
	287
	288	format : str
	289	Format string for text file output. Each entry in the array is formatted to text by
	290	first converting it to the closest Python type, and then using "format" % item.
	291
	292	'''
	293
	294	return 1
	295
	296	SRVI_HEADER = numpy.dtype([
	297	('SignatureSRVI1',numpy.str_,4),#
	298	('SizeOfDataBlock1','<i4'),#
	299	('DataBlockTitleSRVI1',numpy.str_,4),#
	300	('SizeOfSRVI1','<i4'),])#
	301
	302	class SRVIHeader(Header):
	303	def __init__(self, SignatureSRVI1=0, SizeOfDataBlock1=0, DataBlockTitleSRVI1=0, SizeOfSRVI1=0):
	304
	305	self.SignatureSRVI1 = SignatureSRVI1
	306	self.SizeOfDataBlock1 = SizeOfDataBlock1
	307	self.DataBlockTitleSRVI1 = DataBlockTitleSRVI1
	308	self.SizeOfSRVI1 = SizeOfSRVI1
	309
	310	self.SRVIHsize=16
	311
	312	def SRVIread(self, fp):
	313
	314	header = numpy.fromfile(fp, SRVI_HEADER,1)
	315
	316	self.SignatureSRVI1 = str(header['SignatureSRVI1'][0])
	317	self.SizeOfDataBlock1 = header['SizeOfDataBlock1'][0]
	318	self.DataBlockTitleSRVI1 = str(header['DataBlockTitleSRVI1'][0])
	319	self.SizeOfSRVI1 = header['SizeOfSRVI1'][0]
	320	#16
	321	print 'Pointer fp SRVIheader', fp.tell()
	322
	323
	324	SRVI_STRUCTURE = numpy.dtype([
	325	('frame_cnt','<u4'),#
	326	('time_t','<u4'), #
	327	('tpow','<f4'), #
	328	('npw1','<f4'), #
	329	('npw2','<f4'), #
	330	('cpw1','<f4'), #
	331	('pcw2','<f4'), #
	332	('ps_err','<u4'), #
	333	('te_err','<u4'), #
	334	('rc_err','<u4'), #
	335	('grs1','<u4'), #
	336	('grs2','<u4'), #
	337	('azipos','<f4'), #
	338	('azivel','<f4'), #
	339	('elvpos','<f4'), #
	340	('elvvel','<f4'), #
	341	('northAngle','<f4'), #
	342	('microsec','<u4'), #
	343	('azisetvel','<f4'), #
	344	('elvsetpos','<f4'), #
	345	('RadarConst','<f4'),]) #
	346
	347
	348
	349
	350	class RecordHeader(Header):
	351
	352
	353	def __init__(self, frame_cnt=0, time_t= 0, tpow=0, npw1=0, npw2=0,
	354	cpw1=0, pcw2=0, ps_err=0, te_err=0, rc_err=0, grs1=0,
	355	grs2=0, azipos=0, azivel=0, elvpos=0, elvvel=0, northangle=0,
	356	microsec=0, azisetvel=0, elvsetpos=0, RadarConst=0 , RecCounter=0, Off2StartNxtRec=0):
	357
	358
	359	self.frame_cnt = frame_cnt
	360	self.dwell = time_t
	361	self.tpow = tpow
	362	self.npw1 = npw1
	363	self.npw2 = npw2
	364	self.cpw1 = cpw1
	365	self.pcw2 = pcw2
	366	self.ps_err = ps_err
	367	self.te_err = te_err
	368	self.rc_err = rc_err
	369	self.grs1 = grs1
	370	self.grs2 = grs2
	371	self.azipos = azipos
	372	self.azivel = azivel
	373	self.elvpos = elvpos
	374	self.elvvel = elvvel
	375	self.northAngle = northangle
	376	self.microsec = microsec
	377	self.azisetvel = azisetvel
	378	self.elvsetpos = elvsetpos
	379	self.RadarConst = RadarConst
	380	self.RHsize=84
	381	self.RecCounter = RecCounter
	382	self.Off2StartNxtRec=Off2StartNxtRec
	383
	384	def RHread(self, fp):
	385
	386	#startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
	387
	388	#OffRHeader= 1180 + self.RecCounter*(self.Off2StartNxtRec)
	389	#startFp.seek(OffRHeader, os.SEEK_SET)
	390
	391	#print 'Posicion del bloque: ',OffRHeader
	392
	393	header = numpy.fromfile(fp,SRVI_STRUCTURE,1)
	394
	395	self.frame_cnt = header['frame_cnt'][0]#
	396	self.time_t = header['time_t'][0] #
	397	self.tpow = header['tpow'][0] #
	398	self.npw1 = header['npw1'][0] #
	399	self.npw2 = header['npw2'][0] #
	400	self.cpw1 = header['cpw1'][0] #
	401	self.pcw2 = header['pcw2'][0] #
	402	self.ps_err = header['ps_err'][0] #
	403	self.te_err = header['te_err'][0] #
	404	self.rc_err = header['rc_err'][0] #
	405	self.grs1 = header['grs1'][0] #
	406	self.grs2 = header['grs2'][0] #
	407	self.azipos = header['azipos'][0] #
	408	self.azivel = header['azivel'][0] #
	409	self.elvpos = header['elvpos'][0] #
	410	self.elvvel = header['elvvel'][0] #
	411	self.northAngle = header['northAngle'][0] #
	412	self.microsec = header['microsec'][0] #
	413	self.azisetvel = header['azisetvel'][0] #
	414	self.elvsetpos = header['elvsetpos'][0] #
	415	self.RadarConst = header['RadarConst'][0] #
	416	#84
	417
	418	#print 'Pointer fp RECheader', fp.tell()
	419
	420	#self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
	421
	422	#self.RHsize = 180+20*self.nChannels
	423	#self.Datasize= self.nProfilesself.nChannelsself.nHeights24
	424	#print 'Datasize',self.Datasize
	425	#endFp = self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
	426
	427	print '=============================================='
	428
	429	print '=============================================='
	430
	431
	432	return 1
	433
	434	class MIRA35CReader (ProcessingUnit,FileHeaderMIRA35c,SRVIHeader,RecordHeader):
	435
	436	path = None
	437	startDate = None
	438	endDate = None
	439	startTime = None
	440	endTime = None
	441	walk = None
	442	isConfig = False
	443
	444
	445	fileList= None
	446
	447	#metadata
	448	TimeZone= None
	449	Interval= None
	450	heightList= None
	451
	452	#data
	453	data= None
	454	utctime= None
	455
	456
	457
	458	def __init__(self, **kwargs):
	459
	460	#Eliminar de la base la herencia
	461	ProcessingUnit.__init__(self, **kwargs)
	462	self.PointerReader = 0
	463	self.FileHeaderFlag = False
	464	self.utc = None
	465	self.ext = ".zspca"
	466	self.optchar = "P"
	467	self.fpFile=None
	468	self.fp = None
	469	self.BlockCounter=0
	470	self.dtype = None
	471	self.fileSizeByHeader = None
	472	self.filenameList = []
	473	self.fileSelector = 0
	474	self.Off2StartNxtRec=0
	475	self.RecCounter=0
	476	self.flagNoMoreFiles = 0
	477	self.data_spc=None
	478	#self.data_cspc=None
	479	self.data_output=None
	480	self.path = None
	481	self.OffsetStartHeader=0
	482	self.Off2StartData=0
	483	self.ipp = 0
	484	self.nFDTdataRecors=0
	485	self.blocksize = 0
	486	self.dataOut = Spectra()
	487	self.profileIndex = 1 #Always
	488	self.dataOut.flagNoData=False
	489	self.dataOut.nRdPairs = 0
	490	self.dataOut.pairsList = []
	491	self.dataOut.data_spc=None
	492
	493	self.dataOut.normFactor=1
	494	self.nextfileflag = True
	495	self.dataOut.RadarConst = 0
	496	self.dataOut.HSDV = []
	497	self.dataOut.NPW = []
	498	self.dataOut.COFA = []
	499	self.dataOut.noise = 0
	500
	501
	502	def Files2Read(self, fp):
	503	'''
	504	Function that indicates the number of .fdt files that exist in the folder to be read.
	505	It also creates an organized list with the names of the files to read.
	506	'''
	507	#self.__checkPath()
	508
	509	ListaData=os.listdir(fp) #Gets the list of files within the fp address
	510	ListaData=sorted(ListaData) #Sort the list of files from least to largest by names
	511	nFiles=0 #File Counter
	512	FileList=[] #A list is created that will contain the .fdt files
	513	for IndexFile in ListaData :
	514	if '.zspca' in IndexFile and '.gz' not in IndexFile:
	515	FileList.append(IndexFile)
	516	nFiles+=1
	517
	518	#print 'Files2Read'
	519	#print 'Existen '+str(nFiles)+' archivos .fdt'
	520
	521	self.filenameList=FileList #List of files from least to largest by names
	522
	523
	524	def run(self, **kwargs):
	525	'''
	526	This method will be the one that will initiate the data entry, will be called constantly.
	527	You should first verify that your Setup () is set up and then continue to acquire
	528	the data to be processed with getData ().
	529	'''
	530	if not self.isConfig:
	531	self.setup(**kwargs)
	532	self.isConfig = True
	533
	534	self.getData()
	535
	536
	537	def setup(self, path=None,
	538	startDate=None,
	539	endDate=None,
	540	startTime=None,
	541	endTime=None,
	542	walk=True,
	543	timezone='utc',
	544	code = None,
	545	online=False,
	546	ReadMode=None, **kwargs):
	547
	548	self.isConfig = True
	549
	550	self.path=path
	551	self.startDate=startDate
	552	self.endDate=endDate
	553	self.startTime=startTime
	554	self.endTime=endTime
	555	self.walk=walk
	556	#self.ReadMode=int(ReadMode)
	557
	558	pass
	559
	560
	561	def getData(self):
	562	'''
	563	Before starting this function, you should check that there is still an unread file,
	564	If there are still blocks to read or if the data block is empty.
	565
	566	You should call the file "read".
	567
	568	'''
	569
	570	if self.flagNoMoreFiles:
	571	self.dataOut.flagNoData = True
	572	print 'NoData se vuelve true'
	573	return 0
	574
	575	self.fp=self.path
	576	self.Files2Read(self.fp)
	577	self.readFile(self.fp)
	578
	579	self.dataOut.data_spc = self.dataOut_spc#self.data_spc.copy()
	580	self.dataOut.RadarConst = self.RadarConst
	581	self.dataOut.data_output=self.data_output
	582	self.dataOut.noise = self.dataOut.getNoise()
	583	#print 'ACAAAAAA', self.dataOut.noise
	584	self.dataOut.data_spc = self.dataOut.data_spc+self.dataOut.noise
	585	#print 'self.dataOut.noise',self.dataOut.noise
	586
	587
	588	return self.dataOut.data_spc
	589
	590
	591	def readFile(self,fp):
	592	'''
	593	You must indicate if you are reading in Online or Offline mode and load the
	594	The parameters for this file reading mode.
	595
	596	Then you must do 2 actions:
	597
	598	1. Get the BLTR FileHeader.
	599	2. Start reading the first block.
	600	'''
	601
	602	#The address of the folder is generated the name of the .fdt file that will be read
	603	print "File: ",self.fileSelector+1
	604
	605	if self.fileSelector < len(self.filenameList):
	606
	607	self.fpFile=str(fp)+'/'+str(self.filenameList[self.fileSelector])
	608
	609	if self.nextfileflag==True:
	610	self.fp = open(self.fpFile,"rb")
	611	self.nextfileflag==False
	612
	613	'''HERE STARTING THE FILE READING'''
	614
	615
	616	self.fheader = FileHeaderMIRA35c()
	617	self.fheader.FHread(self.fp) #Bltr FileHeader Reading
	618
	619
	620	self.SPARrawGate1 = self.fheader.SPARrawGate1
	621	self.SPARrawGate2 = self.fheader.SPARrawGate2
	622	self.Num_Hei = self.SPARrawGate2 - self.SPARrawGate1
	623	self.Num_Bins = self.fheader.PPARsft
	624	self.dataOut.nFFTPoints = self.fheader.PPARsft
	625
	626
	627	self.Num_inCoh = self.fheader.PPARavc
	628	self.dataOut.PRF = self.fheader.PPARprf
	629	self.dataOut.frequency = 34.8510*9
	630	self.Lambda = SPEED_OF_LIGHT/self.dataOut.frequency
	631	self.dataOut.ippSeconds= 1./float(self.dataOut.PRF)
	632
	633	pulse_width = self.fheader.PPARpdr * 10**-9
	634	self.__deltaHeigth = 0.5 * SPEED_OF_LIGHT * pulse_width
	635
	636	self.data_spc = numpy.zeros((self.Num_Hei, self.Num_Bins,2))#
	637	self.dataOut.HSDV = numpy.zeros((self.Num_Hei, 2))
	638
	639	self.Ze = numpy.zeros(self.Num_Hei)
	640	self.ETA = numpy.zeros(([2,self.Num_Hei]))
	641
	642
	643
	644	self.readBlock() #Block reading
	645
	646	else:
	647	print 'readFile FlagNoData becomes true'
	648	self.flagNoMoreFiles=True
	649	self.dataOut.flagNoData = True
	650	self.FileHeaderFlag == True
	651	return 0
	652
	653
	654
	655	def readBlock(self):
	656	'''
	657	It should be checked if the block has data, if it is not passed to the next file.
	658
	659	Then the following is done:
	660
	661	1. Read the RecordHeader
	662	2. Fill the buffer with the current block number.
	663
	664	'''
	665
	666	if self.PointerReader > 1180:
	667	self.fp.seek(self.PointerReader , os.SEEK_SET)
	668	self.FirstPoint = self.PointerReader
	669
	670	else :
	671	self.FirstPoint = 1180
	672
	673
	674
	675	self.srviHeader = SRVIHeader()
	676
	677	self.srviHeader.SRVIread(self.fp) #Se obtiene la cabecera del SRVI
	678
	679	self.blocksize = self.srviHeader.SizeOfDataBlock1 # Se obtiene el tamao del bloque
	680
	681	if self.blocksize == 148:
	682	print 'blocksize == 148 bug'
	683	jump = numpy.fromfile(self.fp,[('jump',numpy.str_,140)] ,1)
	684
	685	self.srviHeader.SRVIread(self.fp) #Se obtiene la cabecera del SRVI
	686
	687	if not self.srviHeader.SizeOfSRVI1:
	688	self.fileSelector+=1
	689	self.nextfileflag==True
	690	self.FileHeaderFlag == True
	691
	692	self.recordheader = RecordHeader()
	693	self.recordheader.RHread(self.fp)
	694	self.RadarConst = self.recordheader.RadarConst
	695	dwell = self.recordheader.time_t
	696	npw1 = self.recordheader.npw1
	697	npw2 = self.recordheader.npw2
	698
	699
	700	self.dataOut.channelList = range(1)
	701	self.dataOut.nIncohInt = self.Num_inCoh
	702	self.dataOut.nProfiles = self.Num_Bins
	703	self.dataOut.nCohInt = 1
	704	self.dataOut.windowOfFilter = 1
	705	self.dataOut.utctime = dwell
	706	self.dataOut.timeZone=0
	707
	708	self.dataOut.outputInterval = self.dataOut.getTimeInterval()
	709	self.dataOut.heightList = self.SPARrawGate1self.__deltaHeigth + numpy.array(range(self.Num_Hei))self.__deltaHeigth
	710
	711
	712
	713	self.HSDVsign = numpy.fromfile( self.fp, [('HSDV',numpy.str_,4)],1)
	714	self.SizeHSDV = numpy.fromfile( self.fp, [('SizeHSDV','<i4')],1)
	715	self.HSDV_Co = numpy.fromfile( self.fp, [('HSDV_Co','<f4')],self.Num_Hei)
	716	self.HSDV_Cx = numpy.fromfile( self.fp, [('HSDV_Cx','<f4')],self.Num_Hei)
	717
	718	self.COFAsign = numpy.fromfile( self.fp, [('COFA',numpy.str_,4)],1)
	719	self.SizeCOFA = numpy.fromfile( self.fp, [('SizeCOFA','<i4')],1)
	720	self.COFA_Co = numpy.fromfile( self.fp, [('COFA_Co','<f4')],self.Num_Hei)
	721	self.COFA_Cx = numpy.fromfile( self.fp, [('COFA_Cx','<f4')],self.Num_Hei)
	722
	723	self.ZSPCsign = numpy.fromfile(self.fp, [('ZSPCsign',numpy.str_,4)],1)
	724	self.SizeZSPC = numpy.fromfile(self.fp, [('SizeZSPC','<i4')],1)
	725
	726	self.dataOut.HSDV[0]=self.HSDV_Co[:][0]
	727	self.dataOut.HSDV[1]=self.HSDV_Cx[:][0]
	728
	729	for irg in range(self.Num_Hei):
	730	nspc = numpy.fromfile(self.fp, [('nspc','int16')],1)[0][0] # Number of spectral sub pieces containing significant power
	731
	732	for k in range(nspc):
	733	binIndex = numpy.fromfile(self.fp, [('binIndex','int16')],1)[0][0] # Index of the spectral bin where the piece is beginning
	734	nbins = numpy.fromfile(self.fp, [('nbins','int16')],1)[0][0] # Number of bins of the piece
	735
	736	#Co_Channel
	737	jbin = numpy.fromfile(self.fp, [('jbin','uint16')],nbins)[0][0] # Spectrum piece to be normaliced
	738	jmax = numpy.fromfile(self.fp, [('jmax','float32')],1)[0][0] # Maximun piece to be normaliced
	739
	740
	741	self.data_spc[irg,binIndex:binIndex+nbins,0] = self.data_spc[irg,binIndex:binIndex+nbins,0]+jbin/65530.*jmax
	742
	743	#Cx_Channel
	744	jbin = numpy.fromfile(self.fp, [('jbin','uint16')],nbins)[0][0]
	745	jmax = numpy.fromfile(self.fp, [('jmax','float32')],1)[0][0]
	746
	747
	748	self.data_spc[irg,binIndex:binIndex+nbins,1] = self.data_spc[irg,binIndex:binIndex+nbins,1]+jbin/65530.*jmax
	749
	750	for bin in range(self.Num_Bins):
	751
	752	self.data_spc[:,bin,0] = self.data_spc[:,bin,0] - self.dataOut.HSDV[:,0]
	753
	754	self.data_spc[:,bin,1] = self.data_spc[:,bin,1] - self.dataOut.HSDV[:,1]
	755
	756
	757	numpy.set_printoptions(threshold='nan')
	758
	759	self.data_spc = numpy.where(self.data_spc > 0. , self.data_spc, 0)
	760
	761	self.dataOut.COFA = numpy.array([self.COFA_Co , self.COFA_Cx])
	762
	763	print ' '
	764	print 'SPC',numpy.shape(self.dataOut.data_spc)
	765	#print 'SPC',self.dataOut.data_spc
	766
	767	noinor1 = 713031680
	768	noinor2 = 30
	769
	770	npw1 = 1#0*(npw1/10) noinor1 * noinor2
	771	npw2 = 1#0*(npw2/10) noinor1 * noinor2
	772	self.dataOut.NPW = numpy.array([npw1, npw2])
	773
	774	print ' '
	775
	776	self.data_spc = numpy.transpose(self.data_spc, (2,1,0))
	777	self.data_spc = numpy.fft.fftshift(self.data_spc, axes = 1)
	778
	779	self.data_spc = numpy.fliplr(self.data_spc)
	780
	781	self.data_spc = numpy.where(self.data_spc > 0. , self.data_spc, 0)
	782	self.dataOut_spc= numpy.ones([1, self.Num_Bins , self.Num_Hei])
	783	self.dataOut_spc[0,:,:] = self.data_spc[0,:,:]
	784	#print 'SHAPE', self.dataOut_spc.shape
	785	#For nyquist correction:
	786	#fix = 20 # ~3m/s
	787	#shift = self.Num_Bins/2 + fix
	788	#self.data_spc = numpy.array([ self.data_spc[: , self.Num_Bins-shift+1: , :] , self.data_spc[: , 0:self.Num_Bins-shift , :]])
	789
	790
	791
	792	'''Block Reading, the Block Data is received and Reshape is used to give it
	793	shape.
	794	'''
	795
	796	self.PointerReader = self.fp.tell()
	797
	798
	799
	800
	801
	802
	803	No newline at end of file

schainpy/model/proc/bltrproc_parameters.py created 644 +403 0

			@@ -0,0 +1,403
		1	'''
		2	Created on Oct 24, 2016
		3
		4	@author: roj- LouVD
		5	'''
		6
		7	import numpy
		8	import copy
		9	import datetime
		10	import time
		11	from time import gmtime
		12
		13	from numpy import transpose
		14
		15	from jroproc_base import ProcessingUnit, Operation
		16	from schainpy.model.data.jrodata import Parameters
		17
		18
		19	class BLTRParametersProc(ProcessingUnit):
		20	'''
		21	Processing unit for BLTR parameters data (winds)
		22
		23	Inputs:
		24	self.dataOut.nmodes - Number of operation modes
		25	self.dataOut.nchannels - Number of channels
		26	self.dataOut.nranges - Number of ranges
		27
		28	self.dataOut.data_SNR - SNR array
		29	self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
		30	self.dataOut.height - Height array (km)
		31	self.dataOut.time - Time array (seconds)
		32
		33	self.dataOut.fileIndex -Index of the file currently read
		34	self.dataOut.lat - Latitude coordinate of BLTR location
		35
		36	self.dataOut.doy - Experiment doy (number of the day in the current year)
		37	self.dataOut.month - Experiment month
		38	self.dataOut.day - Experiment day
		39	self.dataOut.year - Experiment year
		40	'''
		41
		42	def __init__(self, **kwargs):
		43	'''
		44	Inputs: None
		45	'''
		46	ProcessingUnit.__init__(self, **kwargs)
		47	self.dataOut = Parameters()
		48	self.isConfig = False
		49
		50	def setup(self, mode):
		51	'''
		52	'''
		53	self.dataOut.mode = mode
		54
		55	def run(self, mode, snr_threshold=None):
		56	'''
		57	Inputs:
		58	mode = High resolution (0) or Low resolution (1) data
		59	snr_threshold = snr filter value
		60	'''
		61
		62	if not self.isConfig:
		63	self.setup(mode)
		64	self.isConfig = True
		65
		66	if self.dataIn.type == 'Parameters':
		67	self.dataOut.copy(self.dataIn)
		68
		69	self.dataOut.data_output = self.dataOut.data_output[mode]
		70	self.dataOut.heightList = self.dataOut.height[0]
		71	self.dataOut.data_SNR = self.dataOut.data_SNR[mode]
		72
		73	if snr_threshold is not None:
		74	SNRavg = numpy.average(self.dataOut.data_SNR, axis=0)
		75	SNRavgdB = 10*numpy.log10(SNRavg)
		76	for i in range(3):
		77	self.dataOut.data_output[i][SNRavgdB <= snr_threshold] = numpy.nan
		78
		79	# TODO
		80	class OutliersFilter(Operation):
		81
		82	def __init__(self, **kwargs):
		83	'''
		84	'''
		85	Operation.__init__(self, **kwargs)
		86
		87	def run(self, svalue2, method, factor, filter, npoints=9):
		88	'''
		89	Inputs:
		90	svalue - string to select array velocity
		91	svalue2 - string to choose axis filtering
		92	method - 0 for SMOOTH or 1 for MEDIAN
		93	factor - number used to set threshold
		94	filter - 1 for data filtering using the standard deviation criteria else 0
		95	npoints - number of points for mask filter
		96	'''
		97
		98	print ' Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor)
		99
		100
		101	yaxis = self.dataOut.heightList
		102	xaxis = numpy.array([[self.dataOut.utctime]])
		103
		104	# Zonal
		105	value_temp = self.dataOut.data_output[0]
		106
		107	# Zonal
		108	value_temp = self.dataOut.data_output[1]
		109
		110	# Vertical
		111	value_temp = numpy.transpose(self.dataOut.data_output[2])
		112
		113	htemp = yaxis
		114	std = value_temp
		115	for h in range(len(htemp)):
		116	nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
		117	minvalid = npoints
		118
		119	#only if valid values greater than the minimum required (10%)
		120	if nvalues_valid > minvalid:
		121
		122	if method == 0:
		123	#SMOOTH
		124	w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
		125
		126
		127	if method == 1:
		128	#MEDIAN
		129	w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
		130
		131	dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
		132
		133	threshold = dw*factor
		134	value_temp[numpy.where(w > threshold),h] = numpy.nan
		135	value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
		136
		137
		138	#At the end
		139	if svalue2 == 'inHeight':
		140	value_temp = numpy.transpose(value_temp)
		141	output_array[:,m] = value_temp
		142
		143	if svalue == 'zonal':
		144	self.dataOut.data_output[0] = output_array
		145
		146	elif svalue == 'meridional':
		147	self.dataOut.data_output[1] = output_array
		148
		149	elif svalue == 'vertical':
		150	self.dataOut.data_output[2] = output_array
		151
		152	return self.dataOut.data_output
		153
		154
		155	def Median(self,input,width):
		156	'''
		157	Inputs:
		158	input - Velocity array
		159	width - Number of points for mask filter
		160
		161	'''
		162
		163	if numpy.mod(width,2) == 1:
		164	pc = int((width - 1) / 2)
		165	cont = 0
		166	output = []
		167
		168	for i in range(len(input)):
		169	if i >= pc and i < len(input) - pc:
		170	new2 = input[i-pc:i+pc+1]
		171	temp = numpy.where(numpy.isfinite(new2))
		172	new = new2[temp]
		173	value = numpy.median(new)
		174	output.append(value)
		175
		176	output = numpy.array(output)
		177	output = numpy.hstack((input[0:pc],output))
		178	output = numpy.hstack((output,input[-pc:len(input)]))
		179
		180	return output
		181
		182	def Smooth(self,input,width,edge_truncate = None):
		183	'''
		184	Inputs:
		185	input - Velocity array
		186	width - Number of points for mask filter
		187	edge_truncate - 1 for truncate the convolution product else
		188
		189	'''
		190
		191	if numpy.mod(width,2) == 0:
		192	real_width = width + 1
		193	nzeros = width / 2
		194	else:
		195	real_width = width
		196	nzeros = (width - 1) / 2
		197
		198	half_width = int(real_width)/2
		199	length = len(input)
		200
		201	gate = numpy.ones(real_width,dtype='float')
		202	norm_of_gate = numpy.sum(gate)
		203
		204	nan_process = 0
		205	nan_id = numpy.where(numpy.isnan(input))
		206	if len(nan_id[0]) > 0:
		207	nan_process = 1
		208	pb = numpy.zeros(len(input))
		209	pb[nan_id] = 1.
		210	input[nan_id] = 0.
		211
		212	if edge_truncate == True:
		213	output = numpy.convolve(input/norm_of_gate,gate,mode='same')
		214	elif edge_truncate == False or edge_truncate == None:
		215	output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
		216	output = numpy.hstack((input[0:half_width],output))
		217	output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
		218
		219	if nan_process:
		220	pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
		221	pb = numpy.hstack((numpy.zeros(half_width),pb))
		222	pb = numpy.hstack((pb,numpy.zeros(half_width)))
		223	output[numpy.where(pb > 0.9999)] = numpy.nan
		224	input[nan_id] = numpy.nan
		225	return output
		226
		227	def Average(self,aver=0,nhaver=1):
		228	'''
		229	Inputs:
		230	aver - Indicates the time period over which is averaged or consensus data
		231	nhaver - Indicates the decimation factor in heights
		232
		233	'''
		234	nhpoints = 48
		235
		236	lat_piura = -5.17
		237	lat_huancayo = -12.04
		238	lat_porcuya = -5.8
		239
		240	if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
		241	hcm = 3.
		242	if self.dataOut.year == 2003 :
		243	if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
		244	nhpoints = 12
		245
		246	elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
		247	hcm = 3.
		248	if self.dataOut.year == 2003 :
		249	if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
		250	nhpoints = 12
		251
		252
		253	elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
		254	hcm = 5.#2
		255
		256	pdata = 0.2
		257	taver = [1,2,3,4,6,8,12,24]
		258	t0 = 0
		259	tf = 24
		260	ntime =(tf-t0)/taver[aver]
		261	ti = numpy.arange(ntime)
		262	tf = numpy.arange(ntime) + taver[aver]
		263
		264
		265	old_height = self.dataOut.heightList
		266
		267	if nhaver > 1:
		268	num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
		269	deltha = 0.05*nhaver
		270	minhvalid = pdata*nhaver
		271	for im in range(self.dataOut.nmodes):
		272	new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
		273
		274
		275	data_fHeigths_List = []
		276	data_fZonal_List = []
		277	data_fMeridional_List = []
		278	data_fVertical_List = []
		279	startDTList = []
		280
		281
		282	for i in range(ntime):
		283	height = old_height
		284
		285	start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
		286	stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
		287
		288
		289	limit_sec1 = time.mktime(start.timetuple())
		290	limit_sec2 = time.mktime(stop.timetuple())
		291
		292	t1 = numpy.where(self.f_timesec >= limit_sec1)
		293	t2 = numpy.where(self.f_timesec < limit_sec2)
		294	time_select = []
		295	for val_sec in t1[0]:
		296	if val_sec in t2[0]:
		297	time_select.append(val_sec)
		298
		299
		300	time_select = numpy.array(time_select,dtype = 'int')
		301	minvalid = numpy.ceil(pdata*nhpoints)
		302
		303	zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
		304	mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
		305	ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
		306
		307	if nhaver > 1:
		308	new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
		309	new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
		310	new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
		311
		312	if len(time_select) > minvalid:
		313	time_average = self.f_timesec[time_select]
		314
		315	for im in range(self.dataOut.nmodes):
		316
		317	for ih in range(self.dataOut.nranges):
		318	if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
		319	zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
		320
		321	if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
		322	mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
		323
		324	if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
		325	ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
		326
		327	if nhaver > 1:
		328	for ih in range(num_hei):
		329	hvalid = numpy.arange(nhaver) + nhaver*ih
		330
		331	if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
		332	new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) / numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
		333
		334	if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
		335	new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) / numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
		336
		337	if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
		338	new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) / numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
		339	if nhaver > 1:
		340	zon_aver = new_zon_aver
		341	mer_aver = new_mer_aver
		342	ver_aver = new_ver_aver
		343	height = new_height
		344
		345
		346	tstart = time_average[0]
		347	tend = time_average[-1]
		348	startTime = time.gmtime(tstart)
		349
		350	year = startTime.tm_year
		351	month = startTime.tm_mon
		352	day = startTime.tm_mday
		353	hour = startTime.tm_hour
		354	minute = startTime.tm_min
		355	second = startTime.tm_sec
		356
		357	startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
		358
		359
		360	o_height = numpy.array([])
		361	o_zon_aver = numpy.array([])
		362	o_mer_aver = numpy.array([])
		363	o_ver_aver = numpy.array([])
		364	if self.dataOut.nmodes > 1:
		365	for im in range(self.dataOut.nmodes):
		366
		367	if im == 0:
		368	h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
		369	else:
		370	h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
		371
		372
		373	ht = h_select[0]
		374
		375	o_height = numpy.hstack((o_height,height[im,ht]))
		376	o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
		377	o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
		378	o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
		379
		380	data_fHeigths_List.append(o_height)
		381	data_fZonal_List.append(o_zon_aver)
		382	data_fMeridional_List.append(o_mer_aver)
		383	data_fVertical_List.append(o_ver_aver)
		384
		385
		386	else:
		387	h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
		388	ht = h_select[0]
		389	o_height = numpy.hstack((o_height,height[im,ht]))
		390	o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
		391	o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
		392	o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
		393
		394	data_fHeigths_List.append(o_height)
		395	data_fZonal_List.append(o_zon_aver)
		396	data_fMeridional_List.append(o_mer_aver)
		397	data_fVertical_List.append(o_ver_aver)
		398
		399
		400	return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List
		401
		402
		403	No newline at end of file

schainpy/model/data/jrodata.py +9 0

                  noise_estimation = None
+                 GauSPC = None #Fit gaussian SPC
                  def __init__(self):
                      '''
                      else:
                          return self.paramInterval
+                 def setValue(self, value):
+                     print "This property should not be initialized"
+                     return
                  def getNoise(self):
                      return self.spc_noise
                  timeInterval = property(getTimeInterval)
+                 noise = property(getNoise, setValue, "I'm the 'Noise' property.")

schainpy/model/graphics/jroplot_correlation.py +2 -2

                  HEIGHTPROF = None
                  PREFIX = 'corr'
-                 def __init__(self, **kwargs):
-                     Figure.__init__(self, **kwargs)
+                 def __init__(self):
                      self.isConfig = False
                      self.__nsubplots = 1

schainpy/model/graphics/jroplot_heispectra.py +3 -4

                  HEIGHTPROF = None
                  PREFIX = 'spc'
-                 def __init__(self, **kwargs):
+                 def __init__(self):
-                     Figure.__init__(self, **kwargs)
                      self.isConfig = False
                      self.__nsubplots = 1
                  PREFIX = 'rtinoise'
-                 def __init__(self, **kwargs):
-                     Figure.__init__(self, **kwargs)
+                 def __init__(self):
                      self.timerange = 24*60*60
                      self.isConfig = False
                      self.__nsubplots = 1

schainpy/model/graphics/jroplot_parameters.py +218 -9

		@@ -6,6 +6,217 from figure import Figure, isRealtime, isTimeInHourRange
6	6	from plotting_codes import *
7	7
8	8
	9	class FitGauPlot(Figure):
	10
	11	isConfig = None
	12	__nsubplots = None
	13
	14	WIDTHPROF = None
	15	HEIGHTPROF = None
	16	PREFIX = 'fitgau'
	17
	18	def __init__(self, **kwargs):
	19	Figure.__init__(self, **kwargs)
	20	self.isConfig = False
	21	self.__nsubplots = 1
	22
	23	self.WIDTH = 250
	24	self.HEIGHT = 250
	25	self.WIDTHPROF = 120
	26	self.HEIGHTPROF = 0
	27	self.counter_imagwr = 0
	28
	29	self.PLOT_CODE = SPEC_CODE
	30
	31	self.FTP_WEI = None
	32	self.EXP_CODE = None
	33	self.SUB_EXP_CODE = None
	34	self.PLOT_POS = None
	35
	36	self.__xfilter_ena = False
	37	self.__yfilter_ena = False
	38
	39	def getSubplots(self):
	40
	41	ncol = int(numpy.sqrt(self.nplots)+0.9)
	42	nrow = int(self.nplots*1./ncol + 0.9)
	43
	44	return nrow, ncol
	45
	46	def setup(self, id, nplots, wintitle, showprofile=True, show=True):
	47
	48	self.__showprofile = showprofile
	49	self.nplots = nplots
	50
	51	ncolspan = 1
	52	colspan = 1
	53	if showprofile:
	54	ncolspan = 3
	55	colspan = 2
	56	self.__nsubplots = 2
	57
	58	self.createFigure(id = id,
	59	wintitle = wintitle,
	60	widthplot = self.WIDTH + self.WIDTHPROF,
	61	heightplot = self.HEIGHT + self.HEIGHTPROF,
	62	show=show)
	63
	64	nrow, ncol = self.getSubplots()
	65
	66	counter = 0
	67	for y in range(nrow):
	68	for x in range(ncol):
	69
	70	if counter >= self.nplots:
	71	break
	72
	73	self.addAxes(nrow, ncolncolspan, y, xncolspan, colspan, 1)
	74
	75	if showprofile:
	76	self.addAxes(nrow, ncolncolspan, y, xncolspan+colspan, 1, 1)
	77
	78	counter += 1
	79
	80	def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
	81	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
	82	save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
	83	server=None, folder=None, username=None, password=None,
	84	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
	85	xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 1):
	86
	87	"""
	88
	89	Input:
	90	dataOut :
	91	id :
	92	wintitle :
	93	channelList :
	94	showProfile :
	95	xmin : None,
	96	xmax : None,
	97	ymin : None,
	98	ymax : None,
	99	zmin : None,
	100	zmax : None
	101	"""
	102	if realtime:
	103	if not(isRealtime(utcdatatime = dataOut.utctime)):
	104	print 'Skipping this plot function'
	105	return
	106
	107	if channelList == None:
	108	channelIndexList = dataOut.channelIndexList
	109	else:
	110	channelIndexList = []
	111	for channel in channelList:
	112	if channel not in dataOut.channelList:
	113	raise ValueError, "Channel %d is not in dataOut.channelList" %channel
	114	channelIndexList.append(dataOut.channelList.index(channel))
	115
	116	# if normFactor is None:
	117	# factor = dataOut.normFactor
	118	# else:
	119	# factor = normFactor
	120	if xaxis == "frequency":
	121	x = dataOut.spc_range[0]
	122	xlabel = "Frequency (kHz)"
	123
	124	elif xaxis == "time":
	125	x = dataOut.spc_range[1]
	126	xlabel = "Time (ms)"
	127
	128	else:
	129	x = dataOut.spc_range[2]
	130	xlabel = "Velocity (m/s)"
	131
	132	ylabel = "Range (Km)"
	133
	134	y = dataOut.getHeiRange()
	135
	136	z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor
	137	print 'GausSPC', z[0,32,10:40]
	138	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
	139	zdB = 10*numpy.log10(z)
	140
	141	avg = numpy.average(z, axis=1)
	142	avgdB = 10*numpy.log10(avg)
	143
	144	noise = dataOut.spc_noise
	145	noisedB = 10*numpy.log10(noise)
	146
	147	thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
	148	title = wintitle + " Spectra"
	149	if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
	150	title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
	151
	152	if not self.isConfig:
	153
	154	nplots = len(channelIndexList)
	155
	156	self.setup(id=id,
	157	nplots=nplots,
	158	wintitle=wintitle,
	159	showprofile=showprofile,
	160	show=show)
	161
	162	if xmin == None: xmin = numpy.nanmin(x)
	163	if xmax == None: xmax = numpy.nanmax(x)
	164	if ymin == None: ymin = numpy.nanmin(y)
	165	if ymax == None: ymax = numpy.nanmax(y)
	166	if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
	167	if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
	168
	169	self.FTP_WEI = ftp_wei
	170	self.EXP_CODE = exp_code
	171	self.SUB_EXP_CODE = sub_exp_code
	172	self.PLOT_POS = plot_pos
	173
	174	self.isConfig = True
	175
	176	self.setWinTitle(title)
	177
	178	for i in range(self.nplots):
	179	index = channelIndexList[i]
	180	str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
	181	title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
	182	if len(dataOut.beam.codeList) != 0:
	183	title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime)
	184
	185	axes = self.axesList[i*self.__nsubplots]
	186	axes.pcolor(x, y, zdB[index,:,:],
	187	xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
	188	xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap,
	189	ticksize=9, cblabel='')
	190
	191	if self.__showprofile:
	192	axes = self.axesList[i*self.__nsubplots +1]
	193	axes.pline(avgdB[index,:], y,
	194	xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
	195	xlabel='dB', ylabel='', title='',
	196	ytick_visible=False,
	197	grid='x')
	198
	199	noiseline = numpy.repeat(noisedB[index], len(y))
	200	axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
	201
	202	self.draw()
	203
	204	if figfile == None:
	205	str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
	206	name = str_datetime
	207	if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
	208	name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
	209	figfile = self.getFilename(name)
	210
	211	self.save(figpath=figpath,
	212	figfile=figfile,
	213	save=save,
	214	ftp=ftp,
	215	wr_period=wr_period,
	216	thisDatetime=thisDatetime)
	217
	218
	219
9	220	class MomentsPlot(Figure):
10	221
11	222	isConfig = None
		@@ -446,8 +657,7 class WindProfilerPlot(Figure):
446	657	# tmin = None
447	658	# tmax = None
448	659
449
450		x = dataOut.getTimeRange1(dataOut.outputInterval)
	660	x = dataOut.getTimeRange1(dataOut.paramInterval)
451	661	y = dataOut.heightList
452	662	z = dataOut.data_output.copy()
453	663	nplots = z.shape[0] #Number of wind dimensions estimated
		@@ -559,7 +769,7 class WindProfilerPlot(Figure):
559	769	thisDatetime=thisDatetime,
560	770	update_figfile=update_figfile)
561	771
562		if dataOut.ltctime + dataOut.~~output~~Interval >= self.xmax:
	772	if dataOut.ltctime + dataOut.paramInterval >= self.xmax:
563	773	self.counter_imagwr = wr_period
564	774	self.isConfig = False
565	775	update_figfile = True
		@@ -636,12 +846,12 class ParametersPlot(Figure):
636	846
637	847	counter += 1
638	848
639		def run(self, dataOut, id, wintitle="", channelList=None, paramIndex = 0, colormap=~~True~~,
	849	def run(self, dataOut, id, wintitle="", channelList=None, paramIndex = 0, colormap="jet",
640	850	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, timerange=None,
641	851	showSNR=False, SNRthresh = -numpy.inf, SNRmin=None, SNRmax=None,
642	852	save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
643	853	server=None, folder=None, username=None, password=None,
644		ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
	854	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, HEIGHT=None):
645	855	"""
646	856
647	857	Input:
		@@ -658,10 +868,9 class ParametersPlot(Figure):
658	868	zmax : None
659	869	"""
660	870
661		if colormap:
662		colormap="jet"
663		~~else~~:
664		colormap="RdBu_r"
	871	if HEIGHT is not None:
	872	self.HEIGHT = HEIGHT
	873
665	874
666	875	if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
667	876	return

schainpy/model/graphics/jroplot_spectra.py +23 -16

                          save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
                          server=None, folder=None, username=None, password=None,
                          ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
-                         xaxis="velocity", **kwargs):
+                         xaxis="frequency", colormap='jet', normFactor=None):
                      """
                          zmin            :    None,
                          zmax            :    None
                      """
-                     colormap = kwargs.get('colormap','jet')
                      if realtime:
                          if not(isRealtime(utcdatatime = dataOut.utctime)):
                              print 'Skipping this plot function'
                                  raise ValueError, "Channel %d is not in dataOut.channelList" %channel
                              channelIndexList.append(dataOut.channelList.index(channel))
+                     if normFactor is None:
                      factor = dataOut.normFactor
+                     else:
+                         factor = normFactor
                      if xaxis == "frequency":
                          x = dataOut.getFreqRange(1)/1000.
                          xlabel = "Frequency (kHz)"
                          save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
                          power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
                          server=None, folder=None, username=None, password=None,
-                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0,
+                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None,
                          xaxis='frequency'):
                      """
                      if len(pairsIndexList) > 4:
                          pairsIndexList = pairsIndexList[0:4]
+                     if normFactor is None:
                      factor = dataOut.normFactor
+                     else:
+                         factor = normFactor
                      x = dataOut.getVelRange(1)
                      y = dataOut.getHeiRange()
                      z = dataOut.data_spc[:,:,:]/factor
                  def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
                          xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
-                         timerange=None,
+                         timerange=None, colormap='jet',
                          save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
                          server=None, folder=None, username=None, password=None,
-                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, **kwargs):
+                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, HEIGHT=None):
                      """
                          zmax            :    None
                      """
-                     colormap = kwargs.get('colormap', 'jet')
+                     #colormap = kwargs.get('colormap', 'jet')
+                     if HEIGHT is not None:
+                         self.HEIGHT  = HEIGHT
                      if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
                          return
                                  raise ValueError, "Channel %d is not in dataOut.channelList"
                              channelIndexList.append(dataOut.channelList.index(channel))
-                     if hasattr(dataOut, 'normFactor'):
+                     if normFactor is None:
                          factor = dataOut.normFactor
                      else:
-                         factor = 1
+                         factor = normFactor
              #         factor = dataOut.normFactor
                      x = dataOut.getTimeRange()
                      y = dataOut.getHeiRange()
-             #         z = dataOut.data_spc/factor
-             #         z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
-             #         avg = numpy.average(z, axis=1)
-             #         avgdB = 10.*numpy.log10(avg)
-                     avgdB = dataOut.getPower()
+                     z = dataOut.data_spc/factor
+                     z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
+                     avg = numpy.average(z, axis=1)
+                     avgdB = 10.*numpy.log10(avg)
+                     # avgdB = dataOut.getPower()
                      thisDatetime = dataOut.datatime
              #         thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
                  PREFIX = 'noise'
                  def __init__(self, **kwargs):
                      Figure.__init__(self, **kwargs)
                      self.timerange = 24*60*60

schainpy/model/graphics/mpldriver.py +24 -11

                                                      rowspan=rowspan,
                                                      polar=polar)
+                 axes.grid(True)
                  matplotlib.pyplot.ion()
                  return axes
              def pline(iplot, x, y, xlabel='', ylabel='', title=''):
-                 ax = iplot.axes
+                 ax = iplot.get_axes()
                  printLabels(ax, xlabel, ylabel, title)
                  z = numpy.ma.masked_invalid(z)
                  cmap=matplotlib.pyplot.get_cmap(colormap)
-                 cmap.set_bad('black', 1.)
+                 cmap.set_bad('white',1.)
                  imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=cmap)
                  cb =  matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
                  cb.set_label(cblabel)
                      ax.xaxis.set_major_formatter(FuncFormatter(func))
                      ax.xaxis.set_major_locator(LinearLocator(7))
+                 ax.grid(True)
                  matplotlib.pyplot.ion()
                  return imesh
              def pcolor(imesh, z, xlabel='', ylabel='', title=''):
+                 z = numpy.ma.masked_invalid(z)
+                 cmap=matplotlib.pyplot.get_cmap('jet')
+                 cmap.set_bad('white',1.)
                  z = z.T
-                 ax = imesh.axes
+                 ax = imesh.get_axes()
                  printLabels(ax, xlabel, ylabel, title)
                  imesh.set_array(z.ravel())
+                 ax.grid(True)
              def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
                  printLabels(ax, xlabel, ylabel, title)
+                 z = numpy.ma.masked_invalid(z)
+                 cmap=matplotlib.pyplot.get_cmap(colormap)
+                 cmap.set_bad('white',1.)
                  ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
+                 ax.grid(True)
              def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
                  z = numpy.ma.masked_invalid(z)
                  cmap=matplotlib.pyplot.get_cmap(colormap)
-                 cmap.set_bad('black', 1.)
+                 cmap.set_bad('white',1.)
                  ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=cmap)
+                 ax.grid(True)
              def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
                              ticksize=9, xtick_visible=True, ytick_visible=True,
              def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''):
-                 ax = iplot.axes
+                 ax = iplot.get_axes()
                  printLabels(ax, xlabel, ylabel, title)
              def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
-                 ax = iplot.axes
+                 ax = iplot.get_axes()
                  printLabels(ax, xlabel, ylabel, title)
                  for i in range(len(ax.lines)):
              #     ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
              #     ax.text(0, 50, ylabel, rotation='vertical', va ='center', ha = 'left' ,size='11')
              #     ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
-                 ax.yaxis.labelpad = 40
+                 ax.yaxis.labelpad = 230
                  printLabels(ax, xlabel, ylabel, title)
                  iplot = ax.lines[-1]
              def polar(iplot, x, y, xlabel='', ylabel='', title=''):
-                 ax = iplot.axes
+                 ax = iplot.get_axes()
              #     ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
                  printLabels(ax, xlabel, ylabel, title)

schainpy/model/io/__init__.py +7 0

              from jroIO_kamisr import *
              from jroIO_param import *
              from jroIO_hf import *
+             from jroIO_madrigal import *
+             from bltrIO_param import *
+             from jroIO_bltr import *
+             from jroIO_mira35c import *

schainpy/model/io/jroIO_base.py +5 -4

                      print "[Reading] Number of read blocks %04d" %self.nTotalBlocks
                  def printNumberOfBlock(self):
+                     'SPAM!'
-                     if self.flagIsNewBlock:
-                         print "[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
-                                                                   self.processingHeaderObj.dataBlocksPerFile,
-                                                                   self.dataOut.datatime.ctime())
+             #         if self.flagIsNewBlock:
+             #             print "[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
+             #                                                       self.processingHeaderObj.dataBlocksPerFile,
+             #                                                       self.dataOut.datatime.ctime())
                  def printInfo(self):

schainpy/model/io/jroIO_param.py +2 -2

                      print "[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime)
                      print
-                     for i in range(len(filenameList)):
-                         print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
+             #         for i in range(len(filenameList)):
+             #             print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
                      self.filenameList = filenameList
                      self.datetimeList = datetimeList

schainpy/model/io/jroIO_spectra.py +32 -4

              from schainpy.model.data.jrodata import Spectra
              class SpectraReader(JRODataReader, ProcessingUnit):
                  """
                  Esta clase permite leer datos de espectros desde archivos procesados (.pdata). La lectura
                  de los datos siempre se realiza por bloques. Los datos leidos (array de 3 dimensiones)
                                      paresCanalesDiferentes * alturas * perfiles  (Cross Spectra)
                                      canales * alturas                            (DC Channels)
                  Esta clase contiene instancias (objetos) de las clases BasicHeader, SystemHeader,
                  RadarControllerHeader y Spectra. Los tres primeros se usan para almacenar informacion de la
                  cabecera de datos (metadata), y el cuarto (Spectra) para obtener y almacenar un bloque de
                      Inicializador de la clase SpectraReader para la lectura de datos de espectros.
                      Inputs:
                          dataOut    :    Objeto de la clase Spectra. Este objeto sera utilizado para
                                            almacenar un perfil de datos cada vez que se haga un requerimiento
                                            (getData). El perfil sera obtenido a partir del buffer de datos,
                                            bloque de datos.
                                            Si este parametro no es pasado se creara uno internamente.
                      Affected:
                          self.dataOut
                      Return      : None
                      """
                      #Eliminar de la base la herencia
                      ProcessingUnit.__init__(self, **kwargs)
                      for i in range(0, self.processingHeaderObj.totalSpectra*2, 2):
                          if self.processingHeaderObj.spectraComb[i] == self.processingHeaderObj.spectraComb[i+1]:
                              self.nRdChannels = self.nRdChannels + 1 #par de canales iguales
                          else:
                              self.nRdPairs = self.nRdPairs + 1 #par de canales diferentes
                              self.rdPairList.append((self.processingHeaderObj.spectraComb[i], self.processingHeaderObj.spectraComb[i+1]))
                      Variables afectadas:
                          self.flagIsNewFile
                          self.flagIsNewBlock
                          self.nTotalBlocks
                      Exceptions:
                          Si un bloque leido no es un bloque valido
                      """
+                     print ' ======================================================== '
+                     print ' '
+                     print ' '
+                     print self.processingHeaderObj.totalSpectra, 'TotalSpectra', type(self.processingHeaderObj.totalSpectra)
+                     print self.processingHeaderObj.spectraComb, 'SpectraComb', type(self.processingHeaderObj.spectraComb)
+                     print ' '
+                     print ' '
+                     print ' ======================================================== '
                      blockOk_flag = False
                      fpointer = self.fp.tell()
                      self.data_spc = spc
                      if self.processingHeaderObj.flag_cspc:
                          cspc = numpy.transpose( cspc, (0,2,1) )
                          self.data_cspc = cspc['real'] + cspc['imag']*1j
                      else:
                          self.data_cspc = None
                      if self.processingHeaderObj.flag_dc:
                          self.data_dc = dc['real'] + dc['imag']*1j
                      else:
                              self.dataOut.flagNoData = True
                              return 0
                      #data es un numpy array de 3 dmensiones (perfiles, alturas y canales)
                      if self.data_spc is None:
              #    dataOut = None
-                 def __init__(self, **kwargs):
+                 def __init__(self):
                      """
                      Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
                      Affected:
                          self.dataOut
                          self.basicHeaderObj
                          self.systemHeaderObj
                      Return: None
                      """
-                     Operation.__init__(self, **kwargs)
+                     Operation.__init__(self)
                      self.isConfig = False
                      self.data_cspc = None
                      self.data_dc = None
                      self.fp = None
                      return 1
                  def setBlockDimension(self):
                      """
                      Obtiene las formas dimensionales del los subbloques de datos que componen un bloque
                      """
                      Escribe el buffer en el file designado
                      Affected:
                          self.data_spc
                          self.data_cspc
                          data = data.reshape((-1))
                          data.tofile(self.fp)
                      if self.data_dc is not None:
                          data = numpy.zeros( self.shape_dc_Buffer, self.dtype )
                          dc = self.data_dc
              #         if self.data_cspc is not None:
              #             self.data_cspc.fill(0)
                      self.flagIsNewFile = 0
                      self.flagIsNewBlock = 1
                      self.nTotalBlocks += 1
                      """
                      Setea un bloque de datos y luego los escribe en un file
                      Affected:
                          self.data_spc
                          self.data_cspc
                      if self.dataOut.flagDiscontinuousBlock:
                          self.data_spc.fill(0)
-                         if self.dataOut.data_cspc is not None:
-                             self.data_cspc.fill(0)
-                         if self.dataOut.data_dc is not None:
-                             self.data_dc.fill(0)
                          self.setNextFile()
                      return 1
                  def __getBlockSize(self):
                      '''
                      Este metodos determina el cantidad de bytes para un bloque de datos de tipo Spectra
                      self.processingHeaderObj.totalSpectra = self.dataOut.nPairs + self.dataOut.nChannels
                      self.processingHeaderObj.shif_fft = self.dataOut.flagShiftFFT
                      if self.processingHeaderObj.totalSpectra > 0:
                          channelList = []
                          for channel in range(self.dataOut.nChannels):

schainpy/model/io/jroIO_voltage.py +6 0

                      self.getBasicHeader()
+                     #print self.basicHeaderObj.printInfo()
+                     #print self.systemHeaderObj.printInfo()
+                     #print self.radarControllerHeaderObj.printInfo()
+                     #print self.processingHeaderObj.printInfo()
                      self.dataOut.realtime = self.online
                      return self.dataOut.data
              class VoltageWriter(JRODataWriter, Operation):
                  """
                  Esta clase permite escribir datos de voltajes a archivos procesados (.r). La escritura

schainpy/model/proc/__init__.py +1 0

              from jroproc_parameters import *
              from jroproc_spectra_lags import *
              from jroproc_spectra_acf import *
+             from bltrproc_parameters import *

schainpy/model/proc/jroproc_parameters.py +1305 -60

This diff has been collapsed as it changes many lines, (1365 lines changed) Show them Hide them
			@@ -1,15 +1,62
	1	1	import numpy
	2	2	import math
	3	3	from scipy import optimize, interpolate, signal, stats, ndimage
		4	import scipy
	4	5	import re
	5	6	import datetime
	6	7	import copy
	7	8	import sys
	8	9	import importlib
	9	10	import itertools
	10
		11	from multiprocessing import Pool, TimeoutError
		12	from multiprocessing.pool import ThreadPool
		13	import copy_reg
		14	import cPickle
		15	import types
		16	from functools import partial
		17	import time
		18	#from sklearn.cluster import KMeans
		19
		20	import matplotlib.pyplot as plt
		21
		22	from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
	11	23	from jroproc_base import ProcessingUnit, Operation
	12	24	from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
		25	from scipy import asarray as ar,exp
		26	from scipy.optimize import curve_fit
		27
		28	import warnings
		29	from numpy import NaN
		30	from scipy.optimize.optimize import OptimizeWarning
		31	warnings.filterwarnings('ignore')
		32
		33
		34	SPEED_OF_LIGHT = 299792458
		35
		36
		37	'''solving pickling issue'''
		38
		39	def _pickle_method(method):
		40	func_name = method.im_func.__name__
		41	obj = method.im_self
		42	cls = method.im_class
		43	return _unpickle_method, (func_name, obj, cls)
		44
		45	def _unpickle_method(func_name, obj, cls):
		46	for cls in cls.mro():
		47	try:
		48	func = cls.__dict__[func_name]
		49	except KeyError:
		50	pass
		51	else:
		52	break
		53	return func.__get__(obj, cls)
		54
		55
		56
		57
		58
		59
	13	60
	14	61
	15	62	class ParametersProc(ProcessingUnit):
			@@ -57,7 +104,9 class ParametersProc(ProcessingUnit):
	57	104	self.dataOut.timeInterval1 = self.dataIn.timeInterval
	58	105	self.dataOut.heightList = self.dataIn.getHeiRange()
	59	106	self.dataOut.frequency = self.dataIn.frequency
	60		#self.dataOut.noise = self.dataIn.noise
		107	self.dataOut.noise = self.dataIn.noise
		108
		109
	61	110
	62	111	def run(self):
	63	112
			@@ -77,17 +126,36 class ParametersProc(ProcessingUnit):
	77	126	if self.dataIn.type == "Spectra":
	78	127
	79	128	self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)
	80		~~self~~.~~dataOut~~.~~data_spc~~ = self.dataIn.data_spc
	81		self.dataOut.data_cspc = self.dataIn.data_cspc
	82		self.dataOut.nProfiles = self.dataIn.nProfiles
	83		self.dataOut.nIncohInt = self.dataIn.nIncohInt
	84		self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
	85		self.dataOut.ippFactor = self.dataIn.ippFactor
	86		#self.dataOut.normFactor = self.dataIn.getNormFactor()
	87		self.dataOut.pairsList = self.dataIn.pairsList
	88		self.dataOut.groupList = self.dataIn.pairsList
		129	print 'self.dataIn.data_spc', self.dataIn.data_spc.shape
	89	130	self.dataOut.abscissaList = self.dataIn.getVelRange(1)
		131	self.dataOut.spc_noise = self.dataIn.getNoise()
		132	self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1) )
		133
		134	self.dataOut.normFactor = self.dataIn.normFactor
		135	#self.dataOut.outputInterval = self.dataIn.outputInterval
		136	self.dataOut.groupList = self.dataIn.pairsList
	90	137	self.dataOut.flagNoData = False
		138	#print 'datain chandist ',self.dataIn.ChanDist
		139	if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
		140	self.dataOut.ChanDist = self.dataIn.ChanDist
		141	else: self.dataOut.ChanDist = None
		142
		143	print 'datain chandist ',self.dataOut.ChanDist
		144
		145	if hasattr(self.dataIn, 'VelRange'): #Velocities range
		146	self.dataOut.VelRange = self.dataIn.VelRange
		147	else: self.dataOut.VelRange = None
		148
		149	if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
		150	self.dataOut.RadarConst = self.dataIn.RadarConst
		151
		152	if hasattr(self.dataIn, 'NPW'): #NPW
		153	self.dataOut.NPW = self.dataIn.NPW
		154
		155	if hasattr(self.dataIn, 'COFA'): #COFA
		156	self.dataOut.COFA = self.dataIn.COFA
		157
		158
	91	159
	92	160	#---------------------- Correlation Data ---------------------------
	93	161
			@@ -108,7 +176,6 class ParametersProc(ProcessingUnit):
	108	176
	109	177	if self.dataIn.type == "Parameters":
	110	178	self.dataOut.copy(self.dataIn)
	111		self.dataOut.utctimeInit = self.dataIn.utctime
	112	179	self.dataOut.flagNoData = False
	113	180
	114	181	return True
			@@ -119,6 +186,1186 class ParametersProc(ProcessingUnit):
	119	186
	120	187	return
	121	188
		189
		190	def target(tups):
		191
		192	obj, args = tups
		193	#print 'TARGETTT', obj, args
		194	return obj.FitGau(args)
		195
		196	class GaussianFit(Operation):
		197
		198	'''
		199	Function that fit of one and two generalized gaussians (gg) based
		200	on the PSD shape across an "power band" identified from a cumsum of
		201	the measured spectrum - noise.
		202
		203	Input:
		204	self.dataOut.data_pre : SelfSpectra
		205
		206	Output:
		207	self.dataOut.GauSPC : SPC_ch1, SPC_ch2
		208
		209	'''
		210	def __init__(self, **kwargs):
		211	Operation.__init__(self, **kwargs)
		212	self.i=0
		213
		214
		215	def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
		216	"""This routine will find a couple of generalized Gaussians to a power spectrum
		217	input: spc
		218	output:
		219	Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
		220	"""
		221
		222	self.spc = dataOut.data_pre[0].copy()
		223
		224
		225	print 'SelfSpectra Shape', numpy.asarray(self.spc).shape
		226
		227
		228	#plt.figure(50)
		229	#plt.subplot(121)
		230	#plt.plot(self.spc,'k',label='spc(66)')
		231	#plt.plot(xFrec,ySamples[1],'g',label='Ch1')
		232	#plt.plot(xFrec,ySamples[2],'r',label='Ch2')
		233	#plt.plot(xFrec,FitGauss,'yo:',label='fit')
		234	#plt.legend()
		235	#plt.title('DATOS A ALTURA DE 7500 METROS')
		236	#plt.show()
		237
		238	self.Num_Hei = self.spc.shape[2]
		239	#self.Num_Bin = len(self.spc)
		240	self.Num_Bin = self.spc.shape[1]
		241	self.Num_Chn = self.spc.shape[0]
		242
		243	Vrange = dataOut.abscissaList
		244
		245	#print 'self.spc2', numpy.asarray(self.spc).shape
		246
		247	GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei])
		248	SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
		249	SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
		250	SPC_ch1[:] = numpy.NaN
		251	SPC_ch2[:] = numpy.NaN
		252
		253
		254	start_time = time.time()
		255
		256	noise_ = dataOut.spc_noise[0].copy()
		257
		258
		259
		260	pool = Pool(processes=self.Num_Chn)
		261	args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
		262	objs = [self for __ in range(self.Num_Chn)]
		263	attrs = zip(objs, args)
		264	gauSPC = pool.map(target, attrs)
		265	dataOut.GauSPC = numpy.asarray(gauSPC)
		266	# ret = []
		267	# for n in range(self.Num_Chn):
		268	# self.FitGau(args[n])
		269	# dataOut.GauSPC = ret
		270
		271
		272
		273	# for ch in range(self.Num_Chn):
		274	#
		275	# for ht in range(self.Num_Hei):
		276	# #print (numpy.asarray(self.spc).shape)
		277	# spc = numpy.asarray(self.spc)[ch,:,ht]
		278	#
		279	# #############################################
		280	# # normalizing spc and noise
		281	# # This part differs from gg1
		282	# spc_norm_max = max(spc)
		283	# spc = spc / spc_norm_max
		284	# pnoise = pnoise / spc_norm_max
		285	# #############################################
		286	#
		287	# if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's
		288	# fatspectra=1.0
		289	# else:
		290	# fatspectra=0.5
		291	#
		292	# wnoise = noise_ / spc_norm_max
		293	# #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise
		294	# #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
		295	# #if wnoise>1.1*pnoise: # to be tested later
		296	# # wnoise=pnoise
		297	# noisebl=wnoise0.9; noisebh=wnoise1.1
		298	# spc=spc-wnoise
		299	#
		300	# minx=numpy.argmin(spc)
		301	# spcs=numpy.roll(spc,-minx)
		302	# cum=numpy.cumsum(spcs)
		303	# tot_noise=wnoise * self.Num_Bin #64;
		304	# #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
		305	# #snr=tot_signal/tot_noise
		306	# #snr=cum[-1]/tot_noise
		307	#
		308	# #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
		309	#
		310	# snr = sum(spcs)/tot_noise
		311	# snrdB=10.*numpy.log10(snr)
		312	#
		313	# #if snrdB < -9 :
		314	# # snrdB = numpy.NaN
		315	# # continue
		316	#
		317	# #print 'snr',snrdB # , sum(spcs) , tot_noise
		318	#
		319	#
		320	# #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
		321	# # return [None,]4,[None,]4,None,snrdB,None,None,[None,]5,[None,]9,None
		322	#
		323	# cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
		324	# cumlo=cummax*epsi;
		325	# cumhi=cummax*(1-epsi)
		326	# powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
		327	#
		328	# #if len(powerindex)==1:
		329	# ##return [numpy.mod(powerindex[0]+minx,64),None,None,None,],[None,]4,None,snrdB,None,None,[None,]5,[None,]*9,None
		330	# #return [numpy.mod(powerindex[0]+minx, self.Num_Bin ),None,None,None,],[None,]4,None,snrdB,None,None,[None,]5,[None,]*9,None
		331	# #elif len(powerindex)<4*fatspectra:
		332	# #return [None,]4,[None,]4,None,snrdB,None,None,[None,]5,[None,]9,None
		333	#
		334	# if len(powerindex) < 1:# case for powerindex 0
		335	# continue
		336	# powerlo=powerindex[0]
		337	# powerhi=powerindex[-1]
		338	# powerwidth=powerhi-powerlo
		339	#
		340	# firstpeak=powerlo+powerwidth/10.# first gaussian energy location
		341	# secondpeak=powerhi-powerwidth/10.#second gaussian energy location
		342	# midpeak=(firstpeak+secondpeak)/2.
		343	# firstamp=spcs[int(firstpeak)]
		344	# secondamp=spcs[int(secondpeak)]
		345	# midamp=spcs[int(midpeak)]
		346	# #x=numpy.spc.shape[1]
		347	#
		348	# #x=numpy.arange(64)
		349	# x=numpy.arange( self.Num_Bin )
		350	# y_data=spc+wnoise
		351	#
		352	# # single gaussian
		353	# #shift0=numpy.mod(midpeak+minx,64)
		354	# shift0=numpy.mod(midpeak+minx, self.Num_Bin )
		355	# width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
		356	# power0=2.
		357	# amplitude0=midamp
		358	# state0=[shift0,width0,amplitude0,power0,wnoise]
		359	# #bnds=((0,63),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
		360	# bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
		361	# #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(0.1,0.5))
		362	# # bnds = range of fft, power width, amplitude, power, noise
		363	# lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
		364	#
		365	# chiSq1=lsq1[1];
		366	# jack1= self.y_jacobian1(x,lsq1[0])
		367	#
		368	#
		369	# try:
		370	# sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
		371	# except:
		372	# std1=32.; sigmas1=numpy.ones(5)
		373	# else:
		374	# std1=sigmas1[0]
		375	#
		376	#
		377	# if fatspectra<1.0 and powerwidth<4:
		378	# choice=0
		379	# Amplitude0=lsq1[0][2]
		380	# shift0=lsq1[0][0]
		381	# width0=lsq1[0][1]
		382	# p0=lsq1[0][3]
		383	# Amplitude1=0.
		384	# shift1=0.
		385	# width1=0.
		386	# p1=0.
		387	# noise=lsq1[0][4]
		388	# #return (numpy.array([shift0,width0,Amplitude0,p0]),
		389	# # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
		390	#
		391	# # two gaussians
		392	# #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
		393	# shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
		394	# shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
		395	# width0=powerwidth/6.;
		396	# width1=width0
		397	# power0=2.;
		398	# power1=power0
		399	# amplitude0=firstamp;
		400	# amplitude1=secondamp
		401	# state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
		402	# #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
		403	# bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
		404	# #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
		405	#
		406	# lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
		407	#
		408	#
		409	# chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
		410	#
		411	#
		412	# try:
		413	# sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
		414	# except:
		415	# std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
		416	# else:
		417	# std2a=sigmas2[0]; std2b=sigmas2[4]
		418	#
		419	#
		420	#
		421	# oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
		422	#
		423	# if snrdB>-9: # when SNR is strong pick the peak with least shift (LOS velocity) error
		424	# if oneG:
		425	# choice=0
		426	# else:
		427	# w1=lsq2[0][1]; w2=lsq2[0][5]
		428	# a1=lsq2[0][2]; a2=lsq2[0][6]
		429	# p1=lsq2[0][3]; p2=lsq2[0][7]
		430	# s1=(2*(1+1./p1))scipy.special.gamma(1./p1)/p1; s2=(2*(1+1./p2))scipy.special.gamma(1./p2)/p2;
		431	# gp1=a1w1s1; gp2=a2w2s2 # power content of each ggaussian with proper p scaling
		432	#
		433	# if gp1>gp2:
		434	# if a1>0.7*a2:
		435	# choice=1
		436	# else:
		437	# choice=2
		438	# elif gp2>gp1:
		439	# if a2>0.7*a1:
		440	# choice=2
		441	# else:
		442	# choice=1
		443	# else:
		444	# choice=numpy.argmax([a1,a2])+1
		445	# #else:
		446	# #choice=argmin([std2a,std2b])+1
		447	#
		448	# else: # with low SNR go to the most energetic peak
		449	# choice=numpy.argmax([lsq1[0][2]lsq1[0][1],lsq2[0][2]lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
		450	#
		451	# #print 'choice',choice
		452	#
		453	# if choice==0: # pick the single gaussian fit
		454	# Amplitude0=lsq1[0][2]
		455	# shift0=lsq1[0][0]
		456	# width0=lsq1[0][1]
		457	# p0=lsq1[0][3]
		458	# Amplitude1=0.
		459	# shift1=0.
		460	# width1=0.
		461	# p1=0.
		462	# noise=lsq1[0][4]
		463	# elif choice==1: # take the first one of the 2 gaussians fitted
		464	# Amplitude0 = lsq2[0][2]
		465	# shift0 = lsq2[0][0]
		466	# width0 = lsq2[0][1]
		467	# p0 = lsq2[0][3]
		468	# Amplitude1 = lsq2[0][6] # This is 0 in gg1
		469	# shift1 = lsq2[0][4] # This is 0 in gg1
		470	# width1 = lsq2[0][5] # This is 0 in gg1
		471	# p1 = lsq2[0][7] # This is 0 in gg1
		472	# noise = lsq2[0][8]
		473	# else: # the second one
		474	# Amplitude0 = lsq2[0][6]
		475	# shift0 = lsq2[0][4]
		476	# width0 = lsq2[0][5]
		477	# p0 = lsq2[0][7]
		478	# Amplitude1 = lsq2[0][2] # This is 0 in gg1
		479	# shift1 = lsq2[0][0] # This is 0 in gg1
		480	# width1 = lsq2[0][1] # This is 0 in gg1
		481	# p1 = lsq2[0][3] # This is 0 in gg1
		482	# noise = lsq2[0][8]
		483	#
		484	# #print len(noise + Amplitude0numpy.exp(-0.5(abs(x-shift0))/width0)**p0)
		485	# SPC_ch1[:,ht] = noise + Amplitude0numpy.exp(-0.5(abs(x-shift0))/width0)**p0
		486	# SPC_ch2[:,ht] = noise + Amplitude1numpy.exp(-0.5(abs(x-shift1))/width1)**p1
		487	# #print 'SPC_ch1.shape',SPC_ch1.shape
		488	# #print 'SPC_ch2.shape',SPC_ch2.shape
		489	# #dataOut.data_param = SPC_ch1
		490	# GauSPC[0] = SPC_ch1
		491	# GauSPC[1] = SPC_ch2
		492
		493	# #plt.gcf().clear()
		494	# plt.figure(50+self.i)
		495	# self.i=self.i+1
		496	# #plt.subplot(121)
		497	# plt.plot(self.spc,'k')#,label='spc(66)')
		498	# plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1')
		499	# #plt.plot(SPC_ch2,'r')#,label='gg2')
		500	# #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
		501	# #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
		502	# #plt.plot(xFrec,FitGauss,'yo:',label='fit')
		503	# plt.legend()
		504	# plt.title('DATOS A ALTURA DE 7500 METROS')
		505	# plt.show()
		506	# print 'shift0', shift0
		507	# print 'Amplitude0', Amplitude0
		508	# print 'width0', width0
		509	# print 'p0', p0
		510	# print '========================'
		511	# print 'shift1', shift1
		512	# print 'Amplitude1', Amplitude1
		513	# print 'width1', width1
		514	# print 'p1', p1
		515	# print 'noise', noise
		516	# print 's_noise', wnoise
		517
		518	print '========================================================'
		519	print 'total_time: ', time.time()-start_time
		520
		521	# re-normalizing spc and noise
		522	# This part differs from gg1
		523
		524
		525
		526	''' Parameters:
		527	1. Amplitude
		528	2. Shift
		529	3. Width
		530	4. Power
		531	'''
		532
		533
		534	###############################################################################
		535	def FitGau(self, X):
		536
		537	Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
		538	#print 'VARSSSS', ch, pnoise, noise, num_intg
		539
		540	#print 'HEIGHTS', self.Num_Hei
		541
		542	GauSPC = []
		543	SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
		544	SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
		545	SPC_ch1[:] = 0#numpy.NaN
		546	SPC_ch2[:] = 0#numpy.NaN
		547
		548
		549
		550	for ht in range(self.Num_Hei):
		551	#print (numpy.asarray(self.spc).shape)
		552
		553	#print 'TTTTT', ch , ht
		554	#print self.spc.shape
		555
		556
		557	spc = numpy.asarray(self.spc)[ch,:,ht]
		558
		559	#############################################
		560	# normalizing spc and noise
		561	# This part differs from gg1
		562	spc_norm_max = max(spc)
		563	spc = spc / spc_norm_max
		564	pnoise = pnoise / spc_norm_max
		565	#############################################
		566
		567	fatspectra=1.0
		568
		569	wnoise = noise_ / spc_norm_max
		570	#wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
		571	#if wnoise>1.1*pnoise: # to be tested later
		572	# wnoise=pnoise
		573	noisebl=wnoise0.9; noisebh=wnoise1.1
		574	spc=spc-wnoise
		575	# print 'wnoise', noise_[0], spc_norm_max, wnoise
		576	minx=numpy.argmin(spc)
		577	spcs=numpy.roll(spc,-minx)
		578	cum=numpy.cumsum(spcs)
		579	tot_noise=wnoise * self.Num_Bin #64;
		580	#print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
		581	#tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
		582	#snr=tot_signal/tot_noise
		583	#snr=cum[-1]/tot_noise
		584	snr = sum(spcs)/tot_noise
		585	snrdB=10.*numpy.log10(snr)
		586
		587	if snrdB < SNRlimit :
		588	snr = numpy.NaN
		589	SPC_ch1[:,ht] = 0#numpy.NaN
		590	SPC_ch1[:,ht] = 0#numpy.NaN
		591	GauSPC = (SPC_ch1,SPC_ch2)
		592	continue
		593	#print 'snr',snrdB #, sum(spcs) , tot_noise
		594
		595
		596
		597	#if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
		598	# return [None,]4,[None,]4,None,snrdB,None,None,[None,]5,[None,]9,None
		599
		600	cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
		601	cumlo=cummax*epsi;
		602	cumhi=cummax*(1-epsi)
		603	powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
		604
		605
		606	if len(powerindex) < 1:# case for powerindex 0
		607	continue
		608	powerlo=powerindex[0]
		609	powerhi=powerindex[-1]
		610	powerwidth=powerhi-powerlo
		611
		612	firstpeak=powerlo+powerwidth/10.# first gaussian energy location
		613	secondpeak=powerhi-powerwidth/10.#second gaussian energy location
		614	midpeak=(firstpeak+secondpeak)/2.
		615	firstamp=spcs[int(firstpeak)]
		616	secondamp=spcs[int(secondpeak)]
		617	midamp=spcs[int(midpeak)]
		618
		619	x=numpy.arange( self.Num_Bin )
		620	y_data=spc+wnoise
		621
		622	# single gaussian
		623	shift0=numpy.mod(midpeak+minx, self.Num_Bin )
		624	width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
		625	power0=2.
		626	amplitude0=midamp
		627	state0=[shift0,width0,amplitude0,power0,wnoise]
		628	bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
		629	lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
		630
		631	chiSq1=lsq1[1];
		632	jack1= self.y_jacobian1(x,lsq1[0])
		633
		634
		635	try:
		636	sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
		637	except:
		638	std1=32.; sigmas1=numpy.ones(5)
		639	else:
		640	std1=sigmas1[0]
		641
		642
		643	if fatspectra<1.0 and powerwidth<4:
		644	choice=0
		645	Amplitude0=lsq1[0][2]
		646	shift0=lsq1[0][0]
		647	width0=lsq1[0][1]
		648	p0=lsq1[0][3]
		649	Amplitude1=0.
		650	shift1=0.
		651	width1=0.
		652	p1=0.
		653	noise=lsq1[0][4]
		654	#return (numpy.array([shift0,width0,Amplitude0,p0]),
		655	# numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
		656
		657	# two gaussians
		658	#shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
		659	shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
		660	shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
		661	width0=powerwidth/6.;
		662	width1=width0
		663	power0=2.;
		664	power1=power0
		665	amplitude0=firstamp;
		666	amplitude1=secondamp
		667	state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
		668	#bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
		669	bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
		670	#bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
		671
		672	lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
		673
		674
		675	chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
		676
		677
		678	try:
		679	sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
		680	except:
		681	std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
		682	else:
		683	std2a=sigmas2[0]; std2b=sigmas2[4]
		684
		685
		686
		687	oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
		688
		689	if snrdB>-6: # when SNR is strong pick the peak with least shift (LOS velocity) error
		690	if oneG:
		691	choice=0
		692	else:
		693	w1=lsq2[0][1]; w2=lsq2[0][5]
		694	a1=lsq2[0][2]; a2=lsq2[0][6]
		695	p1=lsq2[0][3]; p2=lsq2[0][7]
		696	s1=(2*(1+1./p1))scipy.special.gamma(1./p1)/p1;
		697	s2=(2*(1+1./p2))scipy.special.gamma(1./p2)/p2;
		698	gp1=a1w1s1; gp2=a2w2s2 # power content of each ggaussian with proper p scaling
		699
		700	if gp1>gp2:
		701	if a1>0.7*a2:
		702	choice=1
		703	else:
		704	choice=2
		705	elif gp2>gp1:
		706	if a2>0.7*a1:
		707	choice=2
		708	else:
		709	choice=1
		710	else:
		711	choice=numpy.argmax([a1,a2])+1
		712	#else:
		713	#choice=argmin([std2a,std2b])+1
		714
		715	else: # with low SNR go to the most energetic peak
		716	choice=numpy.argmax([lsq1[0][2]lsq1[0][1],lsq2[0][2]lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
		717
		718
		719	shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
		720	shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
		721
		722	max_vel = 20
		723
		724	#first peak will be 0, second peak will be 1
		725	if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range
		726	shift0=lsq2[0][0]
		727	width0=lsq2[0][1]
		728	Amplitude0=lsq2[0][2]
		729	p0=lsq2[0][3]
		730
		731	shift1=lsq2[0][4]
		732	width1=lsq2[0][5]
		733	Amplitude1=lsq2[0][6]
		734	p1=lsq2[0][7]
		735	noise=lsq2[0][8]
		736	else:
		737	shift1=lsq2[0][0]
		738	width1=lsq2[0][1]
		739	Amplitude1=lsq2[0][2]
		740	p1=lsq2[0][3]
		741
		742	shift0=lsq2[0][4]
		743	width0=lsq2[0][5]
		744	Amplitude0=lsq2[0][6]
		745	p0=lsq2[0][7]
		746	noise=lsq2[0][8]
		747
		748	if Amplitude0<0.1: # in case the peak is noise
		749	shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
		750	if Amplitude1<0.1:
		751	shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
		752
		753
		754	# if choice==0: # pick the single gaussian fit
		755	# Amplitude0=lsq1[0][2]
		756	# shift0=lsq1[0][0]
		757	# width0=lsq1[0][1]
		758	# p0=lsq1[0][3]
		759	# Amplitude1=0.
		760	# shift1=0.
		761	# width1=0.
		762	# p1=0.
		763	# noise=lsq1[0][4]
		764	# elif choice==1: # take the first one of the 2 gaussians fitted
		765	# Amplitude0 = lsq2[0][2]
		766	# shift0 = lsq2[0][0]
		767	# width0 = lsq2[0][1]
		768	# p0 = lsq2[0][3]
		769	# Amplitude1 = lsq2[0][6] # This is 0 in gg1
		770	# shift1 = lsq2[0][4] # This is 0 in gg1
		771	# width1 = lsq2[0][5] # This is 0 in gg1
		772	# p1 = lsq2[0][7] # This is 0 in gg1
		773	# noise = lsq2[0][8]
		774	# else: # the second one
		775	# Amplitude0 = lsq2[0][6]
		776	# shift0 = lsq2[0][4]
		777	# width0 = lsq2[0][5]
		778	# p0 = lsq2[0][7]
		779	# Amplitude1 = lsq2[0][2] # This is 0 in gg1
		780	# shift1 = lsq2[0][0] # This is 0 in gg1
		781	# width1 = lsq2[0][1] # This is 0 in gg1
		782	# p1 = lsq2[0][3] # This is 0 in gg1
		783	# noise = lsq2[0][8]
		784
		785	#print len(noise + Amplitude0numpy.exp(-0.5(abs(x-shift0))/width0)**p0)
		786	SPC_ch1[:,ht] = noise + Amplitude0numpy.exp(-0.5(abs(x-shift0))/width0)**p0
		787	SPC_ch2[:,ht] = noise + Amplitude1numpy.exp(-0.5(abs(x-shift1))/width1)**p1
		788	#print 'SPC_ch1.shape',SPC_ch1.shape
		789	#print 'SPC_ch2.shape',SPC_ch2.shape
		790	#dataOut.data_param = SPC_ch1
		791	GauSPC = (SPC_ch1,SPC_ch2)
		792	#GauSPC[1] = SPC_ch2
		793
		794	# print 'shift0', shift0
		795	# print 'Amplitude0', Amplitude0
		796	# print 'width0', width0
		797	# print 'p0', p0
		798	# print '========================'
		799	# print 'shift1', shift1
		800	# print 'Amplitude1', Amplitude1
		801	# print 'width1', width1
		802	# print 'p1', p1
		803	# print 'noise', noise
		804	# print 's_noise', wnoise
		805
		806	return GauSPC
		807
		808
		809	def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination.
		810	y_model=self.y_model1(x,state)
		811	s0,w0,a0,p0,n=state
		812	e0=((x-s0)/w0)**2;
		813
		814	e0u=((x-s0-self.Num_Bin)/w0)**2;
		815
		816	e0d=((x-s0+self.Num_Bin)/w0)**2
		817	m0=numpy.exp(-0.5e0*(p0/2.));
		818	m0u=numpy.exp(-0.5e0u*(p0/2.));
		819	m0d=numpy.exp(-0.5e0d*(p0/2.))
		820	JA=m0+m0u+m0d
		821	JP=(-1/4.)a0m0e0(p0/2.)numpy.log(e0)+(-1/4.)a0m0ue0u(p0/2.)numpy.log(e0u)+(-1/4.)a0m0de0d(p0/2.)numpy.log(e0d)
		822
		823	JS=(p0/w0/2.)a0m0e0(p0/2.-1)((x-s0)/w0)+(p0/w0/2.)a0m0ue0u(p0/2.-1)((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)a0m0de0d(p0/2.-1)((x-s0+ self.Num_Bin )/w0)
		824
		825	JW=(p0/w0/2.)a0m0e0(p0/2.-1)((x-s0)/w0)*2+(p0/w0/2.)a0m0ue0u*(p0/2.-1)((x-s0- self.Num_Bin )/w0)*2+(p0/w0/2.)a0m0de0d*(p0/2.-1)((x-s0+ self.Num_Bin )/w0)**2
		826	jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model])
		827	return jack1.T
		828
		829	def y_jacobian2(self,x,state):
		830	y_model=self.y_model2(x,state)
		831	s0,w0,a0,p0,s1,w1,a1,p1,n=state
		832	e0=((x-s0)/w0)**2;
		833
		834	e0u=((x-s0- self.Num_Bin )/w0)**2;
		835
		836	e0d=((x-s0+ self.Num_Bin )/w0)**2
		837	e1=((x-s1)/w1)**2;
		838
		839	e1u=((x-s1- self.Num_Bin )/w1)**2;
		840
		841	e1d=((x-s1+ self.Num_Bin )/w1)**2
		842	m0=numpy.exp(-0.5e0*(p0/2.));
		843	m0u=numpy.exp(-0.5e0u*(p0/2.));
		844	m0d=numpy.exp(-0.5e0d*(p0/2.))
		845	m1=numpy.exp(-0.5e1*(p1/2.));
		846	m1u=numpy.exp(-0.5e1u*(p1/2.));
		847	m1d=numpy.exp(-0.5e1d*(p1/2.))
		848	JA=m0+m0u+m0d
		849	JA1=m1+m1u+m1d
		850	JP=(-1/4.)a0m0e0(p0/2.)numpy.log(e0)+(-1/4.)a0m0ue0u(p0/2.)numpy.log(e0u)+(-1/4.)a0m0de0d(p0/2.)numpy.log(e0d)
		851	JP1=(-1/4.)a1m1e1(p1/2.)numpy.log(e1)+(-1/4.)a1m1ue1u(p1/2.)numpy.log(e1u)+(-1/4.)a1m1de1d(p1/2.)numpy.log(e1d)
		852
		853	JS=(p0/w0/2.)a0m0e0(p0/2.-1)((x-s0)/w0)+(p0/w0/2.)a0m0ue0u(p0/2.-1)((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)a0m0de0d(p0/2.-1)((x-s0+ self.Num_Bin )/w0)
		854
		855	JS1=(p1/w1/2.)a1m1e1(p1/2.-1)((x-s1)/w1)+(p1/w1/2.)a1m1ue1u(p1/2.-1)((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)a1m1de1d(p1/2.-1)((x-s1+ self.Num_Bin )/w1)
		856
		857	JW=(p0/w0/2.)a0m0e0(p0/2.-1)((x-s0)/w0)*2+(p0/w0/2.)a0m0ue0u*(p0/2.-1)((x-s0- self.Num_Bin )/w0)*2+(p0/w0/2.)a0m0de0d*(p0/2.-1)((x-s0+ self.Num_Bin )/w0)**2
		858
		859	JW1=(p1/w1/2.)a1m1e1(p1/2.-1)((x-s1)/w1)*2+(p1/w1/2.)a1m1ue1u*(p1/2.-1)((x-s1- self.Num_Bin )/w1)*2+(p1/w1/2.)a1m1de1d*(p1/2.-1)((x-s1+ self.Num_Bin )/w1)**2
		860	jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model])
		861	return jack2.T
		862
		863	def y_model1(self,x,state):
		864	shift0,width0,amplitude0,power0,noise=state
		865	model0=amplitude0numpy.exp(-0.5abs((x-shift0)/width0)**power0)
		866
		867	model0u=amplitude0numpy.exp(-0.5abs((x-shift0- self.Num_Bin )/width0)**power0)
		868
		869	model0d=amplitude0numpy.exp(-0.5abs((x-shift0+ self.Num_Bin )/width0)**power0)
		870	return model0+model0u+model0d+noise
		871
		872	def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
		873	shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
		874	model0=amplitude0numpy.exp(-0.5abs((x-shift0)/width0)**power0)
		875
		876	model0u=amplitude0numpy.exp(-0.5abs((x-shift0- self.Num_Bin )/width0)**power0)
		877
		878	model0d=amplitude0numpy.exp(-0.5abs((x-shift0+ self.Num_Bin )/width0)**power0)
		879	model1=amplitude1numpy.exp(-0.5abs((x-shift1)/width1)**power1)
		880
		881	model1u=amplitude1numpy.exp(-0.5abs((x-shift1- self.Num_Bin )/width1)**power1)
		882
		883	model1d=amplitude1numpy.exp(-0.5abs((x-shift1+ self.Num_Bin )/width1)**power1)
		884	return model0+model0u+model0d+model1+model1u+model1d+noise
		885
		886	def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
		887
		888	return num_intgsum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))*2)#/(64-5.) # /(64-5.) can be commented
		889
		890	def misfit2(self,state,y_data,x,num_intg):
		891	return num_intgsum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))*2)#/(64-9.)
		892
		893
		894	class PrecipitationProc(Operation):
		895
		896	'''
		897	Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R)
		898
		899	Input:
		900	self.dataOut.data_pre : SelfSpectra
		901
		902	Output:
		903
		904	self.dataOut.data_output : Reflectivity factor, rainfall Rate
		905
		906
		907	Parameters affected:
		908	'''
		909
		910
		911	def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None,
		912	tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None):
		913
		914	self.spc = dataOut.data_pre[0].copy()
		915	self.Num_Hei = self.spc.shape[2]
		916	self.Num_Bin = self.spc.shape[1]
		917	self.Num_Chn = self.spc.shape[0]
		918
		919	Velrange = dataOut.abscissaList
		920
		921	if radar == "MIRA35C" :
		922
		923	Ze = self.dBZeMODE2(dataOut)
		924
		925	else:
		926
		927	self.Pt = Pt
		928	self.Gt = Gt
		929	self.Gr = Gr
		930	self.Lambda = Lambda
		931	self.aL = aL
		932	self.tauW = tauW
		933	self.ThetaT = ThetaT
		934	self.ThetaR = ThetaR
		935
		936	RadarConstant = GetRadarConstant()
		937	SPCmean = numpy.mean(self.spc,0)
		938	ETA = numpy.zeros(self.Num_Hei)
		939	Pr = numpy.sum(SPCmean,0)
		940
		941	#for R in range(self.Num_Hei):
		942	# ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
		943
		944	D_range = numpy.zeros(self.Num_Hei)
		945	EqSec = numpy.zeros(self.Num_Hei)
		946	del_V = numpy.zeros(self.Num_Hei)
		947
		948	for R in range(self.Num_Hei):
		949	ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
		950
		951	h = R + Altitude #Range from ground to radar pulse altitude
		952	del_V[R] = 1 + 3.68 * 10*-5 h + 1.71 * 10*-9 h**2 #Density change correction for velocity
		953
		954	D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity
		955	SIGMA[R] = numpy.pi5 / Lambda4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma)
		956
		957	N_dist[R] = ETA[R] / SIGMA[R]
		958
		959	Ze = (ETA * Lambda*4) / (numpy.pi Km)
		960	Z = numpy.sum( N_dist * D_range**6 )
		961	RR = 610-4numpy.pi * numpy.sum( D_range*3 N_dist * Velrange ) #Rainfall rate
		962
		963
		964	RR = (Ze/200)**(1/1.6)
		965	dBRR = 10*numpy.log10(RR)
		966
		967	dBZe = 10*numpy.log10(Ze)
		968	dataOut.data_output = Ze
		969	dataOut.data_param = numpy.ones([2,self.Num_Hei])
		970	dataOut.channelList = [0,1]
		971	print 'channelList', dataOut.channelList
		972	dataOut.data_param[0]=dBZe
		973	dataOut.data_param[1]=dBRR
		974	print 'RR SHAPE', dBRR.shape
		975	print 'Ze SHAPE', dBZe.shape
		976	print 'dataOut.data_param SHAPE', dataOut.data_param.shape
		977
		978
		979	def dBZeMODE2(self, dataOut): # Processing for MIRA35C
		980
		981	NPW = dataOut.NPW
		982	COFA = dataOut.COFA
		983
		984	SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]])
		985	RadarConst = dataOut.RadarConst
		986	#frequency = 34.8510*9
		987
		988	ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei]))
		989	data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
		990
		991	ETA = numpy.sum(SNR,1)
		992	print 'ETA' , ETA
		993	ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
		994
		995	Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
		996
		997	for r in range(self.Num_Hei):
		998
		999	Ze[0,r] = ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2)
		1000	#Ze[1,r] = ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2)
		1001
		1002	return Ze
		1003
		1004	def GetRadarConstant(self):
		1005
		1006	"""
		1007	Constants:
		1008
		1009	Pt: Transmission Power dB
		1010	Gt: Transmission Gain dB
		1011	Gr: Reception Gain dB
		1012	Lambda: Wavelenght m
		1013	aL: Attenuation loses dB
		1014	tauW: Width of transmission pulse s
		1015	ThetaT: Transmission antenna bean angle rad
		1016	ThetaR: Reception antenna beam angle rad
		1017
		1018	"""
		1019	Numerator = ( (4numpy.pi)3 aL*2 16 * numpy.log(2) )
		1020	Denominator = ( Pt * Gt * Gr * Lambda*2 SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
		1021	RadarConstant = Numerator / Denominator
		1022
		1023	return RadarConstant
		1024
		1025
		1026
		1027	class FullSpectralAnalysis(Operation):
		1028
		1029	"""
		1030	Function that implements Full Spectral Analisys technique.
		1031
		1032	Input:
		1033	self.dataOut.data_pre : SelfSpectra and CrossSPectra data
		1034	self.dataOut.groupList : Pairlist of channels
		1035	self.dataOut.ChanDist : Physical distance between receivers
		1036
		1037
		1038	Output:
		1039
		1040	self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind
		1041
		1042
		1043	Parameters affected: Winds, height range, SNR
		1044
		1045	"""
		1046	def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7):
		1047
		1048	spc = dataOut.data_pre[0].copy()
		1049	cspc = dataOut.data_pre[1].copy()
		1050
		1051	nChannel = spc.shape[0]
		1052	nProfiles = spc.shape[1]
		1053	nHeights = spc.shape[2]
		1054
		1055	pairsList = dataOut.groupList
		1056	if dataOut.ChanDist is not None :
		1057	ChanDist = dataOut.ChanDist
		1058	else:
		1059	ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]])
		1060
		1061	#print 'ChanDist', ChanDist
		1062
		1063	if dataOut.VelRange is not None:
		1064	VelRange= dataOut.VelRange
		1065	else:
		1066	VelRange= dataOut.abscissaList
		1067
		1068	ySamples=numpy.ones([nChannel,nProfiles])
		1069	phase=numpy.ones([nChannel,nProfiles])
		1070	CSPCSamples=numpy.ones([nChannel,nProfiles],dtype=numpy.complex_)
		1071	coherence=numpy.ones([nChannel,nProfiles])
		1072	PhaseSlope=numpy.ones(nChannel)
		1073	PhaseInter=numpy.ones(nChannel)
		1074	dataSNR = dataOut.data_SNR
		1075
		1076
		1077
		1078	data = dataOut.data_pre
		1079	noise = dataOut.noise
		1080	print 'noise',noise
		1081	#SNRdB = 10*numpy.log10(dataOut.data_SNR)
		1082
		1083	FirstMoment = numpy.average(dataOut.data_param[:,1,:],0)
		1084	#SNRdBMean = []
		1085
		1086
		1087	#for j in range(nHeights):
		1088	# FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]]))
		1089	# SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
		1090
		1091	data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN
		1092
		1093	velocityX=[]
		1094	velocityY=[]
		1095	velocityV=[]
		1096
		1097	dbSNR = 10*numpy.log10(dataSNR)
		1098	dbSNR = numpy.average(dbSNR,0)
		1099	for Height in range(nHeights):
		1100
		1101	[Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit)
		1102
		1103	if abs(Vzon)<100. and abs(Vzon)> 0.:
		1104	velocityX=numpy.append(velocityX, Vzon)#Vmag
		1105
		1106	else:
		1107	print 'Vzon',Vzon
		1108	velocityX=numpy.append(velocityX, numpy.NaN)
		1109
		1110	if abs(Vmer)<100. and abs(Vmer) > 0.:
		1111	velocityY=numpy.append(velocityY, Vmer)#Vang
		1112
		1113	else:
		1114	print 'Vmer',Vmer
		1115	velocityY=numpy.append(velocityY, numpy.NaN)
		1116
		1117	if dbSNR[Height] > SNRlimit:
		1118	velocityV=numpy.append(velocityV, FirstMoment[Height])
		1119	else:
		1120	velocityV=numpy.append(velocityV, numpy.NaN)
		1121	#FirstMoment[Height]= numpy.NaN
		1122	# if SNRdBMean[Height] <12:
		1123	# FirstMoment[Height] = numpy.NaN
		1124	# velocityX[Height] = numpy.NaN
		1125	# velocityY[Height] = numpy.NaN
		1126
		1127
		1128	data_output[0]=numpy.array(velocityX)
		1129	data_output[1]=numpy.array(velocityY)
		1130	data_output[2]=-velocityV#FirstMoment
		1131
		1132	print ' '
		1133	#print 'FirstMoment'
		1134	#print FirstMoment
		1135	print 'velocityX',data_output[0]
		1136	print ' '
		1137	print 'velocityY',data_output[1]
		1138	#print numpy.array(velocityY)
		1139	print ' '
		1140	#print 'SNR'
		1141	#print 10*numpy.log10(dataOut.data_SNR)
		1142	#print numpy.shape(10*numpy.log10(dataOut.data_SNR))
		1143	print ' '
		1144
		1145
		1146	dataOut.data_output=data_output
		1147	return
		1148
		1149
		1150	def moving_average(self,x, N=2):
		1151	return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
		1152
		1153	def gaus(self,xSamples,a,x0,sigma):
		1154	return anumpy.exp(-(xSamples-x0)2/(2sigma**2))
		1155
		1156	def Find(self,x,value):
		1157	for index in range(len(x)):
		1158	if x[index]==value:
		1159	return index
		1160
		1161	def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit):
		1162
		1163	ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
		1164	phase=numpy.ones([spc.shape[0],spc.shape[1]])
		1165	CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
		1166	coherence=numpy.ones([spc.shape[0],spc.shape[1]])
		1167	PhaseSlope=numpy.ones(spc.shape[0])
		1168	PhaseInter=numpy.ones(spc.shape[0])
		1169	xFrec=VelRange
		1170
		1171	'''Getting Eij and Nij'''
		1172
		1173	E01=ChanDist[0][0]
		1174	N01=ChanDist[0][1]
		1175
		1176	E02=ChanDist[1][0]
		1177	N02=ChanDist[1][1]
		1178
		1179	E12=ChanDist[2][0]
		1180	N12=ChanDist[2][1]
		1181
		1182	z = spc.copy()
		1183	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
		1184
		1185	for i in range(spc.shape[0]):
		1186
		1187	'''**** Line of Data SPC ****'''
		1188	zline=z[i,:,Height]
		1189
		1190	'''**** SPC is normalized ****'''
		1191	FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy())
		1192	FactNorm= FactNorm/numpy.sum(FactNorm)
		1193
		1194	SmoothSPC=self.moving_average(FactNorm,N=3)
		1195
		1196	xSamples = ar(range(len(SmoothSPC)))
		1197	ySamples[i] = SmoothSPC
		1198
		1199	#dbSNR=10*numpy.log10(dataSNR)
		1200	print ' '
		1201	print ' '
		1202	print ' '
		1203
		1204	#print 'dataSNR', dbSNR.shape, dbSNR[0,40:120]
		1205	print 'SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20]
		1206	print 'noise',noise
		1207	print 'zline',zline.shape, zline[0:20]
		1208	print 'FactNorm',FactNorm.shape, FactNorm[0:20]
		1209	print 'FactNorm suma', numpy.sum(FactNorm)
		1210
		1211	for i in range(spc.shape[0]):
		1212
		1213	'''**** Line of Data CSPC ****'''
		1214	cspcLine=cspc[i,:,Height].copy()
		1215
		1216	'''**** CSPC is normalized ****'''
		1217	chan_index0 = pairsList[i][0]
		1218	chan_index1 = pairsList[i][1]
		1219	CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) #
		1220
		1221	CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor)
		1222
		1223	CSPCSamples[i] = CSPCNorm
		1224	coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
		1225
		1226	coherence[i]= self.moving_average(coherence[i],N=2)
		1227
		1228	phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
		1229
		1230	print 'cspcLine', cspcLine.shape, cspcLine[0:20]
		1231	print 'CSPCFactor', CSPCFactor#, CSPCFactor[0:20]
		1232	print numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i]
		1233	print 'CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20]
		1234	print 'CSPCNorm suma', numpy.sum(CSPCNorm)
		1235	print 'CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20]
		1236
		1237	'''**** Getting fij width ****'''
		1238
		1239	yMean=[]
		1240	yMean2=[]
		1241
		1242	for j in range(len(ySamples[1])):
		1243	yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
		1244
		1245	'''***** Getting fitting Gaussian ****'''
		1246	meanGauss=sum(xSamples*yMean) / len(xSamples)
		1247	sigma=sum(yMean(xSamples-meanGauss)*2) / len(xSamples)
		1248
		1249	print '****************************'
		1250	print 'len(xSamples): ',len(xSamples)
		1251	print 'yMean: ', yMean.shape, yMean[0:20]
		1252	print 'ySamples', ySamples.shape, ySamples[0,0:20]
		1253	print 'xSamples: ',xSamples.shape, xSamples[0:20]
		1254
		1255	print 'meanGauss',meanGauss
		1256	print 'sigma',sigma
		1257
		1258	#if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001):
		1259	if dbSNR > SNRlimit :
		1260	try:
		1261	popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
		1262
		1263	if numpy.amax(popt)>numpy.amax(yMean)*0.3:
		1264	FitGauss=self.gaus(xSamples,*popt)
		1265
		1266	else:
		1267	FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
		1268	print 'Verificador: Dentro', Height
		1269	except :#RuntimeError:
		1270	FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
		1271
		1272
		1273	else:
		1274	FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
		1275
		1276	Maximun=numpy.amax(yMean)
		1277	eMinus1=Maximunnumpy.exp(-1)#0.8
		1278
		1279	HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
		1280	HalfWidth= xFrec[HWpos]
		1281	GCpos=self.Find(FitGauss, numpy.amax(FitGauss))
		1282	Vpos=self.Find(FactNorm, numpy.amax(FactNorm))
		1283
		1284	#Vpos=FirstMoment[]
		1285
		1286	'''**** Getting Fij ****'''
		1287
		1288	GaussCenter=xFrec[GCpos]
		1289	if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
		1290	Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
		1291	else:
		1292	Fij=abs(GaussCenter-HalfWidth)+0.0000001
		1293
		1294	'''**** Getting Frecuency range of significant data ****'''
		1295
		1296	Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
		1297
		1298	if Rangpos<GCpos:
		1299	Range=numpy.array([Rangpos,2*GCpos-Rangpos])
		1300	elif Rangpos< ( len(xFrec)- len(xFrec)*0.1):
		1301	Range=numpy.array([2*GCpos-Rangpos,Rangpos])
		1302	else:
		1303	Range = numpy.array([0,0])
		1304
		1305	print ' '
		1306	print 'GCpos',GCpos, ( len(xFrec)- len(xFrec)*0.1)
		1307	print 'Rangpos',Rangpos
		1308	print 'RANGE: ', Range
		1309	FrecRange=xFrec[Range[0]:Range[1]]
		1310
		1311	'''**** Getting SCPC Slope ****'''
		1312
		1313	for i in range(spc.shape[0]):
		1314
		1315	if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.5:
		1316	PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
		1317
		1318	print 'FrecRange', len(FrecRange) , FrecRange
		1319	print 'PhaseRange', len(PhaseRange), PhaseRange
		1320	print ' '
		1321	if len(FrecRange) == len(PhaseRange):
		1322	slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
		1323	PhaseSlope[i]=slope
		1324	PhaseInter[i]=intercept
		1325	else:
		1326	PhaseSlope[i]=0
		1327	PhaseInter[i]=0
		1328	else:
		1329	PhaseSlope[i]=0
		1330	PhaseInter[i]=0
		1331
		1332	'''Getting constant C'''
		1333	cC=(Fijnumpy.pi)*2
		1334
		1335	'''**** Getting constants F and G ****'''
		1336	MijEijNij=numpy.array([[E02,N02], [E12,N12]])
		1337	MijResult0=(-PhaseSlope[1]cC) / (2numpy.pi)
		1338	MijResult1=(-PhaseSlope[2]cC) / (2numpy.pi)
		1339	MijResults=numpy.array([MijResult0,MijResult1])
		1340	(cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
		1341
		1342	'''**** Getting constants A, B and H ****'''
		1343	W01=numpy.amax(coherence[0])
		1344	W02=numpy.amax(coherence[1])
		1345	W12=numpy.amax(coherence[2])
		1346
		1347	WijResult0=((cFE01+cGN01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
		1348	WijResult1=((cFE02+cGN02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
		1349	WijResult2=((cFE12+cGN12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
		1350
		1351	WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
		1352
		1353	WijEijNij=numpy.array([ [E012, N012, 2E01N01] , [E022, N022, 2E02N02] , [E122, N122, 2E12N12] ])
		1354	(cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
		1355
		1356	VxVy=numpy.array([[cA,cH],[cH,cB]])
		1357
		1358	VxVyResults=numpy.array([-cF,-cG])
		1359	(Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
		1360
		1361	Vzon = Vy
		1362	Vmer = Vx
		1363	Vmag=numpy.sqrt(Vzon2+Vmer2)
		1364	Vang=numpy.arctan2(Vmer,Vzon)
		1365	Vver=xFrec[Vpos]
		1366	print 'vzon y vmer', Vzon, Vmer
		1367	return Vzon, Vmer, Vver, GaussCenter
		1368
	122	1369	class SpectralMoments(Operation):
	123	1370
	124	1371	'''
			@@ -162,25 +1409,23 class SpectralMoments(Operation):
	162	1409
	163	1410	dataOut.data_param = data_param[:,1:,:]
	164	1411	dataOut.data_SNR = data_param[:,0]
	165		dataOut.data_DOP = data_param[:,1]
	166		dataOut.data_MEAN = data_param[:,2]
	167		dataOut.data_STD = data_param[:,3]
	168	1412	return
	169	1413
	170		def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
		1414	def __calculateMoments(self, oldspec, oldfreq, n0,
		1415	nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
	171	1416
	172		if (nicoh is None): nicoh = 1
	173		if (graph is None): graph = 0
	174		if (smooth is None): smooth = 0
		1417	if (nicoh == None): nicoh = 1
		1418	if (graph == None): graph = 0
		1419	if (smooth == None): smooth = 0
	175	1420	elif (self.smooth < 3): smooth = 0
	176	1421
	177		if (type1 is None): type1 = 0
	178		if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
	179		if (snrth is None): snrth = -3
	180		if (dc is None): dc = 0
	181		if (aliasing is None): aliasing = 0
	182		if (oldfd is None): oldfd = 0
	183		if (wwauto is None): wwauto = 0
		1422	if (type1 == None): type1 = 0
		1423	if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
		1424	if (snrth == None): snrth = -3
		1425	if (dc == None): dc = 0
		1426	if (aliasing == None): aliasing = 0
		1427	if (oldfd == None): oldfd = 0
		1428	if (wwauto == None): wwauto = 0
	184	1429
	185	1430	if (n0 < 1.e-20): n0 = 1.e-20
	186	1431
			@@ -247,8 +1492,8 class SpectralMoments(Operation):
	247	1492	num_pairs = len(pairslist)
	248	1493
	249	1494	vel = self.dataOut.abscissaList
	250		spectra = self.dataOut.data_pre[0]
	251		cspectra = self.data~~Out~~.data_~~pre~~[1]
		1495	spectra = self.dataOut.data_pre
		1496	cspectra = self.dataIn.data_cspc
	252	1497	delta_v = vel[1] - vel[0]
	253	1498
	254	1499	#Calculating the power spectrum
			@@ -480,7 +1725,7 class SpectralFitting(Operation):
	480	1725	error1 = p0*numpy.nan
	481	1726
	482	1727	#Save
	483		if self.dataOut.data_param is None:
		1728	if self.dataOut.data_param == None:
	484	1729	self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
	485	1730	self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
	486	1731
			@@ -526,6 +1771,9 class WindProfiler(Operation):
	526	1771
	527	1772	n = None
	528	1773
		1774	def __init__(self):
		1775	Operation.__init__(self)
		1776
	529	1777	def __calculateCosDir(self, elev, azim):
	530	1778	zen = (90 - elev)*numpy.pi/180
	531	1779	azim = azim*numpy.pi/180
			@@ -816,7 +2064,7 class WindProfiler(Operation):
	816	2064	self.__dataReady = True
	817	2065	return
	818	2066
	819		def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax, ~~binkm~~=2):
		2067	def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
	820	2068	'''
	821	2069	Function that implements winds estimation technique with detected meteors.
	822	2070
			@@ -828,7 +2076,7 class WindProfiler(Operation):
	828	2076	'''
	829	2077	# print arrayMeteor.shape
	830	2078	#Settings
	831		nInt = (heightMax - heightMin)/~~binkm~~
		2079	nInt = (heightMax - heightMin)/2
	832	2080	# print nInt
	833	2081	nInt = int(nInt)
	834	2082	# print nInt
			@@ -1130,11 +2378,6 class WindProfiler(Operation):
	1130	2378	hmax = kwargs['hmax']
	1131	2379	else: hmax = 110
	1132	2380
	1133		if kwargs.has_key('BinKm'):
	1134		binkm = kwargs['BinKm']
	1135		else:
	1136		binkm = 2
	1137
	1138	2381	dataOut.outputInterval = nHours*3600
	1139	2382
	1140	2383	if self.__isConfig == False:
			@@ -1145,7 +2388,7 class WindProfiler(Operation):
	1145	2388
	1146	2389	self.__isConfig = True
	1147	2390
	1148		if self.__buffer is None:
		2391	if self.__buffer == None:
	1149	2392	self.__buffer = dataOut.data_param
	1150	2393	self.__firstdata = copy.copy(dataOut)
	1151	2394
			@@ -1159,7 +2402,7 class WindProfiler(Operation):
	1159	2402
	1160	2403	self.__initime += dataOut.outputInterval #to erase time offset
	1161	2404
	1162		dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax, ~~binkm~~)
		2405	dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
	1163	2406	dataOut.flagNoData = False
	1164	2407	self.__buffer = None
	1165	2408
			@@ -1187,7 +2430,7 class WindProfiler(Operation):
	1187	2430	else: mode = 'SA'
	1188	2431
	1189	2432	#Borrar luego esto
	1190		if dataOut.groupList is None:
		2433	if dataOut.groupList == None:
	1191	2434	dataOut.groupList = [(0,1),(0,2),(1,2)]
	1192	2435	groupList = dataOut.groupList
	1193	2436	C = 3e8
			@@ -1209,7 +2452,7 class WindProfiler(Operation):
	1209	2452
	1210	2453	self.__isConfig = True
	1211	2454
	1212		if self.__buffer is None:
		2455	if self.__buffer == None:
	1213	2456	self.__buffer = dataOut.data_param
	1214	2457	self.__firstdata = copy.copy(dataOut)
	1215	2458
			@@ -1235,6 +2478,8 class WindProfiler(Operation):
	1235	2478
	1236	2479	class EWDriftsEstimation(Operation):
	1237	2480
		2481	def __init__(self):
		2482	Operation.__init__(self)
	1238	2483
	1239	2484	def __correctValues(self, heiRang, phi, velRadial, SNR):
	1240	2485	listPhi = phi.tolist()
			@@ -1569,7 +2814,7 class SMDetection(Operation):
	1569	2814
	1570	2815
	1571	2816	#Getting Pairslist
	1572		if channelPositions is None:
		2817	if channelPositions == None:
	1573	2818	# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
	1574	2819	channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
	1575	2820	meteorOps = SMOperations()
			@@ -1700,7 +2945,7 class SMDetection(Operation):
	1700	2945	# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
	1701	2946	dataOut.data_param = arrayParameters
	1702	2947
	1703		if arrayParameters is None:
		2948	if arrayParameters == None:
	1704	2949	dataOut.flagNoData = True
	1705	2950	else:
	1706	2951	dataOut.flagNoData = True
			@@ -1894,7 +3139,7 class SMDetection(Operation):
	1894	3139
	1895	3140	if (indDNthresh.size > 0):
	1896	3141	indEnd = indDNthresh[0] - 1
	1897		indInit = indUPthresh[j] if isinstance(indUPthresh[j], (int, float)) else indUPthresh[j][0] ##CHECK!!!!
		3142	indInit = indUPthresh[j]
	1898	3143
	1899	3144	meteor = powerAux[indInit:indEnd + 1]
	1900	3145	indPeak = meteor.argmax() + indInit
			@@ -1949,19 +3194,19 class SMDetection(Operation):
	1949	3194	indSides = pairsarray[:,1]
	1950	3195
	1951	3196	pairslist1 = list(pairslist)
	1952		pairslist1.append((0,4))
	1953		pairslist1.append((1,3))
		3197	pairslist1.append((0,1))
		3198	pairslist1.append((3,4))
	1954	3199
	1955	3200	listMeteors1 = []
	1956	3201	listPowerSeries = []
	1957	3202	listVoltageSeries = []
	1958	3203	#volts has the war data
	1959	3204
	1960		if frequency == 30~~.175~~e6:
		3205	if frequency == 30e6:
	1961	3206	timeLag = 4510*-3
	1962	3207	else:
	1963	3208	timeLag = 1510*-3
	1964		lag = ~~int~~(numpy.ceil(timeLag/timeInterval))
		3209	lag = numpy.ceil(timeLag/timeInterval)
	1965	3210
	1966	3211	for i in range(len(listMeteors)):
	1967	3212
			@@ -1969,10 +3214,10 class SMDetection(Operation):
	1969	3214	meteorAux = numpy.zeros(16)
	1970	3215
	1971	3216	#Loading meteor Data (mHeight, mStart, mPeak, mEnd)
	1972		mHeight = ~~int~~(listMeteors[i][0])
	1973		mStart = ~~int~~(listMeteors[i][1])
	1974		mPeak = ~~int~~(listMeteors[i][2])
	1975		mEnd = ~~int~~(listMeteors[i][3])
		3217	mHeight = listMeteors[i][0]
		3218	mStart = listMeteors[i][1]
		3219	mPeak = listMeteors[i][2]
		3220	mEnd = listMeteors[i][3]
	1976	3221
	1977	3222	#get the volt data between the start and end times of the meteor
	1978	3223	meteorVolts = volts[:,mStart:mEnd+1,mHeight]
			@@ -2061,11 +3306,11 class SMDetection(Operation):
	2061	3306
	2062	3307	threshError = 10
	2063	3308	#Depending if it is 30 or 50 MHz
	2064		if frequency == 30~~.175~~e6:
		3309	if frequency == 30e6:
	2065	3310	timeLag = 4510*-3
	2066	3311	else:
	2067	3312	timeLag = 1510*-3
	2068		lag = ~~int~~(numpy.ceil(timeLag/timeInterval))
		3313	lag = numpy.ceil(timeLag/timeInterval)
	2069	3314
	2070	3315	listMeteors1 = []
	2071	3316
			@@ -2121,14 +3366,14 class SMDetection(Operation):
	2121	3366	def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
	2122	3367
	2123	3368	pairslist1 = list(pairslist)
	2124		pairslist1.append((0,4))
	2125		pairslist1.append((1,3))
		3369	pairslist1.append((0,1))
		3370	pairslist1.append((3,4))
	2126	3371	numPairs = len(pairslist1)
	2127	3372	#Time Lag
	2128	3373	timeLag = 4510*-3
	2129	3374	c = 3e8
	2130	3375	lag = numpy.ceil(timeLag/timeInterval)
	2131		freq = 30~~.175~~e6
		3376	freq = 30e6
	2132	3377
	2133	3378	listMeteors1 = []
	2134	3379
			@@ -2149,7 +3394,7 class SMDetection(Operation):
	2149	3394	#Method 2
	2150	3395	slopes = numpy.zeros(numPairs)
	2151	3396	time = numpy.array([-2,-1,1,2])*timeInterval
	2152		angAllCCF = numpy.angle(allCCFs[:,[0,4,2,3],0])
		3397	angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
	2153	3398
	2154	3399	#Correct phases
	2155	3400	derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
			@@ -2236,7 +3481,7 class CorrectSMPhases(Operation):
	2236	3481	arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
	2237	3482
	2238	3483	meteorOps = SMOperations()
	2239		if channelPositions is None:
		3484	if channelPositions == None:
	2240	3485	# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
	2241	3486	channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
	2242	3487
			@@ -2413,7 +3658,7 class SMPhaseCalibration(Operation):
	2413	3658
	2414	3659	self.__isConfig = True
	2415	3660
	2416		if self.__buffer is None:
		3661	if self.__buffer == None:
	2417	3662	self.__buffer = dataOut.data_param.copy()
	2418	3663
	2419	3664	else:
			@@ -2469,7 +3714,6 class SMPhaseCalibration(Operation):
	2469	3714	phasesOff = phasesOff.reshape((1,phasesOff.size))
	2470	3715	dataOut.data_output = -phasesOff
	2471	3716	dataOut.flagNoData = False
	2472		dataOut.channelList = pairslist0
	2473	3717	self.__buffer = None
	2474	3718
	2475	3719
			@@ -2803,3 +4047,4 class SMOperations():
	2803	4047	# error[indInvalid1] = 13
	2804	4048	#
	2805	4049	# return heights, error
		4050	No newline at end of file

schainpy/model/proc/jroproc_spectra.py +1 0

                          dataOut.nIncohInt *= self.n
                          dataOut.utctime = avgdatatime
                          dataOut.flagNoData = False

setup.py +1 -1

                    scripts=['schainpy/gui/schainGUI'],
                    ext_modules=[Extension("cSchain", ["schainpy/model/proc/extensions.c"], include_dirs=[numpy.get_include()])],
                    install_requires=[
-                                   "scipy >= 0.14.0",
+                                   "scipy >= 0.13.0",
                                    "h5py >= 2.2.1",
                                    "matplotlib >= 1.4.2",
                                    "pyfits >= 3.4",

schainpy/trash removed 0 -1

			@@ -1,1 +0,0
	1		You should install "digital_rf_hdf5" module if you want to read USRP data

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