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

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

Fixing prints

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