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

Last change on this file since 2828 was 2738, checked in by lfita, 5 years ago

Making use of `linvals'
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            if gen.IsNumber(lvals[0], 'R') and len(lvals) > 1 and lvals[1] != 'hPa':
163                Sline = line.replace('\n', '').replace('\t', '')
164                linvals = gen.getting_fixedline(Sline,[7,14,21,28,35,42,49,56,63,70],\
165                  ['R', 'R', 'R', 'R', 'R', 'R', 'R', 'R', 'R', 'R', 'R'])
166                if (debug):
167                    print '      getting values at pressure:', np.float(lvals[0])
168                pvals[np.float(lvals[0])] = np.array(linvals[1:], dtype=np.float)
169            else:
170                if (debug):
171                    print '      getting text values ...'
172                if lvals[0] == 'hPa':
173                    if idate == 1:
174                        if (debug):
175                            print '      getting units of the variables'
176                        varu = {}
177                        iv = 0
178                        for vn in sndvarn:
179                            varu[vn] = lvals[iv]
180                            iv = iv + 1
181                else:
182                    txt = ''
183                    # Specific work for 'Observation time:'
184                    if gen.searchInlist(lvals, 'Observation'):
185                        txtvals['Observation time'] = lvals[2]
186                    else:
187                        for iv in range(len(lvals)):
188                            if not gen.IsNumber(lvals[iv], 'R'):
189                                if len(txt) == 0:
190                                    txt = lvals[iv]
191                                else:
192                                    txt = txt + ' ' + lvals[iv]
193                            else:
194                                # Specific work for '1000 hPa to 500 hPa thickness:'
195                                if gen.searchInlist(lvals, 'thickness:') and iv <= 5:
196                                    if len(txt) == 0:
197                                        txt = lvals[iv]
198                                    else:
199                                        txt = txt + ' ' + lvals[iv]
200                                else:
201                                    txtn = txt.replace(':', '')
202                                    if gen.searchInlist(intTXTvals, txtn):
203                                        txtvals[txtn] = int(lvals[iv])
204                                    elif gen.searchInlist(txtTXTvals, txtn):
205                                        txtvals[txtn] = lvals[iv]
206                                    else:
207                                        txtvals[txtn] = np.float(lvals[iv])
208                        if (debug):
209                            print "        text value: '" + txtn + "':", txtvals[txtn]
210
211# Including last sounding
212stationdateS = str(txtvals['Station number']) + '_' + txtvals['Observation time']
213print "   data for sounding: '" + stationdateS + "' _______"
214
215if (debug):
216    for ip in pvals.keys():
217        print ip, ':', pvals[ip]
218    for ic in txtvals.keys():
219        print ic, ':', txtvals[ic]
220
221Tsoundings.append(stationdateS)
222soundings[stationdateS] = [pvals, txtvals]
223if len(soundings.keys()) == 1:
224    presvals = list(pvals.keys())
225else:
226    noinc = list(set(pvals.keys()).difference(set(presvals)))
227    presvals = presvals + noinc
228    presvals.sort(reverse=True)
229
230osnd.close()
231varns = list(sndvarn)
232varns.remove('PRES')
233
234# Getting measured values
235Ntimes = len(Tsoundings)
236Npres = len(presvals)
237Nvals = len(varns)
238
239sndvals = np.ones((Ntimes, Npres, Nvals), dtype=np.float)*gen.fillValueF
240for it in range(Ntimes):
241    snditS = Tsoundings[it]
242    sndvs = soundings[snditS]
243    sndv = sndvs[0]
244    for ip in range(Npres):
245        if sndv.has_key(presvals[ip]):
246            sndpv = sndv[presvals[ip]]
247
248            if len(sndpv) == Nvals:
249                sndvals[it,ip,:] = sndpv
250            else:
251                ivv = 0
252                if len(sndpv) == Nvals - len(NOTlowpres):
253                    for iv in range(Nvals):
254                        if not gen.searchInlist(NOTlowpres, varns[iv]):
255                            sndvals[it,ip,iv] = sndpv[ivv]
256                            ivv = ivv + 1
257                elif len(sndpv) == Nvals - len(NOTlowpres) - 2:
258                    for iv in range(Nvals):
259                        if not gen.searchInlist(NOTlowpres, varns[iv]) and         \
260                          not gen.searchInlist(['DRCT', 'SKNT'], varns[iv]):
261                            sndvals[it,ip,iv] = sndpv[ivv]
262                            ivv = ivv + 1
263                elif len(sndpv) == 1:
264                    sndvals[it,ip,ivv] = sndpv[0]
265
266sndvals = ma.masked_equal(sndvals, gen.fillValueF)
267print '  recupered values shape:', sndvals.shape
268if (debug):
269    print '    values from file _______'
270    print sndvals
271
272# Removing not computed values
273statglobalattr = {}
274txtn = list(txtvals.keys())
275snditS = Tsoundings[0]
276sndvs = soundings[snditS]
277sndc = sndvs[1]
278
279for Sn in intTXTvals:
280    SnS = Sn.replace(' ','_')
281    if gen.searchInlist(txtn, Sn):
282        if sndc.has_key(Sn): statglobalattr[SnS] = int(sndc[Sn])
283        else: statglobalattr[SnS] = '-'
284        txtn.remove(Sn)
285    else:
286        if sndc.has_key(Sn): statglobalattr[SnS] = int(sndc[Sn])
287        else: statglobalattr[SnS] = '-'
288
289for Sn in txtTXTvals:
290    SnS = Sn.replace(' ','_')
291    if gen.searchInlist(txtn, Sn):
292        if sndc.has_key(Sn): statglobalattr[SnS] = sndc[Sn]
293        else: statglobalattr[SnS] = '-'
294        txtn.remove(Sn)
295    else:
296        if sndc.has_key(Sn): statglobalattr[SnS] = int(sndc[Sn])
297        else: statglobalattr[SnS] = '-'
298
299for Sn in floatTXTvals:
300    SnS = Sn.replace(' ','_')
301    if gen.searchInlist(txtn, Sn): 
302        if sndc.has_key(Sn): statglobalattr[SnS] = np.float(sndc[Sn])
303        else: statglobalattr[Sn] = '-'
304        txtn.remove(Sn)
305    else:
306        if sndc.has_key(Sn): statglobalattr[SnS] = int(sndc[Sn])
307        else: statglobalattr[SnS] = '-'
308
309Ncomp = len(txtn)
310compvals = np.ones((Ntimes, Ncomp), dtype=np.float)*gen.fillValueF
311
312tvals = []
313for it in range(Ntimes):
314    snditS = Tsoundings[it]
315    sndvs = soundings[snditS]
316    sndc = sndvs[1]
317    tvals.append(sndc['Observation time'])
318    for ic in range(len(txtn)):
319        if sndc.has_key(txtn[ic]):
320            compvals[it,ic] = sndc[txtn[ic]]
321
322# netCDF file creation
323##
324ofilen = 'UWyoming_snd_' + str(txtvals['Station number']) + '.nc'
325
326onewnc = NetCDFFile(ofilen, 'w')
327
328# Creation of dimensions
329onewnc.createDimension('pres', Npres)
330onewnc.createDimension('time', None)
331#onewnc.createDimension('cvals', Ncomp)
332#onewnc.createDimension('Lstring', Lstring)
333
334# Creation of variable-dimensions
335newvar = onewnc.createVariable('pres', 'f8', ('pres'))
336newvar[:] = presvals
337ncvar.basicvardef(newvar, 'pres', 'pressure', varu['PRES'])
338
339# No more computedvalues matrix !
340#newvar = onewnc.createVariable('cvals', 'c', ('cvals', 'Lstring'))
341#ncvar.writing_str_nc(newvar, txtn, Lstring)
342#ncvar.basicvardef(newvar, 'cvals', 'computed values from sounding data', 'hPa')
343
344newvar = onewnc.createVariable('time', 'f8', ('time'))
345timeT = []
346atimeT = np.zeros((len(tvals),6), dtype=int)
347iit = 0
348for it in tvals:
349    dateT = time.strptime(it, tfmt)
350    timeT.append(dateT)
351    atimeT[iit,:] = np.array([dateT.tm_year, dateT.tm_mon, dateT.tm_mday,            \
352      dateT.tm_hour, dateT.tm_min, dateT.tm_sec])
353    iit = iit + 1
354
355CFtimes = gen.realdatetime_CFcompilant(atimeT, TrefYmdHMS, tunits)
356newvar[:] = CFtimes
357ncvar.basicvardef(newvar, 'time', 'time', tunits + ' since ' + TrefS)
358ncvar.set_attribute(newvar, 'calendar', 'gregorian')
359
360# Filling with sounding values
361iv = 0
362for varn in varns:
363    CFvals = gen.variables_values(varn)
364    newvar = onewnc.createVariable(CFvals[0], 'f', ('time', 'pres'),                 \
365      fill_value=gen.fillValueF)
366    newvar[:] = sndvals[:,:,iv]
367    ncvar.basicvardef(newvar, varn, CFvals[4].replace('|', ' '), varu[varn])
368    iv = iv + 1
369onewnc.sync()
370
371# No more computedvalues matrix !
372#newvar = onewnc.createVariable('computedvals', 'f', ('time', 'cvals'),               \
373#  fill_value=gen.fillValueF)
374#newvar[:] = compvals
375#ncvar.basicvardef(newvar, 'computedvals', 'values computed from sounding data', '-')
376#onewnc.sync()
377
378# Getting specific 1D values
379for ic in range(len(txtn)):
380    if len(txtn[ic]) < 2: continue
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.