source: lmdz_wrf/trunk/tools/create_OBSnetcdf.py @ 574

Last change on this file since 574 was 553, checked in by lfita, 10 years ago

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