Context Navigation

source: lmdz_wrf/trunk/tools/UWyoming_snd_nc.py @ 2424

Last change on this file since 2424 was 2060, checked in by lfita, 7 years ago
Adding reference and name of the station
File size: 13.5 KB

Rev	Line
[1926]	1	# Python tools to transform from U. Wyoming sounding file to netCDF
	2	# http://weather.uwyo.edu/upperair/sounding.html
	3	# L. Fita, CIMA August 2017
	4	##
	5	## e.g. # UWyoming_snd_nc.py -f snd_CORDOBA.txt
	6
	7	### ASCII files from U. Wyoming can have more than one sounding per file, but from
	8	### the same station.
	9	### Script spect to find
	10	# [num] [Text] Observations at [time]Z [Date]
	11	#
	12	#-----------------------------------------------------------------------------
	13	# PRES HGHT TEMP DWPT RELH MIXR DRCT SKNT THTA THTE THTV
	14	# hPa m C C % g/kg deg knot K K K
	15	#-----------------------------------------------------------------------------
	16	# lines of data
	17	# Station information and sounding indices
	18
	19	# Station identifier: [Value]
	20	# Station number: [Value]
	21	# Observation time: [date]/[time]
	22	# (...)
	23	#
	24	# [num] [Text] Observations at [time]Z [Date]
	25	# (...)
	26
	27	from optparse import OptionParser
	28	import numpy as np
	29	from netCDF4 import Dataset as NetCDFFile
	30	import os
	31	import re
	32	import numpy.ma as ma
	33	# Importing generic tools file 'nc_var_tools.py'
	34	import nc_var_tools as ncvar
	35	# Importing generic tools file 'generic_tools.py'
	36	import generic_tools as gen
	37	import subprocess as sub
	38	import time
	39
	40	main = 'UWyoming_snd_nc.py'
	41
	42	###### ###### ##### #### ### ## #
	43
	44	# Time reference
	45	TrefS = '1949-12-01 00:00:00'
	46	TrefYmdHMS = '19491201000000'
	47	tunits = 'minutes'
	48	tfmt = '%y%m%d/%H%M'
	49
	50	# Integer text values (all the rest as float)
	51	intTXTvals = ['Station number']
	52
	53	# Text values (all the rest as float)
	54	txtTXTvals = ['Station identifier', 'Observation time']
	55
[1941]	56	# Float values for global attributes:
	57	floatTXTvals = ['Station elevation', 'Station longitude', 'Station latitude']
	58
[1926]	59	# Not lower pressure variables
	60	NOTlowpres = ['DWPT', 'RELH', 'MIXR', 'THTE']
	61
	62	Lstring = 256
	63
	64	# Arguments
	65	##
	66	parser = OptionParser()
	67	parser.add_option("-D", "--Debug", dest="debug",
	68	help="debug prints", metavar="BOOL")
	69	parser.add_option("-f", "--snd_file", dest="sndfile",
	70	help="U. Wyoming sounding (http://weather.uwyo.edu/upperair/sounding.html) file to use",
	71	metavar="FILE")
	72
	73	(opts, args) = parser.parse_args()
	74
	75	####### #######
	76	## MAIN
	77	#######
	78
[1941]	79	# Global attributes for the station
	80	stngattrk = intTXTvals + txtTXTvals + floatTXTvals
	81
[1926]	82	if opts.debug is None:
	83	print gen.infmsg
	84	print ' ' + main + ": no debug value provided!!"
	85	print ' assuming:', False
	86	debug = False
	87	else:
	88	debug = gen.Str_Bool(opts.debug)
	89
	90	if not os.path.isfile(opts.sndfile):
	91	print gen.errormsg
	92	print ' ' + main + ": sounding file '" + opts.sndfile + "' does not exist !!"
	93	quit(-1)
	94
	95	# Reading sounding file
	96	osnd = open(opts.sndfile, 'r')
	97
	98	# Processed sounding
	99	soundings = {}
	100
	101	# Processed dates
	102	idate = 0
	103
[1938]	104	# List with the dates of the soundings, to keep them consecutive!!
	105	Tsoundings = []
	106
[1926]	107	for line in osnd:
	108	if line[0:1] != '#' and line[0:1] != '<' and line[0:1] != '-' and len(line) > 1:
	109	lvals = gen.values_line(line, ' ', ['\t'])
	110	newsnd = lvals.count('Observations')
	111	# Processing a new date
	112	if newsnd != 0:
	113	if idate == 0:
	114	print ' reading first sounding!'
	115	pvals = {}
	116	txtvals = {}
	117	idate = idate + 1
[2060]	118	strefn = lvals[1]
	119	idobsS = lvals.index('Observations')
	120	stnameS = ' '.join(lvals[2:idobsS])
[1926]	121	else:
	122	stationdateS = str(txtvals['Station number']) + '_' + txtvals['Observation time']
	123	print " data for sounding: '" + stationdateS + "' _______"
	124
	125	if (debug):
	126	for ip in pvals.keys():
	127	print ip, ':', pvals[ip]
	128	for ic in txtvals.keys():
	129	print ic, ':', txtvals[ic]
	130
[1938]	131	Tsoundings.append(stationdateS)
[1926]	132	soundings[stationdateS] = [pvals, txtvals]
	133	if len(soundings.keys()) == 1:
	134	presvals = list(pvals.keys())
	135	else:
	136	noinc = list(set(pvals.keys()).difference(set(presvals)))
	137	presvals = presvals + noinc
	138	presvals.sort(reverse=True)
	139	pvals = {}
	140	txtvals = {}
	141	idate = idate + 1
	142	elif lvals[0] == 'hPa':
	143	if (debug):
	144	print ' getting text values ...'
	145	if lvals[0] == 'hPa':
	146	if idate == 1:
	147	if (debug):
	148	print ' getting units of the variables'
	149	varu = {}
	150	iv = 0
	151	for vn in sndvarn:
	152	varu[vn] = lvals[iv]
	153	iv = iv + 1
	154	elif lvals[0] == '-':
	155	print ' '
	156	# elif lvals[0].index() ==
	157	elif lvals[0] == 'PRES':
	158	# doing nothing!
	159	sndvarn = list(lvals)
	160	else:
	161	# Numeric values
	162
	163	if gen.IsNumber(lvals[0], 'R') and lvals[1] != 'hPa':
	164	Sline = line.replace('\n', '').replace('\t', '')
	165	linvals = gen.getting_fixedline(Sline,[7,14,21,28,35,42,49,56,63,70],\
	166	['R', 'R', 'R', 'R', 'R', 'R', 'R', 'R', 'R', 'R', 'R'])
	167	if (debug):
	168	print ' getting values at pressure:', np.float(lvals[0])
	169	pvals[np.float(lvals[0])] = np.array(lvals[1:], dtype=np.float)
	170	else:
	171	if (debug):
	172	print ' getting text values ...'
	173	if lvals[0] == 'hPa':
	174	if idate == 1:
	175	if (debug):
	176	print ' getting units of the variables'
	177	varu = {}
	178	iv = 0
	179	for vn in sndvarn:
	180	varu[vn] = lvals[iv]
	181	iv = iv + 1
	182	else:
	183	txt = ''
	184	# Specific work for 'Observation time:'
	185	if gen.searchInlist(lvals, 'Observation'):
	186	txtvals['Observation time'] = lvals[2]
	187	else:
	188	for iv in range(len(lvals)):
	189	if not gen.IsNumber(lvals[iv], 'R'):
	190	if len(txt) == 0:
	191	txt = lvals[iv]
	192	else:
	193	txt = txt + ' ' + lvals[iv]
	194	else:
	195	# Specific work for '1000 hPa to 500 hPa thickness:'
	196	if gen.searchInlist(lvals, 'thickness:') and iv <= 5:
	197	if len(txt) == 0:
	198	txt = lvals[iv]
	199	else:
	200	txt = txt + ' ' + lvals[iv]
	201	else:
	202	txtn = txt.replace(':', '')
	203	if gen.searchInlist(intTXTvals, txtn):
	204	txtvals[txtn] = int(lvals[iv])
	205	elif gen.searchInlist(txtTXTvals, txtn):
	206	txtvals[txtn] = lvals[iv]
	207	else:
	208	txtvals[txtn] = np.float(lvals[iv])
	209	if (debug):
	210	print " text value: '" + txtn + "':", txtvals[txtn]
	211
	212	# Including last sounding
	213	stationdateS = str(txtvals['Station number']) + '_' + txtvals['Observation time']
	214	print " data for sounding: '" + stationdateS + "' _______"
	215
	216	if (debug):
	217	for ip in pvals.keys():
	218	print ip, ':', pvals[ip]
	219	for ic in txtvals.keys():
	220	print ic, ':', txtvals[ic]
	221
[1941]	222	Tsoundings.append(stationdateS)
[1926]	223	soundings[stationdateS] = [pvals, txtvals]
	224	if len(soundings.keys()) == 1:
	225	presvals = list(pvals.keys())
	226	else:
	227	noinc = list(set(pvals.keys()).difference(set(presvals)))
	228	presvals = presvals + noinc
	229	presvals.sort(reverse=True)
	230
	231	osnd.close()
	232	varns = list(sndvarn)
	233	varns.remove('PRES')
	234
	235	# Getting measured values
[1938]	236	Ntimes = len(Tsoundings)
[1926]	237	Npres = len(presvals)
	238	Nvals = len(varns)
	239
	240	sndvals = np.ones((Ntimes, Npres, Nvals), dtype=np.float)*gen.fillValueF
	241	for it in range(Ntimes):
[1938]	242	snditS = Tsoundings[it]
[1926]	243	sndvs = soundings[snditS]
	244	sndv = sndvs[0]
	245	for ip in range(Npres):
	246	if sndv.has_key(presvals[ip]):
	247	sndpv = sndv[presvals[ip]]
	248
	249	if len(sndpv) == Nvals:
	250	sndvals[it,ip,:] = sndpv
	251	else:
	252	ivv = 0
	253	if len(sndpv) == Nvals - len(NOTlowpres):
	254	for iv in range(Nvals):
	255	if not gen.searchInlist(NOTlowpres, varns[iv]):
	256	sndvals[it,ip,iv] = sndpv[ivv]
	257	ivv = ivv + 1
	258	elif len(sndpv) == Nvals - len(NOTlowpres) - 2:
	259	for iv in range(Nvals):
	260	if not gen.searchInlist(NOTlowpres, varns[iv]) and \
	261	not gen.searchInlist(['DRCT', 'SKNT'], varns[iv]):
	262	sndvals[it,ip,iv] = sndpv[ivv]
	263	ivv = ivv + 1
	264	elif len(sndpv) == 1:
	265	sndvals[it,ip,ivv] = sndpv[0]
	266
	267	sndvals = ma.masked_equal(sndvals, gen.fillValueF)
	268	print ' recupered values shape:', sndvals.shape
	269	if (debug):
	270	print ' values from file _______'
	271	print sndvals
	272
[1941]	273	# Removing not computed values
	274	statglobalattr = {}
[1926]	275	txtn = list(txtvals.keys())
[1941]	276	snditS = Tsoundings[0]
	277	sndvs = soundings[snditS]
	278	sndc = sndvs[1]
	279
[1926]	280	for Sn in intTXTvals:
[1941]	281	SnS = Sn.replace(' ','_')
	282	if gen.searchInlist(txtn, Sn):
	283	if sndc.has_key(Sn): statglobalattr[SnS] = int(sndc[Sn])
	284	else: statglobalattr[SnS] = '-'
	285	txtn.remove(Sn)
	286	else:
	287	if sndc.has_key(Sn): statglobalattr[SnS] = int(sndc[Sn])
	288	else: statglobalattr[SnS] = '-'
	289
[1926]	290	for Sn in txtTXTvals:
[1941]	291	SnS = Sn.replace(' ','_')
	292	if gen.searchInlist(txtn, Sn):
	293	if sndc.has_key(Sn): statglobalattr[SnS] = sndc[Sn]
	294	else: statglobalattr[SnS] = '-'
	295	txtn.remove(Sn)
	296	else:
	297	if sndc.has_key(Sn): statglobalattr[SnS] = int(sndc[Sn])
	298	else: statglobalattr[SnS] = '-'
[1926]	299
[1941]	300	for Sn in floatTXTvals:
	301	SnS = Sn.replace(' ','_')
	302	if gen.searchInlist(txtn, Sn):
	303	if sndc.has_key(Sn): statglobalattr[SnS] = np.float(sndc[Sn])
	304	else: statglobalattr[Sn] = '-'
	305	txtn.remove(Sn)
	306	else:
	307	if sndc.has_key(Sn): statglobalattr[SnS] = int(sndc[Sn])
	308	else: statglobalattr[SnS] = '-'
	309
[1926]	310	Ncomp = len(txtn)
	311	compvals = np.ones((Ntimes, Ncomp), dtype=np.float)*gen.fillValueF
	312
	313	tvals = []
	314	for it in range(Ntimes):
[1938]	315	snditS = Tsoundings[it]
[1926]	316	sndvs = soundings[snditS]
	317	sndc = sndvs[1]
	318	tvals.append(sndc['Observation time'])
	319	for ic in range(len(txtn)):
	320	if sndc.has_key(txtn[ic]):
	321	compvals[it,ic] = sndc[txtn[ic]]
	322
	323	# netCDF file creation
	324	##
	325	ofilen = 'UWyoming_snd_' + str(txtvals['Station number']) + '.nc'
	326
	327	onewnc = NetCDFFile(ofilen, 'w')
	328
	329	# Creation of dimensions
	330	onewnc.createDimension('pres', Npres)
	331	onewnc.createDimension('time', None)
[1941]	332	#onewnc.createDimension('cvals', Ncomp)
	333	#onewnc.createDimension('Lstring', Lstring)
[1926]	334
	335	# Creation of variable-dimensions
	336	newvar = onewnc.createVariable('pres', 'f8', ('pres'))
	337	newvar[:] = presvals
	338	ncvar.basicvardef(newvar, 'pres', 'pressure', varu['PRES'])
	339
[1941]	340	# No more computedvalues matrix !
	341	#newvar = onewnc.createVariable('cvals', 'c', ('cvals', 'Lstring'))
	342	#ncvar.writing_str_nc(newvar, txtn, Lstring)
	343	#ncvar.basicvardef(newvar, 'cvals', 'computed values from sounding data', 'hPa')
[1926]	344
	345	newvar = onewnc.createVariable('time', 'f8', ('time'))
	346	timeT = []
	347	atimeT = np.zeros((len(tvals),6), dtype=int)
	348	iit = 0
	349	for it in tvals:
	350	dateT = time.strptime(it, tfmt)
	351	timeT.append(dateT)
	352	atimeT[iit,:] = np.array([dateT.tm_year, dateT.tm_mon, dateT.tm_mday, \
	353	dateT.tm_hour, dateT.tm_min, dateT.tm_sec])
	354	iit = iit + 1
	355
	356	CFtimes = gen.realdatetime_CFcompilant(atimeT, TrefYmdHMS, tunits)
	357	newvar[:] = CFtimes
	358	ncvar.basicvardef(newvar, 'time', 'time', tunits + ' since ' + TrefS)
	359	ncvar.set_attribute(newvar, 'calendar', 'gregorian')
	360
	361	# Filling with sounding values
	362	iv = 0
	363	for varn in varns:
	364	CFvals = gen.variables_values(varn)
	365	newvar = onewnc.createVariable(CFvals[0], 'f', ('time', 'pres'), \
	366	fill_value=gen.fillValueF)
	367	newvar[:] = sndvals[:,:,iv]
	368	ncvar.basicvardef(newvar, varn, CFvals[4].replace('\|', ' '), varu[varn])
	369	iv = iv + 1
	370	onewnc.sync()
	371
[1941]	372	# No more computedvalues matrix !
	373	#newvar = onewnc.createVariable('computedvals', 'f', ('time', 'cvals'), \
	374	# fill_value=gen.fillValueF)
	375	#newvar[:] = compvals
	376	#ncvar.basicvardef(newvar, 'computedvals', 'values computed from sounding data', '-')
	377	#onewnc.sync()
[1926]	378
[1941]	379	# Getting specific 1D values
	380	for ic in range(len(txtn)):
	381	Sn = txtn[ic]
	382	CFvalues = gen.variables_values(Sn)
	383	newvar=onewnc.createVariable(CFvalues[0],'f', ('time'), fill_value=gen.fillValueF)
	384	newvar[:] = compvals[:,ic]
	385	ncvar.basicvardef(newvar, CFvalues[1], CFvalues[4].replace('\|',' '), CFvalues[5])
	386	onewnc.sync()
	387
[1926]	388	# Global attributes
[1941]	389	ncvar.add_global_PyNCplot(onewnc, main, '', '0.2')
[2060]	390	onewnc.setncattr('Station_ref', strefn)
	391	onewnc.setncattr('Station_name', stnameS)
	392
[1941]	393	for atn in stngattrk:
	394	atnS = atn.replace(' ','_')
	395	onewnc.setncattr(atnS, statglobalattr[atnS])
[1926]	396
	397	onewnc.close()
	398	print main + ": succesful writing of sounding file '" + ofilen + "' !!"
	399

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

Download in other formats: