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

Last change on this file since 1030 was 928, checked in by lfita, 8 years ago

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