source: lmdz_wrf/branches/LMDZ_WRFmeas/tools/create_OBSnetcdf.py @ 1639

Last change on this file since 1639 was 414, checked in by lfita, 10 years ago

Removing WRFV3 folder and keeping only the essential
File size: 32.8 KB
Line 
1# -*- coding: iso-8859-15 -*-
2# Generic program to transfrom ASCII observational data in columns to a netcdf
3# L. Fita, LMD February 2015
4## e.g. # create_OBSnetcdf.py -c '#' -d ACAR/description.dat -e space -f ACAR/2012/10/ACAR_121018.asc -g true -t 19491201000000,seconds -k trajectory
5
6import numpy as np
7from netCDF4 import Dataset as NetCDFFile
8import os
9import re
10from optparse import OptionParser
11
12# version
13version=1.1
14
15# Filling values for floats, integer and string
16fillValueF = 1.e20
17fillValueI = -99999
18fillValueS = '---'
19
20# Length of the string variables
21StringLength = 200
22
23# Size of the map for the complementary variables/maps
24Ndim2D = 100
25
26main = 'create_OBSnetcdf.py'
27errormsg = 'ERROR -- error -- ERROR -- error'
28warnmsg = 'WARNING -- warning -- WARNING -- warning'
29
30fillValue = 1.e20
31
32def searchInlist(listname, nameFind):
33    """ Function to search a value within a list
34    listname = list
35    nameFind = value to find
36    >>> searInlist(['1', '2', '3', '5'], '5')
37    True
38    """
39    for x in listname:
40      if x == nameFind:
41        return True
42    return False
43
44def set_attribute(ncvar, attrname, attrvalue):
45    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
46    ncvar = object netcdf variable
47    attrname = name of the attribute
48    attrvalue = value of the attribute
49    """
50    import numpy as np
51    from netCDF4 import Dataset as NetCDFFile
52
53    attvar = ncvar.ncattrs()
54    if searchInlist(attvar, attrname):
55        attr = ncvar.delncattr(attrname)
56
57    attr = ncvar.setncattr(attrname, attrvalue)
58
59    return ncvar
60
61def basicvardef(varobj, vstname, vlname, vunits):
62    """ Function to give the basic attributes to a variable
63    varobj= netCDF variable object
64    vstname= standard name of the variable
65    vlname= long name of the variable
66    vunits= units of the variable
67    """
68    attr = varobj.setncattr('standard_name', vstname)
69    attr = varobj.setncattr('long_name', vlname)
70    attr = varobj.setncattr('units', vunits)
71
72def remove_NONascii(string):
73    """ Function to remove that characters which are not in the standard 127 ASCII
74      string= string to transform
75    >>> remove_NONascii('Lluís')
76    Lluis
77    """
78    fname = 'remove_NONascii'
79
80    newstring = string
81
82    RTFchar= ['á', 'é', 'í', 'ó', 'ú', 'à', 'Ú', 'ì', 'ò', 'ù', 'â', 'ê', 'î', 'ÃŽ',  \
83      'û', 'À', 'ë', 'ï', 'ö', 'ÃŒ', 'ç', 'ñ','Ê', 'œ', 'Á', 'É', 'Í', 'Ó', 'Ú', 'À', \
84      'È', 'Ì', 'Ò', 'Ù', 'Â', 'Ê', 'Î', 'Ô', 'Û', 'Ä', 'Ë', 'Ï', 'Ö', 'Ü', 'Ç', 'Ñ',\
85      'Æ', 'Œ', '\n', '\t']
86    ASCchar= ['a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o',  \
87      'u', 'a', 'e', 'i', 'o', 'u', 'c', 'n','ae', 'oe', 'A', 'E', 'I', 'O', 'U', 'A', \
88      'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'C', 'N',\
89      'AE', 'OE', '', ' ']
90
91    Nchars = len(RTFchar)
92    for ichar in range(Nchars):
93        foundchar = string.find(RTFchar[ichar])
94        if foundchar != 0:
95            newstring = newstring.replace(RTFchar[ichar], ASCchar[ichar])
96
97    return newstring
98
99def read_description(fdobs, dbg):
100    """ reads the description file of the observational data-set
101    fdobs= descriptive observational data-set
102    dbg= boolean argument for debugging
103    * Each station should have a 'description.dat' file with:
104      institution=Institution who creates the data
105      department=Department within the institution
106      scientists=names of the data producers
107      contact=contact of the data producers
108      description=description of the observations
109      acknowledgement=sentence of acknowlegement
110      comment=comment for the measurements
111
112      MissingValue='|' list of ASCII values for missing values within the data
113        (as they appear!)
114      comment=comments
115
116      varN='|' list of variable names
117      varLN='|' list of long variable names
118      varU='|' list units of the variables
119      varBUFR='|' list BUFR code of the variables
120      varTYPE='|' list of variable types ('D', 'F', 'I', 'I64', 'S')
121
122      NAMElon=name of the variable with the longitude (x position)
123      NAMElat=name of the variable with the latitude (y position)
124      NAMEheight=ind_alt
125      NAMEtime=name of the varibale with the time
126      FMTtime=format of the time (as in 'C', 'CFtime' for already CF-like time)
127
128    """
129    fname = 'read_description'
130
131    descobj = open(fdobs, 'r')
132    desc = {}
133
134    namevalues = []
135
136    for line in descobj:
137        if line[0:1] != '#' and len(line) > 1 :
138            descn = remove_NONascii(line.split('=')[0])
139            descv = remove_NONascii(line.split('=')[1])
140            namevalues.append(descn)
141            if descn[0:3] != 'var':
142                if descn != 'MissingValue':
143                    desc[descn] = descv
144                else:
145                    desc[descn] = []
146                    for dn in descv.split('|'): desc[descn].append(dn)
147                    print '  ' + fname + ': missing values found:',desc[descn]
148            elif descn[0:3] == 'var':
149                desc[descn] = descv.split('|')
150            elif descn[0:4] == 'NAME':
151                desc[descn] = descv
152            elif descn[0:3] == 'FMT':
153                desc[descn] = descv
154
155    if dbg:
156        Nvars = len(desc['varN'])
157        print '  ' + fname + ": description content of '" + fdobs + "'______________"
158        for varn in namevalues:
159            if varn[0:3] != 'var':
160                print '    ' + varn + ':',desc[varn]
161            elif varn == 'varN':
162                print '    * Variables:'
163                for ivar in range(Nvars):
164                    varname = desc['varN'][ivar]
165                    varLname = desc['varLN'][ivar]
166                    varunits = desc['varU'][ivar]
167                    if desc.has_key('varBUFR'): 
168                        varbufr=[ivar]
169                        print '      ', ivar, varname + ':',varLname,'[',varunits,   \
170                          ']','bufr code:',varbufr
171                    else:
172                        print '      ', ivar, varname + ':',varLname,'[',varunits,']'
173
174    descobj.close()
175
176    return desc
177
178def value_fmt(val, miss, fmt):
179    """ Function to transform an ASCII value to a given format
180      val= value to transform
181      miss= list of possible missing values
182      fmt= format to take:
183        'D': float double precission
184        'F': float
185        'I': integer
186        'I64': 64-bits integer
187        'S': string
188    >>> value_fmt('9876.12', '-999', 'F')
189    9876.12
190    """
191
192    fname = 'value_fmt'
193
194    fmts = ['D', 'F', 'I', 'I64', 'S']
195    Nfmts = len(fmts)
196
197    if not searchInlist(miss,val):
198        if not searchInlist(fmts, fmt):
199            print errormsg
200            print '  ' + fname + ": format '" + fmt + "' not ready !!"
201            quit(-1)
202        else:
203            if fmt == 'D':
204                newval = np.float32(val)
205            elif fmt == 'F':
206                newval = np.float(val)
207            elif fmt == 'I':
208                newval = int(val)
209            elif fmt == 'I64':
210                newval = np.int64(val)
211            elif fmt == 'S':
212                newval = val
213    else:
214        newval = None
215
216    return newval
217
218def read_datavalues(dataf, comchar, colchar, fmt, miss, varns, dbg):
219    """ Function to read from an ASCII file values in column
220      dataf= data file
221      comchar= list of the characters indicating comment in the file
222      colchar= character which indicate end of value in the column
223      dbg= debug mode or not
224      fmt= list of kind of values to be found
225      miss= missing value
226      varns= list of name of the variables to find
227    """
228
229    fname = 'read_datavalues'
230
231    ofile = open(dataf, 'r')
232    Nvals = len(fmt)
233
234    finalvalues = {}
235
236    iline = 0
237    for line in ofile:
238        line = line.replace('\n','').replace(chr(13),'')
239        if not searchInlist(comchar,line[0:1]) and len(line) > 1:
240            values0 = line.split(colchar)
241# Removing no-value columns
242            values = []
243            for iv in values0:
244                if len(iv) > 0: values.append(iv)
245
246            Nvalues = len(values)
247# Checkings for wierd characters at the end of lines (use it to check)
248#            if values[Nvalues-1][0:4] == '-999':
249#                print line,'last v:',values[Nvalues-1],'len:',len(values[Nvalues-1])
250#                for ic in range(len(values[Nvalues-1])):
251#                    print ic,ord(values[Nvalues-1][ic:ic+1])
252#                quit()
253
254            if len(values[Nvalues-1]) == 0:
255                Nvalues = Nvalues - 1
256               
257            if Nvalues != Nvals:
258                print warnmsg
259                print '  ' + fname + ': number of formats:',Nvals,' and number of ', \
260                  'values:',Nvalues,' with split character *' + colchar +            \
261                  '* does not coincide!!'
262                print '    * what is found is ________'
263                if Nvalues > Nvals:
264                    Nshow = Nvals
265                    for ivar in range(Nshow):
266                        print '    ',varns[ivar],'fmt:',fmt[ivar],'value:',values[ivar]
267                    print '      missing formats for:',values[Nshow:Nvalues+1]
268                    print '      values not considered, continue'
269                else:
270                    Nshow = Nvalues
271                    for ivar in range(Nshow):
272                        print '   ',varns[ivar],'fmt:',fmt[ivar],'value:',values[ivar]
273                    print '      excess of formats:',fmt[Nshow:Nvals+1]
274                    quit(-1)
275
276# Reading and transforming values
277            if dbg: print '  ' + fname + ': values found _________'
278
279            for ivar in range(Nvals):
280                if dbg: 
281                    print iline, varns[ivar],'value:',values[ivar],miss,fmt[ivar]
282
283                if iline == 0:
284                    listvals = []
285                    listvals.append(value_fmt(values[ivar], miss, fmt[ivar]))
286                    finalvalues[varns[ivar]] = listvals
287                else:
288                    listvals = finalvalues[varns[ivar]]
289                    listvals.append(value_fmt(values[ivar], miss, fmt[ivar]))
290                    finalvalues[varns[ivar]] = listvals
291        else:
292# First line without values
293            if iline == 0: iline = -1
294   
295        iline = iline + 1
296
297    ofile.close()
298
299    return finalvalues
300
301def writing_str_nc(varo, values, Lchar):
302    """ Function to write string values in a netCDF variable as a chain of 1char values
303    varo= netCDF variable object
304    values = list of values to introduce
305    Lchar = length of the string in the netCDF file
306    """
307
308    Nvals = len(values)
309    for iv in range(Nvals):   
310        stringv=values[iv] 
311        charvals = np.chararray(Lchar)
312        Lstr = len(stringv)
313        charvals[Lstr:Lchar] = ''
314
315        for ich in range(Lstr):
316            charvals[ich] = stringv[ich:ich+1]
317
318        varo[iv,:] = charvals
319
320    return
321
322def Stringtimes_CF(tvals, fmt, Srefdate, tunits, dbg):
323    """ Function to transform a given data in String formats to a CF date
324      tvals= string temporal values
325      fmt= format of the the time values
326      Srefdate= reference date in [YYYY][MM][DD][HH][MI][SS] format
327      tunits= units to use ('weeks', 'days', 'hours', 'minutes', 'seconds')
328      dbg= debug
329    >>> Stringtimes_CF(['19760217082712','20150213101837'], '%Y%m%d%H%M%S',
330      '19491201000000', 'hours', False)
331    [229784.45333333  571570.31027778]
332    """
333    import datetime as dt
334
335    fname = 'Stringtimes'
336
337    dimt = len(tvals)
338
339    yrref = int(Srefdate[0:4])
340    monref = int(Srefdate[4:6])
341    dayref = int(Srefdate[6:8])
342    horref = int(Srefdate[8:10])
343    minref = int(Srefdate[10:12])
344    secref = int(Srefdate[12:14])
345    refdate = dt.datetime( yrref, monref, dayref, horref, minref, secref)
346
347    cftimes = np.zeros((dimt), dtype=np.float)
348
349    Nfmt=len(fmt.split('%'))
350
351    if dbg: print '  ' + fname + ': fmt=',fmt,'refdate:',Srefdate,'uits:',tunits,    \
352      'date dt_days dt_time deltaseconds CFtime _______'
353    for it in range(dimt):
354
355# Removing excess of mili-seconds (up to 6 decimals)
356        if fmt.split('%')[Nfmt-1] == 'f':
357            tpoints = tvals[it].split('.')
358            if len(tpoints[len(tpoints)-1]) > 6:
359                milisec = '{0:.6f}'.format(np.float('0.'+tpoints[len(tpoints)-1]))[0:7]
360                newtval = ''
361                for ipt in range(len(tpoints)-1):
362                    newtval = newtval + tpoints[ipt] + '.'
363                newtval = newtval + str(milisec)[2:len(milisec)+1]
364            else:
365                newtval = tvals[it]
366            tval = dt.datetime.strptime(newtval, fmt)
367        else:
368            tval = dt.datetime.strptime(tvals[it], fmt)
369
370        deltatime = tval - refdate
371        deltaseconds = deltatime.days*24.*3600. + deltatime.seconds +                \
372          deltatime.microseconds/100000.
373        if tunits == 'weeks':
374            deltat = 7.*24.*3600.
375        elif tunits == 'days':
376            deltat = 24.*3600.
377        elif tunits == 'hours':
378            deltat = 3600.
379        elif tunits == 'minutes':
380            deltat = 60.
381        elif tunits == 'seconds':
382            deltat = 1.
383        else:
384            print errormsg
385            print '  ' + fname + ": time units '" + tunits + "' not ready !!"
386            quit(-1)
387
388        cftimes[it] = deltaseconds / deltat
389        if dbg:
390            print '  ' + tvals[it], deltatime, deltaseconds, cftimes[it]
391
392    return cftimes
393
394def adding_complementary(onc, dscn, okind, dvalues, tvals, refCFt, CFtu, Nd, dbg):
395    """ Function to add complementary variables as function of the observational type
396      onc= netcdf objecto file to add the variables
397      dscn= description dictionary
398      okind= observational kind
399      dvalues= values
400      tvals= CF time values
401      refCFt= reference time of CF time (in [YYYY][MM][DD][HH][MI][SS])
402      CFtu= CF time units
403      Nd= size of the domain
404      dbg= debugging flag
405    """
406    import numpy.ma as ma
407
408    fname = 'adding_complementary'
409 
410# Kind of observations which require de integer lon/lat (for the 2D map)
411    map2D=['multi-points', 'trajectory']
412
413    SrefCFt = refCFt[0:4] +'-'+ refCFt[4:6] +'-'+ refCFt[6:8] + ' ' + refCFt[8:10] + \
414      ':'+ refCFt[10:12] +':'+ refCFt[12:14]
415
416    if dscn['NAMElon'] == '-' or dscn['NAMElat'] == '-':
417        print errormsg
418        print '  ' + fname + ": to complement a '" + okind + "' observation kind " + \
419          " a given longitude ('NAMElon':",dscn['NAMElon'],") and latitude ('" +     \
420          "'NAMElat:'", dscn['NAMElat'],') from the data has to be provided and ' +  \
421          'any are given !!'
422        quit(-1)
423
424    if not dvalues.has_key(dscn['NAMElon']) or not dvalues.has_key(dscn['NAMElat']):
425        print errormsg
426        print '  ' + fname + ": observations do not have 'NAMElon':",                \
427          dscn['NAMElon'],"and/or 'NAMElat:'", dscn['NAMElat'],' !!'
428        print '    available data:',dvalues.keys()
429        quit(-1)
430
431    if okind == 'trajectory':
432        if dscn['NAMEheight'] == '-':
433            print warnmsg
434            print '  ' + fname + ": to complement a '" + okind + "' observation " +  \
435              "kind a given height ('NAMEheight':",dscn['NAMEheight'],"') might " +  \
436              'be provided and any is given !!'
437            quit(-1)
438
439        if not dvalues.has_key(dscn['NAMEheight']):
440            print errormsg
441            print '  ' + fname + ": observations do not have 'NAMEtime':",           \
442              dscn['NAMEtime'],' !!'
443            print '    available data:',dvalues.keys()
444            quit(-1)
445
446    if searchInlist(map2D, okind):
447# A new 2D map with the number of observation will be added for that 'NAMElon'
448#   and 'NAMElat' are necessary. A NdxNd domain space size will be used.
449        objfile.createDimension('lon2D',Nd)
450        objfile.createDimension('lat2D',Nd)
451        lonvals = ma.masked_equal(dvalues[dscn['NAMElon']], None)
452        latvals = ma.masked_equal(dvalues[dscn['NAMElat']], None)
453
454        minlon = min(lonvals)
455        maxlon = max(lonvals)
456        minlat = min(latvals)
457        maxlat = max(latvals)
458
459        blon = (maxlon - minlon)/(Nd-1)
460        blat = (maxlat - minlat)/(Nd-1)
461
462        newvar = onc.createVariable( 'lon2D', 'f8', ('lon2D'))
463        basicvardef(newvar, 'longitude', 'longitude map observations','degrees_East')
464        newvar[:] = minlon + np.arange(Nd)*blon
465        newattr = set_attribute(newvar, 'axis', 'X')
466
467        newvar = onc.createVariable( 'lat2D', 'f8', ('lat2D'))
468        basicvardef(newvar, 'latitude', 'latitude map observations', 'degrees_North')
469        newvar[:] = minlat + np.arange(Nd)*blat
470        newattr = set_attribute(newvar, 'axis', 'Y')
471
472        if dbg: 
473            print '  ' + fname + ': minlon=',minlon,'maxlon=',maxlon
474            print '  ' + fname + ': minlat=',minlat,'maxlat=',maxlat
475            print '  ' + fname + ': precission on x-axis=', blon*(Nd-1), 'y-axis=',  \
476              blat*(Nd-1)
477       
478    if okind == 'multi-points':
479        map2D = np.ones((Nd, Nd), dtype=np.float)*fillValueI
480
481        dimt = len(tvals)
482        Nlost = 0
483        for it in range(dimt):
484            lon = dvalues[dscn['NAMElon']][it]
485            lat = dvalues[dscn['NAMElat']][it]
486            if lon is not None and lat is not None:
487                ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
488                ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
489
490                if map2D[ilat,ilon] == fillValueI: 
491                    map2D[ilat,ilon] = 1
492                else:
493                    map2D[ilat,ilon] = map2D[ilat,ilon] + 1
494                if dbg: print it, lon, lat, ilon, ilat, map2D[ilat,ilon]
495
496        newvar = onc.createVariable( 'mapobs', 'f4', ('lat2D', 'lon2D'),             \
497          fill_value = fillValueI)
498        basicvardef(newvar, 'map_observations', 'number of observations', '-')
499        newvar[:] = map2D
500        newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D')
501
502    elif okind == 'trajectory':
503# A new 2D map with the trajectory 'NAMElon' and 'NAMElat' and maybe 'NAMEheight'
504#   are necessary. A NdxNdxNd domain space size will be used. Using time as
505#   reference variable
506        if dscn['NAMEheight'] == '-':
507# No height
508            map2D = np.ones((Nd, Nd), dtype=np.float)*fillValueI
509
510            dimt = len(tvals)
511            Nlost = 0
512            if dbg: print ' time-step lon lat ix iy passes _______'
513            for it in range(dimt):
514                lon = dvalues[dscn['NAMElon']][it]
515                lat = dvalues[dscn['NAMElat']][it]
516                if lon is  not None and lat is not None:
517                    ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
518                    ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
519
520                    if map2D[ilat,ilon] == fillValueI: 
521                        map2D[ilat,ilon] = 1
522                    else:
523                        map2D[ilat,ilon] = map2D[ilat,ilon] + 1
524                    if dbg: print it, lon, lat, ilon, ilat, map2D[ilat,ilon]
525
526            newvar = onc.createVariable( 'trjobs', 'i', ('lat2D', 'lon2D'),          \
527              fill_value = fillValueI)
528            basicvardef(newvar, 'trajectory_observations', 'number of passes', '-' )
529            newvar[:] = map2D
530            newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D')
531
532        else:
533            ivn = 0
534            for vn in dscn['varN']:
535                if vn == dscn['NAMEheight']: 
536                    zu = dscn['varU'][ivn]
537                    break
538                ivn = ivn + 1
539                   
540            objfile.createDimension('z2D',Nd)
541            zvals = ma.masked_equal(dvalues[dscn['NAMEheight']], None)
542            minz = min(zvals)
543            maxz = max(zvals)
544
545            bz = (maxz - minz)/(Nd-1)
546
547            newvar = onc.createVariable( 'z2D', 'f8', ('z2D'))
548            basicvardef(newvar, 'z2D', 'z-coordinate map observations', zu)
549            newvar[:] = minz + np.arange(Nd)*bz
550            newattr = set_attribute(newvar, 'axis', 'Z')
551
552            if dbg:
553                print '  ' + fname + ': zmin=',minz,zu,'zmax=',maxz,zu
554                print '  ' + fname + ': precission on z-axis=', bz*(Nd-1), zu
555
556            map3D = np.ones((Nd, Nd, Nd), dtype=int)*fillValueI
557            dimt = len(tvals)
558            Nlost = 0
559            if dbg: print ' time-step lon lat z ix iy iz passes _______'
560            for it in range(dimt):
561                lon = dvalues[dscn['NAMElon']][it]
562                lat = dvalues[dscn['NAMElat']][it]
563                z = dvalues[dscn['NAMEheight']][it]
564                if lon is  not None and lat is not None and z is not None:
565                    ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
566                    ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
567                    iz = int((Nd-1)*(z - minz)/(maxz - minz))
568
569                    if map3D[iz,ilat,ilon] == fillValueI: 
570                        map3D[iz,ilat,ilon] = 1
571                    else:
572                        map3D[iz,ilat,ilon] = map3D[iz,ilat,ilon] + 1
573                    if dbg: print it, lon, lat, z, ilon, ilat, iz, map3D[iz,ilat,ilon]
574
575            newvar = onc.createVariable( 'trjobs', 'i', ('z2D', 'lat2D', 'lon2D'),   \
576              fill_value = fillValueI)
577            basicvardef(newvar, 'trajectory_observations', 'number of passes', '-')
578            newvar[:] = map3D
579            newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D z2D')
580
581    onc.sync()
582    return
583
584def adding_station_desc(onc,stdesc):
585    """ Function to add a station description in a netCDF file
586      onc= netCDF object
587      stdesc= station description lon, lat, height
588    """
589    fname = 'adding_station_desc'
590
591    newdim = onc.createDimension('nst',1)
592
593    if onc.variables.has_key('lon'):
594        print warnmsg
595        print '  ' + fname + ": variable 'lon' already exist !!"
596        print "    renaming it as 'lonst'"
597        lonname = 'lonst'
598    else:
599        lonname = 'lon'
600
601    newvar = objfile.createVariable( lonname, 'f4', ('nst'))
602    basicvardef(newvar, lonname, 'longitude', 'degrees_West' )
603    newvar[:] = stdesc[0]
604
605    if onc.variables.has_key('lat'):
606        print warnmsg
607        print '  ' + fname + ": variable 'lat' already exist !!"
608        print "    renaming it as 'latst'"
609        latname = 'latst'
610    else:
611        latname = 'lat'
612
613    newvar = objfile.createVariable( latname, 'f4', ('nst'))
614    basicvardef(newvar, lonname, 'latitude', 'degrees_North' )
615    newvar[:] = stdesc[1]
616
617    if onc.variables.has_key('height'):
618        print warnmsg
619        print '  ' + fname + ": variable 'height' already exist !!"
620        print "    renaming it as 'heightst'"
621        heightname = 'heightst'
622    else:
623        heightname = 'height'
624
625    newvar = objfile.createVariable( heightname, 'f4', ('nst'))
626    basicvardef(newvar, heightname, 'height above sea level', 'm' )
627    newvar[:] = stdesc[2]
628
629    return
630
631####### ###### ##### #### ### ## #
632
633strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
634  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
635
636kindobs=['multi-points', 'single-station', 'trajectory']
637strkObs="kind of observations; 'multi-points': multiple individual punctual obs " +  \
638  "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
639  "'trajectory': following a trajectory"
640
641parser = OptionParser()
642parser.add_option("-c", "--comments", dest="charcom",
643  help="':', list of characters used for comments", metavar="VALUES")
644parser.add_option("-d", "--descriptionfile", dest="fdesc",
645  help="description file to use", metavar="FILE")
646parser.add_option("-e", "--end_column", dest="endcol",
647  help="character to indicate end of the column ('space', for ' ')", metavar="VALUE")
648parser.add_option("-f", "--file", dest="obsfile",
649  help="observational file to use", metavar="FILE")
650parser.add_option("-g", "--debug", dest="debug",
651  help="whther debug is required ('false', 'true')", metavar="VALUE")
652parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
653  help=strkObs, metavar="FILE")
654parser.add_option("-s", "--stationLocation", dest="stloc", 
655  help="longitude, latitude and height of the station (only for 'single-station')", 
656  metavar="FILE")
657parser.add_option("-t", "--CFtime", dest="CFtime", help=strCFt, metavar="VALUE")
658
659(opts, args) = parser.parse_args()
660
661#######    #######
662## MAIN
663    #######
664
665ofile='OBSnetcdf.nc'
666
667if opts.charcom is None:
668    print warnmsg
669    print '  ' + main + ': No list of comment characters provided!!'
670    print '    assuming no need!'
671    charcomments = []
672else:
673    charcomments = opts.charcom.split(':')   
674
675if opts.fdesc is None:
676    print errormsg
677    print '  ' + main + ': No description file for the observtional data provided!!'
678    quit(-1)
679
680if opts.endcol is None:
681    print warnmsg
682    print '  ' + main + ': No list of comment characters provided!!'
683    print "    assuming 'space'"
684    endcol = ' '
685else:
686    if opts.endcol == 'space':
687        endcol = ' '
688    else:
689        endcol = opts.endcol
690
691if opts.obsfile is None:
692    print errormsg
693    print '  ' + main + ': No observtions file provided!!'
694    quit(-1)
695
696if opts.debug is None:
697    print warnmsg
698    print '  ' + main + ': No debug provided!!'
699    print "    assuming 'False'"
700    debug = False
701else:
702    if opts.debug == 'true':
703        debug = True
704    else:
705        debug = False
706
707if not os.path.isfile(opts.fdesc):
708    print errormsg
709    print '   ' + main + ": description file '" + opts.fdesc + "' does not exist !!"
710    quit(-1)
711
712if not os.path.isfile(opts.obsfile):
713    print errormsg
714    print '   ' + main + ": observational file '" + opts.obsfile + "' does not exist !!"
715    quit(-1)
716
717if opts.CFtime is None:
718    print warnmsg
719    print '  ' + main + ': No CFtime criteria are provided !!'
720    print "    either time is already in CF-format ('timeFMT=CFtime') in '" +        \
721      opts.fdesc + "'"
722    print "    or assuming refdate: '19491201000000' and time units: 'hours'"
723    referencedate = '19491201000000'
724    timeunits = 'hours'
725else:
726    referencedate = opts.CFtime.split(',')[0]
727    timeunits = opts.CFtime.split(',')[1]
728
729if opts.obskind is None:
730    print warnmsg
731    print '  ' + main + ': No kind of observations provided !!'
732    print "    assuming 'single-station': single station on a fixed position at 0,0,0"
733    obskind = 'single-station'
734    stationdesc = [0.0, 0.0, 0.0]
735else:
736    obskind = opts.obskind
737    if obskind == 'single-station':
738        if opts.stloc is None:
739            print errornmsg
740            print '  ' + main + ': No station location provided !!'
741            quit(-1)
742        else:
743            stationdesc = [np.float(opts.stloc.split(',')[0]),                       \
744              np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2])]
745    else:
746        obskind = opts.obskind
747
748# Reading description file
749##
750description = read_description(opts.fdesc, debug)
751
752Nvariables=len(description['varN'])
753formats = description['varTYPE']
754
755if len(formats) != Nvariables: 
756    print errormsg
757    print '  ' + main + ': number of formats:',len(formats),' and number of ' +     \
758      'variables', Nvariables,' does not coincide!!'
759    print '    * what is found is _______'
760    if Nvariables > len(formats):
761        Nshow = len(formats)
762        for ivar in range(Nshow):
763            print '      ',description['varN'][ivar],':', description['varLN'][ivar],\
764               '[', description['varU'][ivar], '] fmt:', formats[ivar]
765        print '      missing values for:', description['varN'][Nshow:Nvariables+1]
766    else:
767        Nshow = Nvariables
768        for ivar in range(Nshow):
769            print '      ',description['varN'][ivar],':', description['varLN'][ivar],\
770               '[', description['varU'][ivar], '] fmt:', formats[ivar]
771        print '      excess of formats for:', formats[Nshow:len(formats)+1]
772
773    quit(-1)
774
775# Reading values
776##
777datavalues = read_datavalues(opts.obsfile, charcomments, endcol, formats, 
778  description['MissingValue'], description['varN'], debug)
779
780# Total number of values
781Ntvalues = len(datavalues[description['varN'][0]])
782print main + ': total temporal values found:',Ntvalues
783
784objfile = NetCDFFile(ofile, 'w')
785 
786# Creation of dimensions
787##
788objfile.createDimension('time',None)
789objfile.createDimension('StrLength',StringLength)
790
791# Creation of variables
792##
793for ivar in range(Nvariables):
794    varn = description['varN'][ivar]
795    print "  including: '" + varn + "' ... .. ."
796
797    if formats[ivar] == 'D':
798        newvar = objfile.createVariable(varn, 'f32', ('time'), fill_value=fillValueF)
799        basicvardef(newvar, varn, description['varLN'][ivar],                        \
800          description['varU'][ivar])
801        newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn])
802    elif formats[ivar] == 'F':
803        newvar = objfile.createVariable(varn, 'f', ('time'), fill_value=fillValueF)
804        basicvardef(newvar, varn, description['varLN'][ivar],                        \
805          description['varU'][ivar])
806        newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn])
807    elif formats[ivar] == 'I':
808        newvar = objfile.createVariable(varn, 'i', ('time'), fill_value=fillValueI)
809        basicvardef(newvar, varn, description['varLN'][ivar],                        \
810          description['varU'][ivar])
811# Why is not wotking with integers?
812        vals = np.array(datavalues[varn])
813        for iv in range(Ntvalues):
814            if vals[iv] is None: vals[iv] = fillValueI
815        newvar[:] = vals
816    elif formats[ivar] == 'S':
817        newvar = objfile.createVariable(varn, 'c', ('time','StrLength'))
818        basicvardef(newvar, varn, description['varLN'][ivar],                        \
819          description['varU'][ivar])
820        writing_str_nc(newvar, datavalues[varn], StringLength)
821
822# Extra variable descriptions/attributes
823    if description.has_key('varBUFR'): 
824        set_attribute(newvar,'bufr_code',description['varBUFR'][ivar])
825
826    objfile.sync()
827
828# Time variable in CF format
829##
830print description['FMTtime'], 'CFtime'
831if description['FMTtime'] == 'CFtime':
832    timevals = datavalues[description['NAMEtime']]
833    iv = 0
834    for ivn in description['varN']:
835        if ivn == description['NAMEtime']:
836            tunits = description['varU'][iv]
837            break
838        iv = iv + 1
839else:
840    # Time as a composition of different columns
841    tcomposite = description['NAMEtime'].find('@')
842    if tcomposite != 0:
843        timevars = description['NAMEtime'].split('@')
844
845        print warnmsg
846        print '  ' + main + ': time values as combination of different columns!'
847        print '    combining:',timevars,' with a final format: ',description['FMTtime']
848
849        timeSvals = []
850        if debug: print '  ' + main + ': creating times _______'
851        for it in range(Ntvalues):
852            tSvals = ''
853            for tvar in timevars:
854                tSvals = tSvals + datavalues[tvar][it] + ' '
855
856            timeSvals.append(tSvals[0:len(tSvals)-1])
857            if debug: print it, '*' + timeSvals[it] + '*'
858
859        timevals = Stringtimes_CF(timeSvals, description['FMTtime'], referencedate,      \
860          timeunits, debug)
861    else:
862        timevals = Stringtimes_CF(datavalues[description['NAMEtime']],                   \
863          description['FMTtime'], referencedate, timeunits, debug)
864
865    CFtimeRef = referencedate[0:4] +'-'+ referencedate[4:6] +'-'+ referencedate[6:8] +   \
866      ' ' + referencedate[8:10] +':'+ referencedate[10:12] +':'+ referencedate[12:14]
867    tunits = timeunits + ' since ' + CFtimeRef
868
869if objfile.variables.has_key('time'):
870    print warnmsg
871    print '  ' + main + ": variable 'time' already exist !!"
872    print "    renaming it as 'CFtime'"
873    timeCFname = 'CFtime'
874    newdim = objfile.renameDimension('time','CFtime')
875    newvar = objfile.createVariable( timeCFname, 'f8', ('CFtime'))
876    basicvardef(newvar, timeCFname, 'time', tunits )
877else:
878    timeCFname = 'time'
879    newvar = objfile.createVariable( timeCFname, 'f8', ('time'))
880    basicvardef(newvar, timeCFname, 'time', tunits )
881
882set_attribute(newvar, 'calendar', 'standard')
883newvar[:] = timevals
884
885# Global attributes
886##
887for descn in description.keys():
888    if descn[0:3] != 'var' and descn[0:4] != 'NAME' and descn[0:3] != 'FMT':
889        set_attribute(objfile, descn, description[descn])
890
891set_attribute(objfile,'author_nc','Lluis Fita')
892set_attribute(objfile,'institution_nc','Laboratoire de Meteorology Dynamique, ' +    \
893  'LMD-Jussieu, UPMC, Paris')
894set_attribute(objfile,'country_nc','France')
895set_attribute(objfile,'script_nc','create_OBSnetcdf.py')
896set_attribute(objfile,'version_script',version)
897set_attribute(objfile,'information',                                                 \
898  'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
899
900objfile.sync()
901
902# Adding new variables as function of the observational type
903## 'multi-points', 'single-station', 'trajectory'
904
905if obskind != 'single-station':
906    adding_complementary(objfile, description, obskind, datavalues, timevals,        \
907      referencedate, timeunits, Ndim2D, debug)
908else:
909# Adding three variables with the station location, longitude, latitude and height
910    adding_station_desc(objfile,stationdesc)
911
912objfile.sync()
913objfile.close()
914
915print main + ": Successfull generation of netcdf observational file '" + ofile + "' !!"
Note: See TracBrowser for help on using the repository browser.