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

Last change on this file since 2816 was 2816, checked in by lfita, 5 years ago

Adding new debug print and better description of the code
File size: 66.5 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.2
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
44
45def typemod(value, typeval):
46    """ Function to give back a value in a given dtype
47    >>> print(typemod(8.2223, 'np.float64'))
48    <type 'numpy.float64'>
49    >>> print(typemod(8.2223, 'tuple'))
50    <type 'tuple'>
51    """
52    fname='typemod'
53
54    if typeval == 'int' or typeval == 'I':
55        return int(value)
56    elif typeval == 'long':
57        return long(value)
58    elif typeval == 'float' or typeval == 'F' or typeval == 'R':
59        return float(value)
60    elif typeval == 'complex':
61        return complex(value)
62    elif typeval == 'str' or typeval == 'S':
63        return str(value)
64    elif typeval == 'bool':
65        return bool(value)
66    elif typeval == 'B':
67        return Str_Bool(value)
68    elif typeval == 'list':
69        newval = []
70        newval.append(value)
71        return newval
72    elif typeval == 'dic':
73        newval = {}
74        newval[value] = value
75        return newval
76    elif typeval == 'tuple':
77        newv = []
78        newv.append(value)
79        newval = tuple(newv)
80        return newval
81    elif typeval == 'np.int8':
82        return np.int8(value)
83    elif typeval == 'np.int16':
84        return np.int16(value)
85    elif typeval == 'np.int32':
86        return np.int32(value)
87    elif typeval == 'np.int64':
88        return np.int64(value)
89    elif typeval == 'np.uint8':
90        return np.uint8(value)
91    elif typeval == 'np.uint16':
92        return np.uint16(value)
93    elif typeval == 'np.np.uint32':
94        return np.uint32(value)
95    elif typeval == 'np.uint64':
96        return np.uint64(value)
97    elif typeval == 'np.float' or typeval == 'R':
98        return np.float(value)
99    elif typeval == 'np.float16':
100        return np.float16(value)
101    elif typeval == 'np.float32':
102        return np.float32(value)
103    elif typeval == 'float32':
104        return np.float32(value)
105    elif typeval == 'np.float64' or typeval == 'D':
106        return np.float64(value)
107    elif typeval == 'float64':
108        return np.float64(value)
109    elif typeval == 'np.complex':
110        return np.complex(value)
111    elif typeval == 'np.complex64':
112        return np.complex64(value)
113    elif typeval == 'np.complex128':
114        return np.complex128(value)
115    else:
116        print errormsg
117        print fname + ':  data type "' + typeval + '" is not ready !'
118        print errormsg
119        quit(-1)
120
121    return
122
123def str_list_k(string, cdiv, kind):
124    """ Function to obtain a list of types of values from a string giving a split character
125      string= String from which to obtain a list ('None' for None)
126      cdiv= character to use to split the string
127      kind= kind of desired value (as string like: 'np.float', 'int', 'np.float64', ....)
128    >>> str_list_k('1:@#:$:56', ':', 'S')
129    ['1', '@#', '$', '56']
130    >>> str_list_k('1:3.4:12.3', ':', 'np.float64')
131    [1.0, 3.3999999999999999, 12.300000000000001]
132    """
133    fname = 'str_list'
134
135    if string.find(cdiv) != -1:
136        listv = string.split(cdiv)
137    else:
138        if string == 'None':
139            listv = None
140        else:
141            listv = [string]
142
143    if listv is not None:
144        finalist = []
145        for lv in listv:
146            finalist.append(typemod(lv, kind))
147    else:
148        finalist = None
149
150    return finalist
151
152def stringS_dictvar(stringv, Dc=',', DVs=':'):
153    """ Function to provide a dictionary from a String which list of DVs separated
154          [key] [value] couples
155      stringv: String with a dictionary content as [key1][DVs][val1][Dc][key2][DVs][Val2][Dc]...
156      Sc: character to separate key values
157      DVs: character to separate key-value couples
158    >>> stringS_dictvar('i:1,iv:4,vii:7',',',':')
159    {'i': '1', 'vii': '7', 'iv': '4'}
160    """
161    fname = 'stringS_dictvar'
162
163    if stringv.count(Dc) != 0:
164        strvv = stringv.split(Dc)
165    else:
166        strvv = [stringv]
167
168    dictv = {}
169    for Svals in strvv:
170        keyn = Svals.split(DVs)[0]
171        valn = Svals.split(DVs)[1]
172        dictv[keyn] = valn
173
174    return dictv
175
176def set_attribute(ncvar, attrname, attrvalue):
177    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
178    ncvar = object netcdf variable
179    attrname = name of the attribute
180    attrvalue = value of the attribute
181    """
182    import numpy as np
183    from netCDF4 import Dataset as NetCDFFile
184
185    attvar = ncvar.ncattrs()
186    if searchInlist(attvar, attrname):
187        attr = ncvar.delncattr(attrname)
188
189    attr = ncvar.setncattr(attrname, attrvalue)
190
191    return ncvar
192
193def set_attributek(ncv, attrname, attrval, attrkind):
194    """ Sets a value of an attribute of a netCDF variable with a kind. Removes previous attribute value if exists
195    ncvar = object netcdf variable
196    attrname = name of the attribute
197    attrvalue = value of the attribute
198    atrtrkind = kind of attribute: 'S', string ('!' as spaces); 'U', unicode ('!' as spaces); 'I', integer;
199      'Inp32', numpy integer 32; 'R', ot 'F' real; ' npfloat', np.float; 'D', np.float64
200    """
201    fname = 'set_attributek'
202    validk = {'S': 'string', 'U': 'unicode', 'I': 'integer',                         \
203      'Inp32': 'integer long (np.int32)', 'F': 'float', 'R': 'float',                \
204      'npfloat': 'np.float', 'D': 'float long (np.float64)'}
205
206    if type(attrkind) == type('s'):
207        if attrkind == 'S':
208            attrvalue = str(attrval.replace('!', ' '))
209        elif attrkind == 'U':
210            attrvalue = unicode(attrval.replace('!',' '))
211        elif attrkind == 'I':
212            attrvalue = np.int(attrval)
213        elif attrkind == 'R' or attrkind == 'F' :
214            attrvalue = float(attrval)
215        elif attrkind == 'npfloat':
216            attrvalue = np.float(attrval)
217        elif attrkind == 'D':
218            attrvalue = np.float64(attrval)
219        else:
220            print errormsg
221            print '  ' + fname + ": '" + attrkind + "' kind of attribute is not ready!"
222            print '  valid values: _______'
223            for key in validk.keys():
224                print '    ', key,':', validk[key]
225            quit(-1)
226    else:
227        if attrkind == type(str('a')):
228            attrvalue = str(attrval.replace('!', ' '))
229        elif attrkind == type(unicode('a')):
230            attrvalue = unicode(attrval.replace('!',' '))
231        elif attrkind == type(np.int(1)):
232            attrvalue = np.int(attrval)
233        elif attrkind == np.dtype('i'):
234            attrvalue = np.int32(attrval)
235        elif attrkind == type(float(1.)):
236            attrvalue = float(attrval)
237        elif attrkind == type(np.float(1.)):
238            attrvalue = np.float(attrval)
239        elif attrkind == np.dtype('float32'):
240            attrvalue = np.array(attrval, dtype='f')
241        elif attrkind == type(np.float32(1.)):
242            attrvalue = np.float32(attrval)
243        elif attrkind == type(np.float64(1.)):
244            attrvalue = np.float64(attrval)
245        elif attrkind == type(np.array([1.,2.])):
246            attrvalue = attrval
247        else:
248            print errormsg
249            print '  ' + fname + ": '" + attrkind + "' kind of attribute is not ready!"
250            print '  valid values: _______'
251            for key in validk.keys():
252                print '    ', key,':', validk[key]
253            quit(-1)
254
255    attvar = ncv.ncattrs()
256    if searchInlist(attvar, attrname):
257        attr = ncv.delncattr(attrname)
258
259    if attrname == 'original_subroutines_author': 
260        attrvalue = 'Cindy Bruyere'
261    attr = ncv.setncattr(attrname, attrvalue)
262       
263    return attr
264
265
266def basicvardef(varobj, vstname, vlname, vunits):
267    """ Function to give the basic attributes to a variable
268    varobj= netCDF variable object
269    vstname= standard name of the variable
270    vlname= long name of the variable
271    vunits= units of the variable
272    """
273    attr = varobj.setncattr('standard_name', vstname)
274    attr = varobj.setncattr('long_name', vlname)
275    attr = varobj.setncattr('units', vunits)
276
277    return
278
279def remove_NONascii(string):
280    """ Function to remove that characters which are not in the standard 127 ASCII
281      string= string to transform
282    >>> remove_NONascii('Lluís')
283    Lluis
284    """
285    fname = 'remove_NONascii'
286
287    newstring = string
288
289    RTFchar= ['á', 'é', 'í', 'ó', 'ú', 'à', 'Ú', 'ì', 'ò', 'ù', 'â', 'ê', 'î', 'ÃŽ',  \
290      'û', 'À', 'ë', 'ï', 'ö', 'ÃŒ', 'ç', 'ñ','Ê', 'œ', 'Á', 'É', 'Í', 'Ó', 'Ú', 'À', \
291      'È', 'Ì', 'Ò', 'Ù', 'Â', 'Ê', 'Î', 'Ô', 'Û', 'Ä', 'Ë', 'Ï', 'Ö', 'Ü', 'Ç', 'Ñ',\
292      'Æ', 'Œ', '\n', '\t']
293    ASCchar= ['a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o',  \
294      'u', 'a', 'e', 'i', 'o', 'u', 'c', 'n','ae', 'oe', 'A', 'E', 'I', 'O', 'U', 'A', \
295      'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'C', 'N',\
296      'AE', 'OE', '', ' ']
297
298    Nchars = len(RTFchar)
299    for ichar in range(Nchars):
300        foundchar = string.find(RTFchar[ichar])
301        if foundchar != -1:
302            newstring = newstring.replace(RTFchar[ichar], ASCchar[ichar])
303
304    return newstring
305
306def read_description(fdobs, dbg):
307    """ reads the description file of the observational data-set
308    fdobs= descriptive observational data-set
309    dbg= boolean argument for debugging
310    * Each station should have a 'description.dat' file with:
311      institution=Institution who creates the data
312      department=Department within the institution
313      scientists=names of the data producers
314      contact=contact of the data producers
315      description=description of the observations
316      acknowledgement=sentence of acknowlegement
317      comment=comment for the measurements
318
319      MissingValue='|' list of ASCII values for missing values within the data
320        (as they appear!, 'empty' for no value at all)
321      comment=comments
322
323      varN='|' list of variable names
324      varLN='|' list of long variable names
325      varU='|' list units of the variables
326      varBUFR='|' list BUFR code of the variables
327      varTYPE='|' list of variable types ('D', 'F', 'I', 'I64', 'S')
328      varOPER='|' list of operations to do to the variables to meet their units ([oper],[val])
329        [oper]:
330          -, nothing
331          sumc, add [val]
332          subc, rest [val]
333          mulc, multiply by [val]
334          divc, divide by [val]
335          rmchar,[val],[pos], remove [val] characters from [pos]='B', beginning, 'E', end
336          repl,[char1]:[val1];...;[charM]:[valM], replace a given character [chari]
337            by a value [vali]
338      NAMElon=name of the variable with the longitude (x position)
339      NAMElat=name of the variable with the latitude (y position)
340      NAMEheight=ind_alt
341      NAMEtime=name of the varibale with the time
342      FMTtime=format of the time (as in 'C', 'CFtime' for already CF-like time)
343
344    """
345    fname = 'read_description'
346
347    descobj = open(fdobs, 'r')
348    desc = {}
349
350    namevalues = []
351
352    for line in descobj:
353        if line[0:1] != '#' and len(line) > 1 :
354            descn = remove_NONascii(line.split('=')[0])
355            descv = remove_NONascii(line.split('=')[1])
356            namevalues.append(descn)
357            if descn[0:3] != 'var':
358                if descn != 'MissingValue':
359                    desc[descn] = descv
360                else:
361                    desc[descn] = []
362                    for dn in descv.split('|'): desc[descn].append(dn)
363                    print '  ' + fname + ': missing values found:',desc[descn]
364            elif descn[0:3] == 'var':
365                desc[descn] = descv.split('|')
366            elif descn[0:4] == 'NAME':
367                desc[descn] = descv
368            elif descn[0:3] == 'FMT':
369                desc[descn] = descv
370    if not desc.has_key('varOPER'): desc['varOPER'] = None
371
372    if dbg:
373        Nvars = len(desc['varN'])
374        print '  ' + fname + ": description content of '" + fdobs + "'______________"
375        for varn in namevalues:
376            if varn[0:3] != 'var':
377                print '    ' + varn + ':',desc[varn]
378            elif varn == 'varN':
379                print '    * Variables:'
380                for ivar in range(Nvars):
381                    varname = desc['varN'][ivar]
382                    varLname = desc['varLN'][ivar]
383                    varunits = desc['varU'][ivar]
384                    if desc.has_key('varBUFR'): 
385                        varbufr = desc['varBUFR'][ivar]
386                    else:
387                        varbufr = None
388                    if desc['varOPER'] is not None:
389                        opv = desc['varOPER'][ivar]
390                    else:
391                        opv = None
392                    print '      ', ivar, varname+':',varLname,'[',varunits, \
393                       ']','bufr code:',varbufr,'oper:',opv
394
395    descobj.close()
396
397    return desc
398
399def value_fmt(val, miss, op, fmt):
400    """ Function to transform an ASCII value to a given format
401      val= value to transform
402      miss= list of possible missing values
403      op= operation to perform to the value
404      fmt= format to take:
405        'D': float double precission
406        'F': float
407        'I': integer
408        'I64': 64-bits integer
409        'S': string
410    >>> value_fmt('9876.12', '-999', 'F')
411    9876.12
412    """
413
414    fname = 'value_fmt'
415
416    aopers = ['sumc','subc','mulc','divc', 'rmchar', 'repl']
417
418    fmts = ['D', 'F', 'I', 'I64', 'S']
419    Nfmts = len(fmts)
420
421    if not searchInlist(miss,val):
422        if searchInlist(miss,'empty') and len(val) == 0: 
423            newval = None
424        else:
425            if op != '-':
426                opern = op.split(',')[0]
427                operv = op.split(',')[1]
428
429                if not searchInlist(aopers,opern):
430                    print errormsg
431                    print '  ' + fname + ": operation '" + opern + "' not ready!!"
432                    print '    availables:',aopers
433                    quit(-1)
434            else:
435                opern = 'sumc'
436                operv = '0'
437
438            if not searchInlist(fmts, fmt):
439                print errormsg
440                print '  ' + fname + ": format '" + fmt + "' not ready !!"
441                quit(-1)
442            else:
443                if fmt == 'D':
444                    opv = np.float32(operv)
445                    if opern == 'sumc': newval = np.float32(val) + opv
446                    elif opern == 'subc': newval = np.float32(val) - opv
447                    elif opern == 'mulc': newval = np.float32(val) * opv
448                    elif opern == 'divc': newval = np.float32(val) / opv
449                    elif opern == 'rmchar': 
450                        opv = int(operv)
451                        Lval = len(val)
452                        if op.split(',')[2] == 'B':
453                            newval = np.float32(val[opv:Lval+1])
454                        elif op.split(',')[2] == 'E':
455                            newval = np.float32(val[Lval-opv:Lval])
456                        else:
457                            print errormsg
458                            print '  ' + fname + ": operation '" + opern + "' not " +\
459                              " work with '" + op.split(',')[2] + "' !!"
460                            quit(-1)
461                    elif opern == 'repl':
462                        replvals = stringS_dictvar(operv, Dc=';', DVs=':')
463                        replaced = False
464                        for repls in replvals.keys():
465                            if val == repls: 
466                                newval = np.float32(replvals[val])
467                                replaced = True
468                                break
469                        if not replaced: newval = np.float32(val)
470
471                elif fmt == 'F':
472                    opv = np.float(operv)
473                    if opern == 'sumc': newval = np.float(val) + opv
474                    elif opern == 'subc': newval = np.float(val) - opv
475                    elif opern == 'mulc': newval = np.float(val) * opv
476                    elif opern == 'divc': newval = np.float(val) / opv
477                    elif opern == 'rmchar': 
478                        opv = int(operv)
479                        Lval = len(val)
480                        if op.split(',')[2] == 'B':
481                            newval = np.float(val[opv:Lval+1])
482                        elif op.split(',')[2] == 'E':
483                            newval = np.float(val[0:Lval-opv])
484                        else:
485                            print errormsg
486                            print '  ' + fname + ": operation '" + opern + "' not " +\
487                              " work with '" + op.split(',')[2] + "' !!"
488                            quit(-1)
489                    elif opern == 'repl':
490                        replvals = stringS_dictvar(operv, Dc=';', DVs=':')
491                        replaced = False
492                        for repls in replvals.keys():
493                            if val == repls: 
494                                newval =  np.float(replvals[val])
495                                replaced = True
496                                break
497                        if not replaced: newval = np.float(val)
498
499                elif fmt == 'I':
500                    if opern == 'sumc': newval = int(val) + int(operv)
501                    elif opern == 'subc': newval = int(val) - int(operv)
502                    elif opern == 'mulc': newval = int(val) * int(operv)
503                    elif opern == 'divc': newval = int(val) / int(operv)
504                    elif opern == 'rmchar': 
505                        opv = int(operv)
506                        Lval = len(val)
507                        if op.split(',')[2] == 'B':
508                            newval = int(val[opv:Lval+1])
509                        elif op.split(',')[2] == 'E':
510                            newval = int(val[0:Lval-opv])
511                        else:
512                            print errormsg
513                            print '  ' + fname + ": operation '" + opern + "' not " +\
514                              " work with '" + op.split(',')[2] + "' !!"
515                            quit(-1)
516                    elif opern == 'repl':
517                        replvals = stringS_dictvar(operv, Dc=';', DVs=':')
518                        replaced = False
519                        for repls in replvals.keys():
520                            if val == repls: 
521                                newval = int(replvals[val])
522                                replaced = True
523                                break
524                        if not replaced: newval = int(val)
525                       
526                elif fmt == 'I64':
527                    opv = np.int64(operv)
528                    if opern == 'sumc': newval = np.int64(val) + opv
529                    elif opern == 'subc': newval = np.int64(val) - opv
530                    elif opern == 'mulc': newval = np.int64(val) * opv
531                    elif opern == 'divc': newval = np.int64(val) / opv
532                    elif opern == 'rmchar': 
533                        opv = int(operv)
534                        Lval = len(val)
535                        if op.split(',')[2] == 'B':
536                            newval = np.int64(val[opv:Lval+1])
537                        elif op.split(',')[2] == 'E':
538                            newval = np.int64(val[0:Lval-opv])
539                        else:
540                            print errormsg
541                            print '  ' + fname + ": operation '" + opern + "' not " +\
542                              " work with '" + op.split(',')[2] + "' !!"
543                            quit(-1)
544                    elif opern == 'repl':
545                        replvals = stringS_dictvar(operv, Dc=';', DVs=':')
546                        replaced = False
547                        for repls in replvals.keys():
548                            if val == repls: 
549                                newval = np.int64(replvals[val])
550                                replaced = True
551                                break
552                        if not replaced: newval = np.int64(val)
553
554                elif fmt == 'S':
555                    if opern == 'rmchar': 
556                        opv = int(operv)
557                        Lval = len(val)
558                        if op.split(',')[2] == 'B':
559                            newval = val[opv:Lval+1]
560                        elif op.split(',')[2] == 'E':
561                            newval = val[0:Lval-opv]
562                        else:
563                            print errormsg
564                            print '  ' + fname + ": operation '" + opern + "' not " +\
565                              " work with '" + op.split(',')[2] + "' !!"
566                            quit(-1)
567                    elif opern == 'repl':
568                        replvals = stringS_dictvar(operv, Dc=';', DVs=':')
569                        replaced = False
570                        for repls in replvals.keys():
571                            if val == repls: 
572                                newval = replvals[val]
573                                replaced = True
574                                break
575                        if not replaced: newval = val
576
577                    else:
578                        newval = val
579
580    else:
581        newval = None
582
583    return newval
584
585def getting_fixedline(line, cuts, types, dbg=False):
586    """ Function to get the values from a line of text with fixed lenght of different values
587      line: line with values
588      cuts: character number where a value ends
589      types: consecutive type of values
590        'I': integer
591        'R': real
592        'D': float64
593        'S': string
594        'B': boolean
595      dbg: debug mode (default False)
596    >>> Sline='   87007 03012015  25.6   6.4             9.4    5   15'
597    >>> getting_fixedline(Sline, [8, 17, 23, 29, 36, 40, 45, 50], ['I', 'R', 'R', 'R', 'I', 'R', 'R', 'R'])
598    [87007, 3012015.0, 25.6, 6.4, -99999, -99999, 9.4, 5.0, 15.0]
599    """
600    fname = 'getting_fixedline'
601
602    if len(cuts) + 1 != len(types):
603        print errormsg
604        print '  ' + fname + ': The number of types :', len(types), 'must be +1',    \
605          'number of cuts:', len(cuts)
606        quit(-1)
607
608    values = []
609    val = line[0:cuts[0]]
610    if len(val.replace(' ','')) >= 1:
611        values.append(typemod(val, types[0]))
612    else:
613        if types[0] == 'I': values.append(fillValueI)
614        elif types[0] == 'R': values.append(fillValueF)
615        elif types[0] == 'D': values.append(fillValueF)
616        elif types[0] == 'S': values.append(fillValueS)
617        elif types[0] == 'B': values.append(fillValueB)
618
619    Ncuts = len(cuts)
620    for ic in range(1,Ncuts):
621        val = line[cuts[ic-1]:cuts[ic]]
622        if dbg: print ic, ':', val, '-->', types[ic]
623        if len(val.replace(' ','')) >= 1:
624            values.append(typemod(val, types[ic]))
625        else:
626            if types[ic] == 'I': values.append(fillValueI)
627            elif types[ic] == 'R': values.append(fillValueF)
628            elif types[ic] == 'D': values.append(fillValueF)
629            elif types[ic] == 'S': values.append(fillValueS)
630            elif types[ic] == 'B': values.append(fillValueB)
631
632    # Last value
633    Lline = len(line)
634    val = line[cuts[Ncuts-1]:Lline]
635    if len(val.replace(' ','')) >= 1:
636        values.append(typemod(val, types[Ncuts]))
637    else:
638        if types[Ncuts] == 'I': values.append(fillValueI)
639        elif types[Ncuts] == 'R': values.append(fillValueF)
640        elif types[Ncuts] == 'D': values.append(fillValueF)
641        elif types[Ncuts] == 'S': values.append(fillValueS)
642        elif types[Ncuts] == 'B': values.append(fillValueB)
643
644    return values
645
646
647def getting_fixedline_NOk(line, cuts, miss, dbg=False):
648    """ Function to get the values from a line of text with fixed lenght of
649        different values without types
650      line: line with values
651      cuts: character number where a value ends
652      miss: character form issed values
653      dbg: debug mode (default False)
654    >>> Sline='   87007 03012015  25.6   6.4             9.4    5   15'
655    >>> getting_fixedline_NOk(Sline, [8, 17, 23, 29, 36, 40, 45, 50], '#')
656    ['87007', '03012015', '25.6', '6.4', '#', '#', '9.4', '5', '15']
657    """
658    fname = 'getting_fixedline_NOk'
659
660    values = []
661    val = line[0:cuts[0]]
662    if len(val.replace(' ','')) >= 1:
663        values.append(val.replace(' ',''))
664    else:
665        values.append(miss)
666
667    Ncuts = len(cuts)
668    for ic in range(1,Ncuts):
669        val = line[cuts[ic-1]:cuts[ic]]
670        if dbg: print ic, ':', val
671        if len(val.replace(' ','')) >= 1:
672            values.append(val.replace(' ',''))
673        else:
674            values.append(miss)
675
676    # Last value
677    Lline = len(line)
678    val = line[cuts[Ncuts-1]:Lline]
679    if len(val.replace(' ','')) >= 1:
680        values.append(val.replace(' ',''))
681    else:
682        values.append(miss)
683
684    return values
685
686def read_datavalues(dataf, comchar, colchar, fmt, jBl, oper, miss, varns, dbg):
687    """ Function to read from an ASCII file values in column
688      dataf= data file
689      comchar= list of the characters indicating comment in the file
690      colchar= character which indicate end of value in the column
691      dbg= debug mode or not
692      fmt= list of kind of values to be found
693      jBl= number of lines to jump from the beginning of file
694      oper= list of operations to perform
695      miss= missing value
696      varns= list of name of the variables to find
697    """
698
699    fname = 'read_datavalues'
700
701    ofile = open(dataf, 'r')
702    Nvals = len(fmt)
703
704    if oper is None:
705        opers = []
706        for ioper in range(Nvals):
707            opers.append('-')
708    else:
709        opers = oper
710
711    finalvalues = {}
712
713    iline = 0
714    for line in ofile:
715        line = line.replace('\n','').replace(chr(13),'')
716        if not searchInlist(comchar,line[0:1]) and len(line) > 1 and iline > jBl-1:
717            if dbg: print line
718            # Getting values
719            if colchar[0:4] != 'None': 
720                values0 = line.split(colchar)
721            else:
722                ltypes=str_list_k(colchar[5:len(colchar)+1],',','I')
723                values0 = getting_fixedline_NOk(line, ltypes, miss[0], dbg=dbg)
724# Removing no-value columns
725            values = []
726            for iv in values0:
727                if len(iv) > 0: values.append(iv)
728
729            Nvalues = len(values)
730# Checkings for wierd characters at the end of lines (use it to check)
731#            if values[Nvalues-1][0:4] == '-999':
732#                print line,'last v:',values[Nvalues-1],'len:',len(values[Nvalues-1])
733#                for ic in range(len(values[Nvalues-1])):
734#                    print ic,ord(values[Nvalues-1][ic:ic+1])
735#                quit()
736
737            if len(values[Nvalues-1]) == 0:
738                Nvalues = Nvalues - 1
739               
740            if Nvalues != Nvals:
741                print warnmsg
742                print '  ' + fname + ': number of formats:',Nvals,' and number of ', \
743                  'values:',Nvalues,' with split character *' + colchar +            \
744                  '* does not coincide!!'
745                print '    * what is found is ________'
746                if Nvalues > Nvals:
747                    Nshow = Nvals
748                    for ivar in range(Nshow):
749                        print '    ',varns[ivar],'fmt:',fmt[ivar],'value:',values[ivar]
750                    print '      missing formats for:',values[Nshow:Nvalues+1]
751                    print '      values not considered, continue'
752                else:
753                    Nshow = Nvalues
754                    for ivar in range(Nshow):
755                        print '   ',varns[ivar],'fmt:',fmt[ivar],'value:',values[ivar]
756                    print '      excess of formats:',fmt[Nshow:Nvals+1]
757                    quit(-1)
758
759# Reading and transforming values
760            if dbg: print '  ' + fname + ': values found _________'
761
762            for ivar in range(Nvals):
763                if dbg: 
764                    print iline, ',', ivar, '  ', varns[ivar],'value:',values[ivar], \
765                      miss,opers[ivar], fmt[ivar]
766
767                if iline == 0:
768                    listvals = []
769                    listvals.append(value_fmt(values[ivar], miss, opers[ivar],       \
770                      fmt[ivar]))
771                    finalvalues[varns[ivar]] = listvals
772                else:
773                    listvals = finalvalues[varns[ivar]]
774                    listvals.append(value_fmt(values[ivar], miss, opers[ivar],       \
775                      fmt[ivar]))
776                    finalvalues[varns[ivar]] = listvals
777        else:
778# First line without values
779            if iline == 0: iline = -1
780   
781        iline = iline + 1
782
783    ofile.close()
784
785    return finalvalues
786
787def writing_str_nc(varo, values, Lchar):
788    """ Function to write string values in a netCDF variable as a chain of 1char values
789    varo= netCDF variable object
790    values = list of values to introduce
791    Lchar = length of the string in the netCDF file
792    """
793
794    Nvals = len(values)
795
796    for iv in range(Nvals):   
797        stringv=values[iv] 
798        charvals = np.chararray(Lchar)
799        Lstr = len(stringv)
800        charvals[Lstr:Lchar] = ''
801
802        for ich in range(Lstr):
803            charvals[ich] = stringv[ich:ich+1]
804
805        varo[iv,:] = charvals
806
807    return
808
809def Stringtimes_CF(tvals, fmt, Srefdate, tunits, dbg):
810    """ Function to transform a given data in String formats to a CF date
811      tvals= string temporal values
812      fmt= format of the the time values
813      Srefdate= reference date in [YYYY][MM][DD][HH][MI][SS] format
814      tunits= units to use ('weeks', 'days', 'hours', 'minutes', 'seconds')
815      dbg= debug
816    >>> Stringtimes_CF(['19760217082712','20150213101837'], '%Y%m%d%H%M%S',
817      '19491201000000', 'hours', False)
818    [229784.45333333  571570.31027778]
819    """
820    import datetime as dt
821
822    fname = 'Stringtimes'
823
824    dimt = len(tvals)
825
826    yrref = int(Srefdate[0:4])
827    monref = int(Srefdate[4:6])
828    dayref = int(Srefdate[6:8])
829    horref = int(Srefdate[8:10])
830    minref = int(Srefdate[10:12])
831    secref = int(Srefdate[12:14])
832    refdate = dt.datetime( yrref, monref, dayref, horref, minref, secref)
833
834    cftimes = np.zeros((dimt), dtype=np.float)
835
836    Nfmt=len(fmt.split('%'))
837
838    if dbg: print '  ' + fname + ': fmt=',fmt,'refdate:',Srefdate,'uits:',tunits,    \
839      'date dt_days dt_time deltaseconds CFtime _______'
840    for it in range(dimt):
841
842# Removing excess of mili-seconds (up to 6 decimals)
843        if fmt.split('%')[Nfmt-1] == 'f':
844            tpoints = tvals[it].split('.')
845            if len(tpoints[len(tpoints)-1]) > 6:
846                milisec = '{0:.6f}'.format(np.float('0.'+tpoints[len(tpoints)-1]))[0:7]
847                newtval = ''
848                for ipt in range(len(tpoints)-1):
849                    newtval = newtval + tpoints[ipt] + '.'
850                newtval = newtval + str(milisec)[2:len(milisec)+1]
851            else:
852                newtval = tvals[it]
853            tval = dt.datetime.strptime(newtval, fmt)
854        else:
855            tval = dt.datetime.strptime(tvals[it], fmt)
856
857        deltatime = tval - refdate
858        deltaseconds = deltatime.days*24.*3600. + deltatime.seconds +                \
859          deltatime.microseconds/100000.
860        if tunits == 'weeks':
861            deltat = 7.*24.*3600.
862        elif tunits == 'days':
863            deltat = 24.*3600.
864        elif tunits == 'hours':
865            deltat = 3600.
866        elif tunits == 'minutes':
867            deltat = 60.
868        elif tunits == 'seconds':
869            deltat = 1.
870        else:
871            print errormsg
872            print '  ' + fname + ": time units '" + tunits + "' not ready !!"
873            quit(-1)
874
875        cftimes[it] = deltaseconds / deltat
876        if dbg:
877            print '  ' + tvals[it], deltatime, deltaseconds, cftimes[it]
878
879    return cftimes
880
881def adding_complementary(onc, dscn, okind, dvalues, tvals, refCFt, CFtu, Nd, dbg):
882    """ Function to add complementary variables as function of the observational type
883      onc= netcdf objecto file to add the variables
884      dscn= description dictionary
885      okind= observational kind
886      dvalues= values
887      tvals= CF time values
888      refCFt= reference time of CF time (in [YYYY][MM][DD][HH][MI][SS])
889      CFtu= CF time units
890      Nd= size of the domain
891      dbg= debugging flag
892    """
893    import numpy.ma as ma
894
895    fname = 'adding_complementary'
896 
897# Kind of observations which require de integer lon/lat (for the 2D map)
898    map2D=['multi-points', 'trajectory']
899
900    SrefCFt = refCFt[0:4] +'-'+ refCFt[4:6] +'-'+ refCFt[6:8] + ' ' + refCFt[8:10] + \
901      ':'+ refCFt[10:12] +':'+ refCFt[12:14]
902
903    if dscn['NAMElon'] == '-' or dscn['NAMElat'] == '-':
904        print errormsg
905        print '  ' + fname + ": to complement a '" + okind + "' observation kind " + \
906          " a given longitude ('NAMElon':",dscn['NAMElon'],") and latitude ('" +     \
907          "'NAMElat:'", dscn['NAMElat'],') from the data has to be provided and ' +  \
908          'any are given !!'
909        quit(-1)
910
911    if not dvalues.has_key(dscn['NAMElon']) or not dvalues.has_key(dscn['NAMElat']):
912        print errormsg
913        print '  ' + fname + ": observations do not have 'NAMElon':",                \
914          dscn['NAMElon'],"and/or 'NAMElat:'", dscn['NAMElat'],' !!'
915        print '    available data:',dvalues.keys()
916        quit(-1)
917
918    if okind == 'trajectory':
919        if dscn['NAMEheight'] == '-':
920            print warnmsg
921            print '  ' + fname + ": to complement a '" + okind + "' observation " +  \
922              "kind a given height ('NAMEheight':",dscn['NAMEheight'],"') might " +  \
923              'be provided and any is given !!'
924            quit(-1)
925
926        if not dvalues.has_key(dscn['NAMEheight']):
927            print errormsg
928            print '  ' + fname + ": observations do not have 'NAMEtime':",           \
929              dscn['NAMEtime'],' !!'
930            print '    available data:',dvalues.keys()
931            quit(-1)
932
933    if searchInlist(map2D, okind):
934# A new 2D map with the number of observation will be added for that 'NAMElon'
935#   and 'NAMElat' are necessary. A NdxNd domain space size will be used.
936        objfile.createDimension('lon2D',Nd)
937        objfile.createDimension('lat2D',Nd)
938        lonvals = ma.masked_equal(dvalues[dscn['NAMElon']], None)
939        latvals = ma.masked_equal(dvalues[dscn['NAMElat']], None)
940
941        minlon = min(lonvals)
942        maxlon = max(lonvals)
943        minlat = min(latvals)
944        maxlat = max(latvals)
945
946        blon = (maxlon - minlon)/(Nd-1)
947        blat = (maxlat - minlat)/(Nd-1)
948
949        newvar = onc.createVariable( 'lon2D', 'f8', ('lon2D'))
950        basicvardef(newvar, 'longitude', 'longitude map observations','degrees_East')
951        newvar[:] = minlon + np.arange(Nd)*blon
952        newattr = set_attribute(newvar, 'axis', 'X')
953
954        newvar = onc.createVariable( 'lat2D', 'f8', ('lat2D'))
955        basicvardef(newvar, 'latitude', 'latitude map observations', 'degrees_North')
956        newvar[:] = minlat + np.arange(Nd)*blat
957        newattr = set_attribute(newvar, 'axis', 'Y')
958
959        if dbg: 
960            print '  ' + fname + ': minlon=',minlon,'maxlon=',maxlon
961            print '  ' + fname + ': minlat=',minlat,'maxlat=',maxlat
962            print '  ' + fname + ': precission on x-axis=', blon*(Nd-1), 'y-axis=',  \
963              blat*(Nd-1)
964       
965    if okind == 'multi-points':
966        map2D = np.ones((Nd, Nd), dtype=np.float)*fillValueI
967
968        dimt = len(tvals)
969        Nlost = 0
970        for it in range(dimt):
971            lon = dvalues[dscn['NAMElon']][it]
972            lat = dvalues[dscn['NAMElat']][it]
973            if lon is not None and lat is not None:
974                ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
975                ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
976
977                if map2D[ilat,ilon] == fillValueI: 
978                    map2D[ilat,ilon] = 1
979                else:
980                    map2D[ilat,ilon] = map2D[ilat,ilon] + 1
981                if dbg: print it, lon, lat, ilon, ilat, map2D[ilat,ilon]
982
983        newvar = onc.createVariable( 'mapobs', 'f4', ('lat2D', 'lon2D'),             \
984          fill_value = fillValueI)
985        basicvardef(newvar, 'map_observations', 'number of observations', '-')
986        newvar[:] = map2D
987        newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D')
988
989    elif okind == 'trajectory':
990# A new 2D map with the trajectory 'NAMElon' and 'NAMElat' and maybe 'NAMEheight'
991#   are necessary. A NdxNdxNd domain space size will be used. Using time as
992#   reference variable
993        if dscn['NAMEheight'] == '-':
994# No height
995            map2D = np.ones((Nd, Nd), dtype=np.float)*fillValueI
996
997            dimt = len(tvals)
998            Nlost = 0
999            if dbg: print ' time-step lon lat ix iy passes _______'
1000            for it in range(dimt):
1001                lon = dvalues[dscn['NAMElon']][it]
1002                lat = dvalues[dscn['NAMElat']][it]
1003                if lon is  not None and lat is not None:
1004                    ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
1005                    ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
1006
1007                    if map2D[ilat,ilon] == fillValueI: 
1008                        map2D[ilat,ilon] = 1
1009                    else:
1010                        map2D[ilat,ilon] = map2D[ilat,ilon] + 1
1011                    if dbg: print it, lon, lat, ilon, ilat, map2D[ilat,ilon]
1012
1013            newvar = onc.createVariable( 'trjobs', 'i', ('lat2D', 'lon2D'),          \
1014              fill_value = fillValueI)
1015            basicvardef(newvar, 'trajectory_observations', 'number of passes', '-' )
1016            newvar[:] = map2D
1017            newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D')
1018
1019        else:
1020            ivn = 0
1021            for vn in dscn['varN']:
1022                if vn == dscn['NAMEheight']: 
1023                    zu = dscn['varU'][ivn]
1024                    break
1025                ivn = ivn + 1
1026                   
1027            objfile.createDimension('z2D',Nd)
1028            zvals = ma.masked_equal(dvalues[dscn['NAMEheight']], None)
1029            minz = min(zvals)
1030            maxz = max(zvals)
1031
1032            bz = (maxz - minz)/(Nd-1)
1033
1034            newvar = onc.createVariable( 'z2D', 'f8', ('z2D'))
1035            basicvardef(newvar, 'z2D', 'z-coordinate map observations', zu)
1036            newvar[:] = minz + np.arange(Nd)*bz
1037            newattr = set_attribute(newvar, 'axis', 'Z')
1038
1039            if dbg:
1040                print '  ' + fname + ': zmin=',minz,zu,'zmax=',maxz,zu
1041                print '  ' + fname + ': precission on z-axis=', bz*(Nd-1), zu
1042
1043            map3D = np.ones((Nd, Nd, Nd), dtype=int)*fillValueI
1044            dimt = len(tvals)
1045            Nlost = 0
1046            if dbg: print ' time-step lon lat z ix iy iz passes _______'
1047            for it in range(dimt):
1048                lon = dvalues[dscn['NAMElon']][it]
1049                lat = dvalues[dscn['NAMElat']][it]
1050                z = dvalues[dscn['NAMEheight']][it]
1051                if lon is  not None and lat is not None and z is not None:
1052                    ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon))
1053                    ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat))
1054                    iz = int((Nd-1)*(z - minz)/(maxz - minz))
1055
1056                    if map3D[iz,ilat,ilon] == fillValueI: 
1057                        map3D[iz,ilat,ilon] = 1
1058                    else:
1059                        map3D[iz,ilat,ilon] = map3D[iz,ilat,ilon] + 1
1060                    if dbg: print it, lon, lat, z, ilon, ilat, iz, map3D[iz,ilat,ilon]
1061
1062            newvar = onc.createVariable( 'trjobs', 'i', ('z2D', 'lat2D', 'lon2D'),   \
1063              fill_value = fillValueI)
1064            basicvardef(newvar, 'trajectory_observations', 'number of passes', '-')
1065            newvar[:] = map3D
1066            newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D z2D')
1067
1068    onc.sync()
1069    return
1070
1071def adding_station_desc(onc,stdesc):
1072    """ Function to add a station description in a netCDF file
1073      onc= netCDF object
1074      stdesc= station description name, lon, lat, height
1075    """
1076    fname = 'adding_station_desc'
1077
1078    newdim = onc.createDimension('nst',1)
1079
1080    newvar = objfile.createVariable( 'station', 'c', ('nst','StrLength'))
1081    writing_str_nc(newvar, [stdesc[0].replace('!', ' ')], StringLength)
1082
1083    newvar = objfile.createVariable( 'lonstGDM', 'c', ('nst','StrLength'))
1084    Gv = int(stdesc[1])
1085    Dv = int((stdesc[1] - Gv)*60.)
1086    Mv = int((stdesc[1] - Gv - Dv/60.)*3600.)
1087    writing_str_nc(newvar, [str(Gv)+"d" + str(Dv)+"m" + str(Mv)+'s'], StringLength)
1088
1089    if onc.variables.has_key('lon'):
1090        print warnmsg
1091        print '  ' + fname + ": variable 'lon' already exist !!"
1092        print "    renaming it as 'lonst'"
1093        lonname = 'lonst'
1094    else:
1095        lonname = 'lon'
1096
1097    newvar = objfile.createVariable( lonname, 'f4', ('nst'))
1098    basicvardef(newvar, lonname, 'longitude', 'degrees_West' )
1099    newvar[:] = stdesc[1]
1100
1101    newvar = objfile.createVariable( 'latstGDM', 'c', ('nst','StrLength'))
1102    Gv = int(stdesc[2])
1103    Dv = int((stdesc[2] - Gv)*60.)
1104    Mv = int((stdesc[2] - Gv - Dv/60.)*3600.)
1105    writing_str_nc(newvar, [str(Gv)+"d" + str(Dv)+"m" + str(Mv)+'s'], StringLength)
1106
1107    if onc.variables.has_key('lat'):
1108        print warnmsg
1109        print '  ' + fname + ": variable 'lat' already exist !!"
1110        print "    renaming it as 'latst'"
1111        latname = 'latst'
1112    else:
1113        latname = 'lat'
1114
1115    newvar = objfile.createVariable( latname, 'f4', ('nst'))
1116    basicvardef(newvar, lonname, 'latitude', 'degrees_North' )
1117    newvar[:] = stdesc[2]
1118
1119    if onc.variables.has_key('height'):
1120        print warnmsg
1121        print '  ' + fname + ": variable 'height' already exist !!"
1122        print "    renaming it as 'heightst'"
1123        heightname = 'heightst'
1124    else:
1125        heightname = 'height'
1126
1127    newvar = objfile.createVariable( heightname, 'f4', ('nst'))
1128    basicvardef(newvar, heightname, 'height above sea level', 'm' )
1129    newvar[:] = stdesc[3]
1130
1131    return
1132
1133def oper_values(dvals, opers):
1134    """ Function to operate the values according to given parameters
1135      dvals= datavalues
1136      opers= operations
1137    """
1138    fname = 'oper_values'
1139
1140    newdvals = {}
1141    varnames = dvals.keys()
1142
1143    aopers = ['sumc','subc','mulc','divc']
1144
1145    Nopers = len(opers)
1146    for iop in range(Nopers):
1147        vn = varnames[iop]
1148        print vn,'oper:',opers[iop]
1149        if opers[iop] != '-':
1150            opern = opers[iop].split(',')[0]
1151            operv = np.float(opers[iop].split(',')[1])
1152
1153            vvals = np.array(dvals[vn])
1154
1155            if opern == 'sumc': 
1156              newvals = np.where(vvals is None, None, vvals+operv)
1157            elif opern == 'subc': 
1158              newvals = np.where(vvals is None, None, vvals-operv)
1159            elif opern == 'mulc': 
1160              newvals = np.where(vvals is None, None, vvals*operv)
1161            elif opern == 'divc': 
1162              newvals = np.where(vvals is None, None, vvals/operv)
1163            else: 
1164                print errormsg
1165                print '  ' + fname + ": operation '" + opern + "' not ready!!"
1166                print '    availables:',aopers
1167                quit(-1)
1168
1169            newdvals[vn] = list(newvals)
1170        else:
1171            newdvals[vn] = dvals[vn]
1172
1173    return newdvals
1174
1175def WMOcodevar(unitsn, onc, Lstrcode=1024):
1176    """ Function to add a variabe providing description of a units based on WMO code
1177      unitsn= name of the units (all WMO codes derived units must be labelled as
1178        'wmo_code_[num]' associated to a file valled wmo_[num].code)
1179      onc= netCDF object file to add the description
1180      Lstrcode= length of description of codes
1181
1182      wmo_[num].code must have the structure: ('#' for comments)
1183        reference|web page, document reference with the code
1184        short_description|main description of the code
1185        long_description|long description of the code
1186        codeTYPE|type of the code (andy of 'D', 'F', 'I', 'S')
1187        @| Values line giving the start of the values
1188        [val1]|[meaning of first value]
1189        (...)
1190        [valN]|[meaning of last value]
1191    """
1192    fname = 'WMOcodevar'
1193
1194# From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python
1195    folder = os.path.dirname(os.path.realpath(__file__))
1196
1197    code = unitsn.split('_')[2]
1198
1199    infile = folder + '/wmo_' + code +'.code'
1200
1201    if not os.path.isfile(infile):
1202        print warnmsg
1203        print '  ' + fname + ": WMO code file '" + infile + "' does not exist !!"
1204        return
1205    # main expected values
1206    descvals = ['wmo_code', 'reference', 'short_description', 'long_description',    \
1207      'codeTYPE', '@']
1208
1209    availcodetype = ['D', 'F', 'I', 'S']
1210
1211    codevalues = {}
1212    ocode = open(infile, 'r')
1213    inivals = False
1214    codvals = []
1215    codmeanings = []
1216    for line in ocode:
1217        if len(line) > 1 and line[0:1] != '#':
1218            linev = line.replace('\n','').replace('\t','    ').replace('\r','')
1219            if not inivals:
1220                Sv = linev.split('|')[0]
1221                Vv = linev.split('|')[1]
1222                if searchInlist(descvals, Sv): 
1223                    if Sv != '@': codevalues[Sv] = Vv
1224                    else: inivals = True
1225            else:
1226                Svv = linev.split('|')[0]
1227                Vvv = linev.split('|')[1]
1228                codvals.append(Svv)
1229                codmeanings.append(Vvv)
1230                       
1231    # Creating variable
1232    if not searchInlist(onc.dimensions, 'Lstringcode'):
1233        print '  ' + fname + ": Adding string length dimension 'Lstringcode' for " + \
1234          " code descriptions"
1235        newdim = onc.createDimension('Lstringcode', Lstrcode)
1236    Ncodes = len(codvals)
1237    codedimn = 'wmo_code_' + str(code)
1238    if not searchInlist(onc.dimensions, codedimn):
1239        print '  ' + fname + ": Adding '" + codedimn + "' dimension for code " +     \
1240          "description"
1241        newdim = onc.createDimension(codedimn, Ncodes)
1242    onc.sync()
1243
1244    # Variable with the value of the code
1245    if not onc.variables.has_key(codedimn):
1246        if codevalues['codeTYPE'] == 'D':
1247            newvar = onc.createVariable(codedimn, 'f8', (codedimn))
1248            for iv in range(Ncodes): newvar[iv] = np.float64(codvals[iv])
1249        elif codevalues['codeTYPE'] == 'F':
1250            newvar = onc.createVariable(codedimn, 'f', (codedimn))
1251            for iv in range(Ncodes): newvar[iv] = np.float(codvals[iv])
1252        elif codevalues['codeTYPE'] == 'I':
1253            newvar = onc.createVariable(codedimn, 'i', (codedimn))
1254            for iv in range(Ncodes): newvar[iv] = int(codvals[iv])
1255        elif codevalues['codeTYPE'] == 'S':
1256            newvar = onc.createVariable(codedimn, 'c', (codedimn, 'Lstringcode'))
1257            writing_str_nc(newvar, codvals, Lstrcode)
1258        else:
1259            print errormsg
1260            print '  ' + fname + ": codeTYPE '" + codevalues['codeTYPE'] + "' not" + \
1261              " ready !!"
1262            print '   available ones:', availcodetype
1263            quit(-1)
1264
1265        for descv in descvals:
1266            if descv != '@' and descv != 'codeTYPE':
1267                newvar.setncattr(descv, codevalues[descv])
1268    # Variable with the meaning of the code
1269    if not onc.variables.has_key(codedimn+'_meaning'):
1270        print '  '+fname+": Adding '" + codedimn + "_meaning' variable for code " +  \
1271          "description"
1272        newvar = onc.createVariable(codedimn+'_meaning','c',(codedimn,'Lstringcode'))
1273        writing_str_nc(newvar, codmeanings, Lstrcode)
1274        newvar.setncattr('description', 'meaning of WMO code ' + str(code))
1275
1276    onc.sync()
1277
1278    return
1279
1280def EXTRAcodevar(unitsn, onc, Lstrcode=1024):
1281    """ Function to add a variabe providing description of a units based on an extra
1282        code
1283      unitsn= name of the units (all codes derived units must be labelled as
1284        'extra_code_[ref]' associated to a file valled extra_[ref].code)
1285      onc= netCDF object file to add the description
1286      Lstrcode= length of description of codes
1287
1288      extra_[ref].code must have the structure: ('#' for comments)
1289        reference|web page, document reference with the code
1290        short_description|main description of the code
1291        long_description|long description of the code
1292        codeTYPE|type of the code (andy of 'D', 'F', 'I', 'S')
1293        @| Values line giving the start of the values
1294        [val1]|[meaning of first value]
1295        (...)
1296        [valN]|[meaning of last value]
1297    """
1298    fname = 'EXTRAcodevar'
1299
1300# From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python
1301    folder = os.path.dirname(os.path.realpath(__file__))
1302
1303    code = unitsn.split('_')[2]
1304
1305    infile = folder + '/extra_' + code +'.code'
1306
1307    if not os.path.isfile(infile):
1308        print warnmsg
1309        print '  ' + fname + ": EXTRA code file '" + infile + "' does not exist !!"
1310        return
1311    # main expected values
1312    descvals = ['reference', 'short_description', 'long_description',    \
1313      'codeTYPE', '@']
1314
1315    availcodetype = ['D', 'F', 'I', 'S']
1316
1317    codevalues = {}
1318    ocode = open(infile, 'r')
1319    inivals = False
1320    codvals = []
1321    codmeanings = []
1322    for line in ocode:
1323        if len(line) > 1 and line[0:1] != '#':
1324            linev = line.replace('\n','').replace('\t','    ').replace('\r','')
1325            if not inivals:
1326                Sv = linev.split('|')[0]
1327                Vv = linev.split('|')[1]
1328                if searchInlist(descvals, Sv): 
1329                    if Sv != '@': codevalues[Sv] = Vv
1330                    else: inivals = True
1331            else:
1332                Svv = linev.split('|')[0]
1333                Vvv = linev.split('|')[1]
1334                codvals.append(Svv)
1335                codmeanings.append(Vvv)
1336                       
1337    # Creating variable
1338    if not searchInlist(onc.dimensions, 'Lstringcode'):
1339        print '  ' + fname + ": Adding string length dimension 'Lstringcode' for " + \
1340          " code descriptions"
1341        newdim = onc.createDimension('Lstringcode', Lstrcode)
1342    Ncodes = len(codvals)
1343    codedimn = 'extra_code_' + str(code)
1344    if not searchInlist(onc.dimensions, codedimn):
1345        print '  ' + fname + ": Adding '" + codedimn + "' dimension for code " +     \
1346          "description"
1347        newdim = onc.createDimension(codedimn, Ncodes)
1348    onc.sync()
1349
1350    # Variable with the value of the code
1351    if not onc.variables.has_key(codedimn):
1352        if codevalues['codeTYPE'] == 'D':
1353            newvar = onc.createVariable(codedimn, 'f8', (codedimn))
1354            for iv in range(Ncodes): newvar[iv] = np.float64(codvals[iv])
1355        elif codevalues['codeTYPE'] == 'F':
1356            newvar = onc.createVariable(codedimn, 'f', (codedimn))
1357            for iv in range(Ncodes): newvar[iv] = np.float(codvals[iv])
1358        elif codevalues['codeTYPE'] == 'I':
1359            newvar = onc.createVariable(codedimn, 'i', (codedimn))
1360            for iv in range(Ncodes): newvar[iv] = int(codvals[iv])
1361        elif codevalues['codeTYPE'] == 'S':
1362            newvar = onc.createVariable(codedimn, 'c', (codedimn, 'Lstringcode'))
1363            writing_str_nc(newvar, codvals, Lstrcode)
1364        else:
1365            print errormsg
1366            print '  ' + fname + ": codeTYPE '" + codevalues['codeTYPE'] + "' not" + \
1367              " ready !!"
1368            print '   available ones:', availcodetype
1369            quit(-1)
1370
1371        for descv in descvals:
1372            if descv != '@' and descv != 'codeTYPE':
1373                newvar.setncattr(descv, codevalues[descv])
1374    # Variable with the meaning of the code
1375    if not onc.variables.has_key(codedimn+'_meaning'):
1376        print '  '+fname+": Adding '" + codedimn + "_meaning' variable for code " +  \
1377          "description"
1378        newvar = onc.createVariable(codedimn+'_meaning','c',(codedimn,'Lstringcode'))
1379        writing_str_nc(newvar, codmeanings, Lstrcode)
1380        newvar.setncattr('description', 'meaning of EXTRA code ' + str(code))
1381
1382    onc.sync()
1383
1384    return
1385
1386def add_global_PyNCplot(ObjFile, pyscript, funcname, version):
1387    """ Function to add global attributes from 'PyNCplot' to a given netCDF
1388      ObjFile= netCDF file object to which add the global attributes
1389      funcname= name of the function from which file comes from
1390      version= version of the function
1391    """
1392    fname = 'add_global_PyNCplot'
1393
1394    # Global values
1395    ObjFile.setncattr('file_creation', '--------')
1396    ObjFile.setncattr('file_author', 'L. Fita')
1397    newattr = set_attributek(ObjFile, 'file_a_institution', unicode('Centro de ' +   \
1398     'Investigaciones del Mar y la Atm' + unichr(243) + 'sfera (CIMA)'), 'U')
1399    newattr = set_attributek(ObjFile, 'file_a_institution2', unicode('Instituto ' +  \
1400       'Franco-Argentino sobre Estudios de Clima y sus Impactos (CNRS, UMI-3351-' +  \
1401       'IFAECI'), 'U')
1402    newattr = set_attributek(ObjFile, 'center', unicode('Consejo Nacional de ' +     \
1403      'Investigaciones Cient' + unichr(237) + 'ficas y T' + unichr(233) +            \
1404      'cnicas (CONICET)'), 'U')
1405    ObjFile.setncattr('university', 'Universidad de Buenos Aires (UBA)')
1406    ObjFile.setncattr('city', 'C. A. Buenos Aires')
1407    ObjFile.setncattr('country', 'Argentina')
1408    ObjFile.setncattr('tool', 'PyNCplot')
1409    ObjFile.setncattr('url', 'http://www.xn--llusfb-5va.cat/python/PyNCplot')
1410    ObjFile.setncattr('script', pyscript)
1411    if funcname is not None:
1412        ObjFile.setncattr('function', funcname)
1413    ObjFile.setncattr('version', version)
1414
1415    ObjFile.sync() 
1416
1417    return
1418
1419####### ###### ##### #### ### ## #
1420
1421strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
1422  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
1423
1424kindobs=['stations-map','multi-points', 'single-station', 'trajectory']
1425strkObs="kind of observations; 'multi-points': multiple individual punctual obs " +  \
1426  "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
1427  "'trajectory': following a trajectory"
1428
1429parser = OptionParser()
1430parser.add_option("-c", "--comments", dest="charcom",
1431  help="':', list of characters used for comments", metavar="VALUES")
1432parser.add_option("-d", "--descriptionfile", dest="fdesc",
1433  help="description file to use" + read_description.__doc__, metavar="FILE")
1434parser.add_option("-e", "--end_column", dest="endcol",
1435  help="character to indicate end of the column ('space', for ' ', or 'None,[fixcoljumps]' for fixed size columns with [fixcoljumps]: ',' separated list of positions of column ends within line)", metavar="VALUE")
1436parser.add_option("-f", "--file", dest="obsfile",
1437  help="observational file to use", metavar="FILE")
1438parser.add_option("-j", "--jumpBlines", dest="jumpBlines",
1439  help="number of lines to jump from the beggining of file", metavar="VALUE")
1440parser.add_option("-g", "--debug", dest="debug",
1441  help="whether debug is required ('false', 'true')", metavar="VALUE")
1442parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
1443  help=strkObs, metavar="VALUE")
1444parser.add_option("-s", "--stationLocation", dest="stloc", 
1445  help="name ('!' for spaces), longitude, latitude and height of the station (only for 'single-station')", 
1446  metavar="FILE")
1447parser.add_option("-t", "--CFtime", dest="CFtime", help=strCFt, metavar="VALUE")
1448
1449(opts, args) = parser.parse_args()
1450
1451#######    #######
1452## MAIN
1453    #######
1454
1455ofile='OBSnetcdf.nc'
1456
1457if opts.charcom is None:
1458    print warnmsg
1459    print '  ' + main + ': No list of comment characters provided!!'
1460    print '    assuming no need!'
1461    charcomments = []
1462else:
1463    charcomments = opts.charcom.split(':')   
1464
1465if opts.jumpBlines is None:
1466    print warnmsg
1467    print '  ' + main + ': No number of lines to jump from beggining of file provided!!'
1468    print '    assuming no need!'
1469    jumpBlines = 0
1470else:
1471    jumpBlines = int(opts.jumpBlines)
1472
1473
1474if opts.fdesc is None:
1475    print errormsg
1476    print '  ' + main + ': No description file for the observtional data provided!!'
1477    quit(-1)
1478
1479if opts.endcol is None:
1480    print warnmsg
1481    print '  ' + main + ': No list of comment characters provided!!'
1482    print "    assuming 'space'"
1483    endcol = ' '
1484else:
1485    if opts.endcol == 'space':
1486        endcol = ' '
1487    else:
1488        endcol = opts.endcol
1489
1490if opts.obsfile is None:
1491    print errormsg
1492    print '  ' + main + ': No observations file provided!!'
1493    quit(-1)
1494
1495if opts.debug is None:
1496    print warnmsg
1497    print '  ' + main + ': No debug provided!!'
1498    print "    assuming 'False'"
1499    debug = False
1500else:
1501    if opts.debug == 'true':
1502        debug = True
1503    else:
1504        debug = False
1505
1506if not os.path.isfile(opts.fdesc):
1507    print errormsg
1508    print '   ' + main + ": description file '" + opts.fdesc + "' does not exist !!"
1509    quit(-1)
1510
1511if not os.path.isfile(opts.obsfile):
1512    print errormsg
1513    print '   ' + main + ": observational file '" + opts.obsfile + "' does not exist !!"
1514    quit(-1)
1515
1516if opts.CFtime is None:
1517    print warnmsg
1518    print '  ' + main + ': No CFtime criteria are provided !!'
1519    print "    either time is already in CF-format ('timeFMT=CFtime') in '" +        \
1520      opts.fdesc + "'"
1521    print "    or assuming refdate: '19491201000000' and time units: 'hours'"
1522    referencedate = '19491201000000'
1523    timeunits = 'hours'
1524else:
1525    referencedate = opts.CFtime.split(',')[0]
1526    timeunits = opts.CFtime.split(',')[1]
1527
1528if opts.obskind is None:
1529    print warnmsg
1530    print '  ' + main + ': No kind of observations provided !!'
1531    print "    assuming 'single-station': single station on a fixed position at 0,0,0"
1532    obskind = 'single-station'
1533    stationdesc = [0.0, 0.0, 0.0, 0.0]
1534else:
1535    obskind = opts.obskind
1536    if obskind == 'single-station':
1537        if opts.stloc is None:
1538            print errormsg
1539            print '  ' + main + ': No station location provided !!'
1540            quit(-1)
1541        else:
1542            stvals = opts.stloc.split(',')
1543            stationdesc = [stvals[0], np.float(stvals[1]), np.float(stvals[2]),      \
1544              np.float(stvals[3])]
1545    else:
1546        obskind = opts.obskind
1547
1548# Reading description file
1549##
1550description = read_description(opts.fdesc, debug)
1551
1552Nvariables=len(description['varN'])
1553NlongN=len(description['varLN'])
1554NvarU=len(description['varU'])
1555formats = description['varTYPE']
1556
1557if Nvariables != NlongN:
1558    print errormsg
1559    print '  ' + main + ': number of variables:', Nvariables,' and number of ' +     \
1560      'long names', NlongN,' does not coincide!!'
1561    print '    * what is found is _______'
1562    if Nvariables > NlongN:
1563        Nshow = NlongN
1564        for ivar in range(Nshow):
1565            print ivar,'      ',description['varN'][ivar],':', description['varLN'][ivar],\
1566               '[', description['varU'][ivar], '] fmt:', formats[ivar]
1567        print '      missing values for:', description['varN'][Nshow:Nvariables+1]
1568    else:
1569        Nshow = Nvariables
1570        for ivar in range(Nshow):
1571            print ivar,'      ',description['varN'][ivar],':', description['varLN'][ivar],\
1572               '[', description['varU'][ivar], '] fmt:', formats[ivar]
1573        print '      excess of long names for :', description['varLN'][Nshow:NlongN+1]
1574
1575    quit(-1)
1576
1577if Nvariables != NvarU:
1578    print errormsg
1579    print '  ' + main + ': number of variables:', Nvariables,' and number of ' +     \
1580      'units', NvarU,' does not coincide!!'
1581    print '    * what is found is _______'
1582    if Nvariables > NvarU:
1583        Nshow = NvarU
1584        for ivar in range(Nshow):
1585            print ivar,'      ',description['varN'][ivar],':', description['varLN'][ivar],\
1586               '[', description['varU'][ivar], '] fmt:', formats[ivar]
1587        print '      missing values for:', description['varN'][Nshow:Nvariables+1]
1588    else:
1589        Nshow = Nvariables
1590        for ivar in range(Nshow):
1591            print ivar,'      ',description['varN'][ivar],':', description['varLN'][ivar],\
1592               '[', description['varU'][ivar], '] fmt:', formats[ivar]
1593        print '      excess of units for :', description['varU'][Nshow:NvarU+1]
1594
1595    quit(-1)
1596
1597print 'Number of variables', Nvariables
1598print 'Number of long names', len(description['varLN'])
1599
1600if len(formats) != Nvariables: 
1601    print errormsg
1602    print '  ' + main + ': number of formats:',len(formats),' and number of ' +     \
1603      'variables', Nvariables,' does not coincide!!'
1604    print '    * what is found is _______'
1605    if Nvariables > len(formats):
1606        Nshow = len(formats)
1607        for ivar in range(Nshow):
1608            print ivar,'      ',description['varN'][ivar],':', description['varLN'][ivar],\
1609               '[', description['varU'][ivar], '] fmt:', formats[ivar]
1610        print '      missing values for:', description['varN'][Nshow:Nvariables+1]
1611    else:
1612        Nshow = Nvariables
1613        for ivar in range(Nshow):
1614            print ivar,'      ',description['varN'][ivar],':', description['varLN'][ivar],\
1615               '[', description['varU'][ivar], '] fmt:', formats[ivar]
1616        print '      excess of formats for:', formats[Nshow:len(formats)+1]
1617
1618#    quit(-1)
1619
1620# Reading values
1621##
1622datavalues = read_datavalues(opts.obsfile, charcomments, endcol, formats, jumpBlines, 
1623  description['varOPER'], description['MissingValue'], description['varN'], debug)
1624
1625# Total number of values
1626Ntvalues = len(datavalues[description['varN'][0]])
1627if obskind == 'stations-map':
1628    print main + ': total values found:',Ntvalues
1629else:
1630    print main + ': total temporal values found:',Ntvalues
1631
1632objfile = NetCDFFile(ofile, 'w')
1633 
1634# Creation of dimensions
1635##
1636if obskind == 'stations-map':
1637    rowsdim = 'Npoints'
1638    dimlength = Ntvalues
1639else:
1640    rowsdim = 'time'
1641    dimlength = None
1642
1643objfile.createDimension(rowsdim,dimlength)
1644objfile.createDimension('StrLength',StringLength)
1645
1646# Creation of variables
1647##
1648for ivar in range(Nvariables):
1649    varn = description['varN'][ivar]
1650    print "  including: '" + varn + "' ... .. ."
1651
1652    if formats[ivar] == 'D':
1653        newvar = objfile.createVariable(varn, 'f32', (rowsdim), fill_value=fillValueF)
1654        basicvardef(newvar, varn, description['varLN'][ivar],                        \
1655          description['varU'][ivar])
1656        newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn])
1657    elif formats[ivar] == 'F':
1658        newvar = objfile.createVariable(varn, 'f', (rowsdim), fill_value=fillValueF)
1659        basicvardef(newvar, varn, description['varLN'][ivar],                        \
1660          description['varU'][ivar])
1661        newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn])
1662    elif formats[ivar] == 'I':
1663        newvar = objfile.createVariable(varn, 'i', (rowsdim), fill_value=fillValueI)
1664        basicvardef(newvar, varn, description['varLN'][ivar],                        \
1665          description['varU'][ivar])
1666# Why is not wotking with integers?
1667        vals = np.array(datavalues[varn])
1668        for iv in range(Ntvalues):
1669            if vals[iv] is None: vals[iv] = fillValueI
1670        newvar[:] = vals
1671    elif formats[ivar] == 'S':
1672        newvar = objfile.createVariable(varn, 'c', (rowsdim,'StrLength'))
1673        basicvardef(newvar, varn, description['varLN'][ivar],                        \
1674          description['varU'][ivar])
1675        vals = datavalues[varn]
1676        for iv in range(Ntvalues):
1677            if vals[iv] is None: vals[iv] = fillValueS
1678        writing_str_nc(newvar, vals, StringLength)
1679
1680    # Getting new variables to describe certain units as codeWMO_[num] from an
1681    #   external file
1682    if description['varU'][ivar][0:8] == 'wmo_code':
1683        WMOcodevar(description['varU'][ivar], objfile)
1684    # Getting new variables to describe certain units as codeEXTRA_[ref] from an
1685    #   external file
1686    print 'desc varU', description['varU'][ivar][0:10]
1687    if description['varU'][ivar][0:10] == 'extra_code':
1688        EXTRAcodevar(description['varU'][ivar], objfile)
1689
1690# Extra variable descriptions/attributes
1691    if description.has_key('varBUFR'): 
1692        set_attribute(newvar,'bufr_code',description['varBUFR'][ivar])
1693
1694    objfile.sync()
1695
1696# Time variable in CF format
1697##
1698if description['FMTtime'] == 'CFtime':
1699    timevals = datavalues[description['NAMEtime']]
1700    iv = 0
1701    for ivn in description['varN']:
1702        if ivn == description['NAMEtime']:
1703            tunits = description['varU'][iv]
1704            break
1705        iv = iv + 1
1706else:
1707    # Time as a composition of different columns
1708    tcomposite = description['NAMEtime'].find('@')
1709    if tcomposite != -1:
1710        timevars = description['NAMEtime'].split('@')
1711
1712        print warnmsg
1713        print '  ' + main + ': time values as combination of different columns!'
1714        print '    combining:',timevars,' with a final format: ',description['FMTtime']
1715
1716        timeSvals = []
1717        if debug: print '  ' + main + ': creating times _______'
1718        for it in range(Ntvalues):
1719            tSvals = ''
1720            for tvar in timevars:
1721                tSvals = tSvals + datavalues[tvar][it] + ' '
1722
1723            timeSvals.append(tSvals[0:len(tSvals)-1])
1724            if debug: print it, '*' + timeSvals[it] + '*'
1725
1726        timevals = Stringtimes_CF(timeSvals, description['FMTtime'].replace('@',' '),\
1727          referencedate, timeunits, debug)
1728    else:
1729        timevals = Stringtimes_CF(datavalues[description['NAMEtime']],                   \
1730          description['FMTtime'], referencedate, timeunits, debug)
1731
1732    CFtimeRef = referencedate[0:4] +'-'+ referencedate[4:6] +'-'+ referencedate[6:8] +   \
1733      ' ' + referencedate[8:10] +':'+ referencedate[10:12] +':'+ referencedate[12:14]
1734    tunits = timeunits + ' since ' + CFtimeRef
1735
1736if objfile.variables.has_key('time'):
1737    print warnmsg
1738    print '  ' + main + ": variable 'time' already exist !!"
1739    print "    renaming it as 'CFtime'"
1740    timeCFname = 'CFtime'
1741    newdim = objfile.renameDimension('time','CFtime')
1742    newvar = objfile.createVariable( timeCFname, 'f8', ('CFtime'))
1743    basicvardef(newvar, timeCFname, 'time', tunits )
1744else:
1745    if not searchInlist(objfile.dimensions, 'time'):
1746        newdim = objfile.createDimension('time',None)
1747    timeCFname = 'time'
1748    newvar = objfile.createVariable( timeCFname, 'f8', ('time'))
1749    newvar[:] = np.zeros(timevals.shape[0])
1750    basicvardef(newvar, timeCFname, 'time', tunits )
1751
1752set_attribute(newvar, 'calendar', 'standard')
1753if obskind == 'stations-map':
1754    newvar[:] = timevals[0]
1755else:
1756    newvar[:] = timevals
1757
1758# Global attributes
1759##
1760for descn in description.keys():
1761    if descn[0:3] != 'var' and descn[0:4] != 'NAME' and descn[0:3] != 'FMT':
1762        set_attribute(objfile, descn, description[descn])
1763
1764add_global_PyNCplot(objfile, main, 'main', version)
1765
1766objfile.sync()
1767
1768# Adding new variables as function of the observational type
1769## 'multi-points', 'single-station', 'trajectory'
1770
1771if obskind != 'single-station':
1772    adding_complementary(objfile, description, obskind, datavalues, timevals,        \
1773      referencedate, timeunits, Ndim2D, debug)
1774else:
1775# Adding three variables with the station location, longitude, latitude and height
1776    adding_station_desc(objfile,stationdesc)
1777
1778objfile.sync()
1779objfile.close()
1780
1781print main + ": Successfull generation of netcdf observational file '" + ofile + "' !!"
Note: See TracBrowser for help on using the repository browser.