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

Last change on this file since 2424 was 2060, checked in by lfita, 6 years ago

Adding reference and name of the station
File size: 13.5 KB
Line 
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 
27from optparse import OptionParser
28import numpy as np
29from netCDF4 import Dataset as NetCDFFile
30import os
31import re
32import numpy.ma as ma
33# Importing generic tools file 'nc_var_tools.py'
34import nc_var_tools as ncvar
35# Importing generic tools file 'generic_tools.py'
36import generic_tools as gen
37import subprocess as sub
38import time
39
40main = 'UWyoming_snd_nc.py'
41
42###### ###### ##### #### ### ## #
43
44# Time reference
45TrefS = '1949-12-01 00:00:00'
46TrefYmdHMS = '19491201000000'
47tunits = 'minutes'
48tfmt = '%y%m%d/%H%M'
49
50# Integer text values (all the rest as float)
51intTXTvals = ['Station number']
52
53# Text values (all the rest as float)
54txtTXTvals = ['Station identifier', 'Observation time']
55
56# Float values for global attributes:
57floatTXTvals = ['Station elevation', 'Station longitude', 'Station latitude']
58
59# Not lower pressure variables
60NOTlowpres = ['DWPT', 'RELH', 'MIXR', 'THTE']
61
62Lstring = 256
63
64# Arguments
65##
66parser = OptionParser()
67parser.add_option("-D", "--Debug", dest="debug", 
68  help="debug prints",  metavar="BOOL")
69parser.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
79# Global attributes for the station
80stngattrk = intTXTvals + txtTXTvals + floatTXTvals
81
82if opts.debug is None:
83    print gen.infmsg
84    print '  ' + main + ": no debug value provided!!"
85    print '    assuming:', False
86    debug = False
87else:
88    debug = gen.Str_Bool(opts.debug)
89
90if 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
96osnd = open(opts.sndfile, 'r')
97
98# Processed sounding
99soundings = {}
100
101# Processed dates
102idate = 0
103
104# List with the dates of  the soundings, to keep them consecutive!!
105Tsoundings = []
106
107for 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
118                strefn = lvals[1]
119                idobsS = lvals.index('Observations')
120                stnameS = ' '.join(lvals[2:idobsS])
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
131                Tsoundings.append(stationdateS)
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
213stationdateS = str(txtvals['Station number']) + '_' + txtvals['Observation time']
214print "   data for sounding: '" + stationdateS + "' _______"
215
216if (debug):
217    for ip in pvals.keys():
218        print ip, ':', pvals[ip]
219    for ic in txtvals.keys():
220        print ic, ':', txtvals[ic]
221
222Tsoundings.append(stationdateS)
223soundings[stationdateS] = [pvals, txtvals]
224if len(soundings.keys()) == 1:
225    presvals = list(pvals.keys())
226else:
227    noinc = list(set(pvals.keys()).difference(set(presvals)))
228    presvals = presvals + noinc
229    presvals.sort(reverse=True)
230
231osnd.close()
232varns = list(sndvarn)
233varns.remove('PRES')
234
235# Getting measured values
236Ntimes = len(Tsoundings)
237Npres = len(presvals)
238Nvals = len(varns)
239
240sndvals = np.ones((Ntimes, Npres, Nvals), dtype=np.float)*gen.fillValueF
241for it in range(Ntimes):
242    snditS = Tsoundings[it]
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
267sndvals = ma.masked_equal(sndvals, gen.fillValueF)
268print '  recupered values shape:', sndvals.shape
269if (debug):
270    print '    values from file _______'
271    print sndvals
272
273# Removing not computed values
274statglobalattr = {}
275txtn = list(txtvals.keys())
276snditS = Tsoundings[0]
277sndvs = soundings[snditS]
278sndc = sndvs[1]
279
280for Sn in intTXTvals:
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
290for Sn in txtTXTvals:
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] = '-'
299
300for 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
310Ncomp = len(txtn)
311compvals = np.ones((Ntimes, Ncomp), dtype=np.float)*gen.fillValueF
312
313tvals = []
314for it in range(Ntimes):
315    snditS = Tsoundings[it]
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##
325ofilen = 'UWyoming_snd_' + str(txtvals['Station number']) + '.nc'
326
327onewnc = NetCDFFile(ofilen, 'w')
328
329# Creation of dimensions
330onewnc.createDimension('pres', Npres)
331onewnc.createDimension('time', None)
332#onewnc.createDimension('cvals', Ncomp)
333#onewnc.createDimension('Lstring', Lstring)
334
335# Creation of variable-dimensions
336newvar = onewnc.createVariable('pres', 'f8', ('pres'))
337newvar[:] = presvals
338ncvar.basicvardef(newvar, 'pres', 'pressure', varu['PRES'])
339
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')
344
345newvar = onewnc.createVariable('time', 'f8', ('time'))
346timeT = []
347atimeT = np.zeros((len(tvals),6), dtype=int)
348iit = 0
349for 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
356CFtimes = gen.realdatetime_CFcompilant(atimeT, TrefYmdHMS, tunits)
357newvar[:] = CFtimes
358ncvar.basicvardef(newvar, 'time', 'time', tunits + ' since ' + TrefS)
359ncvar.set_attribute(newvar, 'calendar', 'gregorian')
360
361# Filling with sounding values
362iv = 0
363for 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
370onewnc.sync()
371
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()
378
379# Getting specific 1D values
380for 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
388# Global attributes
389ncvar.add_global_PyNCplot(onewnc, main, '', '0.2')
390onewnc.setncattr('Station_ref', strefn)
391onewnc.setncattr('Station_name', stnameS)
392
393for atn in stngattrk:
394    atnS = atn.replace(' ','_')
395    onewnc.setncattr(atnS, statglobalattr[atnS])
396
397onewnc.close()
398print main + ": succesful writing of sounding file '" + ofilen + "' !!"
399
Note: See TracBrowser for help on using the repository browser.