source: lmdz_wrf/trunk/tools/TS_ASCII_netCDF.py @ 1906

Last change on this file since 1906 was 1219, checked in by lfita, 8 years ago

Changing url of web page to remove the accent
Adding warranty and license to TS_ASCII_netCDF.py
File size: 14.5 KB
Line 
1# Python script to transfomr ASCII Time-series WRF outputs to netCDF
2# From L. Fita work in different places: LMD (France)
3# More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot
4#
5# pyNCplot and its component TS_ASCII_netCDF.py comes with ABSOLUTELY NO WARRANTY.
6# This work is licendes under a Creative Commons
7#   Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0)
8#
9
10## g.e. # TS_ASCII_netCDF.py -f //media/ExtDiskD/bkup/ciclad/etudes/WL_HyMeX/iop15/wrf/run/control/stations_20121018000000-20121022000000/h0001.d01.TS -s 20121018000000
11
12import os
13from optparse import OptionParser
14import numpy as np
15from netCDF4 import Dataset as NetCDFFile
16import re
17
18main = 'TS_ASCII_netCDF.py'
19errormsg='ERROR -- error -- ERROR -- error'
20warnmsg='WARNING -- warning -- WARNING -- warning'
21
22fillValue = 1.e20
23
24def searchInlist(listname, nameFind):
25    """ Function to search a value within a list
26    listname = list
27    nameFind = value to find
28    >>> searInlist(['1', '2', '3', '5'], '5')
29    True
30    """
31    for x in listname:
32      if x == nameFind:
33        return True
34    return False
35
36def reduce_spaces(string):
37    """ Function to give words of a line of text removing any extra space
38    """
39    values = string.replace('\n','').split(' ')
40    vals = []
41    for val in values:
42         if len(val) > 0:
43             vals.append(val)
44
45    return vals
46
47def reduce_last_spaces(string):
48    """ Function to reduce the last right spaces from a string
49      string= string to remove the spaces at the righ side of the string
50    >>> reduce_last_spaces('Hola Lluis      ')
51    'Hola Lluis'
52    """
53    fname = 'reduce_last_spaces'
54
55    Lstring = len(string)
56
57    firstspace = -1
58    for istr in range(Lstring-1,0,-1):
59        if string[istr:istr+1] == ' ':
60            firstspace = istr
61        else:
62            break
63
64    if firstspace == -1:
65        print errormsg
66        print '  ' + fname + ": string '" + string + "' is not ended by spaces!!"
67        print '    char. number of last character:',ord(string[Lstring:Lstring+1])
68        quit(-1)
69
70    newstr = string[0:firstspace]
71
72    return newstr
73
74def set_attribute(ncv, attrname, attrvalue):
75    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
76    ncv = object netcdf variable
77    attrname = name of the attribute
78    attrvalue = value of the attribute
79    """
80
81    attvar = ncv.ncattrs()
82    if searchInlist(attvar, attrname):
83        attr = ncv.delncattr(attrname)
84
85    attr = ncv.setncattr(attrname, attrvalue)
86
87    return attr
88
89def rmNOnum(string):
90    """ Removing from a string all that characters which are not numbers
91    # From: http://stackoverflow.com/questions/4289331/python-extract-numbers-from-a-string
92    """
93    fname = 'rmNOnum'
94
95    newstr = str(re.findall("[-+]?\d+[\.]?\d*", string)[0])
96
97    return newstr
98
99def basicvardef(varobj, vstname, vlname, vunits):
100    """ Function to give the basic attributes to a variable
101    varobj= netCDF variable object
102    vstname= standard name of the variable
103    vlname= long name of the variable
104    vunits= units of the variable
105    """
106    attr = varobj.setncattr('standard_name', vstname)
107    attr = varobj.setncattr('long_name', vlname)
108    attr = varobj.setncattr('units', vunits)
109
110    return
111
112def values_fortran_fmt(lin,fFMT):
113    """ Function to give back the values of an ASCII line according to its fortran printing format
114      lin= ASCII line
115      fFMT= list with the fortran FORMAT formats
116    forline = 'Natchitoches (RGNL)        1 11 0011  ( 31.733, -93.100) (  28, 157) ( 31.761, -93.113)   41.2 meters'
117    formats = ['A26', 'I2', 'I3', 'A6', 'A2', 'F7.3', 'A1', 'F8.3', 'A3', 'I4', 'A1', 'I4',
118      'A3', 'F7.3', 'A1', 'F8.3', 'A2', 'F6.1', 'A7']
119    >>> values_fortran_fmt(forline, formats)
120    ['Natchitoches (RGNL)        ', 1, 11, ' 0011  ', ' ( ', 31.733, ', ', -93.1, ') ( ', 28, ', ', 157, ')
121      ( ', 31.761, ', ', -93.113, ')  ', 41.2, ' meters']
122    """
123    fname = 'values_fortran_fmt'
124
125    afmts = ['A', 'D', 'F', 'I']
126
127    if lin == 'h':
128        print fname + '_____________________________________________________________'
129        print values_fortran_fmt.__doc__
130        quit()
131
132    fvalues = []
133    ichar=0
134    ival = 0
135    for ft in fFMT:
136        Lft = len(ft)
137        if ft[0:1] == 'A' or ft[0:1] == 'a':
138            Lfmt = int(ft[1:Lft+1])
139            fvalues.append(lin[ichar:ichar+Lfmt+1])
140            ichar = ichar + Lfmt
141        elif ft[0:1] == 'F' or ft[0:1] == 'f':
142            if ft.find('.') != -1:
143                Lft = len(ft.split('.')[0])
144            Lfmt = int(ft[1:Lft])
145            fvalues.append(np.float(rmNOnum(lin[ichar:ichar+Lfmt+1])))
146            ichar = ichar + Lfmt
147        elif ft[0:1] == 'D' or ft[0:1] == 'd':
148            if ft.find('.') != -1:
149                Lft = len(ft.split('.')[0])
150            Lfmt = int(ft[1:Lft])
151            fvalues.append(np.float64(rmNOnum(lin[ichar:ichar+Lfmt+1])))
152            ichar = ichar + Lfmt
153        elif ft[0:1] == 'I' or ft[0:1] == 'i':
154            Lfmt = int(ft[1:Lft+1])
155            fvalues.append(int(rmNOnum(lin[ichar:ichar+Lfmt+1])))
156            ichar = ichar + Lfmt
157        elif ft.find('X') != -1 or ft.find('x') != -1:
158            print '    ' + fname + ': skipping space'
159            ichar = ichar + int(rmNOnum(ft))
160        else:
161            print errormsg
162            print '  ' + fname + ": format '" + ft[0:1] + "' not ready!!"
163            print '    Available formats:',afmts
164            quit(-1)
165
166        ival = ival + 1
167
168    return fvalues
169
170
171def ts_header(ln):
172    """ Function to get the values of the header of the *.TS files
173      line=ASCII lines with the header of the TS file
174      getting the line format from WRFV3.3 'EMCORE' in file 'share/wrf_timeseries.F'
175    """
176    fname = 'ts_header'
177
178    fmts=['A26', 'I2', 'I3', 'A6', 'A2', 'F7.3', 'A1', 'F8.3', 'A3', 'I4', 'A1', 'I4',\
179      'A3', 'F7.3', 'A1', 'F8.3', 'A2', 'F6.1', 'A7']
180
181    headervalues = values_fortran_fmt(ln,fmts)
182
183    return headervalues
184
185def variables_values(varName):
186    """ Function to provide values to plot the different variables values from ASCII file
187      'variables_values.dat'
188    variables_values(varName)
189      [varName]= name of the variable
190        return: [var name], [std name], [minimum], [maximum],
191          [long name]('|' for spaces), [units], [color palette] (following:
192          http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
193     [varn]: original name of the variable
194       NOTE: It might be better doing it with an external ASII file. But then we
195         got an extra dependency...
196    >>> variables_values('WRFght')
197    ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow']
198    """
199    import subprocess as sub
200
201    fname='variables_values'
202
203    if varName == 'h':
204        print fname + '_____________________________________________________________'
205        print variables_values.__doc__
206        quit()
207
208# This does not work....
209#    folderins = sub.Popen(["pwd"], stdout=sub.PIPE)
210#    folder = list(folderins.communicate())[0].replace('\n','')
211# From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python
212    folder = os.path.dirname(os.path.realpath(__file__))
213
214    infile = folder + '/variables_values.dat'
215
216    if not os.path.isfile(infile):
217        print errormsg
218        print '  ' + fname + ": File '" + infile + "' does not exist !!"
219        quit(-1)
220
221# Variable name might come with a statistical surname...
222    stats=['min','max','mean','stdv', 'sum']
223
224# Variables with a statistical section on their name...
225    NOstatsvars = ['zmaxth', 'zmax_th', 'lmax_th', 'lmaxth']
226
227    ifst = False
228    if not searchInlist(NOstatsvars, varName.lower()):
229        for st in stats:
230            if varName.find(st) > -1:
231                print '    '+ fname + ": varibale '" + varName + "' with a " +       \
232                  "statistical surname: '",st,"' !!"
233                Lst = len(st)
234                LvarName = len(varName)
235                varn = varName[0:LvarName - Lst]
236                ifst = True
237                break
238    if not ifst:
239        varn = varName
240
241    ncf = open(infile, 'r')
242
243    for line in ncf:
244        if line[0:1] != '#':
245            values = line.replace('\n','').split(',')
246            if len(values) != 8:
247                print errormsg
248                print "problem in varibale:'", values[0],                            \
249                  'it should have 8 values and it has',len(values)
250                quit(-1)
251
252            if varn[0:6] == 'varDIM': 
253# Variable from a dimension (all with 'varDIM' prefix)
254                Lvarn = len(varn)
255                varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1.,                 \
256                  "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1',  \
257                   'rainbow']
258            else:
259                varvals = [values[1].replace(' ',''), values[2].replace(' ',''),     \
260                  np.float(values[3]), np.float(values[4]),values[5].replace(' ',''),\
261                  values[6].replace(' ',''), values[7].replace(' ','')]
262            if values[0] == varn:
263                ncf.close()
264                return varvals
265                break
266
267    print errormsg
268    print '  ' + fname + ": variable '" + varn + "' not defined !!!"
269    ncf.close()
270    quit(-1)
271
272    return 
273
274####### ###### ##### #### ### ## #
275
276parser = OptionParser()
277parser.add_option("-f", "--TS_file", dest="lfile",
278  help="Time Series ASCII text file to use", metavar="FILE")
279parser.add_option("-s", "--SimulationStartTime", dest="stime",
280  help="Starting time of the simulation ([YYYY][MM][DD][HH][MI][SS] format)", metavar="DATE")
281
282(opts, args) = parser.parse_args()
283
284
285tsvn = ['t', 'q', 'u', 'v', 'psfc', 'glw', 'gsw', 'hfx', 'lh', 'tsk', 'tslb1', 'rainc', 'rainnc', 'clw']
286
287tsvln = ['2 m Temperature', '2 m vapor mixing ratio', '10 m U wind (earth-relative)', '10 m V wind (earth-relative)', 'surface pressure', 'downward longwave radiation flux at the ground (downward is positive)', 'net shortwave radiation flux at the ground (downward is positive)', 'surface sensible heat flux (upward is positive)', 'surface latent heat flux (upward is positive)', 'skin temperature', 'top soil layer temperature', 'rainfall from a cumulus scheme', 'rainfall from an explicit scheme', 'total column-integrated water vapor and cloud variables']
288
289tsvu = ['K', 'kg/kg', 'm/s', 'm/s', 'Pa', 'W/m2', 'W/m2', 'W/m2', 'W/m2', 'K', 'K', 'mm', 'mm', '1']
290
291
292#######    #######
293## MAIN
294    #######
295
296ofile = 'ts.nc'
297Ntsvariables = len(tsvn)
298
299if not os.path.isfile(opts.lfile):
300    print errormsg
301    print '  ' + main + ': Time-Series ASCII text file "' + opts.lfile +             \
302      '" does not exist !!'
303    print errormsg
304    quit()
305
306if opts.stime is None:
307    print errormsg
308    print '  ' + main + ': No initial date/time of the simulation is provided!'
309    quit(-1)
310else:
311    stime = opts.stime
312    refdate = stime[0:4] + '-' + stime[4:6] + '-' + stime[6:8] + ' ' + stime[8:10] + \
313      ':' + stime[10:12] + ':' + stime[12:14]
314
315objlfile = open(opts.lfile, 'r')
316
317objofile = NetCDFFile(ofile, 'w')
318
319# Creation of dimensions
320##
321objofile.createDimension('time',None)
322
323set_attribute(objofile, 'author', 'Lluis Fita Borrell')
324set_attribute(objofile, 'institution', 'Laboratoire Meteorologique Dynamique')
325set_attribute(objofile, 'university', 'University Pierre et Marie Curie')
326set_attribute(objofile, 'center', 'Centre national de la recherche scientifique')
327set_attribute(objofile, 'country', 'France')
328set_attribute(objofile, 'city', 'Paris')
329set_attribute(objofile, 'script', 'TS_ASCII_netCFD.py')
330set_attribute(objofile, 'version', '1.0')
331
332time_step = []
333psfc = []
334rainc = []
335rainnc = []
336drydens = []
337
338tsvals = {}
339
340iline=0
341itz = 0
342for line in objlfile:
343    values = reduce_spaces(line)
344#    print iline, values[0], dimz, Searchdimz
345# Writting general information
346    if iline == 0:
347        newvar = objofile.createVariable('station','c')
348        valueshead = ts_header(line)
349
350        set_attribute(newvar, 'name', reduce_last_spaces(valueshead[0]))
351        set_attribute(newvar, 'acronym',valueshead[3].replace(' ',''))
352
353        set_attribute(newvar, 'real_lon', valueshead[5])
354        set_attribute(newvar, 'real_lat', valueshead[7])
355
356        set_attribute(newvar, 'x_grid_point', valueshead[9])
357        set_attribute(newvar, 'y_grid_point', valueshead[11])
358
359        set_attribute(newvar, 'model_lon', valueshead[13])
360        set_attribute(newvar, 'model_lat', valueshead[15])
361        set_attribute(newvar, 'model_height', valueshead[17])
362        simstarttime = refdate
363    else:
364        tsvals[itz] = values
365        time_step.append(np.float(values[1]))
366        itz = itz + 1
367    iline = iline + 1
368
369dimt = len(time_step)
370
371print '  Found:',dimt,'time steps'
372objlfile.close()
373
374time_stepv = np.zeros((dimt), dtype=np.float)
375tsvaluesv = np.zeros( (dimt,Ntsvariables), dtype= np.float)
376
377pracc = np.zeros((dimt), dtype=np.float)
378
379itz = 0
380for it in range(dimt):
381    time_stepv[it] = np.float(time_step[it])
382
383    for iv in range(Ntsvariables):
384        tsvaluesv[it,iv] = np.float(tsvals[itz][iv+5])
385
386    pracc[it] = np.float(tsvals[it][16]) + np.float(tsvals[it][17])
387    itz = itz + 1
388
389# time
390newvar = objofile.createVariable('time','f8',('time'))
391newvar[:] = time_stepv*3600.
392newattr = basicvardef(newvar, 'time', 'time', 'seconds since ' +               \
393  simstarttime.replace('_',' '))
394set_attribute(newvar, 'calendar', 'standard')
395
396dt = time_stepv[1] - time_stepv[0]
397
398# time-series variables
399for iv in range(Ntsvariables):
400    if tsvn[iv] == 't' or tsvn[iv] == 'u' or tsvn[iv] == 'v' or tsvn[iv] == 'q': 
401      varname, stdname, minvar, maxvar, longname, unitsvar, cbarvar =                \
402      variables_values('TS' + tsvn[iv])
403      tsu = unitsvar
404    else:
405      varname, stdname, minvar, maxvar, longname, unitsvar, cbarvar =                \
406      variables_values(tsvn[iv])
407      tsu = tsvu[iv]
408
409    newvar = objofile.createVariable(varname, 'f4', ('time'))
410    newvar[:] = tsvaluesv[:,iv]
411
412    newattr = basicvardef(newvar, stdname, longname.replace('|',' '), tsu)
413    newattr = set_attribute(newvar, 'wrfTSname', tsvn[iv])
414    newattr = set_attribute(newvar, 'wrfTSdesc', tsvln[iv])
415
416# Extra vars
417
418# pr
419varvals = np.zeros((dimt), dtype=np.float)
420varvals[1:dimt] = pracc[1:dimt] - pracc[0:dimt-1]
421varname, stdname, minvar, maxvar, longname, unitsvar, cbarvar =                      \
422      variables_values('RAINTOT')
423
424newvar = objofile.createVariable(varname, 'f4', ('time'))
425newvar[:] = varvals / dt
426newattr = basicvardef(newvar, stdname, longname.replace('|',' '), unitsvar )
427
428objofile.sync()
429objofile.close()
430
431print 'Successfull generation of Time-Series netCDF file "' + ofile + '" !!!!!'
Note: See TracBrowser for help on using the repository browser.