Context Navigation

source: lmdz_wrf/trunk/tools/create_OBSnetcdf.py @ 1949

Last change on this file since 1949 was 1949, checked in by lfita, 7 years ago

Adding description of WMO codes into `create_OBSnetcdf.py'

wmo_0509.code: WMO code of description og high clouds

File size: 52.5 KB

Rev	Line
[259]	1	# -- coding: iso-8859-15 --
	2	# Generic program to transfrom ASCII observational data in columns to a netcdf
	3	# L. Fita, LMD February 2015
	4	## e.g. # create_OBSnetcdf.py -c '#' -d ACAR/description.dat -e space -f ACAR/2012/10/ACAR_121018.asc -g true -t 19491201000000,seconds -k trajectory
	5
	6	import numpy as np
	7	from netCDF4 import Dataset as NetCDFFile
	8	import os
	9	import re
	10	from optparse import OptionParser
	11
	12	# version
	13	version=1.1
	14
	15	# Filling values for floats, integer and string
	16	fillValueF = 1.e20
	17	fillValueI = -99999
	18	fillValueS = '---'
	19
	20	# Length of the string variables
	21	StringLength = 200
	22
	23	# Size of the map for the complementary variables/maps
	24	Ndim2D = 100
	25
	26	main = 'create_OBSnetcdf.py'
	27	errormsg = 'ERROR -- error -- ERROR -- error'
	28	warnmsg = 'WARNING -- warning -- WARNING -- warning'
	29
	30	fillValue = 1.e20
	31
	32	def searchInlist(listname, nameFind):
	33	""" Function to search a value within a list
	34	listname = list
	35	nameFind = value to find
	36	>>> searInlist(['1', '2', '3', '5'], '5')
	37	True
	38	"""
	39	for x in listname:
	40	if x == nameFind:
	41	return True
	42	return False
	43
[1948]	44
	45	def typemod(value, typeval):
	46	""" Function to give back a value in a given dtype
	47	>>> print(typemod(8.2223, 'np.float64'))
	48	<type 'numpy.float64'>
	49	>>> print(typemod(8.2223, 'tuple'))
	50	<type 'tuple'>
	51	"""
	52	fname='typemod'
	53
	54	if typeval == 'int' or typeval == 'I':
	55	return int(value)
	56	elif typeval == 'long':
	57	return long(value)
	58	elif typeval == 'float' or typeval == 'F' or typeval == 'R':
	59	return float(value)
	60	elif typeval == 'complex':
	61	return complex(value)
	62	elif typeval == 'str' or typeval == 'S':
	63	return str(value)
	64	elif typeval == 'bool':
	65	return bool(value)
	66	elif typeval == 'B':
	67	return Str_Bool(value)
	68	elif typeval == 'list':
	69	newval = []
	70	newval.append(value)
	71	return newval
	72	elif typeval == 'dic':
	73	newval = {}
	74	newval[value] = value
	75	return newval
	76	elif typeval == 'tuple':
	77	newv = []
	78	newv.append(value)
	79	newval = tuple(newv)
	80	return newval
	81	elif typeval == 'np.int8':
	82	return np.int8(value)
	83	elif typeval == 'np.int16':
	84	return np.int16(value)
	85	elif typeval == 'np.int32':
	86	return np.int32(value)
	87	elif typeval == 'np.int64':
	88	return np.int64(value)
	89	elif typeval == 'np.uint8':
	90	return np.uint8(value)
	91	elif typeval == 'np.uint16':
	92	return np.uint16(value)
	93	elif typeval == 'np.np.uint32':
	94	return np.uint32(value)
	95	elif typeval == 'np.uint64':
	96	return np.uint64(value)
	97	elif typeval == 'np.float' or typeval == 'R':
	98	return np.float(value)
	99	elif typeval == 'np.float16':
	100	return np.float16(value)
	101	elif typeval == 'np.float32':
	102	return np.float32(value)
	103	elif typeval == 'float32':
	104	return np.float32(value)
	105	elif typeval == 'np.float64' or typeval == 'D':
	106	return np.float64(value)
	107	elif typeval == 'float64':
	108	return np.float64(value)
	109	elif typeval == 'np.complex':
	110	return np.complex(value)
	111	elif typeval == 'np.complex64':
	112	return np.complex64(value)
	113	elif typeval == 'np.complex128':
	114	return np.complex128(value)
	115	else:
	116	print errormsg
	117	print fname + ': data type "' + typeval + '" is not ready !'
	118	print errormsg
	119	quit(-1)
	120
	121	return
	122
	123	def str_list_k(string, cdiv, kind):
	124	""" Function to obtain a list of types of values from a string giving a split character
	125	string= String from which to obtain a list ('None' for None)
	126	cdiv= character to use to split the string
	127	kind= kind of desired value (as string like: 'np.float', 'int', 'np.float64', ....)
	128	>>> str_list_k('1:@#:$:56', ':', 'S')
	129	['1', '@#', '$', '56']
	130	>>> str_list_k('1:3.4:12.3', ':', 'np.float64')
	131	[1.0, 3.3999999999999999, 12.300000000000001]
	132	"""
	133	fname = 'str_list'
	134
	135	if string.find(cdiv) != -1:
	136	listv = string.split(cdiv)
	137	else:
	138	if string == 'None':
	139	listv = None
	140	else:
	141	listv = [string]
	142
	143	if listv is not None:
	144	finalist = []
	145	for lv in listv:
	146	finalist.append(typemod(lv, kind))
	147	else:
	148	finalist = None
	149
	150	return finalist
	151
	152
[259]	153	def set_attribute(ncvar, attrname, attrvalue):
	154	""" Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
	155	ncvar = object netcdf variable
	156	attrname = name of the attribute
	157	attrvalue = value of the attribute
	158	"""
	159	import numpy as np
	160	from netCDF4 import Dataset as NetCDFFile
	161
	162	attvar = ncvar.ncattrs()
	163	if searchInlist(attvar, attrname):
	164	attr = ncvar.delncattr(attrname)
	165
	166	attr = ncvar.setncattr(attrname, attrvalue)
	167
	168	return ncvar
	169
	170	def basicvardef(varobj, vstname, vlname, vunits):
	171	""" Function to give the basic attributes to a variable
	172	varobj= netCDF variable object
	173	vstname= standard name of the variable
	174	vlname= long name of the variable
	175	vunits= units of the variable
	176	"""
	177	attr = varobj.setncattr('standard_name', vstname)
	178	attr = varobj.setncattr('long_name', vlname)
	179	attr = varobj.setncattr('units', vunits)
	180
[332]	181	return
	182
[259]	183	def remove_NONascii(string):
	184	""" Function to remove that characters which are not in the standard 127 ASCII
	185	string= string to transform
	186	>>> remove_NONascii('LluÃs')
	187	Lluis
	188	"""
	189	fname = 'remove_NONascii'
	190
	191	newstring = string
	192
	193	RTFchar= ['Ã¡', 'Ã©', 'Ã', 'Ã³', 'Ãº', 'Ã ', 'Ãš', 'Ã¬', 'Ã²', 'Ã¹', 'Ã¢', 'Ãª', 'Ã®', 'ÃŽ', \
	194	'Ã»', 'Ã€', 'Ã«', 'Ã¯', 'Ã¶', 'ÃŒ', 'Ã§', 'Ã±','ÃŠ', 'Å', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', \
	195	'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã',\
	196	'Ã', 'Å', '\n', '\t']
	197	ASCchar= ['a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o', \
	198	'u', 'a', 'e', 'i', 'o', 'u', 'c', 'n','ae', 'oe', 'A', 'E', 'I', 'O', 'U', 'A', \
	199	'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'C', 'N',\
	200	'AE', 'OE', '', ' ']
	201
	202	Nchars = len(RTFchar)
	203	for ichar in range(Nchars):
	204	foundchar = string.find(RTFchar[ichar])
[332]	205	if foundchar != -1:
[259]	206	newstring = newstring.replace(RTFchar[ichar], ASCchar[ichar])
	207
	208	return newstring
	209
	210	def read_description(fdobs, dbg):
	211	""" reads the description file of the observational data-set
	212	fdobs= descriptive observational data-set
	213	dbg= boolean argument for debugging
	214	* Each station should have a 'description.dat' file with:
	215	institution=Institution who creates the data
	216	department=Department within the institution
	217	scientists=names of the data producers
	218	contact=contact of the data producers
	219	description=description of the observations
	220	acknowledgement=sentence of acknowlegement
	221	comment=comment for the measurements
	222
	223	MissingValue='\|' list of ASCII values for missing values within the data
[553]	224	(as they appear!, 'empty' for no value at all)
[259]	225	comment=comments
	226
	227	varN='\|' list of variable names
	228	varLN='\|' list of long variable names
	229	varU='\|' list units of the variables
	230	varBUFR='\|' list BUFR code of the variables
	231	varTYPE='\|' list of variable types ('D', 'F', 'I', 'I64', 'S')
[553]	232	varOPER='\|' list of operations to do to the variables to meet their units ([oper],[val])
[587]	233	[oper]:
	234	-, nothing
	235	sumc, add [val]
	236	subc, rest [val]
	237	mulc, multiply by [val]
	238	divc, divide by [val]
	239	rmchar,[val],[pos], remove [val] characters from [pos]='B', beginning, 'E', end
[259]	240	NAMElon=name of the variable with the longitude (x position)
	241	NAMElat=name of the variable with the latitude (y position)
	242	NAMEheight=ind_alt
	243	NAMEtime=name of the varibale with the time
	244	FMTtime=format of the time (as in 'C', 'CFtime' for already CF-like time)
	245
	246	"""
	247	fname = 'read_description'
	248
	249	descobj = open(fdobs, 'r')
	250	desc = {}
	251
	252	namevalues = []
	253
	254	for line in descobj:
	255	if line[0:1] != '#' and len(line) > 1 :
	256	descn = remove_NONascii(line.split('=')[0])
	257	descv = remove_NONascii(line.split('=')[1])
	258	namevalues.append(descn)
	259	if descn[0:3] != 'var':
	260	if descn != 'MissingValue':
	261	desc[descn] = descv
	262	else:
	263	desc[descn] = []
	264	for dn in descv.split('\|'): desc[descn].append(dn)
	265	print ' ' + fname + ': missing values found:',desc[descn]
	266	elif descn[0:3] == 'var':
	267	desc[descn] = descv.split('\|')
	268	elif descn[0:4] == 'NAME':
	269	desc[descn] = descv
	270	elif descn[0:3] == 'FMT':
	271	desc[descn] = descv
[553]	272	if not desc.has_key('varOPER'): desc['varOPER'] = None
[259]	273
	274	if dbg:
	275	Nvars = len(desc['varN'])
	276	print ' ' + fname + ": description content of '" + fdobs + "'______________"
	277	for varn in namevalues:
	278	if varn[0:3] != 'var':
	279	print ' ' + varn + ':',desc[varn]
	280	elif varn == 'varN':
	281	print ' * Variables:'
	282	for ivar in range(Nvars):
	283	varname = desc['varN'][ivar]
	284	varLname = desc['varLN'][ivar]
	285	varunits = desc['varU'][ivar]
	286	if desc.has_key('varBUFR'):
[553]	287	varbufr = desc['varBUFR'][ivar]
[259]	288	else:
[553]	289	varbufr = None
	290	if desc['varOPER'] is not None:
	291	opv = desc['varOPER'][ivar]
	292	else:
	293	opv = None
	294	print ' ', ivar, varname+':',varLname,'[',varunits, \
	295	']','bufr code:',varbufr,'oper:',opv
[259]	296
	297	descobj.close()
	298
	299	return desc
	300
[553]	301	def value_fmt(val, miss, op, fmt):
[259]	302	""" Function to transform an ASCII value to a given format
	303	val= value to transform
	304	miss= list of possible missing values
[553]	305	op= operation to perform to the value
[259]	306	fmt= format to take:
	307	'D': float double precission
	308	'F': float
	309	'I': integer
	310	'I64': 64-bits integer
	311	'S': string
	312	>>> value_fmt('9876.12', '-999', 'F')
	313	9876.12
	314	"""
	315
	316	fname = 'value_fmt'
	317
[587]	318	aopers = ['sumc','subc','mulc','divc', 'rmchar']
[553]	319
[259]	320	fmts = ['D', 'F', 'I', 'I64', 'S']
	321	Nfmts = len(fmts)
	322
	323	if not searchInlist(miss,val):
[553]	324	if searchInlist(miss,'empty') and len(val) == 0:
	325	newval = None
[259]	326	else:
[553]	327	if op != '-':
	328	opern = op.split(',')[0]
[587]	329	operv = op.split(',')[1]
[553]	330
	331	if not searchInlist(aopers,opern):
	332	print errormsg
	333	print ' ' + fname + ": operation '" + opern + "' not ready!!"
	334	print ' availables:',aopers
	335	quit(-1)
	336	else:
	337	opern = 'sumc'
[1948]	338	operv = '0'
[553]	339
	340	if not searchInlist(fmts, fmt):
	341	print errormsg
	342	print ' ' + fname + ": format '" + fmt + "' not ready !!"
	343	quit(-1)
	344	else:
	345	if fmt == 'D':
	346	opv = np.float32(operv)
	347	if opern == 'sumc': newval = np.float32(val) + opv
	348	elif opern == 'subc': newval = np.float32(val) - opv
	349	elif opern == 'mulc': newval = np.float32(val) * opv
	350	elif opern == 'divc': newval = np.float32(val) / opv
[587]	351	elif opern == 'rmchar':
	352	opv = int(operv)
	353	Lval = len(val)
	354	if op.split(',')[2] == 'B':
	355	newval = np.float32(val[opv:Lval+1])
	356	elif op.split(',')[2] == 'E':
	357	newval = np.float32(val[Lval-opv:Lval])
	358	else:
	359	print errormsg
	360	print ' ' + fname + ": operation '" + opern + "' not " +\
	361	" work with '" + op.split(',')[2] + "' !!"
	362	quit(-1)
[553]	363	elif fmt == 'F':
	364	opv = np.float(operv)
	365	if opern == 'sumc': newval = np.float(val) + opv
	366	elif opern == 'subc': newval = np.float(val) - opv
	367	elif opern == 'mulc': newval = np.float(val) * opv
	368	elif opern == 'divc': newval = np.float(val) / opv
[587]	369	elif opern == 'rmchar':
	370	opv = int(operv)
	371	Lval = len(val)
	372	if op.split(',')[2] == 'B':
	373	newval = np.float(val[opv:Lval+1])
	374	elif op.split(',')[2] == 'E':
	375	newval = np.float(val[0:Lval-opv])
	376	else:
	377	print errormsg
	378	print ' ' + fname + ": operation '" + opern + "' not " +\
	379	" work with '" + op.split(',')[2] + "' !!"
	380	quit(-1)
	381
[553]	382	elif fmt == 'I':
	383	opv = int(operv)
	384	if opern == 'sumc': newval = int(val) + opv
	385	elif opern == 'subc': newval = int(val) - opv
	386	elif opern == 'mulc': newval = int(val) * opv
	387	elif opern == 'divc': newval = int(val) / opv
[587]	388	elif opern == 'rmchar':
	389	opv = int(operv)
	390	Lval = len(val)
	391	if op.split(',')[2] == 'B':
	392	newval = int(val[opv:Lval+1])
	393	elif op.split(',')[2] == 'E':
	394	newval = int(val[0:Lval-opv])
	395	else:
	396	print errormsg
	397	print ' ' + fname + ": operation '" + opern + "' not " +\
	398	" work with '" + op.split(',')[2] + "' !!"
	399	quit(-1)
[553]	400	elif fmt == 'I64':
	401	opv = np.int64(operv)
	402	if opern == 'sumc': newval = np.int64(val) + opv
	403	elif opern == 'subc': newval = np.int64(val) - opv
	404	elif opern == 'mulc': newval = np.int64(val) * opv
	405	elif opern == 'divc': newval = np.int64(val) / opv
[587]	406	elif opern == 'rmchar':
	407	opv = int(operv)
	408	Lval = len(val)
	409	if op.split(',')[2] == 'B':
	410	newval = np.int64(val[opv:Lval+1])
	411	elif op.split(',')[2] == 'E':
	412	newval = np.int64(val[0:Lval-opv])
	413	else:
	414	print errormsg
	415	print ' ' + fname + ": operation '" + opern + "' not " +\
	416	" work with '" + op.split(',')[2] + "' !!"
	417	quit(-1)
[553]	418	elif fmt == 'S':
[587]	419	if opern == 'rmchar':
	420	opv = int(operv)
	421	Lval = len(val)
	422	if op.split(',')[2] == 'B':
	423	newval = val[opv:Lval+1]
	424	elif op.split(',')[2] == 'E':
	425	newval = val[0:Lval-opv]
	426	else:
	427	print errormsg
	428	print ' ' + fname + ": operation '" + opern + "' not " +\
	429	" work with '" + op.split(',')[2] + "' !!"
	430	quit(-1)
	431	else:
	432	newval = val
	433
[259]	434	else:
	435	newval = None
	436
	437	return newval
	438
[1948]	439	def getting_fixedline(line, cuts, types, dbg=False):
	440	""" Function to get the values from a line of text with fixed lenght of different values
	441	line: line with values
	442	cuts: character number where a value ends
	443	types: consecutive type of values
	444	'I': integer
	445	'R': real
	446	'D': float64
	447	'S': string
	448	'B': boolean
	449	dbg: debug mode (default False)
	450	>>> Sline=' 87007 03012015 25.6 6.4 9.4 5 15'
	451	>>> getting_fixedline(Sline, [8, 17, 23, 29, 36, 40, 45, 50], ['I', 'R', 'R', 'R', 'I', 'R', 'R', 'R'])
	452	[87007, 3012015.0, 25.6, 6.4, -99999, -99999, 9.4, 5.0, 15.0]
	453	"""
	454	fname = 'getting_fixedline'
	455
	456	if len(cuts) + 1 != len(types):
	457	print errormsg
	458	print ' ' + fname + ': The number of types :', len(types), 'must be +1', \
	459	'number of cuts:', len(cuts)
	460	quit(-1)
	461
	462	values = []
	463	val = line[0:cuts[0]]
	464	if len(val.replace(' ','')) >= 1:
	465	values.append(typemod(val, types[0]))
	466	else:
	467	if types[0] == 'I': values.append(fillValueI)
	468	elif types[0] == 'R': values.append(fillValueF)
	469	elif types[0] == 'D': values.append(fillValueF)
	470	elif types[0] == 'S': values.append(fillValueS)
	471	elif types[0] == 'B': values.append(fillValueB)
	472
	473	Ncuts = len(cuts)
	474	for ic in range(1,Ncuts):
	475	val = line[cuts[ic-1]:cuts[ic]]
	476	if dbg: print ic, ':', val, '-->', types[ic]
	477	if len(val.replace(' ','')) >= 1:
	478	values.append(typemod(val, types[ic]))
	479	else:
	480	if types[ic] == 'I': values.append(fillValueI)
	481	elif types[ic] == 'R': values.append(fillValueF)
	482	elif types[ic] == 'D': values.append(fillValueF)
	483	elif types[ic] == 'S': values.append(fillValueS)
	484	elif types[ic] == 'B': values.append(fillValueB)
	485
	486	# Last value
	487	Lline = len(line)
	488	val = line[cuts[Ncuts-1]:Lline]
	489	if len(val.replace(' ','')) >= 1:
	490	values.append(typemod(val, types[Ncuts]))
	491	else:
	492	if types[Ncuts] == 'I': values.append(fillValueI)
	493	elif types[Ncuts] == 'R': values.append(fillValueF)
	494	elif types[Ncuts] == 'D': values.append(fillValueF)
	495	elif types[Ncuts] == 'S': values.append(fillValueS)
	496	elif types[Ncuts] == 'B': values.append(fillValueB)
	497
	498	return values
	499
	500
	501	def getting_fixedline_NOk(line, cuts, miss, dbg=False):
	502	""" Function to get the values from a line of text with fixed lenght of
	503	different values without types
	504	line: line with values
	505	cuts: character number where a value ends
	506	miss: character form issed values
	507	dbg: debug mode (default False)
	508	>>> Sline=' 87007 03012015 25.6 6.4 9.4 5 15'
	509	>>> getting_fixedline_NOk(Sline, [8, 17, 23, 29, 36, 40, 45, 50], '#')
	510	['87007', '03012015', '25.6', '6.4', '#', '#', '9.4', '5', '15']
	511	"""
	512	fname = 'getting_fixedline_NOk'
	513
	514	values = []
	515	val = line[0:cuts[0]]
	516	if len(val.replace(' ','')) >= 1:
	517	values.append(val.replace(' ',''))
	518	else:
	519	values.append(miss)
	520
	521	Ncuts = len(cuts)
	522	for ic in range(1,Ncuts):
	523	val = line[cuts[ic-1]:cuts[ic]]
	524	if dbg: print ic, ':', val
	525	if len(val.replace(' ','')) >= 1:
	526	values.append(val.replace(' ',''))
	527	else:
	528	values.append(miss)
	529
	530	# Last value
	531	Lline = len(line)
	532	val = line[cuts[Ncuts-1]:Lline]
	533	if len(val.replace(' ','')) >= 1:
	534	values.append(val.replace(' ',''))
	535	else:
	536	values.append(miss)
	537
	538	return values
	539
[553]	540	def read_datavalues(dataf, comchar, colchar, fmt, oper, miss, varns, dbg):
[259]	541	""" Function to read from an ASCII file values in column
	542	dataf= data file
	543	comchar= list of the characters indicating comment in the file
	544	colchar= character which indicate end of value in the column
	545	dbg= debug mode or not
	546	fmt= list of kind of values to be found
[553]	547	oper= list of operations to perform
[259]	548	miss= missing value
	549	varns= list of name of the variables to find
	550	"""
	551
	552	fname = 'read_datavalues'
	553
	554	ofile = open(dataf, 'r')
	555	Nvals = len(fmt)
	556
[553]	557	if oper is None:
	558	opers = []
	559	for ioper in range(Nvals):
	560	opers.append('-')
	561	else:
	562	opers = oper
	563
[259]	564	finalvalues = {}
	565
	566	iline = 0
	567	for line in ofile:
	568	line = line.replace('\n','').replace(chr(13),'')
	569	if not searchInlist(comchar,line[0:1]) and len(line) > 1:
[1948]	570	# Getting values
	571	if colchar[0:4] != 'None':
	572	values0 = line.split(colchar)
	573	else:
	574	ltypes=str_list_k(colchar[5:len(colchar)+1],',','I')
	575	values0 = getting_fixedline_NOk(line, ltypes, miss[0], dbg=dbg)
[259]	576	# Removing no-value columns
	577	values = []
	578	for iv in values0:
	579	if len(iv) > 0: values.append(iv)
	580
	581	Nvalues = len(values)
	582	# Checkings for wierd characters at the end of lines (use it to check)
	583	# if values[Nvalues-1][0:4] == '-999':
	584	# print line,'last v:',values[Nvalues-1],'len:',len(values[Nvalues-1])
	585	# for ic in range(len(values[Nvalues-1])):
	586	# print ic,ord(values[Nvalues-1][ic:ic+1])
	587	# quit()
	588
	589	if len(values[Nvalues-1]) == 0:
	590	Nvalues = Nvalues - 1
	591
	592	if Nvalues != Nvals:
	593	print warnmsg
	594	print ' ' + fname + ': number of formats:',Nvals,' and number of ', \
	595	'values:',Nvalues,' with split character *' + colchar + \
	596	'* does not coincide!!'
	597	print ' * what is found is ________'
	598	if Nvalues > Nvals:
	599	Nshow = Nvals
	600	for ivar in range(Nshow):
	601	print ' ',varns[ivar],'fmt:',fmt[ivar],'value:',values[ivar]
	602	print ' missing formats for:',values[Nshow:Nvalues+1]
	603	print ' values not considered, continue'
	604	else:
	605	Nshow = Nvalues
	606	for ivar in range(Nshow):
	607	print ' ',varns[ivar],'fmt:',fmt[ivar],'value:',values[ivar]
	608	print ' excess of formats:',fmt[Nshow:Nvals+1]
	609	quit(-1)
	610
	611	# Reading and transforming values
	612	if dbg: print ' ' + fname + ': values found _________'
	613
	614	for ivar in range(Nvals):
	615	if dbg:
[553]	616	print iline, varns[ivar],'value:',values[ivar],miss,opers[ivar], \
	617	fmt[ivar]
[259]	618
	619	if iline == 0:
	620	listvals = []
[553]	621	listvals.append(value_fmt(values[ivar], miss, opers[ivar], \
	622	fmt[ivar]))
[259]	623	finalvalues[varns[ivar]] = listvals
	624	else:
	625	listvals = finalvalues[varns[ivar]]
[553]	626	listvals.append(value_fmt(values[ivar], miss, opers[ivar], \
	627	fmt[ivar]))
[259]	628	finalvalues[varns[ivar]] = listvals
	629	else:
	630	# First line without values
	631	if iline == 0: iline = -1
	632
	633	iline = iline + 1
	634
	635	ofile.close()
	636
	637	return finalvalues
	638
	639	def writing_str_nc(varo, values, Lchar):
	640	""" Function to write string values in a netCDF variable as a chain of 1char values
	641	varo= netCDF variable object
	642	values = list of values to introduce
	643	Lchar = length of the string in the netCDF file
	644	"""
	645
	646	Nvals = len(values)
[553]	647
[259]	648	for iv in range(Nvals):
	649	stringv=values[iv]
	650	charvals = np.chararray(Lchar)
	651	Lstr = len(stringv)
	652	charvals[Lstr:Lchar] = ''
	653
	654	for ich in range(Lstr):
	655	charvals[ich] = stringv[ich:ich+1]
	656
	657	varo[iv,:] = charvals
	658
	659	return
	660
	661	def Stringtimes_CF(tvals, fmt, Srefdate, tunits, dbg):
	662	""" Function to transform a given data in String formats to a CF date
	663	tvals= string temporal values
	664	fmt= format of the the time values
	665	Srefdate= reference date in [YYYY][MM][DD][HH][MI][SS] format
	666	tunits= units to use ('weeks', 'days', 'hours', 'minutes', 'seconds')
	667	dbg= debug
	668	>>> Stringtimes_CF(['19760217082712','20150213101837'], '%Y%m%d%H%M%S',
	669	'19491201000000', 'hours', False)
	670	[229784.45333333 571570.31027778]
	671	"""
	672	import datetime as dt
	673
	674	fname = 'Stringtimes'
	675
	676	dimt = len(tvals)
	677
	678	yrref = int(Srefdate[0:4])
	679	monref = int(Srefdate[4:6])
	680	dayref = int(Srefdate[6:8])
	681	horref = int(Srefdate[8:10])
	682	minref = int(Srefdate[10:12])
	683	secref = int(Srefdate[12:14])
	684	refdate = dt.datetime( yrref, monref, dayref, horref, minref, secref)
	685
	686	cftimes = np.zeros((dimt), dtype=np.float)
	687
	688	Nfmt=len(fmt.split('%'))
	689
	690	if dbg: print ' ' + fname + ': fmt=',fmt,'refdate:',Srefdate,'uits:',tunits, \
	691	'date dt_days dt_time deltaseconds CFtime _______'
	692	for it in range(dimt):
	693
	694	# Removing excess of mili-seconds (up to 6 decimals)
	695	if fmt.split('%')[Nfmt-1] == 'f':
	696	tpoints = tvals[it].split('.')
	697	if len(tpoints[len(tpoints)-1]) > 6:
	698	milisec = '{0:.6f}'.format(np.float('0.'+tpoints[len(tpoints)-1]))[0:7]
	699	newtval = ''
	700	for ipt in range(len(tpoints)-1):
	701	newtval = newtval + tpoints[ipt] + '.'
	702	newtval = newtval + str(milisec)[2:len(milisec)+1]
	703	else:
	704	newtval = tvals[it]
	705	tval = dt.datetime.strptime(newtval, fmt)
	706	else:
	707	tval = dt.datetime.strptime(tvals[it], fmt)
	708
	709	deltatime = tval - refdate
	710	deltaseconds = deltatime.days24.3600. + deltatime.seconds + \
	711	deltatime.microseconds/100000.
	712	if tunits == 'weeks':
	713	deltat = 7.24.3600.
	714	elif tunits == 'days':
	715	deltat = 24.*3600.
	716	elif tunits == 'hours':
	717	deltat = 3600.
	718	elif tunits == 'minutes':
	719	deltat = 60.
	720	elif tunits == 'seconds':
	721	deltat = 1.
	722	else:
	723	print errormsg
	724	print ' ' + fname + ": time units '" + tunits + "' not ready !!"
	725	quit(-1)
	726
	727	cftimes[it] = deltaseconds / deltat
	728	if dbg:
	729	print ' ' + tvals[it], deltatime, deltaseconds, cftimes[it]
	730
	731	return cftimes
	732
	733	def adding_complementary(onc, dscn, okind, dvalues, tvals, refCFt, CFtu, Nd, dbg):
	734	""" Function to add complementary variables as function of the observational type
	735	onc= netcdf objecto file to add the variables
	736	dscn= description dictionary
	737	okind= observational kind
	738	dvalues= values
	739	tvals= CF time values
	740	refCFt= reference time of CF time (in [YYYY][MM][DD][HH][MI][SS])
	741	CFtu= CF time units
	742	Nd= size of the domain
	743	dbg= debugging flag
	744	"""
	745	import numpy.ma as ma
	746
	747	fname = 'adding_complementary'
	748
	749	# Kind of observations which require de integer lon/lat (for the 2D map)
	750	map2D=['multi-points', 'trajectory']
	751
	752	SrefCFt = refCFt[0:4] +'-'+ refCFt[4:6] +'-'+ refCFt[6:8] + ' ' + refCFt[8:10] + \
	753	':'+ refCFt[10:12] +':'+ refCFt[12:14]
	754
	755	if dscn['NAMElon'] == '-' or dscn['NAMElat'] == '-':
	756	print errormsg
	757	print ' ' + fname + ": to complement a '" + okind + "' observation kind " + \
	758	" a given longitude ('NAMElon':",dscn['NAMElon'],") and latitude ('" + \
	759	"'NAMElat:'", dscn['NAMElat'],') from the data has to be provided and ' + \
	760	'any are given !!'
	761	quit(-1)
	762
	763	if not dvalues.has_key(dscn['NAMElon']) or not dvalues.has_key(dscn['NAMElat']):
	764	print errormsg
	765	print ' ' + fname + ": observations do not have 'NAMElon':", \
	766	dscn['NAMElon'],"and/or 'NAMElat:'", dscn['NAMElat'],' !!'
	767	print ' available data:',dvalues.keys()
	768	quit(-1)
	769
	770	if okind == 'trajectory':
	771	if dscn['NAMEheight'] == '-':
	772	print warnmsg
	773	print ' ' + fname + ": to complement a '" + okind + "' observation " + \
	774	"kind a given height ('NAMEheight':",dscn['NAMEheight'],"') might " + \
	775	'be provided and any is given !!'
	776	quit(-1)
	777
	778	if not dvalues.has_key(dscn['NAMEheight']):
	779	print errormsg
	780	print ' ' + fname + ": observations do not have 'NAMEtime':", \
	781	dscn['NAMEtime'],' !!'
	782	print ' available data:',dvalues.keys()
	783	quit(-1)
	784
	785	if searchInlist(map2D, okind):
	786	# A new 2D map with the number of observation will be added for that 'NAMElon'
	787	# and 'NAMElat' are necessary. A NdxNd domain space size will be used.
	788	objfile.createDimension('lon2D',Nd)
	789	objfile.createDimension('lat2D',Nd)
	790	lonvals = ma.masked_equal(dvalues[dscn['NAMElon']], None)
	791	latvals = ma.masked_equal(dvalues[dscn['NAMElat']], None)
	792
	793	minlon = min(lonvals)
	794	maxlon = max(lonvals)
	795	minlat = min(latvals)
	796	maxlat = max(latvals)
	797
	798	blon = (maxlon - minlon)/(Nd-1)
	799	blat = (maxlat - minlat)/(Nd-1)
	800
	801	newvar = onc.createVariable( 'lon2D', 'f8', ('lon2D'))
	802	basicvardef(newvar, 'longitude', 'longitude map observations','degrees_East')
	803	newvar[:] = minlon + np.arange(Nd)*blon
	804	newattr = set_attribute(newvar, 'axis', 'X')
	805
	806	newvar = onc.createVariable( 'lat2D', 'f8', ('lat2D'))
	807	basicvardef(newvar, 'latitude', 'latitude map observations', 'degrees_North')
	808	newvar[:] = minlat + np.arange(Nd)*blat
	809	newattr = set_attribute(newvar, 'axis', 'Y')
	810
	811	if dbg:
	812	print ' ' + fname + ': minlon=',minlon,'maxlon=',maxlon
	813	print ' ' + fname + ': minlat=',minlat,'maxlat=',maxlat
	814	print ' ' + fname + ': precission on x-axis=', blon*(Nd-1), 'y-axis=', \
	815	blat*(Nd-1)
	816
	817	if okind == 'multi-points':
	818	map2D = np.ones((Nd, Nd), dtype=np.float)*fillValueI
	819
	820	dimt = len(tvals)
	821	Nlost = 0
	822	for it in range(dimt):
	823	lon = dvalues[dscn['NAMElon']][it]
	824	lat = dvalues[dscn['NAMElat']][it]
	825	if lon is not None and lat is not None:
	826	ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
	827	ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
	828
	829	if map2D[ilat,ilon] == fillValueI:
	830	map2D[ilat,ilon] = 1
	831	else:
	832	map2D[ilat,ilon] = map2D[ilat,ilon] + 1
	833	if dbg: print it, lon, lat, ilon, ilat, map2D[ilat,ilon]
	834
	835	newvar = onc.createVariable( 'mapobs', 'f4', ('lat2D', 'lon2D'), \
	836	fill_value = fillValueI)
	837	basicvardef(newvar, 'map_observations', 'number of observations', '-')
	838	newvar[:] = map2D
	839	newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D')
	840
	841	elif okind == 'trajectory':
	842	# A new 2D map with the trajectory 'NAMElon' and 'NAMElat' and maybe 'NAMEheight'
	843	# are necessary. A NdxNdxNd domain space size will be used. Using time as
	844	# reference variable
	845	if dscn['NAMEheight'] == '-':
	846	# No height
	847	map2D = np.ones((Nd, Nd), dtype=np.float)*fillValueI
	848
	849	dimt = len(tvals)
	850	Nlost = 0
	851	if dbg: print ' time-step lon lat ix iy passes _______'
	852	for it in range(dimt):
	853	lon = dvalues[dscn['NAMElon']][it]
	854	lat = dvalues[dscn['NAMElat']][it]
	855	if lon is not None and lat is not None:
	856	ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
	857	ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
	858
	859	if map2D[ilat,ilon] == fillValueI:
	860	map2D[ilat,ilon] = 1
	861	else:
	862	map2D[ilat,ilon] = map2D[ilat,ilon] + 1
	863	if dbg: print it, lon, lat, ilon, ilat, map2D[ilat,ilon]
	864
	865	newvar = onc.createVariable( 'trjobs', 'i', ('lat2D', 'lon2D'), \
	866	fill_value = fillValueI)
	867	basicvardef(newvar, 'trajectory_observations', 'number of passes', '-' )
	868	newvar[:] = map2D
	869	newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D')
	870
	871	else:
	872	ivn = 0
	873	for vn in dscn['varN']:
	874	if vn == dscn['NAMEheight']:
	875	zu = dscn['varU'][ivn]
	876	break
	877	ivn = ivn + 1
	878
	879	objfile.createDimension('z2D',Nd)
	880	zvals = ma.masked_equal(dvalues[dscn['NAMEheight']], None)
	881	minz = min(zvals)
	882	maxz = max(zvals)
	883
	884	bz = (maxz - minz)/(Nd-1)
	885
	886	newvar = onc.createVariable( 'z2D', 'f8', ('z2D'))
	887	basicvardef(newvar, 'z2D', 'z-coordinate map observations', zu)
	888	newvar[:] = minz + np.arange(Nd)*bz
	889	newattr = set_attribute(newvar, 'axis', 'Z')
	890
	891	if dbg:
	892	print ' ' + fname + ': zmin=',minz,zu,'zmax=',maxz,zu
	893	print ' ' + fname + ': precission on z-axis=', bz*(Nd-1), zu
	894
	895	map3D = np.ones((Nd, Nd, Nd), dtype=int)*fillValueI
	896	dimt = len(tvals)
	897	Nlost = 0
	898	if dbg: print ' time-step lon lat z ix iy iz passes _______'
	899	for it in range(dimt):
	900	lon = dvalues[dscn['NAMElon']][it]
	901	lat = dvalues[dscn['NAMElat']][it]
	902	z = dvalues[dscn['NAMEheight']][it]
	903	if lon is not None and lat is not None and z is not None:
	904	ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
	905	ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
	906	iz = int((Nd-1)*(z - minz)/(maxz - minz))
	907
	908	if map3D[iz,ilat,ilon] == fillValueI:
	909	map3D[iz,ilat,ilon] = 1
	910	else:
	911	map3D[iz,ilat,ilon] = map3D[iz,ilat,ilon] + 1
	912	if dbg: print it, lon, lat, z, ilon, ilat, iz, map3D[iz,ilat,ilon]
	913
	914	newvar = onc.createVariable( 'trjobs', 'i', ('z2D', 'lat2D', 'lon2D'), \
	915	fill_value = fillValueI)
	916	basicvardef(newvar, 'trajectory_observations', 'number of passes', '-')
	917	newvar[:] = map3D
	918	newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D z2D')
	919
	920	onc.sync()
	921	return
	922
	923	def adding_station_desc(onc,stdesc):
	924	""" Function to add a station description in a netCDF file
	925	onc= netCDF object
[553]	926	stdesc= station description name, lon, lat, height
[259]	927	"""
	928	fname = 'adding_station_desc'
	929
	930	newdim = onc.createDimension('nst',1)
	931
[553]	932	newvar = objfile.createVariable( 'station', 'c', ('nst','StrLength'))
	933	writing_str_nc(newvar, [stdesc[0].replace('!', ' ')], StringLength)
	934
	935	newvar = objfile.createVariable( 'lonstGDM', 'c', ('nst','StrLength'))
	936	Gv = int(stdesc[1])
	937	Dv = int((stdesc[1] - Gv)*60.)
	938	Mv = int((stdesc[1] - Gv - Dv/60.)*3600.)
	939	writing_str_nc(newvar, [str(Gv)+"d" + str(Dv)+"m" + str(Mv)+'s'], StringLength)
	940
[259]	941	if onc.variables.has_key('lon'):
	942	print warnmsg
	943	print ' ' + fname + ": variable 'lon' already exist !!"
	944	print " renaming it as 'lonst'"
	945	lonname = 'lonst'
	946	else:
	947	lonname = 'lon'
	948
	949	newvar = objfile.createVariable( lonname, 'f4', ('nst'))
	950	basicvardef(newvar, lonname, 'longitude', 'degrees_West' )
[553]	951	newvar[:] = stdesc[1]
[259]	952
[553]	953	newvar = objfile.createVariable( 'latstGDM', 'c', ('nst','StrLength'))
	954	Gv = int(stdesc[2])
	955	Dv = int((stdesc[2] - Gv)*60.)
	956	Mv = int((stdesc[2] - Gv - Dv/60.)*3600.)
	957	writing_str_nc(newvar, [str(Gv)+"d" + str(Dv)+"m" + str(Mv)+'s'], StringLength)
	958
[259]	959	if onc.variables.has_key('lat'):
	960	print warnmsg
	961	print ' ' + fname + ": variable 'lat' already exist !!"
	962	print " renaming it as 'latst'"
	963	latname = 'latst'
	964	else:
	965	latname = 'lat'
	966
	967	newvar = objfile.createVariable( latname, 'f4', ('nst'))
	968	basicvardef(newvar, lonname, 'latitude', 'degrees_North' )
[553]	969	newvar[:] = stdesc[2]
[259]	970
	971	if onc.variables.has_key('height'):
	972	print warnmsg
	973	print ' ' + fname + ": variable 'height' already exist !!"
	974	print " renaming it as 'heightst'"
	975	heightname = 'heightst'
	976	else:
	977	heightname = 'height'
	978
	979	newvar = objfile.createVariable( heightname, 'f4', ('nst'))
	980	basicvardef(newvar, heightname, 'height above sea level', 'm' )
[553]	981	newvar[:] = stdesc[3]
[259]	982
	983	return
	984
[553]	985	def oper_values(dvals, opers):
	986	""" Function to operate the values according to given parameters
	987	dvals= datavalues
	988	opers= operations
	989	"""
	990	fname = 'oper_values'
	991
	992	newdvals = {}
	993	varnames = dvals.keys()
	994
	995	aopers = ['sumc','subc','mulc','divc']
	996
	997	Nopers = len(opers)
	998	for iop in range(Nopers):
	999	vn = varnames[iop]
	1000	print vn,'oper:',opers[iop]
	1001	if opers[iop] != '-':
	1002	opern = opers[iop].split(',')[0]
	1003	operv = np.float(opers[iop].split(',')[1])
	1004
	1005	vvals = np.array(dvals[vn])
	1006
	1007	if opern == 'sumc':
	1008	newvals = np.where(vvals is None, None, vvals+operv)
	1009	elif opern == 'subc':
	1010	newvals = np.where(vvals is None, None, vvals-operv)
	1011	elif opern == 'mulc':
	1012	newvals = np.where(vvals is None, None, vvals*operv)
	1013	elif opern == 'divc':
	1014	newvals = np.where(vvals is None, None, vvals/operv)
	1015	else:
	1016	print errormsg
	1017	print ' ' + fname + ": operation '" + opern + "' not ready!!"
	1018	print ' availables:',aopers
	1019	quit(-1)
	1020
	1021	newdvals[vn] = list(newvals)
	1022	else:
	1023	newdvals[vn] = dvals[vn]
	1024
	1025	return newdvals
	1026
[1949]	1027	def WMOcodevar(unitsn, onc, Lstrcode=1024):
	1028	""" Function to add a variabe providing description of a units based on WMO code
	1029	unitsn= name of the units (all WMO codes derived units must be labelled as
	1030	'codeWMO_[num]' associated to a file valled wmo_[num].code)
	1031	onc= netCDF object file to add the description
	1032	Lstrcode= length of description of codes
	1033
	1034	wmo_[num].code must have the structure: ('#' for comments)
	1035	reference:web page, document reference with the code
	1036	short_description:main description of the code
	1037	long_description:long description of the code
	1038	codeTYPE:type of the code (andy of 'D', 'F', 'I', 'S')
	1039	@: Values line giving the start of the values
	1040	[val1]=[meaning of first value]
	1041	(...)
	1042	[valN]=[meaning of last value]
	1043	"""
	1044	fname = 'WMOcodevar'
	1045
	1046	# From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python
	1047	folder = os.path.dirname(os.path.realpath(__file__))
	1048
	1049	code = unitsn.split('_')[1]
	1050
	1051	infile = folder + '/wmo_' + code +'.code'
	1052
	1053	if not os.path.isfile(infile):
	1054	print warnmsg
	1055	print ' ' + fname + ": WMO code file '" + infile + "' does not exist !!"
	1056	return
	1057	# main expected values
	1058	descvals = ['reference', 'short_description', 'long_description', 'codeTYPE', \
	1059	'@']
	1060
	1061	availcodetype = ['D', 'F', 'I', 'S']
	1062
	1063	codevalues = {}
	1064	ocode = open(infile, 'r')
	1065	inivals = False
	1066	codvals = []
	1067	codmeanings = []
	1068	for line in ocode:
	1069	if len(line) > 1 and line[0:1] != '#':
	1070	linev = line.replace('\n','').replace('\t',' ').replace('\r','')
	1071	print inivals, '<>', linev
	1072	if not inivals:
	1073	Sv = linev.split(':')[0]
	1074	Vv = linev.split(':')[1]
	1075	if searchInlist(descvals, Sv):
	1076	if Sv != '@': codevalues[Sv] = Vv
	1077	else: inivals = True
	1078	else:
	1079	Svv = linev.split('=')[0]
	1080	Vvv = linev.split('=')[1]
	1081	print 'Lluis Sv:', Svv, 'Vv:', Vvv
	1082	codvals.append(Svv)
	1083	codmeanings.append(Vvv)
	1084
	1085	# Creating variable
	1086	if not searchInlist(onc.dimensions, 'Lstringcode'):
	1087	print ' ' + fname + ": Adding string length dimension 'Lstringcode' for " + \
	1088	" code descriptions"
	1089	newdim = onc.createDimension('Lstringcode', Lstrcode)
	1090	print 'Lluis Svv:', Svv
	1091	Ncodes = len(codvals)
	1092	codedimn = 'wmo_code_' + str(code)
	1093	if not searchInlist(onc.dimensions, codedimn):
	1094	print ' ' + fname + ": Adding '" + codedimn + "' dimension for code " + \
	1095	"description"
	1096	newdim = onc.createDimension(codedimn, Ncodes)
	1097	onc.sync()
	1098
	1099	# Variable with the value of the code
	1100	if not onc.variables.has_key(codedimn):
	1101	if codevalues['codeTYPE'] == 'D':
	1102	newvar = onc.createVariable(codedimn, 'f8', (codedimn))
	1103	for iv in range(Ncodes): newvar[iv] = np.float64(codvals[iv])
	1104	elif codevalues['codeTYPE'] == 'F':
	1105	newvar = onc.createVariable(codedimn, 'f', (codedimn))
	1106	for iv in range(Ncodes): newvar[iv] = np.float(codvals[iv])
	1107	elif codevalues['codeTYPE'] == 'I':
	1108	newvar = onc.createVariable(codedimn, 'i', (codedimn))
	1109	for iv in range(Ncodes): newvar[iv] = integer(codvals[iv])
	1110	elif codevalues['codeTYPE'] == 'S':
	1111	newvar = onc.createVariable(codedimn, 'c', (codedimn, 'Lstringcode'))
	1112	writing_str_nc(newvar, codvals, Lstrcode)
	1113	else:
	1114	print errormsg
	1115	print ' ' + fname + ": codeTYPE '" + codevalues['codeTYPE'] + "' not" + \
	1116	" ready !!"
	1117	print ' available ones:', availcodetype
	1118	quit(-1)
	1119
	1120	for descv in descvals:
	1121	if descv != '@' and descv != 'codeTYPE':
	1122	newvar.setncattr(descv, codevalues[descv])
	1123	# Variable with the meaning of the code
	1124	if not onc.variables.has_key(codedimn+'_meaning'):
	1125	newvar = onc.createVariable(codedimn+'_meaning','c',(codedimn,'Lstringcode'))
	1126	print 'Lluis len Vc:', len(Vvv), 'Ncodes:', Ncodes
	1127	writing_str_nc(newvar, codmeanings, Lstrcode)
	1128	newvar.setncattr('description', 'meaning of WMO code ' + str(code))
	1129
	1130	onc.sync()
	1131
	1132	return
	1133
[259]	1134	####### ###### ##### #### ### ## #
	1135
	1136	strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " + \
	1137	" and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
	1138
[648]	1139	kindobs=['stations-map','multi-points', 'single-station', 'trajectory']
[259]	1140	strkObs="kind of observations; 'multi-points': multiple individual punctual obs " + \
	1141	"(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
	1142	"'trajectory': following a trajectory"
	1143
	1144	parser = OptionParser()
	1145	parser.add_option("-c", "--comments", dest="charcom",
	1146	help="':', list of characters used for comments", metavar="VALUES")
	1147	parser.add_option("-d", "--descriptionfile", dest="fdesc",
[553]	1148	help="description file to use" + read_description.__doc__, metavar="FILE")
[259]	1149	parser.add_option("-e", "--end_column", dest="endcol",
[1948]	1150	help="character to indicate end of the column ('space', for ' ', or 'None,[fixcoljumps]' for fixed size columns with [fixcoljumps]: ',' separated list of positions of column ends within line)", metavar="VALUE")
[259]	1151	parser.add_option("-f", "--file", dest="obsfile",
	1152	help="observational file to use", metavar="FILE")
	1153	parser.add_option("-g", "--debug", dest="debug",
	1154	help="whther debug is required ('false', 'true')", metavar="VALUE")
	1155	parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs,
[553]	1156	help=strkObs, metavar="VALUE")
[259]	1157	parser.add_option("-s", "--stationLocation", dest="stloc",
[553]	1158	help="name ('!' for spaces), longitude, latitude and height of the station (only for 'single-station')",
[259]	1159	metavar="FILE")
	1160	parser.add_option("-t", "--CFtime", dest="CFtime", help=strCFt, metavar="VALUE")
	1161
	1162	(opts, args) = parser.parse_args()
	1163
	1164	####### #######
	1165	## MAIN
	1166	#######
	1167
	1168	ofile='OBSnetcdf.nc'
	1169
	1170	if opts.charcom is None:
	1171	print warnmsg
	1172	print ' ' + main + ': No list of comment characters provided!!'
	1173	print ' assuming no need!'
	1174	charcomments = []
	1175	else:
	1176	charcomments = opts.charcom.split(':')
	1177
	1178	if opts.fdesc is None:
	1179	print errormsg
	1180	print ' ' + main + ': No description file for the observtional data provided!!'
	1181	quit(-1)
	1182
	1183	if opts.endcol is None:
	1184	print warnmsg
	1185	print ' ' + main + ': No list of comment characters provided!!'
	1186	print " assuming 'space'"
	1187	endcol = ' '
	1188	else:
	1189	if opts.endcol == 'space':
	1190	endcol = ' '
	1191	else:
	1192	endcol = opts.endcol
	1193
	1194	if opts.obsfile is None:
	1195	print errormsg
[332]	1196	print ' ' + main + ': No observations file provided!!'
[259]	1197	quit(-1)
	1198
	1199	if opts.debug is None:
	1200	print warnmsg
	1201	print ' ' + main + ': No debug provided!!'
	1202	print " assuming 'False'"
	1203	debug = False
	1204	else:
	1205	if opts.debug == 'true':
	1206	debug = True
	1207	else:
	1208	debug = False
	1209
	1210	if not os.path.isfile(opts.fdesc):
	1211	print errormsg
	1212	print ' ' + main + ": description file '" + opts.fdesc + "' does not exist !!"
	1213	quit(-1)
	1214
	1215	if not os.path.isfile(opts.obsfile):
	1216	print errormsg
	1217	print ' ' + main + ": observational file '" + opts.obsfile + "' does not exist !!"
	1218	quit(-1)
	1219
	1220	if opts.CFtime is None:
	1221	print warnmsg
	1222	print ' ' + main + ': No CFtime criteria are provided !!'
	1223	print " either time is already in CF-format ('timeFMT=CFtime') in '" + \
	1224	opts.fdesc + "'"
	1225	print " or assuming refdate: '19491201000000' and time units: 'hours'"
	1226	referencedate = '19491201000000'
	1227	timeunits = 'hours'
	1228	else:
	1229	referencedate = opts.CFtime.split(',')[0]
	1230	timeunits = opts.CFtime.split(',')[1]
	1231
	1232	if opts.obskind is None:
	1233	print warnmsg
	1234	print ' ' + main + ': No kind of observations provided !!'
	1235	print " assuming 'single-station': single station on a fixed position at 0,0,0"
	1236	obskind = 'single-station'
[553]	1237	stationdesc = [0.0, 0.0, 0.0, 0.0]
[259]	1238	else:
	1239	obskind = opts.obskind
	1240	if obskind == 'single-station':
	1241	if opts.stloc is None:
	1242	print errornmsg
	1243	print ' ' + main + ': No station location provided !!'
	1244	quit(-1)
	1245	else:
[553]	1246	stvals = opts.stloc.split(',')
	1247	stationdesc = [stvals[0], np.float(stvals[1]), np.float(stvals[2]), \
	1248	np.float(stvals[3])]
[259]	1249	else:
	1250	obskind = opts.obskind
	1251
	1252	# Reading description file
	1253	##
	1254	description = read_description(opts.fdesc, debug)
	1255
	1256	Nvariables=len(description['varN'])
	1257	formats = description['varTYPE']
	1258
	1259	if len(formats) != Nvariables:
	1260	print errormsg
	1261	print ' ' + main + ': number of formats:',len(formats),' and number of ' + \
	1262	'variables', Nvariables,' does not coincide!!'
	1263	print ' * what is found is _______'
	1264	if Nvariables > len(formats):
	1265	Nshow = len(formats)
	1266	for ivar in range(Nshow):
	1267	print ' ',description['varN'][ivar],':', description['varLN'][ivar],\
	1268	'[', description['varU'][ivar], '] fmt:', formats[ivar]
	1269	print ' missing values for:', description['varN'][Nshow:Nvariables+1]
	1270	else:
	1271	Nshow = Nvariables
	1272	for ivar in range(Nshow):
	1273	print ' ',description['varN'][ivar],':', description['varLN'][ivar],\
	1274	'[', description['varU'][ivar], '] fmt:', formats[ivar]
	1275	print ' excess of formats for:', formats[Nshow:len(formats)+1]
	1276
[587]	1277	# quit(-1)
[259]	1278
	1279	# Reading values
	1280	##
	1281	datavalues = read_datavalues(opts.obsfile, charcomments, endcol, formats,
[553]	1282	description['varOPER'], description['MissingValue'], description['varN'], debug)
[259]	1283
	1284	# Total number of values
	1285	Ntvalues = len(datavalues[description['varN'][0]])
[648]	1286	if obskind == 'stations-map':
	1287	print main + ': total values found:',Ntvalues
	1288	else:
	1289	print main + ': total temporal values found:',Ntvalues
[259]	1290
	1291	objfile = NetCDFFile(ofile, 'w')
	1292
	1293	# Creation of dimensions
	1294	##
[648]	1295	if obskind == 'stations-map':
	1296	rowsdim = 'Npoints'
	1297	dimlength = Ntvalues
	1298	else:
	1299	rowsdim = 'time'
	1300	dimlength = None
	1301
	1302	objfile.createDimension(rowsdim,dimlength)
[259]	1303	objfile.createDimension('StrLength',StringLength)
	1304
	1305	# Creation of variables
	1306	##
	1307	for ivar in range(Nvariables):
	1308	varn = description['varN'][ivar]
	1309	print " including: '" + varn + "' ... .. ."
	1310
	1311	if formats[ivar] == 'D':
[648]	1312	newvar = objfile.createVariable(varn, 'f32', (rowsdim), fill_value=fillValueF)
[259]	1313	basicvardef(newvar, varn, description['varLN'][ivar], \
	1314	description['varU'][ivar])
	1315	newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn])
	1316	elif formats[ivar] == 'F':
[648]	1317	newvar = objfile.createVariable(varn, 'f', (rowsdim), fill_value=fillValueF)
[259]	1318	basicvardef(newvar, varn, description['varLN'][ivar], \
	1319	description['varU'][ivar])
	1320	newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn])
	1321	elif formats[ivar] == 'I':
[648]	1322	newvar = objfile.createVariable(varn, 'i', (rowsdim), fill_value=fillValueI)
[259]	1323	basicvardef(newvar, varn, description['varLN'][ivar], \
	1324	description['varU'][ivar])
	1325	# Why is not wotking with integers?
	1326	vals = np.array(datavalues[varn])
	1327	for iv in range(Ntvalues):
	1328	if vals[iv] is None: vals[iv] = fillValueI
	1329	newvar[:] = vals
	1330	elif formats[ivar] == 'S':
[648]	1331	newvar = objfile.createVariable(varn, 'c', (rowsdim,'StrLength'))
[259]	1332	basicvardef(newvar, varn, description['varLN'][ivar], \
	1333	description['varU'][ivar])
[1949]	1334	vals = datavalues[varn]
	1335	for iv in range(Ntvalues):
	1336	if vals[iv] is None: vals[iv] = fillValueS
	1337	writing_str_nc(newvar, vals, StringLength)
[259]	1338
[1949]	1339	# Getting new variables to describe certain units as codeWMO_[num] from an
	1340	# external file
	1341	if description['varU'][ivar][0:7] == 'codeWMO':
	1342	WMOcodevar(description['varU'][ivar], objfile)
	1343
[259]	1344	# Extra variable descriptions/attributes
	1345	if description.has_key('varBUFR'):
	1346	set_attribute(newvar,'bufr_code',description['varBUFR'][ivar])
	1347
	1348	objfile.sync()
	1349
	1350	# Time variable in CF format
	1351	##
	1352	if description['FMTtime'] == 'CFtime':
	1353	timevals = datavalues[description['NAMEtime']]
	1354	iv = 0
	1355	for ivn in description['varN']:
	1356	if ivn == description['NAMEtime']:
	1357	tunits = description['varU'][iv]
	1358	break
	1359	iv = iv + 1
	1360	else:
	1361	# Time as a composition of different columns
	1362	tcomposite = description['NAMEtime'].find('@')
[332]	1363	if tcomposite != -1:
[259]	1364	timevars = description['NAMEtime'].split('@')
	1365
	1366	print warnmsg
	1367	print ' ' + main + ': time values as combination of different columns!'
	1368	print ' combining:',timevars,' with a final format: ',description['FMTtime']
	1369
	1370	timeSvals = []
	1371	if debug: print ' ' + main + ': creating times _______'
	1372	for it in range(Ntvalues):
	1373	tSvals = ''
	1374	for tvar in timevars:
	1375	tSvals = tSvals + datavalues[tvar][it] + ' '
	1376
	1377	timeSvals.append(tSvals[0:len(tSvals)-1])
	1378	if debug: print it, '' + timeSvals[it] + ''
	1379
	1380	timevals = Stringtimes_CF(timeSvals, description['FMTtime'], referencedate, \
	1381	timeunits, debug)
	1382	else:
	1383	timevals = Stringtimes_CF(datavalues[description['NAMEtime']], \
	1384	description['FMTtime'], referencedate, timeunits, debug)
	1385
	1386	CFtimeRef = referencedate[0:4] +'-'+ referencedate[4:6] +'-'+ referencedate[6:8] + \
	1387	' ' + referencedate[8:10] +':'+ referencedate[10:12] +':'+ referencedate[12:14]
	1388	tunits = timeunits + ' since ' + CFtimeRef
	1389
	1390	if objfile.variables.has_key('time'):
	1391	print warnmsg
	1392	print ' ' + main + ": variable 'time' already exist !!"
	1393	print " renaming it as 'CFtime'"
	1394	timeCFname = 'CFtime'
	1395	newdim = objfile.renameDimension('time','CFtime')
	1396	newvar = objfile.createVariable( timeCFname, 'f8', ('CFtime'))
	1397	basicvardef(newvar, timeCFname, 'time', tunits )
	1398	else:
[1949]	1399	if not searchInlist(objfile.dimensions, 'time'):
	1400	newdim = objfile.createDimension('time',None)
[259]	1401	timeCFname = 'time'
[1949]	1402	newvar = objfile.createVariable( timeCFname, 'f8', ('time'))
[1948]	1403	newvar[:] = np.zeros(timevals.shape[0])
[259]	1404	basicvardef(newvar, timeCFname, 'time', tunits )
	1405
	1406	set_attribute(newvar, 'calendar', 'standard')
[648]	1407	if obskind == 'stations-map':
	1408	newvar[:] = timevals[0]
	1409	else:
	1410	newvar[:] = timevals
[259]	1411
	1412	# Global attributes
	1413	##
	1414	for descn in description.keys():
	1415	if descn[0:3] != 'var' and descn[0:4] != 'NAME' and descn[0:3] != 'FMT':
	1416	set_attribute(objfile, descn, description[descn])
	1417
	1418	set_attribute(objfile,'author_nc','Lluis Fita')
	1419	set_attribute(objfile,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \
	1420	'LMD-Jussieu, UPMC, Paris')
	1421	set_attribute(objfile,'country_nc','France')
	1422	set_attribute(objfile,'script_nc','create_OBSnetcdf.py')
	1423	set_attribute(objfile,'version_script',version)
	1424	set_attribute(objfile,'information', \
	1425	'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
	1426
	1427	objfile.sync()
	1428
	1429	# Adding new variables as function of the observational type
	1430	## 'multi-points', 'single-station', 'trajectory'
	1431
	1432	if obskind != 'single-station':
	1433	adding_complementary(objfile, description, obskind, datavalues, timevals, \
	1434	referencedate, timeunits, Ndim2D, debug)
	1435	else:
	1436	# Adding three variables with the station location, longitude, latitude and height
	1437	adding_station_desc(objfile,stationdesc)
	1438
	1439	objfile.sync()
	1440	objfile.close()
	1441
	1442	print main + ": Successfull generation of netcdf observational file '" + ofile + "' !!"

Note: See TracBrowser for help on using the repository browser.

Download in other formats: