schain-jc Commit - r1032:c711a14430f7

merge with BLTR and Madrigal modules

Juan C. Espinoza -

r1032:c711a14430f7

parent child

Context file:

r1032:c711a14430f7

Expand all files

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

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

                 nAvg = None
                 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 +47 -47

             from figure import Figure, isRealtime
             class CorrelationPlot(Figure):
                 isConfig = None
                 __nsubplots = None
                 WIDTHPROF = None
                 HEIGHTPROF = None
                 PREFIX = 'corr'
-                def __init__(self, **kwargs):
+                def __init__(self):
-                    Figure.__init__(self, **kwargs)
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.WIDTH = 280
                     self.HEIGHT = 250
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.counter_imagwr = 0
                     self.PLOT_CODE = 1
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                 def getSubplots(self):
                     ncol = int(numpy.sqrt(self.nplots)+0.9)
                     nrow = int(self.nplots*1./ncol + 0.9)
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
                     showprofile = False
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 1
                     colspan = 1
                     if showprofile:
                         ncolspan = 3
                         colspan = 2
                         self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH + self.WIDTHPROF,
                                       heightplot = self.HEIGHT + self.HEIGHTPROF,
                                       show=show)
                     nrow, ncol = self.getSubplots()
                     counter = 0
                     for y in range(nrow):
                         for x in range(ncol):
                             if counter >= self.nplots:
                                 break
                             self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
                             if showprofile:
                                 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
                             counter += 1
                 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
                         xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
                         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):
                     """
                     Input:
                         dataOut         :
                         id        :
                         zmin            :    None,
                         zmax            :    None
                     """
                     if dataOut.flagNoData:
                         return None
                     if realtime:
                         if not(isRealtime(utcdatatime = dataOut.utctime)):
                             print 'Skipping this plot function'
                             return
                     if channelList == None:
                         channelIndexList = dataOut.channelIndexList
                     else:
                             if channel not in dataOut.channelList:
                                 raise ValueError, "Channel %d is not in dataOut.channelList"
                             channelIndexList.append(dataOut.channelList.index(channel))
                     factor = dataOut.normFactor
                     lenfactor = factor.shape[1]
                     x = dataOut.getLagTRange(1)
                     y = dataOut.getHeiRange()
                     z = copy.copy(dataOut.data_corr[:,:,0,:])
                     for i in range(dataOut.data_corr.shape[0]):
                         z[i,:,:] = z[i,:,:]/factor[i,:]
                     zdB = numpy.abs(z)
                     avg = numpy.average(z, axis=1)
             #         avg = numpy.nanmean(z, axis=1)
             #         noise = dataOut.noise/factor
                     #thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
                     title = wintitle + " Correlation"
                     xlabel = "Lag T (s)"
                     ylabel = "Range (Km)"
                     if not self.isConfig:
                         nplots = dataOut.data_corr.shape[0]
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         if zmin == None: zmin = 0
                         if zmax == None: zmax = 1
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.isConfig = True
                     self.setWinTitle(title)
                     for i in range(self.nplots):
                         str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
                         title = "Channel %d and %d: : %s" %(dataOut.pairsList[i][0],dataOut.pairsList[i][1] , str_datetime)
                                     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
                                     xlabel=xlabel, ylabel=ylabel, title=title,
                                     ticksize=9, cblabel='')
             #             if self.__showprofile:
             #                 axes = self.axesList[i*self.__nsubplots +1]
             #                 axes.pline(avgdB[i], y,
             #                         xlabel='dB', ylabel='', title='',
             #                         ytick_visible=False,
             #                         grid='x')
+            #
             #                 noiseline = numpy.repeat(noisedB[i], len(y))
             #                 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
                     self.draw()
                     self.save(figpath=figpath,
                               figfile=figfile,
                               save=save,
                               ftp=ftp,
                               wr_period=wr_period,
                               thisDatetime=thisDatetime)

schainpy/model/graphics/jroplot_heispectra.py +81 -82

             from plotting_codes import *
             class SpectraHeisScope(Figure):
                 isConfig = None
                 __nsubplots = None
                 WIDTHPROF = None
                 HEIGHTPROF = None
                 PREFIX = 'spc'
-                def __init__(self, **kwargs):
+                def __init__(self):
-                    Figure.__init__(self, **kwargs)
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.WIDTH = 230
                     self.HEIGHT = 250
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.counter_imagwr = 0
                     self.PLOT_CODE = SPEC_CODE
                 def getSubplots(self):
                     ncol = int(numpy.sqrt(self.nplots)+0.9)
                     nrow = int(self.nplots*1./ncol + 0.9)
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, show):
                     showprofile = False
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 1
                     colspan = 1
                     if showprofile:
                         ncolspan = 3
                         colspan = 2
                         self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH + self.WIDTHPROF,
                                       heightplot = self.HEIGHT + self.HEIGHTPROF,
                                       show = show)
                     nrow, ncol = self.getSubplots()
                     counter = 0
                     for y in range(nrow):
                         for x in range(ncol):
                             if counter >= self.nplots:
                                 break
                             self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
                             if showprofile:
                                 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
                             counter += 1
                 def run(self, dataOut, id, wintitle="", channelList=None,
                         xmin=None, xmax=None, ymin=None, ymax=None, save=False,
                         figpath='./', 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):
                     """
                     Input:
                         dataOut         :
                         id        :
                         ymin            :    None,
                         ymax            :    None,
                     """
                     if dataOut.realtime:
                         if not(isRealtime(utcdatatime = dataOut.utctime)):
                             print 'Skipping this plot function'
                             return
                     if channelList == None:
                         channelIndexList = dataOut.channelIndexList
                     else:
                             if channel not in dataOut.channelList:
                                 raise ValueError, "Channel %d is not in dataOut.channelList"
                             channelIndexList.append(dataOut.channelList.index(channel))
             #        x = dataOut.heightList
                     c   =  3E8
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     #deberia cambiar para el caso de 1Mhz y 100KHz
                     x = numpy.arange(-1*dataOut.nHeights/2.,dataOut.nHeights/2.)*(c/(2*deltaHeight*dataOut.nHeights*1000))
                     data = dataOut.data_spc / factor
                     datadB = 10.*numpy.log10(data)
                     y = datadB
                     #thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
                     title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
                     #para 1Mhz descomentar la siguiente linea
                     #xlabel = "Frequency x 10000"
                     ylabel = "Intensity (dB)"
                     if not self.isConfig:
                         nplots = len(channelIndexList)
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    show=show)
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.isConfig = True
                     self.setWinTitle(title)
                     for i in range(len(self.axesList)):
                         ychannel = y[i,:]
                         str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
                         axes.pline(x, ychannel,
                                     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
                                     xlabel=xlabel, ylabel=ylabel, title=title, grid='both')
                     self.draw()
                     self.save(figpath=figpath,
                               figfile=figfile,
                               save=save,
                               thisDatetime=thisDatetime)
             class RTIfromSpectraHeis(Figure):
                 isConfig = None
                 __nsubplots = None
                 PREFIX = 'rtinoise'
-                def __init__(self, **kwargs):
+                def __init__(self):
-                    Figure.__init__(self, **kwargs)
                     self.timerange = 24*60*60
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.WIDTH = 820
                     self.HEIGHT = 200
                     self.WIDTHPROF = 120
                     self.xdata = None
                     self.ydata = None
                     self.figfile = None
                     self.PLOT_CODE = RTI_CODE
                 def getSubplots(self):
                     ncol = 1
                     nrow = 1
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 7
                     colspan = 6
                     self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH+self.WIDTHPROF,
                                       heightplot = self.HEIGHT+self.HEIGHTPROF,
                                       show = show)
                     nrow, ncol = self.getSubplots()
                     self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
                 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
                         xmin=None, xmax=None, ymin=None, ymax=None,
                         timerange=None,
                         save=False, figpath='./', 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):
                     if channelList == None:
                         channelIndexList = dataOut.channelIndexList
                         channelList = dataOut.channelList
                             if channel not in dataOut.channelList:
                                 raise ValueError, "Channel %d is not in dataOut.channelList"
                             channelIndexList.append(dataOut.channelList.index(channel))
                     if timerange != None:
                         self.timerange = timerange
                     x = dataOut.getTimeRange()
                     y = dataOut.getHeiRange()
                     factor = dataOut.normFactor
                     data = dataOut.data_spc / factor
                     data = numpy.average(data,axis=1)
                     datadB = 10*numpy.log10(data)
             #        factor = dataOut.normFactor
             #        noise = dataOut.getNoise()/factor
             #        noisedB = 10*numpy.log10(noise)
                     #thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
                     title = wintitle + " RTI: %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     xlabel = "Local Time"
                     ylabel = "Intensity (dB)"
                     if not self.isConfig:
                         nplots = 1
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         self.tmin, self.tmax = self.getTimeLim(x, xmin, xmax)
                         if ymin == None: ymin = numpy.nanmin(datadB)
                         if ymax == None: ymax = numpy.nanmax(datadB)
                         self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.isConfig = True
                         self.figfile = figfile
                         self.xdata = numpy.array([])
                         self.ydata = numpy.array([])
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                     self.setWinTitle(title)
             #        title = "RTI %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     title = "RTI - %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
                     legendlabels = ["channel %d"%idchannel for idchannel in channelList]
                     axes = self.axesList[0]
                     self.xdata = numpy.hstack((self.xdata, x[0:1]))
                     if len(self.ydata)==0:
                         self.ydata = datadB[channelIndexList].reshape(-1,1)
                     else:
                         self.ydata = numpy.hstack((self.ydata, datadB[channelIndexList].reshape(-1,1)))
                     axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
                                 xmin=self.tmin, xmax=self.tmax, ymin=ymin, ymax=ymax,
                                 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='.', markersize=8, linestyle="solid", grid='both',
                                 XAxisAsTime=True
                                 )
                     self.draw()
                     update_figfile = False
                     if dataOut.ltctime >= self.tmax:
                         self.counter_imagwr = wr_period
                         self.isConfig = False
                         update_figfile = True
                     self.save(figpath=figpath,
                               figfile=figfile,
                               save=save,

schainpy/model/graphics/jroplot_parameters.py +222 -13

@@ -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,10 +657,9 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)	661	y = dataOut.heightList
451	y = dataOut.heightList	662	z = dataOut.data_output.copy()
452	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
454	nplotsw = nplots	664	nplotsw = nplots
455		665
@@ -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:
@@ -657,12 +867,11 class ParametersPlot(Figure):
657	zmin : None,	867	zmin : None,
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"	874
665
666	if not isTimeInHourRange(dataOut.datatime, xmin, xmax):	875	if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
667	return	876	return
668		877

schainpy/model/graphics/jroplot_spectra.py +26 -19

                         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))
-                    factor = dataOut.normFactor
+                    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]
-                    factor = dataOut.normFactor
+                    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 +25 -12

                                                     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'):
                 ax.collections.remove(ax.collections[0])
                 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 +200 -172

             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:
+                    Affected:
                         self.dataOut
                     Return      : None
                     """
                     #Eliminar de la base la herencia
                     ProcessingUnit.__init__(self, **kwargs)
             #         self.isConfig = False
                     self.pts2read_SelfSpectra = 0
                     self.pts2read_CrossSpectra = 0
                     self.pts2read_DCchannels = 0
                     self.datablock = None
                     self.utc = None
                     self.ext = ".pdata"
                     self.optchar = "P"
                     self.basicHeaderObj = BasicHeader(LOCALTIME)
                     self.systemHeaderObj = SystemHeader()
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.processingHeaderObj = ProcessingHeader()
                     self.online = 0
                     self.fp = None
                     self.idFile = None
                     self.dtype = None
                     self.fileSizeByHeader = None
                     self.filenameList = []
                     self.filename = None
                     self.fileSize = None
                     self.firstHeaderSize = 0
                     self.basicHeaderSize = 24
                     self.pathList = []
                     self.lastUTTime = 0
                     self.maxTimeStep = 30
                     self.flagNoMoreFiles = 0
                     self.set = 0
                     self.path = None
                     self.delay  = 60   #seconds
                     self.nTries = 3  #quantity tries
                     self.nFiles = 3   #number of files for searching
                     self.nReadBlocks = 0
                     self.flagIsNewFile = 1
                     self.__isFirstTimeOnline = 1
             #         self.ippSeconds = 0
                     self.flagDiscontinuousBlock = 0
                     self.flagIsNewBlock = 0
                     self.nTotalBlocks = 0
                     self.blocksize = 0
                     self.dataOut = self.createObjByDefault()
                     self.profileIndex = 1 #Always
                 def createObjByDefault(self):
                     dataObj = Spectra()
                     return dataObj
                 def __hasNotDataInBuffer(self):
                     return 1
                 def getBlockDimension(self):
                     """
                     Obtiene la cantidad de puntos a leer por cada bloque de datos
                     Affected:
                         self.nRdChannels
                         self.nRdPairs
                     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]))
                     pts2read = self.processingHeaderObj.nHeights * self.processingHeaderObj.profilesPerBlock
                     self.pts2read_SelfSpectra = int(self.nRdChannels * pts2read)
                     self.blocksize = self.pts2read_SelfSpectra
                     if self.processingHeaderObj.flag_cspc:
                         self.pts2read_CrossSpectra = int(self.nRdPairs * pts2read)
                         self.blocksize += self.pts2read_CrossSpectra
                     if self.processingHeaderObj.flag_dc:
                         self.pts2read_DCchannels = int(self.systemHeaderObj.nChannels * self.processingHeaderObj.nHeights)
                         self.blocksize += self.pts2read_DCchannels
             #        self.blocksize = self.pts2read_SelfSpectra + self.pts2read_CrossSpectra + self.pts2read_DCchannels
                 def readBlock(self):
                     """
                     Lee el bloque de datos desde la posicion actual del puntero del archivo
                     (self.fp) y actualiza todos los parametros relacionados al bloque de datos
                     (metadata + data). La data leida es almacenada en el buffer y el contador del buffer
                     es seteado a 0
                     Return: None
                     Variables afectadas:
                         self.flagIsNewFile
                         self.flagIsNewBlock
                         self.data_cspc
                         self.data_dc
                     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()
                     if self.processingHeaderObj.flag_cspc:
                         cspc = numpy.fromfile( self.fp, self.dtype, self.pts2read_CrossSpectra )
                         cspc = cspc.reshape( (self.nRdPairs, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D
                     if self.processingHeaderObj.flag_dc:
                         dc = numpy.fromfile( self.fp, self.dtype, self.pts2read_DCchannels ) #int(self.processingHeaderObj.nHeights*self.systemHeaderObj.nChannels) )
                         dc = dc.reshape( (self.systemHeaderObj.nChannels, self.processingHeaderObj.nHeights) ) #transforma a un arreglo 2D
                     if not(self.processingHeaderObj.shif_fft):
                         #desplaza a la derecha en el eje 2 determinadas posiciones
                         shift = int(self.processingHeaderObj.profilesPerBlock/2)
                         spc = numpy.roll( spc, shift , axis=2 )
                         if self.processingHeaderObj.flag_cspc:
                             #desplaza a la derecha en el eje 2 determinadas posiciones
                             cspc = numpy.roll( cspc, shift, axis=2 )
                     #Dimensions : nChannels, nProfiles, nSamples
                     spc = numpy.transpose( spc, (0,2,1) )
                     self.data_spc = spc
+                    if self.processingHeaderObj.flag_cspc:
-                    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
                     self.nReadBlocks += 1
                     return 1
                 def getFirstHeader(self):
                     self.getBasicHeader()
                     self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
                     self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
             #         self.dataOut.ippSeconds = self.ippSeconds
             #         self.dataOut.timeInterval = self.radarControllerHeaderObj.ippSeconds * self.processingHeaderObj.nCohInt * self.processingHeaderObj.nIncohInt * self.processingHeaderObj.profilesPerBlock
                     self.dataOut.dtype = self.dtype
             #         self.dataOut.nPairs = self.nPairs
                     self.dataOut.pairsList = self.rdPairList
                     self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
                     self.dataOut.nFFTPoints = self.processingHeaderObj.profilesPerBlock
                     self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
                     self.dataOut.nIncohInt = self.processingHeaderObj.nIncohInt
                     xf = self.processingHeaderObj.firstHeight + self.processingHeaderObj.nHeights*self.processingHeaderObj.deltaHeight
                     self.dataOut.heightList = numpy.arange(self.processingHeaderObj.firstHeight, xf, self.processingHeaderObj.deltaHeight)
                     self.dataOut.channelList = range(self.systemHeaderObj.nChannels)
                     self.dataOut.flagShiftFFT = True    #Data is always shifted
                     self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode #asumo q la data no esta decodificada
                     self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip #asumo q la data esta sin flip
                 def getData(self):
                     """
                     First method to execute before "RUN" is called.
                     Copia el buffer de lectura a la clase "Spectra",
                     con todos los parametros asociados a este (metadata). cuando no hay datos en el buffer de
                     lectura es necesario hacer una nueva lectura de los bloques de datos usando "readNextBlock"
                     Return:
 :    Si no hay mas archivos disponibles
 :    Si hizo una buena copia del buffer
                     Affected:
                         self.dataOut
                         self.flagDiscontinuousBlock
                         self.flagIsNewBlock
                     """
                         self.dataOut.flagNoData = True
                         print 'Process finished'
                         return 0
                     self.flagDiscontinuousBlock = 0
                     self.flagIsNewBlock = 0
                     if self.__hasNotDataInBuffer():
                         if not( self.readNextBlock() ):
                             self.dataOut.flagNoData = True
                             return 0
                     #data es un numpy array de 3 dmensiones (perfiles, alturas y canales)
                     if self.data_spc is None:
                         self.dataOut.flagNoData = True
                         return 0
                     self.getBasicHeader()
                     self.getFirstHeader()
                     self.dataOut.data_spc = self.data_spc
                     self.dataOut.data_cspc = self.data_cspc
                     self.dataOut.data_dc = self.data_dc
                     self.dataOut.flagNoData = False
                     self.dataOut.realtime = self.online
                     return self.dataOut.data_spc
             class SpectraWriter(JRODataWriter, Operation):
                 """
                 Esta clase permite escribir datos de espectros a archivos procesados (.pdata). La escritura
                 de los datos siempre se realiza por bloques.
                 """
                 ext = ".pdata"
                 optchar = "P"
                 shape_spc_Buffer = None
                 shape_cspc_Buffer = None
                 shape_dc_Buffer = None
                 data_spc = None
                 data_cspc = None
                 data_dc = None
             #    dataOut = None
-                def __init__(self, **kwargs):
+                def __init__(self):
                     """
                     Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
+                    Affected:
-                    Affected:
                         self.dataOut
                         self.basicHeaderObj
                         self.systemHeaderObj
                     Return: None
                     """
-                    Operation.__init__(self, **kwargs)
+                    Operation.__init__(self)
                     self.isConfig = False
                     self.nTotalBlocks = 0
                     self.data_spc = None
                     self.data_cspc = None
                     self.data_dc = None
                     self.fp = None
                     self.flagIsNewFile = 1
                     self.nTotalBlocks = 0
                     self.flagIsNewBlock = 0
                     self.setFile = None
                     self.dtype = None
                     self.path = None
                     self.noMoreFiles = 0
                     self.filename = None
                     self.basicHeaderObj = BasicHeader(LOCALTIME)
                     self.systemHeaderObj = SystemHeader()
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.processingHeaderObj = ProcessingHeader()
                 def hasAllDataInBuffer(self):
                     return 1
                 def setBlockDimension(self):
                     """
                     self.shape_cspc_Buffer = (self.dataOut.nPairs,
                                               self.processingHeaderObj.nHeights,
                                               self.processingHeaderObj.profilesPerBlock)
                     self.shape_dc_Buffer = (self.dataOut.nChannels,
                                             self.processingHeaderObj.nHeights)
                 def writeBlock(self):
                     """
                     Escribe el buffer en el file designado
                     Affected:
                         self.data_spc
                         self.flagIsNewFile
                         self.flagIsNewBlock
                         self.nTotalBlocks
                         self.nWriteBlocks
                     Return: None
                     """
                     spc = numpy.transpose( self.data_spc, (0,2,1) )
                     if not( self.processingHeaderObj.shif_fft ):
                         spc = numpy.roll( spc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
                         data['imag'] = cspc.imag
                         data = data.reshape((-1))
                         data.tofile(self.fp)
                     if self.data_dc is not None:
                         data = numpy.zeros( self.shape_dc_Buffer, self.dtype )
                         data.tofile(self.fp)
             #         self.data_spc.fill(0)
+            #
             #         if self.data_dc is not None:
             #             self.data_dc.fill(0)
+            #
             #         if self.data_cspc is not None:
             #             self.data_cspc.fill(0)
                     self.flagIsNewFile = 0
                     self.flagIsNewBlock = 1
                     self.nTotalBlocks += 1
                     self.nWriteBlocks += 1
                     self.blockIndex += 1
             #         print "[Writing] Block = %d04" %self.blockIndex
                 def putData(self):
                     """
                     Setea un bloque de datos y luego los escribe en un file
                     Affected:
                         self.data_spc
                         self.data_cspc
                         self.data_dc
                     Return:
 :    Si no hay data o no hay mas files que puedan escribirse
 :    Si se escribio la data de un bloque en un file
                     """
                     if self.dataOut.flagNoData:
                         return 0
                     self.flagIsNewBlock = 0
                     if self.dataOut.flagDiscontinuousBlock:
                         self.data_spc.fill(0)
-                        if self.dataOut.data_cspc is not None:
+                        self.data_cspc.fill(0)
-                            self.data_cspc.fill(0)
+                        self.data_dc.fill(0)
-                        if self.dataOut.data_dc is not None:
-                            self.data_dc.fill(0)
                         self.setNextFile()
                     if self.flagIsNewFile == 0:
                         self.setBasicHeader()
                     self.data_spc = self.dataOut.data_spc.copy()
                     if self.dataOut.data_cspc is not None:
                         self.data_cspc = self.dataOut.data_cspc.copy()
                     if self.dataOut.data_dc is not None:
                         self.data_dc = self.dataOut.data_dc.copy()
                     # #self.processingHeaderObj.dataBlocksPerFile)
                     if self.hasAllDataInBuffer():
             #            self.setFirstHeader()
                         self.writeNextBlock()
                     return 1
                 def __getBlockSize(self):
                     '''
                     Este metodos determina el cantidad de bytes para un bloque de datos de tipo Spectra
                     '''
                     dtype_width = self.getDtypeWidth()
                     pts2write = self.dataOut.nHeights * self.dataOut.nFFTPoints
                     pts2write_SelfSpectra = int(self.dataOut.nChannels * pts2write)
                     blocksize = (pts2write_SelfSpectra*dtype_width)
                     if self.dataOut.data_cspc is not None:
                         pts2write_CrossSpectra = int(self.dataOut.nPairs * pts2write)
                         blocksize += (pts2write_CrossSpectra*dtype_width*2)
                     if self.dataOut.data_dc is not None:
                         pts2write_DCchannels = int(self.dataOut.nChannels * self.dataOut.nHeights)
                         blocksize += (pts2write_DCchannels*dtype_width*2)
             #         blocksize = blocksize #* datatypeValue * 2 #CORREGIR ESTO
                     return blocksize
                 def setFirstHeader(self):
                     """
                     Obtiene una copia del First Header
                     Affected:
                         self.systemHeaderObj
                         self.radarControllerHeaderObj
                         self.dtype
                     Return:
                         None
                     """
                     self.systemHeaderObj = self.dataOut.systemHeaderObj.copy()
                     self.systemHeaderObj.nChannels = self.dataOut.nChannels
                     self.radarControllerHeaderObj = self.dataOut.radarControllerHeaderObj.copy()
                     self.processingHeaderObj.dtype = 1 # Spectra
                     self.processingHeaderObj.blockSize = self.__getBlockSize()
                     self.processingHeaderObj.profilesPerBlock = self.dataOut.nFFTPoints
                     self.processingHeaderObj.dataBlocksPerFile = self.blocksPerFile
                     self.processingHeaderObj.nWindows = 1 #podria ser 1 o self.dataOut.processingHeaderObj.nWindows
                     self.processingHeaderObj.nCohInt = self.dataOut.nCohInt# Se requiere para determinar el valor de timeInterval
                     self.processingHeaderObj.nIncohInt = self.dataOut.nIncohInt
                     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):
                             channelList.append(channel)
                             channelList.append(channel)
                         pairsList = []
                         if self.dataOut.nPairs > 0:
                             for pair in self.dataOut.pairsList:
                                 pairsList.append(pair[0])
                                 pairsList.append(pair[1])
                         spectraComb = channelList + pairsList
                         spectraComb = numpy.array(spectraComb, dtype="u1")
                         self.processingHeaderObj.spectraComb = spectraComb
                     if self.dataOut.code is not None:
                         self.processingHeaderObj.code = self.dataOut.code
                         self.processingHeaderObj.nCode = self.dataOut.nCode
                         self.processingHeaderObj.nBaud = self.dataOut.nBaud
                     if self.processingHeaderObj.nWindows != 0:
                         self.processingHeaderObj.firstHeight = self.dataOut.heightList[0]
                         self.processingHeaderObj.deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
                         self.processingHeaderObj.nHeights = self.dataOut.nHeights
                         self.processingHeaderObj.samplesWin = self.dataOut.nHeights
                     self.processingHeaderObj.processFlags = self.getProcessFlags()
                     self.setBasicHeader()

schainpy/model/io/jroIO_voltage.py +7 -1

                     self.dataOut.flagNoData = False
                     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):
                 """

schainpy/model/proc/__init__.py +2 -1

             from jroproc_correlation import *
             from jroproc_parameters import *
             from jroproc_spectra_lags import *
-            from jroproc_spectra_acf import *
  No newline at end of file
+            from jroproc_spectra_acf import *
+            from bltrproc_parameters import *

schainpy/model/proc/jroproc_parameters.py 0 0

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

schainpy/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