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

Last change on this file since 1784 was 648, checked in by lfita, 9 years ago

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