schain-jc Commit - r1021:2aa9b10a9a70

Operation MAD2Writer done Task #343

Juan C. Espinoza -

r1021:2aa9b10a9a70

parent child

Context file:

r1021:2aa9b10a9a70

Expand all files

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

             '''
             $Author: murco $
             $Id: JRODataIO.py 169 2012-11-19 21:57:03Z murco $
             '''
             from jroIO_voltage import *
             from jroIO_spectra import *
             from jroIO_heispectra import *
             from jroIO_usrp import *
             from jroIO_kamisr import *
             from jroIO_param import *
             from jroIO_hf import *
+            from jroIO_madrigal import *
             from bltrIO_param import *
             from jroIO_bltr import *
             from jroIO_mira35c import *

schainpy/model/io/bltrIO_param.py +1 -3

             '''
             Created on Nov 9, 2016
             @author: roj- LouVD
             '''
             import os
             import sys
             import time
             import glob
             import datetime
             import numpy
             from schainpy.model.proc.jroproc_base import ProcessingUnit
             from schainpy.model.data.jrodata import Parameters
             from schainpy.model.io.jroIO_base import JRODataReader, isNumber
             FILE_HEADER_STRUCTURE = numpy.dtype([
                 ('FMN', '<u4'),
                 ('nrec', '<u4'),
                 ('fr_offset', '<u4'),
                 ('id', '<u4'),
                 ('site', 'u1', (32,))
             ])
             REC_HEADER_STRUCTURE = numpy.dtype([
                 ('rmn', '<u4'),
                 ('rcounter', '<u4'),
                 ('nr_offset', '<u4'),
                 ('tr_offset', '<u4'),
                 ('time', '<u4'),
                 ('time_msec', '<u4'),
                 ('tag', 'u1', (32,)),
                 ('comments', 'u1', (32,)),
                 ('lat', '<f4'),
                 ('lon', '<f4'),
                 ('gps_status', '<u4'),
                 ('freq', '<u4'),
                 ('freq0', '<u4'),
                 ('nchan', '<u4'),
                 ('delta_r', '<u4'),
                 ('nranges', '<u4'),
                 ('r0', '<u4'),
                 ('prf', '<u4'),
                 ('ncoh', '<u4'),
                 ('npoints', '<u4'),
                 ('polarization', '<i4'),
                 ('rx_filter', '<u4'),
                 ('nmodes', '<u4'),
                 ('dmode_index', '<u4'),
                 ('dmode_rngcorr', '<u4'),
                 ('nrxs', '<u4'),
                 ('acf_length', '<u4'),
                 ('acf_lags', '<u4'),
                 ('sea_to_atmos', '<f4'),
                 ('sea_notch', '<u4'),
                 ('lh_sea', '<u4'),
                 ('hh_sea', '<u4'),
                 ('nbins_sea', '<u4'),
                 ('min_snr', '<f4'),
                 ('min_cc', '<f4'),
                 ('max_time_diff', '<f4')
             ])
             DATA_STRUCTURE = numpy.dtype([
                 ('range', '<u4'),
                 ('status', '<u4'),
                 ('zonal', '<f4'),
                 ('meridional', '<f4'),
                 ('vertical', '<f4'),
                 ('zonal_a', '<f4'),
                 ('meridional_a', '<f4'),
                 ('corrected_fading', '<f4'),  # seconds
                 ('uncorrected_fading', '<f4'),  # seconds
                 ('time_diff', '<f4'),
                 ('major_axis', '<f4'),
                 ('axial_ratio', '<f4'),
                 ('orientation', '<f4'),
                 ('sea_power', '<u4'),
                 ('sea_algorithm', '<u4')
             ])
             class BLTRParamReader(JRODataReader, ProcessingUnit):
                 '''
                 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR from *.sswma files
                 '''
                 ext = '.sswma'
                 def __init__(self, **kwargs):
                     ProcessingUnit.__init__(self , **kwargs)
                     self.dataOut = Parameters()
                     self.counter_records = 0
                     self.flagNoMoreFiles = 0
                     self.isConfig = False
                     self.filename = None
                 def setup(self,
                           path=None,
                           startDate=None,
                           endDate=None,
                           ext=None,
                           startTime=datetime.time(0, 0, 0),
                           endTime=datetime.time(23, 59, 59),
                           timezone=0,
                           status_value=0,
                           **kwargs):
                     self.path = path
                     self.startTime = startTime
                     self.endTime = endTime
                     self.status_value = status_value
                     if self.path is None:
                         raise ValueError, "The path is not valid"
                     if ext is None:
                         ext = self.ext
                     self.search_files(self.path, startDate, endDate, ext)
                     self.timezone = timezone
                     self.fileIndex = 0
                     if not self.fileList:
                         raise  Warning, "There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' "%(path)
                     self.setNextFile()
                 def search_files(self, path, startDate, endDate, ext):
                     '''
                      Searching for BLTR rawdata file in path
                      Creating a list of file to proces included in [startDate,endDate]
                      Input:
                          path - Path to find BLTR rawdata files
                          startDate - Select file from this date
                          enDate - Select file until this date
                          ext - Extension of the file to read
                     '''
                     print 'Searching file in %s ' % (path)
                     foldercounter = 0
                     fileList0 = glob.glob1(path, "*%s" % ext)
                     fileList0.sort()
                     self.fileList = []
                     self.dateFileList = []
                     for thisFile in fileList0:
                         year = thisFile[-14:-10]
                         if not isNumber(year):
                             continue
                         month = thisFile[-10:-8]
                         if not isNumber(month):
                             continue
                         day = thisFile[-8:-6]
                         if not isNumber(day):
                             continue
                         year, month, day = int(year), int(month), int(day)
                         dateFile = datetime.date(year, month, day)
                         if (startDate > dateFile) or (endDate < dateFile):
                             continue
                         self.fileList.append(thisFile)
                         self.dateFileList.append(dateFile)
                     return
                 def setNextFile(self):
                     file_id = self.fileIndex
                     if file_id == len(self.fileList):
                         print '\nNo more files in the folder'
                         print 'Total number of file(s) read : {}'.format(self.fileIndex + 1)
                         self.flagNoMoreFiles = 1
                         return 0
                     print '\n[Setting file] (%s) ...' % self.fileList[file_id]
                     filename = os.path.join(self.path, self.fileList[file_id])
                     dirname, name = os.path.split(filename)
                     self.siteFile = name.split('.')[0]  # 'peru2' ---> Piura  -   'peru1' ---> Huancayo or Porcuya
                     if self.filename is not None:
                         self.fp.close()
                     self.filename = filename
                     self.fp = open(self.filename, 'rb')
                     self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
                     self.nrecords = self.header_file['nrec'][0]
                     self.sizeOfFile = os.path.getsize(self.filename)
                     self.counter_records = 0
                     self.flagIsNewFile = 0
                     self.fileIndex += 1
                     return 1
                 def readNextBlock(self):
                     while True:
                         if self.counter_records == self.nrecords:
                             self.flagIsNewFile = 1
                             if not self.setNextFile():
                                 return 0
                         self.readBlock()
                         if (self.datatime.time() < self.startTime) or (self.datatime.time() > self.endTime):
                             print "[Reading] Record No. %d/%d -> %s [Skipping]" %(
                                 self.counter_records,
                                 self.nrecords,
                                 self.datatime.ctime())
                             continue
                         break
                     print "[Reading] Record No. %d/%d -> %s" %(
                         self.counter_records,
                         self.nrecords,
                         self.datatime.ctime())
                     return 1
                 def readBlock(self):
                     pointer = self.fp.tell()
                     header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
                     self.nchannels = header_rec['nchan'][0]/2
                     self.kchan = header_rec['nrxs'][0]
                     self.nmodes = header_rec['nmodes'][0]
                     self.nranges = header_rec['nranges'][0]
                     self.fp.seek(pointer)
                     self.height = numpy.empty((self.nmodes, self.nranges))
                     self.snr = numpy.empty((self.nmodes, self.nchannels, self.nranges))
                     self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
                     for mode in range(self.nmodes):
                         self.readHeader()
                         data = self.readData()
                         self.height[mode] = (data[0] - self.correction) / 1000.
                         self.buffer[mode] = data[1]
                         self.snr[mode] = data[2]
                     self.counter_records = self.counter_records + self.nmodes
                     return
                 def readHeader(self):
                     '''
                     RecordHeader of BLTR rawdata file
                     '''
                     header_structure = numpy.dtype(
                         REC_HEADER_STRUCTURE.descr + [
                             ('antenna_coord', 'f4', (2, self.nchannels)),
                             ('rx_gains', 'u4', (self.nchannels,)),
                             ('rx_analysis', 'u4', (self.nchannels,))
                         ]
                     )
                     self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
                     self.lat = self.header_rec['lat'][0]
                     self.lon = self.header_rec['lon'][0]
                     self.delta = self.header_rec['delta_r'][0]
                     self.correction = self.header_rec['dmode_rngcorr'][0]
                     self.imode = self.header_rec['dmode_index'][0]
                     self.antenna = self.header_rec['antenna_coord']
                     self.rx_gains = self.header_rec['rx_gains']
                     self.time = self.header_rec['time'][0]
                     tseconds = self.header_rec['time'][0]
                     local_t1 = time.localtime(tseconds)
                     self.year = local_t1.tm_year
                     self.month = local_t1.tm_mon
                     self.day = local_t1.tm_mday
                     self.t = datetime.datetime(self.year, self.month, self.day)
                     self.datatime = datetime.datetime.utcfromtimestamp(self.time)
                 def readData(self):
                     '''
                     Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
                     Input:
                         status_value - Array data is set to NAN for values that are not equal to status_value
                     '''
                     data_structure = numpy.dtype(
                         DATA_STRUCTURE.descr + [
                             ('rx_saturation', 'u4', (self.nchannels,)),
                             ('chan_offset', 'u4', (2 * self.nchannels,)),
                             ('rx_amp', 'u4', (self.nchannels,)),
                             ('rx_snr', 'f4', (self.nchannels,)),
                             ('cross_snr', 'f4', (self.kchan,)),
                             ('sea_power_relative', 'f4', (self.kchan,))]
                     )
                     data = numpy.fromfile(self.fp, data_structure, self.nranges)
                     height = data['range']
                     winds = numpy.array((data['zonal'], data['meridional'], data['vertical']))
                     snr = data['rx_snr'].T
                     winds[numpy.where(winds == -9999.)] = numpy.nan
                     winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
                     snr[numpy.where(snr == -9999.)] = numpy.nan
                     snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
                     snr = numpy.power(10, snr / 10)
                     return height, winds, snr
                 def set_output(self):
                     '''
                     Storing data from databuffer to dataOut object
                     '''
                     self.dataOut.data_SNR = self.snr
                     self.dataOut.height = self.height
                     self.dataOut.data_output = self.buffer
                     self.dataOut.utctimeInit = self.time
                     self.dataOut.utctime = self.dataOut.utctimeInit
-                    self.dataOut.counter_records = self.counter_records
-                    self.dataOut.nrecords = self.nrecords
                     self.dataOut.useLocalTime = False
                     self.dataOut.paramInterval = 157
                     self.dataOut.timezone = self.timezone
                     self.dataOut.site = self.siteFile
-                    self.dataOut.nrecords = self.nrecords
+                    self.dataOut.nrecords = self.nrecords/self.nmodes
                     self.dataOut.sizeOfFile = self.sizeOfFile
                     self.dataOut.lat = self.lat
                     self.dataOut.lon = self.lon
                     self.dataOut.channelList = range(self.nchannels)
                     self.dataOut.kchan = self.kchan
                     # self.dataOut.nHeights = self.nranges
                     self.dataOut.delta = self.delta
                     self.dataOut.correction = self.correction
                     self.dataOut.nmodes = self.nmodes
                     self.dataOut.imode = self.imode
                     self.dataOut.antenna = self.antenna
                     self.dataOut.rx_gains = self.rx_gains
                     self.dataOut.flagNoData = False
                 def getData(self):
                     '''
                     Storing data from databuffer to dataOut object
                     '''
                     if self.flagNoMoreFiles:
                         self.dataOut.flagNoData = True
                         print 'No file left to process'
                         return 0
                     if not  self.readNextBlock():
                         self.dataOut.flagNoData = True
                         return 0
                     self.set_output()
                     return 1

schainpy/model/io/jroIO_madrigal.py +168 -300

@@ -1,375 +1,243
1	'''	1	'''
2	Created on Aug 1, 2017	2	Created on Aug 1, 2017
3		3
4	@author: Juan C. Espinoza	4	@author: Juan C. Espinoza
5	'''	5	'''
6		6
7	import os	7	import os
8	import sys	8	import sys
9	import time	9	import time
		10	import json
10	import datetime	11	import datetime
11		12
12	import numpy	13	import numpy
13		14
14	from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
15	from schainpy.model.data.jrodata import Parameters
16	from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
17	from schainpy.model.graphics.jroplot_parameters import WindProfilerPlot
18	from schainpy.model.io.jroIO_base import *
19
20	try:	15	try:
21	import madrigal	16	import madrigal
22	import madrigal.cedar	17	import madrigal.cedar
23	from madrigal.cedar import MadrigalCatalogRecord
24	except:	18	except:
25	print 'You should install "madrigal library" module if you want to read/write Madrigal data'	19	print 'You should install "madrigal library" module if you want to read/write Madrigal data'
26		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
27		61
28	class MADWriter(Operation):	62	class MAD2Writer(Operation):
29		63
30	def __init__(self):	64	def __init__(self, **kwargs):
31		65
32	Operation.__init__(self)	66	Operation.__init__(self, **kwargs)
33	self.dataOut = Parameters()	67	self.dataOut = Parameters()
34	self.path = None	68	self.path = None
35	self.dataOut = None	69	self.dataOut = None
36	self.flagIsNewFile=1	70	self.ext = '.dat'
37	self.ext = ".hdf5"
38		71
39	return	72	return
40		73
41	def run(self, dataOut, path , ~~modetowrite~~,**kwargs):	74	def run(self, dataOut, path, oneDList, twoDParam='', twoDList='{}', metadata='{}', **kwargs):
42		75	'''
43	if self.flagIsNewFile:	76	Inputs:
44	flagdata = self.setup(dataOut, path, modetowrite)	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
45		102
46	self.putData()	103	self.putData()
47	return	104	return
48		105
49	def setup(self, dataOut, path, ~~modetowrite~~):	106	def setup(self, dataOut, path, oneDList, twoDParam, twoDList, metadata, **kwargs):
50	'''	107	'''
51	Recovering data to write in new *.hdf5 file	108	Configure Operation
52	Inputs:
53	modew -- mode to write (1 or 2)
54	path -- destination path
55
56	'''	109	'''
57		110
58	self.im = modetowrite-1
59	if self.im!=0 and self.im!=1:
60	raise ValueError, 'Check "modetowrite" value. Must be egual to 1 or 2, "{}" is not valid. '.format(modetowrite)
61
62	self.dataOut = dataOut	111	self.dataOut = dataOut
63	self.nmodes = self.dataOut.nmodes	112	self.nmodes = self.dataOut.nmodes
64	self.nchannels = self.dataOut.nchannels
65	self.lat = self.dataOut.lat
66	self.lon = self.dataOut.lon
67	self.hcm = 3
68	self.thisDate = self.dataOut.utctimeInit
69	self.year = self.dataOut.year
70	self.month = self.dataOut.month
71	self.day = self.dataOut.day
72	self.path = path	113	self.path = path
73		114	self.blocks = kwargs.get('blocks', None)
74	self.flagIsNewFile = 0	115	self.counter = 0
75		116	self.oneDList = load_json(oneDList)
76	return 1	117	self.twoDList = load_json(twoDList)
77		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
78	def setFile(self):	127	def setFile(self):
79	'''	128	'''
80	- Determining the file name for each mode of operation	129	Create new cedar file object
81	kinst - Kind of Instrument (mnemotic)	130	'''
82	kindat - Kind of Data (mnemotic)	131
83		132	self.mnemonic = MNEMONICS[self.kinst] #TODO get mnemonic from madrigal
84	- Creating a cedarObject	133	date = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
85		134
86	'''	135	filename = '%s%s_%s%s' % (self.mnemonic,
87	lat_piura = -5.17	136	date.strftime('%Y%m%d_%H%M%S'),
88	lat_huancayo = -12.04	137	self.dataOut.mode,
89	lat_porcuya = -5.8	138	self.ext)
90
91	if '%2.2f' % self.lat == '%2.2f' % lat_piura:
92	self.instMnemonic = 'pbr'
93
94	elif '%2.2f' % self.lat == '%2.2f' % lat_huancayo:
95	self.instMnemonic = 'hbr'
96
97	elif '%2.2f' % self.lat == '%2.2f' % lat_porcuya:
98	self.instMnemonic = 'obr'
99	else: raise Warning, "The site of file read doesn't match any site known. Only file from Huancayo, Piura and Porcuya can be processed.\n Check the file "
100
101	mode = ['_mode1','_mode2']
102
103	self.hdf5filename = '%s%4.4d%2.2d%2.2d%s%s' % (self.instMnemonic,
104	self.year,
105	self.month,
106	self.day,
107	mode[self.im],
108	self.ext)
109		139
110	self.fullname=os.path.join(self.path,~~self~~.~~hdf5~~filename)	140	self.fullname = os.path.join(self.path, filename)
111		141
112	if os.path.isfile(self.fullname) :	142	if os.path.isfile(self.fullname) :
113	print "Destination path '%s' already exists. Previous file deleted. " %self.fullname	143	print "Destination path '%s' already exists. Previous file deleted. " %self.fullname
114	os.remove(self.fullname)	144	os.remove(self.fullname)
115
116	# Identify kinst and kindat
117	InstName = self.hdf5filename[0:3]
118	KinstList = [1000, 1001, 1002]
119	KinstId = {'pbr':0, 'hbr':1, 'obr':2} # pbr:piura, hbr:huancayo, obr:porcuya
120	KindatList = [1600, 1601] # mode 1, mode 2
121	self.type = KinstId[InstName]
122	self.kinst = KinstList[self.type]
123	self.kindat = KindatList[self.im]
124		145
125	try:	146	try:
		147	print '[Writing] creating file : %s' % (self.fullname)
126	self.cedarObj = madrigal.cedar.MadrigalCedarFile(self.fullname, True)	148	self.cedarObj = madrigal.cedar.MadrigalCedarFile(self.fullname, True)
127	except ValueError, ~~messag~~e:	149	except ValueError, e:
128	print '[Error]: Impossible to create a cedar object with "madrigal.cedar.MadrigalCedarFile" '	150	print '[Error]: Impossible to create a cedar object with "madrigal.cedar.MadrigalCedarFile" '
129	return	151	return
130		152
131	return 1	153	return 1
132		154
133	def writeBlock(self):	155	def writeBlock(self):
134	'''	156	'''
135	- Selecting mode of operation:	157	Add data records to cedar file taking data from oneDList and twoDList
136		158	attributes.
137	bltr high resolution mode 1 - Low Atmosphere (0 - 3km) // bltr high resolution mode 2 - High Atmosphere (0 - 10km)	159	Allowed parameters in: parcodes.tab
138	msnr - Average Signal Noise Ratio in dB
139	hcm - 3 km
140
141	- Filling the cedarObject by a block: each array data entry is assigned a code that defines the parameter to write to the file
142
143	GDLATR - Reference geod latitude (deg)
144	GDLONR - Reference geographic longitude (deg)
145	GDLAT2 - Geodetic latitude of second inst (deg)
146	GLON2 - Geographic longitude of second inst (deg)
147
148	GDALT - Geodetic altitude (height) (km)
149	SNL - Log10 (signal to noise ratio)
150	VN1P2 - Neutral wind in direction 1 (eastward) (m/s), ie zonal wind
151	VN2P2 - Neutral wind in direction 2 (northward) (m/s), ie meridional wind
152	EL2 - Ending elevation angle (deg), ie vertical wind
153
154	Other parameters: /madrigal3/metadata/parcodes.tab
155
156	'''	160	'''
157		161
158	self.z_zon = self.dataOut.data_output[0,:,:]	162	startTime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
159	self.z_mer =self.dataOut.data_output[1,:,:]	163	endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
160	self.z_ver = self.dataOut.data_output[2,:,:]	164	nrows = len(getattr(self.dataOut, self.twoDParam))
161		165
162	if self.im == 0:	166	rec = madrigal.cedar.MadrigalDataRecord(
163	h_select = numpy.where(numpy.bitwise_and(self.dataOut.height[0, :] >= 0., self.dataOut.height[0, :] <= self.hcm, numpy.isfinite(self.dataOut.height[0, :])))	167	self.kinst,
164	else:	168	self.kindat,
165	h_select = numpy.where(numpy.bitwise_and(self.dataOut.height[0, :] >= 0., self.dataOut.height[0, :] < 20, numpy.isfinite(self.dataOut.height[0, :])))	169	startTime.year,
166		170	startTime.month,
167	ht = h_select[0]	171	startTime.day,
168		172	startTime.hour,
169	self.o_height = self.dataOut.height[self.im, ht]	173	startTime.minute,
170	self.o_zon = self.z_zon[ht, self.im]	174	startTime.second,
171	self.o_mer = self.z_mer[ht, self.im]	175	startTime.microsecond/10000,
172	self.o_ver = self.z_ver[ht, self.im]	176	endTime.year,
173	o_snr = self.dataOut.data_SNR[ :, :, self.im]	177	endTime.month,
174		178	endTime.day,
175	o_snr = o_snr[ht, :]	179	endTime.hour,
176		180	endTime.minute,
177	ndiv = numpy.nansum((numpy.isfinite(o_snr)), 1)	181	endTime.second,
178	ndiv = ndiv.astype(float)	182	endTime.microsecond/10000,
179		183	self.oneDList.keys(),
180	sel_div = numpy.where(ndiv == 0.)	184	self.twoDList.keys(),
181	ndiv[sel_div] = numpy.nan	185	nrows
182		186	)
183	if self.nchannels > 1:
184	msnr = numpy.nansum(o_snr, axis=1)
185	else:
186	msnr = o_snr
187
188	try:
189	self.msnr = 10 * numpy.log10(msnr / ndiv)
190	except ZeroDivisionError:
191	self.msnr = 10 * numpy.log10(msnr /1)
192	print 'Number of division (ndiv) egal to 1 by default. Check SNR'
193
194	time_t = time.gmtime(self.dataOut.time1)
195	year = time_t.tm_year
196	month = time_t.tm_mon
197	day = time_t.tm_mday
198	hour = time_t.tm_hour
199	minute = time_t.tm_min
200	second = time_t.tm_sec
201	timedate_0 = datetime.datetime(year, month, day, hour, minute, second)
202
203	# 1d parameters
204	GDLATR = self.lat
205	GDLONR = self.lon
206	GDLAT2 = self.lat
207	GLON2 = self.lon
208
209	# 2d parameters
210	GDALT = self.o_height
211
212	SNL = self.msnr
213	VN1P2 = self.o_zon
214	VN2P2 = self.o_mer
215	EL2 = self.o_ver
216	NROW = len(self.o_height)
217
218	startTime = timedate_0
219	endTime = startTime
220	self.dataRec = madrigal.cedar.MadrigalDataRecord(self.kinst,
221	self.kindat,
222	startTime.year,
223	startTime.month,
224	startTime.day,
225	startTime.hour,
226	startTime.minute,
227	startTime.second,
228	0,
229	endTime.year,
230	endTime.month,
231	endTime.day,
232	endTime.hour,
233	endTime.minute,
234	endTime.second,
235	0,
236	('gdlatr', 'gdlonr', 'gdlat2', 'glon2'),
237	('gdalt', 'snl', 'vn1p2', 'vn2p2', 'el2'),
238	NROW, ind2DList=['gdalt'])
239		187
240	# Setting 1d values	188	# Setting 1d values
241	self.dataRec.set1D('gdlatr', GDLATR)	189	for key in self.oneDList:
242	self.dataRec.set1D('gdlonr', GDLONR)	190	rec.set1D(key, getattr(self.dataOut, self.oneDList[key]))
243	self.dataRec.set1D('gdlat2', GDLAT2)
244	self.dataRec.set1D('glon2', GLON2)
245		191
246	# Setting 2d values	192	# Setting 2d values
247	for n in range(self.o_height.shape[0]):	193	invalid = numpy.isnan(self.dataOut.data_output)
248	self.dataRec.set2D('gdalt', n, GDALT[n])	194	self.dataOut.data_output[invalid] = MISSING
249	self.dataRec.set2D('snl', n, SNL[n])	195	out = {}
250	self.dataRec.set2D('vn1p2', n, VN1P2[n])	196	for key, value in self.twoDList.items():
251	self.dataRec.set2D('vn2p2', n, VN2P2[n])	197	if isinstance(value, str):
252	self.dataRec.set2D('el2', n, EL2[n])	198	out[key] = getattr(self.dataOut, value)
253		199	elif isinstance(value, tuple):
254	# Appending new data record	200	attr, x = value
255	'''	201	if isinstance(x, (int, float)):
256	[MADRIGAL3]There are two ways to write to a MadrigalCedarFile. Either this method (write) is called after all the	202	out[key] = getattr(self.dataOut, attr)[int(x)]
257	records have been appended to the MadrigalCedarFile, or dump is called after a certain number of records are appended,	203	elif x.lower()=='db':
258	and then at the end dump is called a final time if there were any records not yet dumped, followed by addArray.	204	tmp = getattr(self.dataOut, attr)
259	'''	205	SNRavg = numpy.average(tmp, axis=0)
260		206	out[key] = 10*numpy.log10(SNRavg)
261	self.cedarObj.append(self.dataRec)	207
262	print ' [Writing] records {} (mode {}).'.format(self.dataOut.counter_records,self.im+1)	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)
263	self.cedarObj.dump()	213	self.cedarObj.dump()
264		214	print '[Writing] Record No. {} (mode {}).'.format(
265		215	self.counter,
		216	self.dataOut.mode
		217	)
266		218
267
268	def setHeader(self):	219	def setHeader(self):
269	'''	220	'''
270	- Creating self.catHeadObj	221	Create an add catalog and header to cedar file
271	- Adding information catalog
272	- Writing file header
273
274	'''	222	'''
275	self.catHeadObj = madrigal.cedar.CatalogHeaderCreator(self.fullname)
276	kindatDesc, comments, analyst, history, principleInvestigator = self._info_BLTR()
277
278	self.catHeadObj.createCatalog(principleInvestigator="Jarjar",
279	expPurpose='characterize the atmospheric dynamics in this region where frequently it happens the El Nino',
280	sciRemarks="http://madrigal3.haystack.mit.edu/static/CEDARMadrigalHdf5Format.pdf")
281
282	self.catHeadObj.createHeader(kindatDesc, analyst, comments, history)
283
284	self.catHeadObj.write()
285		223
286	print '[File created] path: %s' % (self.fullname)	224	header = madrigal.cedar.CatalogHeaderCreator(self.fullname)
		225	header.createCatalog(**self.catalog)
		226	header.createHeader(**self.header)
		227	header.write()
287		228
288	def putData(self):	229	def putData(self):
289		230
290	if self.dataOut.flagNoData:	231	if self.dataOut.flagNoData:
291	return 0	232	return 0
292		233
293	if self.~~dataOut~~.counter~~_records~~ == 1:	234	if self.counter == 0:
294	self.setFile()	235	self.setFile()
295	print '[Writing] Setting new hdf5 file for the mode {}'.format(self.im+1)
296		236
297	if self.~~dataOut~~.counter~~_records~~ <= self.dataOut.nrecords:	237	if self.counter <= self.dataOut.nrecords:
298	self.writeBlock()	238	self.writeBlock()
299		239	self.counter += 1
300		240
301	if self.~~dataOut~~.counter~~_records~~ == self.dataOut.nrecords:	241	if self.counter == self.dataOut.nrecords or self.counter == self.blocks:
302	self.cedarObj.addArray()
303
304	self.setHeader()	242	self.setHeader()
305	self.~~flagIsNewFile~~ = 1	243	self.counter = 0
306
307	def _info_BLTR(self):
308
309	kindatDesc = '''--This header is for KINDAT = %d''' % self.kindat
310	history = None
311	analyst = '''Jarjar'''
312	principleInvestigator = '''
313	Jarjar
314	Radio Observatorio de Jicamarca
315	Instituto Geofisico del Peru
316
317	'''
318	if self.type == 1:
319	comments = '''
320
321	--These data are provided by two Boundary Layer and Tropospheric Radar (BLTR) deployed at two different locations at Peru(GMT-5), one of them at Piura(5.17 S, 80.64W) and another located at Huancayo (12.04 S, 75.32 W).
322
323	--The purpose of conducting these observations is to measure wind in the differents levels of height, this radar makes measurements the Zonal(U), Meridional(V) and Vertical(W) wind velocities component in northcoast from Peru. And the main purpose of these mensurations is to characterize the atmospheric dynamics in this region where frequently it happens the 'El Nino Phenomenon'
324
325	--In Kindat = 1600, contains information of wind velocities component since 0 Km to 3 Km.
326
327	--In Kindat = 1601, contains information of wind velocities component since 0 Km to 10 Km.
328
329	--The Huancayo-BLTR is a VHF Profiler Radar System is a 3 channel coherent receiver pulsed radar utilising state-of-the-art software and computing techniques to acquire, decode, and translate signals obtained from partial reflection echoes in the troposphere, lower stratosphere and mesosphere. It uses an array of three horizontal spaced and vertically directed receiving antennas. The data is recorded thirty seconds, averaged to one minute mean values of Height, Zonal, Meridional and Vertical wind.
330
331	--The Huancayo-BLTR was installed in January 2010. This instrument was designed and constructed by Genesis Soft Pty. Ltd. Is constituted by three groups of spaced antennas (distributed) forming an isosceles triangle.
332
333
334	Station _______ Geographic Coord ______ Geomagnetic Coord
335
336	_______________ Latitude _ Longitude __ Latitude _ Longitude
337
338	Huancayo (HUA) __12.04 S ___ 75.32 W _____ -12.05 ____ 352.85
339	Piura (PIU) _____ 5.17 S ___ 80.64 W ______ 5.18 ____ 350.93
340
341	WIND OBSERVATIONS
342
343	--To obtain wind the BLTR uses Spaced Antenna technique (e.g., Briggs 1984). The scatter and reflection it still provided by variations in the refractive index as in the Doppler method(Gage and Basley,1978; Balsley and Gage 1982; Larsen and Rottger 1982), but instead of using the Doppler shift to derive the velocity components, the cross-correlation between signals in an array of three horizontally spaced and vertically directed receiving antennas is used.
344
345	......................................................................
346	For more information, consult the following references:
347	- Balsley, B. B., and K. S. Gage., On the use of radars for operational wind profiling, Bull. Amer. Meteor.Soc.,63, 1009-1018, 1982.
348
349	- Briggs, B. H., The analysis of spaced sensor data by correations techniques, Handbook for MAP, Vol. 13, SCOTEP Secretariat, University of Illinois, Urbana, 166-186, 1984.
350
351	- Gage, K. S., and B.B. Balsley., Doppler radar probing of the clear atmosphere, Bull. Amer. Meteor.Soc., 59, 1074-1093, 1978.
352
353	- Larsen, M. F., The Spaced Antenna Technique for Radar Wind Profiling, Journal of Atm. and Ocean. Technology. , Vol.6, 920-937, 1989.
354
355	- Larsen, M. F., A method for single radar voracity measurements?, Handbook for MAP,SCOSTEP Secretariat, University of the Illinois, Urban, in press, 1989.
356	......................................................................
357
358	ACKNOWLEDGEMENTS:
359
360	--The Piura and Huancayo BLTR are part of the network of instruments operated by the Jicamarca Radio Observatory.
361
362	--The Jicamarca Radio Observatory is a facility of the Instituto Geofisico del Peru operated with support from the NSF Cooperative Agreement ATM-0432565 through Cornell University
363
364	......................................................................
365
366	Further questions and comments should be addressed to:
367	Radio Observatorio de Jicamarca
368	Instituto Geofisico del Peru
369	Lima, Peru
370	Web URL: http://jro.igp.gob.pe
371	......................................................................
372	'''
373
374	return kindatDesc, comments, analyst, history, principleInvestigator
375

schainpy/model/proc/bltrproc_parameters.py +11 -1

             '''
             Created on Oct 24, 2016
             @author: roj- LouVD
             '''
             import numpy
             import copy
             import datetime
             import time
             from time import gmtime
             from numpy import transpose
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Parameters
             class BLTRParametersProc(ProcessingUnit):
                 '''
                 Processing unit for BLTR parameters data (winds)
                 Inputs:
                     self.dataOut.nmodes - Number of operation modes
                     self.dataOut.nchannels - Number of channels
                     self.dataOut.nranges - Number of ranges
                     self.dataOut.data_SNR - SNR array
                     self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
                     self.dataOut.height - Height array (km)
                     self.dataOut.time - Time array (seconds)
                     self.dataOut.fileIndex -Index of the file currently read
                     self.dataOut.lat - Latitude coordinate of BLTR location
                     self.dataOut.doy - Experiment doy (number of the day in the current year)
                     self.dataOut.month - Experiment month
                     self.dataOut.day - Experiment day
                     self.dataOut.year - Experiment year
                 '''
                 def __init__(self, **kwargs):
                     '''
                     Inputs: None
                     '''
                     ProcessingUnit.__init__(self, **kwargs)
                     self.dataOut = Parameters()
+                    self.isConfig = False
-                def run(self, mode, snr_threshold=None):
+                def setup(self, mode):
+                    '''
                     '''
+                    self.dataOut.mode = mode
+                def run(self, mode, snr_threshold=None):
+                    '''
                     Inputs:
                         mode = High resolution (0) or Low resolution (1) data
                         snr_threshold = snr filter value
                     '''
+                    if not self.isConfig:
+                        self.setup(mode)
+                        self.isConfig = True
                     if self.dataIn.type == 'Parameters':
                         self.dataOut.copy(self.dataIn)
                     self.dataOut.data_output = self.dataOut.data_output[mode]
                     self.dataOut.heightList = self.dataOut.height[0]
                     self.dataOut.data_SNR = self.dataOut.data_SNR[mode]
                     if snr_threshold is not None:
                         SNRavg = numpy.average(self.dataOut.data_SNR, axis=0)
                         SNRavgdB = 10*numpy.log10(SNRavg)
                         for i in range(3):
                             self.dataOut.data_output[i][SNRavgdB <= snr_threshold] = numpy.nan
             # TODO
             class OutliersFilter(Operation):
                 def __init__(self, **kwargs):
                     '''
                     '''
                     Operation.__init__(self, **kwargs)
                 def run(self, svalue2, method, factor, filter, npoints=9):
                     '''
                     Inputs:
                         svalue    -  string to select array velocity
                         svalue2    -  string to choose axis filtering
                         method    - 0 for SMOOTH or 1 for MEDIAN
                         factor    - number used to set threshold
                         filter    - 1 for data filtering using the standard deviation criteria else 0
                         npoints    - number of points for mask filter
                     '''
                     print '    Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor)
                     yaxis = self.dataOut.heightList
                     xaxis = numpy.array([[self.dataOut.utctime]])
                     # Zonal
                     value_temp = self.dataOut.data_output[0]
                     # Zonal
                     value_temp = self.dataOut.data_output[1]
                     # Vertical
                     value_temp = numpy.transpose(self.dataOut.data_output[2])
                     htemp = yaxis
                     std = value_temp
                     for h in range(len(htemp)):
                         nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
                         minvalid = npoints
                         #only if valid values greater than the minimum required (10%)
                         if nvalues_valid > minvalid:
                             if method == 0:
                                 #SMOOTH
                                 w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
                             if method == 1:
                                 #MEDIAN
                                 w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
                             dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
                             threshold = dw*factor
                             value_temp[numpy.where(w > threshold),h] = numpy.nan
                             value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
                     #At the end
                     if svalue2 == 'inHeight':
                         value_temp = numpy.transpose(value_temp)
                     output_array[:,m] = value_temp
                     if svalue == 'zonal':
                         self.dataOut.data_output[0] = output_array
                     elif svalue == 'meridional':
                         self.dataOut.data_output[1] = output_array
                     elif svalue == 'vertical':
                         self.dataOut.data_output[2] = output_array
                     return self.dataOut.data_output
                 def Median(self,input,width):
                     '''
                     Inputs:
                         input - Velocity array
                         width - Number of points for mask filter
                     '''
                     if numpy.mod(width,2) == 1:
                         pc = int((width - 1) / 2)
                     cont = 0
                     output = []
                     for i in range(len(input)):
                         if i >= pc and i < len(input) - pc:
                             new2 = input[i-pc:i+pc+1]
                             temp = numpy.where(numpy.isfinite(new2))
                             new = new2[temp]
                             value = numpy.median(new)
                             output.append(value)
                     output = numpy.array(output)
                     output = numpy.hstack((input[0:pc],output))
                     output = numpy.hstack((output,input[-pc:len(input)]))
                     return output
                 def Smooth(self,input,width,edge_truncate = None):
                     '''
                     Inputs:
                         input - Velocity array
                         width - Number of points for mask filter
                         edge_truncate - 1 for truncate the convolution product else
                     '''
                     if numpy.mod(width,2) == 0:
                         real_width = width + 1
                         nzeros = width / 2
                     else:
                         real_width = width
                         nzeros = (width - 1) / 2
                     half_width = int(real_width)/2
                     length = len(input)
                     gate = numpy.ones(real_width,dtype='float')
                     norm_of_gate = numpy.sum(gate)
                     nan_process = 0
                     nan_id = numpy.where(numpy.isnan(input))
                     if len(nan_id[0]) > 0:
                         nan_process = 1
                         pb = numpy.zeros(len(input))
                         pb[nan_id] = 1.
                         input[nan_id] = 0.
                     if edge_truncate == True:
                         output = numpy.convolve(input/norm_of_gate,gate,mode='same')
                     elif edge_truncate == False or edge_truncate == None:
                         output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
                         output = numpy.hstack((input[0:half_width],output))
                         output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
                     if nan_process:
                         pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
                         pb = numpy.hstack((numpy.zeros(half_width),pb))
                         pb = numpy.hstack((pb,numpy.zeros(half_width)))
                         output[numpy.where(pb > 0.9999)] = numpy.nan
                         input[nan_id] = numpy.nan
                     return output
                 def Average(self,aver=0,nhaver=1):
                     '''
                     Inputs:
                         aver - Indicates the time period over which is averaged or consensus data
                         nhaver - Indicates the decimation factor in heights
                     '''
                     nhpoints = 48
                     lat_piura = -5.17
                     lat_huancayo = -12.04
                     lat_porcuya = -5.8
                     if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
                         hcm = 3.
                         if self.dataOut.year == 2003 :
                             if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
                                 nhpoints = 12
                     elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
                         hcm = 3.
                         if self.dataOut.year == 2003 :
                             if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
                                 nhpoints = 12
                     elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
                         hcm = 5.#2
                     pdata = 0.2
                     taver = [1,2,3,4,6,8,12,24]
                     t0 = 0
                     tf = 24
                     ntime =(tf-t0)/taver[aver]
                     ti = numpy.arange(ntime)
                     tf = numpy.arange(ntime) + taver[aver]
                     old_height = self.dataOut.heightList
                     if nhaver > 1:
                         num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
                         deltha = 0.05*nhaver
                         minhvalid = pdata*nhaver
                         for im in range(self.dataOut.nmodes):
                             new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
                     data_fHeigths_List = []
                     data_fZonal_List = []
                     data_fMeridional_List = []
                     data_fVertical_List = []
                     startDTList = []
                     for i in range(ntime):
                         height = old_height
                         start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
                         stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
                         limit_sec1 = time.mktime(start.timetuple())
                         limit_sec2 = time.mktime(stop.timetuple())
                         t1 = numpy.where(self.f_timesec >= limit_sec1)
                         t2 = numpy.where(self.f_timesec < limit_sec2)
                         time_select = []
                         for val_sec in t1[0]:
                             if val_sec in t2[0]:
                                 time_select.append(val_sec)
                         time_select = numpy.array(time_select,dtype = 'int')
                         minvalid = numpy.ceil(pdata*nhpoints)
                         zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
                         mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
                         ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
                         if nhaver > 1:
                             new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
                             new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
                             new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
                         if len(time_select) > minvalid:
                             time_average = self.f_timesec[time_select]
                             for im in range(self.dataOut.nmodes):
                                 for ih in range(self.dataOut.nranges):
                                     if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
                                         zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
                                     if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
                                         mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
                                     if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
                                         ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
                                 if nhaver > 1:
                                     for ih in range(num_hei):
                                         hvalid = numpy.arange(nhaver) + nhaver*ih
                                         if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
                                             new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) /  numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
                                         if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
                                             new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) /  numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
                                         if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
                                             new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) /  numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
                             if nhaver > 1:
                                 zon_aver = new_zon_aver
                                 mer_aver = new_mer_aver
                                 ver_aver = new_ver_aver
                                 height = new_height
                             tstart = time_average[0]
                             tend = time_average[-1]
                             startTime = time.gmtime(tstart)
                             year = startTime.tm_year
                             month = startTime.tm_mon
                             day = startTime.tm_mday
                             hour = startTime.tm_hour
                             minute = startTime.tm_min
                             second = startTime.tm_sec
                             startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
                             o_height = numpy.array([])
                             o_zon_aver = numpy.array([])
                             o_mer_aver = numpy.array([])
                             o_ver_aver = numpy.array([])
                             if self.dataOut.nmodes > 1:
                                 for im in range(self.dataOut.nmodes):
                                     if im == 0:
                                         h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
                                     else:
                                         h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
                                     ht = h_select[0]
                                     o_height = numpy.hstack((o_height,height[im,ht]))
                                     o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
                                     o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
                                     o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
                                 data_fHeigths_List.append(o_height)
                                 data_fZonal_List.append(o_zon_aver)
                                 data_fMeridional_List.append(o_mer_aver)
                                 data_fVertical_List.append(o_ver_aver)
                             else:
                                 h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
                                 ht = h_select[0]
                                 o_height = numpy.hstack((o_height,height[im,ht]))
                                 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
                                 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
                                 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
                                 data_fHeigths_List.append(o_height)
                                 data_fZonal_List.append(o_zon_aver)
                                 data_fMeridional_List.append(o_mer_aver)
                                 data_fVertical_List.append(o_ver_aver)
                     return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List
   No newline at end of file

schainpy/__init__.pyc removed binary 0 0