schain Commit - r1032:c711a14430f7 · Jicamarca Repository

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):
+                def __init__(self):
-                    Figure.__init__(self, **kwargs)
                     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):
+                def __init__(self):
-                    Figure.__init__(self, **kwargs)
                     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	from plotting_codes import *	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	class MomentsPlot(Figure):	220	class MomentsPlot(Figure):
10		221
11	isConfig = None	222	isConfig = None
@@ -446,8 +657,7 class WindProfilerPlot(Figure):
446	# tmin = None	657	# tmin = None
447	# tmax = None	658	# tmax = None
448		659
449		660	x = dataOut.getTimeRange1(dataOut.paramInterval)
450	x = dataOut.getTimeRange1(dataOut.outputInterval)
451	y = dataOut.heightList	661	y = dataOut.heightList
452	z = dataOut.data_output.copy()	662	z = dataOut.data_output.copy()
453	nplots = z.shape[0] #Number of wind dimensions estimated	663	nplots = z.shape[0] #Number of wind dimensions estimated
@@ -559,7 +769,7 class WindProfilerPlot(Figure):
559	thisDatetime=thisDatetime,	769	thisDatetime=thisDatetime,
560	update_figfile=update_figfile)	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	self.counter_imagwr = wr_period	773	self.counter_imagwr = wr_period
564	self.isConfig = False	774	self.isConfig = False
565	update_figfile = True	775	update_figfile = True
@@ -636,12 +846,12 class ParametersPlot(Figure):
636		846
637	counter += 1	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	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, timerange=None,	850	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, timerange=None,
641	showSNR=False, SNRthresh = -numpy.inf, SNRmin=None, SNRmax=None,	851	showSNR=False, SNRthresh = -numpy.inf, SNRmin=None, SNRmax=None,
642	save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,	852	save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
643	server=None, folder=None, username=None, password=None,	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	Input:	857	Input:
@@ -658,10 +868,9 class ParametersPlot(Figure):
658	zmax : None	868	zmax : None
659	"""	869	"""
660		870
661	if colormap:	871	if HEIGHT is not None:
662	colormap="jet"	872	self.HEIGHT = HEIGHT
663	~~else~~:	873
664	colormap="RdBu_r"
665		874
666	if not isTimeInHourRange(dataOut.datatime, xmin, xmax):	875	if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
667	return	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 = dataOut.data_spc/factor
-            #         z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
+                    z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
-            #         avg = numpy.average(z, axis=1)
+                    avg = numpy.average(z, axis=1)
-            #         avgdB = 10.*numpy.log10(avg)
+                    avgdB = 10.*numpy.log10(avg)
-                    avgdB = dataOut.getPower()
+                    # 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

@@ -12,3 +12,10 from jroIO_usrp import *
12	from jroIO_kamisr import *	12	from jroIO_kamisr import *
13	from jroIO_param import *	13	from jroIO_param import *
14	from jroIO_hf import *	14	from jroIO_hf import *
		15
		16	from jroIO_madrigal import *
		17
		18	from bltrIO_param import *
		19	from jroIO_bltr import *
		20	from jroIO_mira35c import *
		21

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

                     print "[Reading] Number of read blocks %04d" %self.nTotalBlocks
                 def printNumberOfBlock(self):
+                    'SPAM!'
-                    if self.flagIsNewBlock:
+            #         if self.flagIsNewBlock:
-                        print "[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
+            #             print "[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
-                                                                  self.processingHeaderObj.dataBlocksPerFile,
+            #                                                       self.processingHeaderObj.dataBlocksPerFile,
-                                                                  self.dataOut.datatime.ctime())
+            #                                                       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)):
+            #         for i in range(len(filenameList)):
-                        print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
+            #             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	import numpy	1	import numpy
	2	import math	2	import math
	3	from scipy import optimize, interpolate, signal, stats, ndimage	3	from scipy import optimize, interpolate, signal, stats, ndimage
			4	import scipy
	4	import re	5	import re
	5	import datetime	6	import datetime
	6	import copy	7	import copy
	7	import sys	8	import sys
	8	import importlib	9	import importlib
	9	import itertools	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	from jroproc_base import ProcessingUnit, Operation	23	from jroproc_base import ProcessingUnit, Operation
	12	from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon	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	class ParametersProc(ProcessingUnit):	62	class ParametersProc(ProcessingUnit):
	@@ -57,7 +104,9 class ParametersProc(ProcessingUnit):
	57	self.dataOut.timeInterval1 = self.dataIn.timeInterval	104	self.dataOut.timeInterval1 = self.dataIn.timeInterval
	58	self.dataOut.heightList = self.dataIn.getHeiRange()	105	self.dataOut.heightList = self.dataIn.getHeiRange()
	59	self.dataOut.frequency = self.dataIn.frequency	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	def run(self):	111	def run(self):
	63		112
	@@ -77,17 +126,36 class ParametersProc(ProcessingUnit):
	77	if self.dataIn.type == "Spectra":	126	if self.dataIn.type == "Spectra":
	78		127
	79	self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)	128	self.dataOut.data_pre = (self.dataIn.data_spc,self.dataIn.data_cspc)
	80	~~self~~.~~dataOut~~.~~data_spc~~ = self.dataIn.data_spc	129	print 'self.dataIn.data_spc', self.dataIn.data_spc.shape
	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
	89	self.dataOut.abscissaList = self.dataIn.getVelRange(1)	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	self.dataOut.flagNoData = False	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	#---------------------- Correlation Data ---------------------------	160	#---------------------- Correlation Data ---------------------------
	93		161
	@@ -108,7 +176,6 class ParametersProc(ProcessingUnit):
	108		176
	109	if self.dataIn.type == "Parameters":	177	if self.dataIn.type == "Parameters":
	110	self.dataOut.copy(self.dataIn)	178	self.dataOut.copy(self.dataIn)
	111	self.dataOut.utctimeInit = self.dataIn.utctime
	112	self.dataOut.flagNoData = False	179	self.dataOut.flagNoData = False
	113		180
	114	return True	181	return True
	@@ -119,6 +186,1186 class ParametersProc(ProcessingUnit):
	119		186
	120	return	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	class SpectralMoments(Operation):	1369	class SpectralMoments(Operation):
	123		1370
	124	'''	1371	'''
	@@ -162,25 +1409,23 class SpectralMoments(Operation):
	162		1409
	163	dataOut.data_param = data_param[:,1:,:]	1410	dataOut.data_param = data_param[:,1:,:]
	164	dataOut.data_SNR = data_param[:,0]	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	return	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	1417	if (nicoh == None): nicoh = 1
	173	if (graph is None): graph = 0	1418	if (graph == None): graph = 0
	174	if (smooth is None): smooth = 0	1419	if (smooth == None): smooth = 0
	175	elif (self.smooth < 3): smooth = 0	1420	elif (self.smooth < 3): smooth = 0
	176		1421
	177	if (type1 is None): type1 = 0	1422	if (type1 == None): type1 = 0
	178	if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1	1423	if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
	179	if (snrth is None): snrth = -3	1424	if (snrth == None): snrth = -3
	180	if (dc is None): dc = 0	1425	if (dc == None): dc = 0
	181	if (aliasing is None): aliasing = 0	1426	if (aliasing == None): aliasing = 0
	182	if (oldfd is None): oldfd = 0	1427	if (oldfd == None): oldfd = 0
	183	if (wwauto is None): wwauto = 0	1428	if (wwauto == None): wwauto = 0
	184		1429
	185	if (n0 < 1.e-20): n0 = 1.e-20	1430	if (n0 < 1.e-20): n0 = 1.e-20
	186		1431
	@@ -247,8 +1492,8 class SpectralMoments(Operation):
	247	num_pairs = len(pairslist)	1492	num_pairs = len(pairslist)
	248		1493
	249	vel = self.dataOut.abscissaList	1494	vel = self.dataOut.abscissaList
	250	spectra = self.dataOut.data_pre[0]	1495	spectra = self.dataOut.data_pre
	251	cspectra = self.data~~Out~~.data_~~pre~~[1]	1496	cspectra = self.dataIn.data_cspc
	252	delta_v = vel[1] - vel[0]	1497	delta_v = vel[1] - vel[0]
	253		1498
	254	#Calculating the power spectrum	1499	#Calculating the power spectrum
	@@ -480,7 +1725,7 class SpectralFitting(Operation):
	480	error1 = p0*numpy.nan	1725	error1 = p0*numpy.nan
	481		1726
	482	#Save	1727	#Save
	483	if self.dataOut.data_param is None:	1728	if self.dataOut.data_param == None:
	484	self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan	1729	self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
	485	self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan	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	n = None	1772	n = None
	528		1773
			1774	def __init__(self):
			1775	Operation.__init__(self)
			1776
	529	def __calculateCosDir(self, elev, azim):	1777	def __calculateCosDir(self, elev, azim):
	530	zen = (90 - elev)*numpy.pi/180	1778	zen = (90 - elev)*numpy.pi/180
	531	azim = azim*numpy.pi/180	1779	azim = azim*numpy.pi/180
	@@ -816,7 +2064,7 class WindProfiler(Operation):
	816	self.__dataReady = True	2064	self.__dataReady = True
	817	return	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	Function that implements winds estimation technique with detected meteors.	2069	Function that implements winds estimation technique with detected meteors.
	822		2070
	@@ -828,7 +2076,7 class WindProfiler(Operation):
	828	'''	2076	'''
	829	# print arrayMeteor.shape	2077	# print arrayMeteor.shape
	830	#Settings	2078	#Settings
	831	nInt = (heightMax - heightMin)/~~binkm~~	2079	nInt = (heightMax - heightMin)/2
	832	# print nInt	2080	# print nInt
	833	nInt = int(nInt)	2081	nInt = int(nInt)
	834	# print nInt	2082	# print nInt
	@@ -1130,11 +2378,6 class WindProfiler(Operation):
	1130	hmax = kwargs['hmax']	2378	hmax = kwargs['hmax']
	1131	else: hmax = 110	2379	else: hmax = 110
	1132		2380
	1133	if kwargs.has_key('BinKm'):
	1134	binkm = kwargs['BinKm']
	1135	else:
	1136	binkm = 2
	1137
	1138	dataOut.outputInterval = nHours*3600	2381	dataOut.outputInterval = nHours*3600
	1139		2382
	1140	if self.__isConfig == False:	2383	if self.__isConfig == False:
	@@ -1145,7 +2388,7 class WindProfiler(Operation):
	1145		2388
	1146	self.__isConfig = True	2389	self.__isConfig = True
	1147		2390
	1148	if self.__buffer is None:	2391	if self.__buffer == None:
	1149	self.__buffer = dataOut.data_param	2392	self.__buffer = dataOut.data_param
	1150	self.__firstdata = copy.copy(dataOut)	2393	self.__firstdata = copy.copy(dataOut)
	1151		2394
	@@ -1159,7 +2402,7 class WindProfiler(Operation):
	1159		2402
	1160	self.__initime += dataOut.outputInterval #to erase time offset	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	dataOut.flagNoData = False	2406	dataOut.flagNoData = False
	1164	self.__buffer = None	2407	self.__buffer = None
	1165		2408
	@@ -1187,7 +2430,7 class WindProfiler(Operation):
	1187	else: mode = 'SA'	2430	else: mode = 'SA'
	1188		2431
	1189	#Borrar luego esto	2432	#Borrar luego esto
	1190	if dataOut.groupList is None:	2433	if dataOut.groupList == None:
	1191	dataOut.groupList = [(0,1),(0,2),(1,2)]	2434	dataOut.groupList = [(0,1),(0,2),(1,2)]
	1192	groupList = dataOut.groupList	2435	groupList = dataOut.groupList
	1193	C = 3e8	2436	C = 3e8
	@@ -1209,7 +2452,7 class WindProfiler(Operation):
	1209		2452
	1210	self.__isConfig = True	2453	self.__isConfig = True
	1211		2454
	1212	if self.__buffer is None:	2455	if self.__buffer == None:
	1213	self.__buffer = dataOut.data_param	2456	self.__buffer = dataOut.data_param
	1214	self.__firstdata = copy.copy(dataOut)	2457	self.__firstdata = copy.copy(dataOut)
	1215		2458
	@@ -1235,6 +2478,8 class WindProfiler(Operation):
	1235		2478
	1236	class EWDriftsEstimation(Operation):	2479	class EWDriftsEstimation(Operation):
	1237		2480
			2481	def __init__(self):
			2482	Operation.__init__(self)
	1238		2483
	1239	def __correctValues(self, heiRang, phi, velRadial, SNR):	2484	def __correctValues(self, heiRang, phi, velRadial, SNR):
	1240	listPhi = phi.tolist()	2485	listPhi = phi.tolist()
	@@ -1569,7 +2814,7 class SMDetection(Operation):
	1569		2814
	1570		2815
	1571	#Getting Pairslist	2816	#Getting Pairslist
	1572	if channelPositions is None:	2817	if channelPositions == None:
	1573	# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T	2818	# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
	1574	channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella	2819	channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
	1575	meteorOps = SMOperations()	2820	meteorOps = SMOperations()
	@@ -1700,7 +2945,7 class SMDetection(Operation):
	1700	# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))	2945	# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
	1701	dataOut.data_param = arrayParameters	2946	dataOut.data_param = arrayParameters
	1702		2947
	1703	if arrayParameters is None:	2948	if arrayParameters == None:
	1704	dataOut.flagNoData = True	2949	dataOut.flagNoData = True
	1705	else:	2950	else:
	1706	dataOut.flagNoData = True	2951	dataOut.flagNoData = True
	@@ -1894,7 +3139,7 class SMDetection(Operation):
	1894		3139
	1895	if (indDNthresh.size > 0):	3140	if (indDNthresh.size > 0):
	1896	indEnd = indDNthresh[0] - 1	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	meteor = powerAux[indInit:indEnd + 1]	3144	meteor = powerAux[indInit:indEnd + 1]
	1900	indPeak = meteor.argmax() + indInit	3145	indPeak = meteor.argmax() + indInit
	@@ -1949,19 +3194,19 class SMDetection(Operation):
	1949	indSides = pairsarray[:,1]	3194	indSides = pairsarray[:,1]
	1950		3195
	1951	pairslist1 = list(pairslist)	3196	pairslist1 = list(pairslist)
	1952	pairslist1.append((0,4))	3197	pairslist1.append((0,1))
	1953	pairslist1.append((1,3))	3198	pairslist1.append((3,4))
	1954		3199
	1955	listMeteors1 = []	3200	listMeteors1 = []
	1956	listPowerSeries = []	3201	listPowerSeries = []
	1957	listVoltageSeries = []	3202	listVoltageSeries = []
	1958	#volts has the war data	3203	#volts has the war data
	1959		3204
	1960	if frequency == 30~~.175~~e6:	3205	if frequency == 30e6:
	1961	timeLag = 4510*-3	3206	timeLag = 4510*-3
	1962	else:	3207	else:
	1963	timeLag = 1510*-3	3208	timeLag = 1510*-3
	1964	lag = ~~int~~(numpy.ceil(timeLag/timeInterval))	3209	lag = numpy.ceil(timeLag/timeInterval)
	1965		3210
	1966	for i in range(len(listMeteors)):	3211	for i in range(len(listMeteors)):
	1967		3212
	@@ -1969,10 +3214,10 class SMDetection(Operation):
	1969	meteorAux = numpy.zeros(16)	3214	meteorAux = numpy.zeros(16)
	1970		3215
	1971	#Loading meteor Data (mHeight, mStart, mPeak, mEnd)	3216	#Loading meteor Data (mHeight, mStart, mPeak, mEnd)
	1972	mHeight = ~~int~~(listMeteors[i][0])	3217	mHeight = listMeteors[i][0]
	1973	mStart = ~~int~~(listMeteors[i][1])	3218	mStart = listMeteors[i][1]
	1974	mPeak = ~~int~~(listMeteors[i][2])	3219	mPeak = listMeteors[i][2]
	1975	mEnd = ~~int~~(listMeteors[i][3])	3220	mEnd = listMeteors[i][3]
	1976		3221
	1977	#get the volt data between the start and end times of the meteor	3222	#get the volt data between the start and end times of the meteor
	1978	meteorVolts = volts[:,mStart:mEnd+1,mHeight]	3223	meteorVolts = volts[:,mStart:mEnd+1,mHeight]
	@@ -2061,11 +3306,11 class SMDetection(Operation):
	2061		3306
	2062	threshError = 10	3307	threshError = 10
	2063	#Depending if it is 30 or 50 MHz	3308	#Depending if it is 30 or 50 MHz
	2064	if frequency == 30~~.175~~e6:	3309	if frequency == 30e6:
	2065	timeLag = 4510*-3	3310	timeLag = 4510*-3
	2066	else:	3311	else:
	2067	timeLag = 1510*-3	3312	timeLag = 1510*-3
	2068	lag = ~~int~~(numpy.ceil(timeLag/timeInterval))	3313	lag = numpy.ceil(timeLag/timeInterval)
	2069		3314
	2070	listMeteors1 = []	3315	listMeteors1 = []
	2071		3316
	@@ -2121,14 +3366,14 class SMDetection(Operation):
	2121	def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):	3366	def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
	2122		3367
	2123	pairslist1 = list(pairslist)	3368	pairslist1 = list(pairslist)
	2124	pairslist1.append((0,4))	3369	pairslist1.append((0,1))
	2125	pairslist1.append((1,3))	3370	pairslist1.append((3,4))
	2126	numPairs = len(pairslist1)	3371	numPairs = len(pairslist1)
	2127	#Time Lag	3372	#Time Lag
	2128	timeLag = 4510*-3	3373	timeLag = 4510*-3
	2129	c = 3e8	3374	c = 3e8
	2130	lag = numpy.ceil(timeLag/timeInterval)	3375	lag = numpy.ceil(timeLag/timeInterval)
	2131	freq = 30~~.175~~e6	3376	freq = 30e6
	2132		3377
	2133	listMeteors1 = []	3378	listMeteors1 = []
	2134		3379
	@@ -2149,7 +3394,7 class SMDetection(Operation):
	2149	#Method 2	3394	#Method 2
	2150	slopes = numpy.zeros(numPairs)	3395	slopes = numpy.zeros(numPairs)
	2151	time = numpy.array([-2,-1,1,2])*timeInterval	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	#Correct phases	3399	#Correct phases
	2155	derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]	3400	derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
	@@ -2236,7 +3481,7 class CorrectSMPhases(Operation):
	2236	arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))	3481	arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
	2237		3482
	2238	meteorOps = SMOperations()	3483	meteorOps = SMOperations()
	2239	if channelPositions is None:	3484	if channelPositions == None:
	2240	# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T	3485	# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
	2241	channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella	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	self.__isConfig = True	3659	self.__isConfig = True
	2415		3660
	2416	if self.__buffer is None:	3661	if self.__buffer == None:
	2417	self.__buffer = dataOut.data_param.copy()	3662	self.__buffer = dataOut.data_param.copy()
	2418		3663
	2419	else:	3664	else:
	@@ -2469,7 +3714,6 class SMPhaseCalibration(Operation):
	2469	phasesOff = phasesOff.reshape((1,phasesOff.size))	3714	phasesOff = phasesOff.reshape((1,phasesOff.size))
	2470	dataOut.data_output = -phasesOff	3715	dataOut.data_output = -phasesOff
	2471	dataOut.flagNoData = False	3716	dataOut.flagNoData = False
	2472	dataOut.channelList = pairslist0
	2473	self.__buffer = None	3717	self.__buffer = None
	2474		3718
	2475		3719
	@@ -2803,3 +4047,4 class SMOperations():
	2803	# error[indInvalid1] = 13	4047	# error[indInvalid1] = 13
	2804	#	4048	#
	2805	# return heights, error	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

NO CONTENT: file was removed

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