Documentation/Forcings/CfStandard: file_rectifier.py

File file_rectifier.py, 55.8 KB (added by llfita, 8 years ago)

Python file to rectify files to be CFstandard-like for XIOS

Line 
1# Python script to rectify ORCHIDEE netCDF files forcings
2# L. Fita, LMD. April 2016
3#
4# More about:
5#   https://forge.ipsl.jussieu.fr/orchidee/wiki/Documentation/Forcings/CfStandard
6#   
7## e.g. # file_rectifier.py -f alb_bg_jrctip.nc -v bg_albedo -d longitude,latitude -r 0.,1. -l longitude,latitude -m maskvarname@mask -o alb_bg_jrctip_XIOS.nc -T None,range -t alb_bg_jrctip_types.txt -I None,time_counter -Y monthly@means,months@since@0001-01-01@00:00:00,standard -S albedo -L background@albedo@from@MODIS -U '1' -A 'dataset:MODIS!albedo!with!the!JRC-TIP!package:U,article:Pinty et al. 2011:U'
8## e.g. # file_rectifier.py -f lai2D.nc -v LAI -d longitude,latitude -r 0.,7. -l longitude,latitude -m maskabove@0. -o lai2D_XIOS.nc -T None,PFT -t PFTmap_types.txt -I None,time_counter -Y monthly@means,months@since@0001-01-01@00:00:00,standard -S lai -L leaf@area@index -U '-' -z 0.
9## e.g. # file_rectifier.py -f PFTmap.nc -v maxvegetfrac -r 0.,1. -d lon,lat -l lon,lat -m masksumrange@0.9999999,1.e-8,2. -o PFTmap_XIOS.nc -T None,veget -I time_counter,time_counter -t PFTmap_types.txt -Y time,years@since@0001-01-01@00:00:00,standard -S maxvegetfrac
10## e.g. # file_rectifier.py -f soils_param.nc -v soilcolor -r 1,8 -d x,y -l nav_lon,nav_lat -m maskabove@1.e-8 -o soils_param_colorXIOS.nc -t colorvalues.txt -A 'dataset:Henderson-Sellers!and!Wilson:S' -S 'soil_color'
11## e.g. # file_rectifier.py -f soils_param.nc -v soiltext -r 0,7 -d x,y -l nav_lon,nav_lat -m maskabove@1.e-8 -o soils_param_Zobler_textXIOS.nc -t soilsparam_text.txt -S soil_texture -A 'article:Zobler!L.;!1986;!NASA:S'
12## e.g. # file_rectifier.py -f PFTmap_IPCC_2000.nc -v maxvegetfrac -r 0.,1. -d lon,lat -l lon,lat -m masksumrange@0.9999999,1.e-8,2. -o PFTmap_IPCC_2000.XIOS.nc -T None,veget -I time_counter,time_counter -t PFTmap_types.txt -Y time,years@since@0001-01-01@00:00:00,standard -S maxvegetfrac
13## e.g. # file_rectifier.py -f cartepente2d_15min.nc -v pente -r 0.,10. -d longitude,latitude -l longitude,latitude -m maskmissing -o cartepente2d_15min_XIOS.nc -S slope -A 'database:NGDC!-!etopo2:S'
14## e.g. # file_rectifier.py -f orchidee_bbg_clim.nc -v flx_co2_bbg_year -r 0.,2. -d lon,lat -l lon,lat -m maskfile@soils_param_Zobler_textXIOS.nc@mask@lon@lat -o orchidee_bbg_clim_XIOS.nc -S flx_co2_bbg_year -z 0. -A 'dataset:D.!Haugustaine!and!G.!van!der!Werf!12-2003:S'
15## e.g. # file_rectifier.py -f orchidee_fertilizer_1995.nc -v N_qt_WRICE_year -r 0.,6.e+06 -d lon,lat -l lon,lat -m maskabove@0. -o orchidee_fertilizer_1995_N_qt_WRICE_yearXIOS.nc -z 0. -S N_qt_WRICE_year -A 'dataset:A.!BOUWMAN!RIVM!07-2003:S'
16## e.g. # file_rectifier.py -f orchidee_fertilizer_1995.nc -v N_qt_OTHER_year -r 0.,3.e+07 -d lon,lat -l lon,lat -m maskabove@0. -o orchidee_fertilizer_1995_N_qt_OTHER_yearXIOS.nc -z 0. -S N_qt_OTHER_year -A 'dataset:A.!BOUWMAN!RIVM!07-2003:S'
17## e.g. # file_rectifier.py -f routing.nc -v topoind -r 1.,9900. -d x,y -l nav_lon,nav_lat -m maskmissing -o routing_topoind_XIOS.nc -S topo_index -A 'associate:nav_lat!na:S'
18## e.g. # file_rectifier.py -f routing.nc -v trip -r 1.,99. -d x,y -l nav_lon,nav_lat -m maskmissing -o routing_trip_XIOS.nc -S trip -A 'associate:nav_lat!na:S'
19## e.g. # file_rectifier.py -f routing.nc -v basins -r 1.,6740. -d x,y -l nav_lon,nav_lat -m maskmissing -o routing_basins_XIOS.nc -S trip -A 'associate:nav_lat!na:S'
20#
21
22from optparse import OptionParser
23import numpy as np
24from netCDF4 import Dataset as NetCDFFile
25import os
26import re
27import numpy.ma as ma
28
29main = 'file_rectifier.py'
30version='1.0'
31date='April 2016'
32errormsg = 'ERROR -- error -- ERROR -- error'
33warnmsg = 'WARNING -- warning -- WARNING -- warning'
34
35def searchInlist(listname, nameFind):
36    """ Function to search a value within a list
37    listname = list
38    nameFind = value to find
39    >>> searInlist(['1', '2', '3', '5'], '5')
40    True
41    """
42    for x in listname:
43      if x == nameFind:
44        return True
45    return False
46
47def printing_class(classobj):
48    """ Function to print all the values of a given class
49    """
50
51    valscls = vars(classobj)
52    for attrcls in valscls:
53        print attrcls, ':', valscls[attrcls]
54
55def count_cond(var, val, con):
56    """ Function to count values of a variable which attain a condition
57      [var]= variable
58      [val]= value
59      [con]= statistics/condition
60        'eq': equal than [val]
61        'neq': not equal than [val]
62        'lt': less than [val]
63        'le': less and equal than [val]
64        'gt': greater than [val]
65        'ge': greater and equal than [val]
66        'in': inside the range [val[0]] <= [var] <= [val[1]] (where [val]=[val1, val2]
67        'out': inside the range [val[0]] > [var]; [val[1]] < [var] (where [val]=[val1, val2]
68      >>> vals = np.arange(100)
69      >>> count_cond(vals,50,'eq')
70      1
71      >>> count_cond(vals,50,'neq')
72      99
73      >>> count_cond(vals,50,'lt')
74      50
75      >>> count_cond(vals,50,'le')
76      51
77      >>> count_cond(vals,50,'gt')
78      49
79      >>> count_cond(vals,50,'ge')
80      50
81      >>> count_cond(vals,[25,75],'in')
82      51
83      >>> count_cond(vals,[25,75],'out')
84      49
85    """
86    fname = 'count_cond'
87
88    cons = ['eq', 'neq', 'lt', 'le', 'gt', 'ge', 'in', 'out']
89
90    if con == 'eq':
91        Nvals = np.sum(var == val)
92    elif con == 'neq':
93        Nvals = np.sum(var != val)
94    elif con == 'lt':
95        Nvals = np.sum(var < val)
96    elif con == 'le':
97        Nvals = np.sum(var <= val)
98    elif con == 'gt':
99        Nvals = np.sum(var > val)
100    elif con == 'ge':
101        Nvals = np.sum(var >= val)
102    elif con == 'in':
103        if len(val) != 2:
104            print errormsg
105            print '  ' + fname + ": for condition '" + con + "' we must have two " + \
106              " values!!" 
107            print '  we have:',  val
108            quit(-1)
109        Nvals = np.sum((var >= val[0])*(var <= val[1]))
110    elif con == 'out':
111        if len(val) != 2:
112            print errormsg
113            print '  ' + fname + ": for condition '" + con + "' we must have two " + \
114              " values!!" 
115            print '  we have:',  val
116            quit(-1)
117        Nvals = np.sum((var < val[0])+(var > val[1]))
118    else:
119        print errormsg
120        print '  ' + fname + ": condition '" + con + "' not ready!!"
121        print '  only ready:', cons
122        quit(-1)
123
124    return Nvals
125
126def index_mat(mat,val):
127    """ Function to provide the coordinates of a given value inside a matrix
128    index_mat(mat,val)
129      mat= matrix with values
130      val= value to search
131    >>> index_mat(np.arange(27).reshape(3,3,3),22)
132    [2 1 1]
133    """
134
135    fname = 'index_mat'
136
137    matshape = mat.shape
138
139    matlist = list(mat.flatten())
140    ifound = matlist.index(val)
141
142    Ndims = len(matshape)
143    valpos = np.zeros((Ndims), dtype=int)
144    baseprevdims = np.zeros((Ndims), dtype=int)
145
146    for dimid in range(Ndims):
147        baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
148        if dimid == 0:
149            alreadyplaced = 0
150        else:
151            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
152        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
153
154    return valpos
155
156def multi_index_mat(mat,val):
157    """ Function to provide the multiple coordinates of a given value inside a matrix
158    index_mat(mat,val)
159      mat= matrix with values
160      val= value to search
161    >>> vals = np.ones((24), dtype=np.float).reshape(2,3,4)
162    vals[:,:,2] = 0.
163    vals[1,:,:] = np.pi
164    vals[:,2,:] = -1.
165    multi_index_mat(vals,1.)
166    [array([0, 0, 0]), array([0, 0, 1]), array([0, 0, 3]), array([0, 1, 0]), array([0, 1, 1]), array([0, 1, 3])]
167    """
168    fname = 'multi_index_mat'
169
170    matshape = mat.shape
171
172    ivalpos = []
173    matlist = list(mat.flatten())
174    Lmatlist = len(matlist)
175   
176    val0 = val - val
177    if val != val0:
178        valdiff = val0
179    else:
180        valdiff = np.ones((1), dtype = type(val))
181   
182    ifound = 0
183    while ifound < Lmatlist:
184        if matlist.count(val) == 0:
185            ifound = Lmatlist + 1
186        else:
187            ifound = matlist.index(val)
188
189            Ndims = len(matshape)
190            valpos = np.zeros((Ndims), dtype=int)
191            baseprevdims = np.zeros((Ndims), dtype=int)
192
193            for dimid in range(Ndims):
194                baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
195                if dimid == 0:
196                    alreadyplaced = 0
197                else:
198                    alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
199                valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
200            matlist[ifound] = valdiff
201            ivalpos.append(valpos)
202
203    return ivalpos
204
205class variable_inf(object):
206    """ Class to provide the information from a given variable
207    var = object netCDF variable
208    self.name: name of the variable
209    self.dtype: type of the variable
210    self.attributes: list with the name of attributes
211    self.FillValue: value of the missing value
212    self.dimns: name of the dimensions of the variable
213    self.dims: dimensions of the variable
214    self.Ndims: number of dimensions
215    self.dimx: length of dimension along x-axis
216    self.dimy: length of dimension along y-axis
217    self.sname: standard name
218    self.lname: long name
219    self.corr: attribute 'coordinates'
220    self.units: units of the variable
221    """
222
223    def __init__(self, var):
224
225        if var is None:
226            self.name = None
227            self.dimns = None
228            self.dims = None
229            self.Ndims = None
230            self.dimx = None
231            self.dimy = None
232            self.sname = None
233            self.lname = None
234            self.corr = None
235            self.units = None
236            self.FillValue = None
237            self.dtype = None
238            self.attributes = None
239        else:
240            self.name = var._name
241            self.dimns = var.dimensions
242            self.dtype = var.dtype
243            self.attributes = var.ncattrs()
244            self.dims = var.shape
245
246            if searchInlist(self.attributes, 'standard_name'):
247                self.sname = var.getncattr('standard_name')
248            else:
249#                print '    variable_inf.classpy: variable "' + self.name + '" does not have attribute "standard_name"'
250                self.sname = None
251
252            if searchInlist(self.attributes, 'long_name'):
253                self.lname = var.getncattr('long_name')
254            else:
255#                print '    variable_inf.classpy: variable "' + self.name + '" does not have attribute "long_name"'
256                self.lname = None
257
258            if searchInlist(self.attributes, 'coordinates'):
259                self.coor = var.getncattr('coordinates')
260            else:
261#                print '    variable_inf.classpy: variable "' + self.name + '" does not have attribute "coordinates"'
262                self.coor = None
263
264            if searchInlist(self.attributes, 'units'):
265                self.units = var.getncattr('units')
266            else:
267#                print '    variable_inf.classpy: variable "' + self.name + '" does not have attribute "units"'
268                self.units = None
269
270            if searchInlist(self.attributes, '_FillValue'):
271                self.FillValue = var.getncattr('_FillValue')
272            else:
273#                print '    variable_inf.classpy: variable "' + self.name + '" does not have attribute "_FillValue"'
274                self.FillValue = None
275             
276            self.Ndims = len(self.dims)
277            if self.Ndims == 1:
278                self.dimx=self.dims[0]
279            if self.Ndims == 2:
280                self.dimy=self.dims[0]
281                self.dimx=self.dims[1]
282            if self.Ndims == 3:
283                self.dimy=self.dims[1]
284                self.dimx=self.dims[2]
285            if self.Ndims == 4:
286                self.dimy=self.dims[2]
287                self.dimx=self.dims[3]
288
289
290def lonlat2D(ncfobj, lonn, latn):
291    """ Function to create 2D matricies of longitudes and latitudes from the file
292    ncfobj= netCDF object file
293    lonn= longitude name
294    latn= latitude name
295    """
296    fname='lonlat2D'
297
298    lonvar = ncfobj.variables[lonn]
299    latvar = ncfobj.variables[latn]
300
301    loninf = variable_inf(lonvar)
302    latinf = variable_inf(latvar)
303
304    if loninf.Ndims == 1:
305        lonv = lonvar[:]
306        latv = latvar[:]
307        dx=loninf.dims[0]
308        dy=latinf.dims[0]
309        lonval = np.zeros((dy, dx), dtype=float)   
310        latval = np.zeros((dy, dx), dtype=float)   
311        for iy in range(dy):
312            lonval[iy,:] = lonv
313        for ix in range(dx):
314            latval[:,ix] = latv
315    elif loninf.Ndims == 2:
316        lonval = lonvar[:]
317        latval = latvar[:]
318    elif loninf.Ndims == 3:
319        lonval = lonvar[0,:,:]
320        latval = latvar[0,:,:]
321    else:
322        print errormsg
323        print '  ' + fname + ' dimension ', lonvar.Ndims, ' of the lon/lat matrices not ready !!!!'
324
325    return lonval, latval
326
327def varval_ind(var, ind):
328    """ Function to provide a value from a variable providing a value for each index
329      var= variable
330      ind= list of indices for each dimension of the variable
331    """
332    fname = 'varval_ind'
333    if len(var.shape) != len(ind):
334        print errormsg
335        print '  ' + fname + ': different shape:', var.shape,'and indices:', ind, '!!'
336        quit(-1)
337
338    varslice=[]
339    for iv in ind:
340        varslice.append(iv)
341
342    val = varval_ind[slice(varslice)]
343
344    return val
345
346def retype(val, vtype):
347    """ Function to transform a value to a given type
348    retype(val, vtype)
349      [val]= value
350      [vtype]= type to transform
351    >>> retype(0, type(True))
352    False
353    """
354
355    fname = 'retype'
356
357    if val == 'h':
358        print fname + '_____________________________________________________________'
359        print retype.__doc__
360        quit()
361
362    if vtype == type(int(1)):
363        newval = int(val)
364    elif vtype == type(float(1)):
365        newval = float(val)
366#    elif vtype == type(float32(1)):
367#        newval = float32(val)
368    elif vtype == type(np.int(1)):
369        typeinf = np.iinfo(np.int)
370        if (abs(np.float64(val)) > typeinf.max):
371            print warnmsg
372            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
373               vtype,' limits!!'
374            print '  Changing value to kind largest allowed value:', typeinf.max
375            newval = abs(np.float64(val))/np.float64(val) * np.int(typeinf.max)
376        else:
377            newval = np.int(val)
378    elif vtype == type(np.int16(1)):
379        typeinf = np.iinfo(np.int16)
380        if (abs(np.float64(val)) > typeinf.max):
381            print warnmsg
382            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
383               vtype,' limits!!'
384            print '  Changing value to kind largest allowed value:', typeinf.max
385            newval = abs(np.float64(val))/np.float64(val) * np.int16(typeinf.max)
386        else:
387            newval = np.int16(val)
388    elif vtype == type(np.int32(1)):
389        typeinf = np.iinfo(np.int32)
390        if (abs(np.float64(val)) > typeinf.max):
391            print warnmsg
392            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
393               vtype,' limits!!'
394            print '  Changing value to kind largest allowed value:', typeinf.max
395            newval = abs(np.float64(val))/np.float64(val) * np.int32(typeinf.max)
396        else:
397            newval = np.int32(val)
398    elif vtype == type(np.int64(1)):
399        newval = np.int64(val)
400    elif vtype == type(np.float(1)):
401        newval = np.float(val)
402#    elif vtype == type(np.float16(1)):
403#        newval = np.float16(val)
404    elif vtype == type(np.float32(1)):
405        newval = np.float32(val)
406    elif vtype == type(np.float64(1)):
407        newval = np.float64(val)
408    elif vtype == type(True):
409        if val == 0:
410            newval = False
411        else:
412            newval = True
413    elif vtype == '|S1':
414        newval = str(val)
415    else:
416        print errormsg
417        print '  ' + fname + ': variable type "', vtype, '" not ready!!!'
418        quit(-1)
419
420    return newval
421
422def fillvalue_kind(vartype, fillval0):
423    """ Function to provide the fiven fill_Value for a kind of variable
424      [vartype]= type of the variable
425      [fillval0]= value to take as fill value, 'std' for the standard value
426           Float = 1.e20
427           Character = '-'
428           Integer = -99999
429           Integer16 = -99999
430           Float64 = 1.e20
431           Integer32 = -99999
432    """
433    fname = 'fillvalue_kind'
434    fillValue = 1.e20
435    fillValueF = 1.e20
436    fillValueC = '-'
437    fillValueI = -99999
438    fillValueI16 = -99999
439    fillValueF64 = 1.e20
440    fillValueI32 = -99999
441
442    if vartype == type('c'):
443        if fillval0 == 'std':
444            fillval = fillvalueC
445        else:
446            fillval = fillval0
447    elif vartype == type(int(1)):
448        if fillval0 == 'std':
449            fillval = fillValueI
450        else:
451            fillval = int(fillval0)
452    elif vartype == type(np.int16(1)):
453        if fillval0 == 'std':
454            fillval = np.int(9999999)
455        else:
456            fillval = np.int16(fillval0)
457    elif vartype == type(np.int32(1)):
458        if fillval0 == 'std':
459            fillval = np.int32(fillValueI)
460        else:
461            fillval = np.int32(fillval0)
462    elif vartype == type(np.float(1.)):
463        if fillval0 == 'std':
464            fillval = np.float(fillValueF)
465        else:
466            fillval = np.float(fillval0)
467    elif vartype == type(np.float32(1.)):
468        if fillval0 == 'std':
469            fillval = np.float32(fillValueF)
470        else:
471            fillval = np.float32(fillval0)
472    elif vartype == type(np.float64(1.)):
473        if fillval0 == 'std':
474            fillval = np.float64(fillValueF)
475        else:
476            fillval = np.float64(fillval0)
477    else:
478        print errormsg
479        print '  ' + fname + ': variable type', vartype, 'not ready !!'
480        quit(-1)
481
482    return fillval
483
484def read_ASCIIlist(filen):
485    """ Function to build a list from an ASCII lines  file
486      filen= name of the ASCII file
487    """
488    import os
489    fname = 'read_ASCIIlist'
490
491    if not os.path.isfile(filen):
492        print errormsg
493        print '  ' + fname + ":  ASCII types file '" + filen + "' does not exist !!"
494        quit(-1)
495
496    of = open(filen, 'r')
497    filevals = []
498    for linef in of:
499        filevals.append(linef.replace('\n',''))
500    of.close()
501
502    return filevals
503
504def set_attributek(ncv, attrname, attrval, attrkind):
505    """ Sets a value of an attribute of a netCDF variable with a kind. Removes previous attribute value if exists
506    ncvar = object netcdf variable
507    attrname = name of the attribute
508    attrvalue = value of the attribute
509    atrtrkind = kind of attribute: 'S', string ('!' as spaces); 'U', unicode ('!' as spaces); 'I', integer;
510      'Inp32', numpy integer 32; 'R', real; ' npfloat', np.float; 'D', np.float64
511    """
512    fname = 'set_attributek'
513    validk = {'S': 'string', 'U': 'unicode', 'I': 'integer',                         \
514      'Inp32': 'integer long (np.int32)', 'R': 'float', 'npfloat': 'np.float',       \
515      'D': 'float long (np.float64)'}
516
517    if type(attrkind) == type('s'):
518        if attrkind == 'S':
519            attrvalue = str(attrval.replace('!', ' '))
520        elif attrkind == 'U':
521            attrvalue = unicode(attrval.replace('!',' '))
522        elif attrkind == 'I':
523            attrvalue = np.int(attrval)
524        elif attrkind == 'R':
525            attrvalue = float(attrval)
526        elif attrkind == 'npfloat':
527            attrvalue = np.float(attrval)
528        elif attrkind == 'D':
529            attrvalue = np.float64(attrval)
530        else:
531            print errormsg
532            print '  ' + fname + ": '" + attrkind + "' kind of attribute is not ready!"
533            print '  valid values: _______'
534            for key in validk.keys():
535                print '    ', key,':', validk[key]
536            quit(-1)
537    else:
538        if attrkind == type(str('a')):
539            attrvalue = str(attrval.replace('!', ' '))
540        elif attrkind == type(unicode('a')):
541            attrvalue = unicode(attrval.replace('!',' '))
542        elif attrkind == type(np.int(1)):
543            attrvalue = np.int(attrval)
544        elif attrkind == np.dtype('i'):
545            attrvalue = np.int32(attrval)
546        elif attrkind == type(float(1.)):
547            attrvalue = float(attrval)
548        elif attrkind == type(np.float(1.)):
549            attrvalue = np.float(attrval)
550        elif attrkind == np.dtype('float32'):
551            attrvalue = np.array(attrval, dtype='f')
552        elif attrkind == type(np.float32(1.)):
553            attrvalue = np.float32(attrval)
554        elif attrkind == type(np.float64(1.)):
555            attrvalue = np.float64(attrval)
556        else:
557            print errormsg
558            print '  ' + fname + ": '" + attrkind + "' kind of attribute is not ready!"
559            print '  valid values: _______'
560            for key in validk.keys():
561                print '    ', key,':', validk[key]
562            quit(-1)
563
564    attvar = ncv.ncattrs()
565    if searchInlist(attvar, attrname):
566        attr = ncv.delncattr(attrname)
567
568    attr = ncv.setncattr(attrname, attrvalue)
569   
570    return attr
571
572def nctype(vartype):
573    """ Function to provide the string for the variable creation in a netCDF file
574      [vartype]= type of the variable
575    """
576    fname = 'nctype'
577
578    if vartype == type('c'):
579         ncs = 'c'
580    elif vartype == type(int(1)):
581         ncs = 'i'
582    elif vartype == type(np.int16(1)):
583         ncs = 'i'
584    elif vartype == type(np.int32(1)):
585         ncs = 'i8'
586    elif vartype == type(np.float(1.)):
587         ncs = 'f'
588    elif vartype == type(np.float32(1.)):
589         ncs = 'f4'
590    elif vartype == type(np.float64(1.)):
591         ncs = 'f8'
592    else:
593        print errormsg
594        print '  ' + fname + ': variable type', vartype, 'not ready !!'
595        quit(-1)
596
597    return ncs
598
599def writing_str_nc(varo, values, Lchar):
600    """ Function to write string values in a netCDF variable as a chain of 1char values
601    varo= netCDF variable object
602    values = list of values to introduce
603    Lchar = length of the string in the netCDF file
604    """
605    fname = 'writing_str_nc'
606
607    Nvals = len(values)
608    for iv in range(Nvals):   
609        if values[iv] is None:
610            stringv = 'None'
611        else:
612            stringv = values[iv] 
613        charvals = np.zeros((Lchar), dtype='string_')
614        Lstr = len(stringv)
615        charvals[Lstr:Lchar] = ''
616
617        for ich in range(Lstr):
618            chrv = stringv[ich:ich+1]
619            charvals[ich] = stringv[ich:ich+1]
620        varo[iv,:] = charvals
621
622    return
623
624####### ###### ##### #### ### ## #
625maskhelp="(option to gnerate the mask: [maskoption]@[values]\
626'maskfile'@[filename]@[maskname]@[lonname]@[latname]: take mask values from another file [filename] with mask in variable [maskname] and [lonname], [latname] lon lat variables;\
627'maskvarname'@[varname]: reference file contents a variable [varname] with its current mask;\
628'maskbelow'@[value]: mask below a given [value];\
629'maskabove'@[value]: mask above a given [value];\
630'maskrange'@[value1],[value2]: mask within a given range [value1] <= [value] <= [value2];\
631'masksumrange'@[value1],[value2],pvalue3] masking taking the sum of all the types at each grid point (SUMij) which [value2] <= SUMij <= [value1] and [value1] < SUMij <= [value3]);\
632'maskmissing': use missing value for the mask"
633extrattrhelp="',' list of extra attributes to be added to the variable as [attrn]:[value]:[attrk] ([attrk], attribute kind: 'S', string ('!' as spaces); 'U', unicode ('!' as spaces); 'I', integer; 'Inp32', numpy integer 32; 'R', real; ' npfloat', np.float; 'D', np.float64)"
634
635parser = OptionParser()
636parser.add_option("-b", "--lonlat_bounds", dest="Llbounds", 
637  help="[OPTIONAL] when is set and it is given a 'True' value, lon/lat bounds are computed", metavar="FILE")
638parser.add_option("-d", "--lonlatdimname", dest="lonlatdimn",
639  help="current names of the longitude and latitude dimensions as londimn,latdimn in the reference file", metavar="LABEL")
640parser.add_option("-f", "--netCDF_file", dest="ncfile", 
641  help="file to use as reference", metavar="FILE")
642parser.add_option("-l", "--lonlatname", dest="lonlatn",
643  help="current names of the longitude and latitude variables as lonn,latn in the reference file", metavar="LABEL")
644parser.add_option("-m", "--maskname", dest="maskn",
645  help=maskhelp, metavar="LABEL")
646parser.add_option("-o", "--outputfile", dest="ofile", 
647  help="output file", metavar="FILE")
648parser.add_option("-r", "--range", dest="range",
649  help="range of the variable as varmin,varmax", metavar="VALUE")
650parser.add_option("-A", "--ExtraAttributes", dest="extrattrs", help=extrattrhelp,  metavar="LABEL")
651parser.add_option("-T", "--typesvarn", dest="typen",
652  help="[IF APPLICABLE] current name of the type variable and dimension in the reference file (as [typen],[typedimn], " + 
653    "'None' for no variable)",  metavar="LABEL")
654parser.add_option("-I", "--timevarn", dest="timen",
655  help="[IF APPLICABLE] current name of the time variable and dimension in the reference file (as [timen],[timedimn])",
656  metavar="LABEL")
657parser.add_option("-t", "--types_file", dest="typesfile", 
658  help="[OPTIONAL] ASCII file with the values (one consecutive line per value, @ for spaces) for the description of " + 
659  "the types of the variable", metavar="FILE")
660parser.add_option("-v", "--variable", dest="varname",
661  help="variable to check", metavar="LABEL")
662parser.add_option("-z", "--ocean_value", dest="oceanval",
663  help="[OPTIONAL] current ocean value to assign as fill_value of the variable", metavar="VALUE")
664parser.add_option("-S", "--varstd_name", dest="svname",
665  help="[OPTIONAL] label as the standard_name for the variable", metavar="LABEL")
666parser.add_option("-L", "--varlong_name", dest="lvname",
667  help="[OPTIONAL] label as the long_name for the variable (@ for spaces)", metavar="LABEL")
668parser.add_option("-U", "--varunits", dest="vunits",
669  help="[OPTIONAL] label as units for the variable", metavar="LABEL")
670parser.add_option("-Y", "--time_attributes", dest="timeattrs",
671  help="[OPTIONAL] long name and units of the time variable (as [long_name],[units],[calendar] @ for spaces)", metavar="LABEL")
672
673(opts, args) = parser.parse_args()
674
675#######    #######
676## MAIN
677    #######
678ofile = opts.ofile
679
680if not os.path.isfile(opts.ncfile):
681    print errormsg
682    print '  ' + main + ":  File '" + opts.ncfile + "' does not exist !!"
683    quit(-1)
684
685ncobj = NetCDFFile(opts.ncfile, 'r')
686newncobj = NetCDFFile(opts.ofile, 'w')
687
688if opts.typen is not None:
689    typen = opts.typen.split(',')[0]
690    typedimn = opts.typen.split(',')[1]
691
692if opts.timen is not None:
693    timen = opts.timen.split(',')[0]
694    timedimn = opts.timen.split(',')[1]
695
696londimn = opts.lonlatdimn.split(',')[0]
697latdimn = opts.lonlatdimn.split(',')[1]
698# Checking presence of mandatory variables
699##
700lonn = opts.lonlatn.split(',')[0]
701latn = opts.lonlatn.split(',')[1]
702varcheck = [lonn, latn, opts.varname]
703
704for varchk in varcheck:
705    if not searchInlist(ncobj.variables, varchk):
706        print errormsg
707        print '  ' + main + ": File '" + opts.ncfile +"' does not have variable '" + \
708          varchk + "' !!"
709        quit(-1)
710
711# Longitude
712olon = ncobj.variables[lonn]
713loninf = variable_inf(olon)
714
715if len(olon.shape) > 1:
716    newdim = newncobj.createDimension('lon',loninf.dims[1])
717    newdim = newncobj.createDimension('lat',loninf.dims[0])
718    vardims = ('lat', 'lon')
719else:
720    newdim = newncobj.createDimension('lon',loninf.dims[0])
721    vardims = ('lon')
722 
723print " Creation of variable 'longitude' ..."
724newlon = newncobj.createVariable('lon', 'f8', vardims)
725newattr = set_attributek(newlon, 'standard_name', 'longitude', 'S')
726newattr = set_attributek(newlon, 'long_name', 'Longitude', 'S')
727newattr = set_attributek(newlon, 'units', 'degrees_East', 'S')
728newattr = set_attributek(newlon, 'axis', 'X', 'S')
729newattr = set_attributek(newlon, '_CoordinateAxisType', 'Lon', 'S')
730newlon[:] = olon[:]
731if opts.Llbounds is not None and opts.Llbounds == 'True':
732    newattr = set_attributek(newlon, 'bounds', 'lon_bnds', 'S')
733
734# Latitude
735olat = ncobj.variables[latn]
736latinf = variable_inf(olat)
737
738if len(olat.shape) > 1:
739    vardims = ('lat', 'lon')
740else:
741    newdim = newncobj.createDimension('lat',latinf.dims[0])
742    vardims = ('lat')
743 
744print " Creation of variable 'latitude' ..."
745newlat = newncobj.createVariable('lat', 'f8', vardims)
746newattr = set_attributek(newlat, 'standard_name', 'latitude', 'S')
747newattr = set_attributek(newlat, 'long_name', 'Latitude', 'S')
748newattr = set_attributek(newlat, 'units', 'degrees_North', 'S')
749newattr = set_attributek(newlat, 'axis', 'Y', 'S')
750newattr = set_attributek(newlat, '_CoordinateAxisType', 'Lat', 'S')
751newlat[:] = olat[:]
752if opts.Llbounds is not None and opts.Llbounds == 'True':
753    newattr = set_attributek(newlat, 'bounds', 'lat_bnds', 'S')
754
755newncobj.sync()
756if searchInlist(ncobj.variables, 'lon_bnds') or searchInlist(ncobj.variables,        \
757  'bounds_lon'):
758    print '  ' + main + ": file contens 'longitude/latitude bounds'!!"
759    print '    writting them in output file'
760    if searchInlist(ncobj.variables, 'lon_bnds'):
761        olbnd = ncobj.variables['lon_bnds']
762        oLbnd = ncobj.variables['lat_bnds']
763    else:
764        olbnd = ncobj.variables['bounds_lon']
765        oLbnd = ncobj.variables['bounds_lat']
766
767    lbnd = olbnd[:]
768    Lbnd = oLbnd[:]
769    nv = len(olbnd.shape)
770
771    newdim = newncobj.createDimension('nv',nv)
772    newvar = newncobj.createVariable('lon_bnds', 'f8', ('lon','nv'))
773    newvar = lbnd[:]
774    newvar = newncobj.createVariable('lat_bnds', 'f8', ('lat','nv'))
775    newvar = Lbnd[:]
776    newattr = set_attributek(newlon, 'bounds', 'lon_bnds', 'S')
777    newattr = set_attributek(newlat, 'bounds', 'lat_bnds', 'S')
778
779 
780if opts.Llbounds is not None and opts.Llbounds == 'True':
781# Only working for 'regular' projections! lon(i) < lon(i+1); lat(j) < lat(j+1)
782# Bounds
783    print " Creation of variable 'lon_bnds', 'lat_bnds' ..."
784    if len(olat.shape) == 1:
785        nv = 2
786        newdim = newncobj.createDimension('nv',nv)
787        lon_bounds = np.zeros((loninf.dims[0],nv), dtype=np.float32)
788        for i in range(loninf.dims[0]-1):
789            dlon = (olon[i] - olon[i+1])*0.5
790            lon_bounds[i,0] = olon[i] - dlon
791            lon_bounds[i,1] = olon[i] + dlon
792        i = loninf.dims[0]
793        dlon = olon[i-1] - olon[i]       
794        lon_bounds[i,0] = olon[i] - dlon
795        lon_bounds[i,1] = olon[i] + dlon
796        newvar = newncobj.createVariable('lon_bnds', 'f8', ('lon','nv'))
797        newvar[:] = lon_bounds[:]
798   
799        lat_bounds = np.zeros((latinf.dims[0],nv), dtype=np.float32)
800        for i in range(latinf.dims[0]-1):
801            dlat = (olat[i] - olat[i+1])*0.5
802            lat_bounds[i,0] = olat[i] - dlat
803            lat_bounds[i,1] = olat[i] + dlat
804        i = latinf.dims[0]
805        dlat = olat[i-1] - olat[i]       
806        lat_bounds[i,0] = olat[i] - dlat
807        lat_bounds[i,1] = olat[i] + dlat
808        newvar = newncobj.createVariable('lat_bnds', 'f8', ('lat','nv'))
809        newvar[:] = lat_bounds[:]
810    elif len(olat.shape) == 2:
811# At this stage assuming contiguos projections
812        nv = 4
813        newdim = newncobj.createDimension('nv',nv)
814        lon_bounds = np.zeros((loninf.dims[0],loninf.dims[1],nv), dtype=np.float32)
815        lat_bounds = np.zeros((loninf.dims[0],loninf.dims[1],nv), dtype=np.float32)
816        for j in range(1,latinf.dims[0]-1):
817            for i in range(1,loninf.dims[1]-1):
818                lon_bounds[j,i,0] = olon[j,i] - (olon[j+1,i-1] - olon[j,i])*0.5 
819                lon_bounds[j,i,1] = olon[j,i] - (olon[j+1,i+1] - olon[j,i])*0.5 
820                lon_bounds[j,i,2] = olon[j,i] - (olon[j-1,i+1] - olon[j,i])*0.5 
821                lon_bounds[j,i,3] = olon[j,i] - (olon[j-1,i-1] - olon[j,i])*0.5 
822                lat_bounds[j,i,0] = olat[j,i] - (olat[j+1,i-1] - olat[j,i])*0.5 
823                lat_bounds[j,i,1] = olat[j,i] - (olat[j+1,i+1] - olat[j,i])*0.5 
824                lat_bounds[j,i,2] = olat[j,i] - (olat[j-1,i+1] - olat[j,i])*0.5 
825                lat_bounds[j,i,3] = olat[j,i] - (olat[j-1,i-1] - olat[j,i])*0.5 
826   
827        i = loninf.dims[1]-1
828        for j in range(1,latinf.dims[0]):
829            lon_bounds[j,i,0] = olon[j,i] - (olon[j+1,i-1] - olon[j,i])*0.5 
830            lon_bounds[j,i,1] = olon[j,i] + (olon[j+1,i-1] - olon[j,i])*0.5 
831            lon_bounds[j,i,2] = olon[j,i] + (olon[j-1,i-1] - olon[j,i])*0.5 
832            lon_bounds[j,i,3] = olon[j,i] - (olon[j-1,i-1] - olon[j,i])*0.5 
833        j = 0
834        lon_bounds[j,i,0] = olon[j,i] - np.abs(olon[j+1,i-1] - olon[j,i])*0.5 
835        lon_bounds[j,i,1] = olon[j,i] + np.abs(olon[j+1,i-1] - olon[j,i])*0.5 
836        lon_bounds[j,i,2] = olon[j,i] + np.abs(olon[j+1,i-1] - olon[j,i])*0.5 
837        lon_bounds[j,i,3] = olon[j,i] - np.abs(olon[j+1,i-1] - olon[j,i])*0.5 
838        j = latinf.dims[0]
839        lon_bounds[j,i,0] = olon[j,i] - np.abs(olon[j-1,i-1] - olon[j,i])*0.5 
840        lon_bounds[j,i,1] = olon[j,i] + np.abs(olon[j-1,i-1] - olon[j,i])*0.5 
841        lon_bounds[j,i,2] = olon[j,i] + np.abs(olon[j-1,i-1] - olon[j,i])*0.5 
842        lon_bounds[j,i,3] = olon[j,i] - np.abs(olon[j-1,i-1] - olon[j,i])*0.5 
843   
844        j = latinf.dims[0]-1
845        for i in range(1,loninf.dims[1]):
846            lat_bounds[j,i,0] = olat[j,i] + (olat[j-1,i-1] - olat[j,i])*0.5 
847            lat_bounds[j,i,1] = olat[j,i] + (olat[j-1,i+1] - olat[j,i])*0.5 
848            lat_bounds[j,i,2] = olat[j,i] - (olat[j-1,i+1] - olat[j,i])*0.5 
849            lat_bounds[j,i,3] = olat[j,i] - (olat[j-1,i-1] - olat[j,i])*0.5 
850        i = 0
851        lat_bounds[j,i,0] = olat[j,i] + np.abs(olat[j-1,i+1] - olat[j,i])*0.5 
852        lat_bounds[j,i,1] = olat[j,i] + np.abs(olat[j-1,i+1] - olat[j,i])*0.5 
853        lat_bounds[j,i,2] = olat[j,i] - np.abs(olat[j-1,i+1] - olat[j,i])*0.5 
854        lat_bounds[j,i,3] = olat[j,i] - np.abs(olat[j-1,i+1] - olat[j,i])*0.5 
855        i = loninf.dims[1]
856        lat_bounds[j,i,0] = olat[j,i] + np.abs(olat[j-1,i-1] - olat[j,i])*0.5 
857        lat_bounds[j,i,1] = olat[j,i] + np.abs(olat[j-1,i-1] - olat[j,i])*0.5 
858        lat_bounds[j,i,2] = olat[j,i] - np.abs(olat[j-1,i-1] - olat[j,i])*0.5 
859        lat_bounds[j,i,3] = olat[j,i] - np.abs(olat[j-1,i-1] - olat[j,i])*0.5 
860        i = 0
861        j = 0
862        lon_bounds[j,i,0] = olon[j,i] - np.abs(olon[j+1,i+1] - olon[j,i])*0.5 
863        lon_bounds[j,i,1] = olon[j,i] + np.abs(olon[j+1,i+1] - olon[j,i])*0.5 
864        lon_bounds[j,i,2] = olon[j,i] + np.abs(olon[j+1,i+1] - olon[j,i])*0.5 
865        lon_bounds[j,i,3] = olon[j,i] - np.abs(olon[j+1,i+1] - olon[j,i])*0.5 
866        lat_bounds[j,i,0] = olat[j,i] + np.abs(olat[j+1,i+1] - olat[j,i])*0.5 
867        lat_bounds[j,i,1] = olat[j,i] + np.abs(olat[j+1,i+1] - olat[j,i])*0.5 
868        lat_bounds[j,i,2] = olat[j,i] - np.abs(olat[j+1,i+1] - olat[j,i])*0.5 
869        lat_bounds[j,i,3] = olat[j,i] - np.abs(olat[j+1,i+1] - olat[j,i])*0.5 
870     
871        newvar = newncobj.createVariable('lon_bnds', 'f8', ('lat', 'lon','nv'))
872        newvar[:] = lon_bounds[:]
873        newvar = newncobj.createVariable('lat_bnds', 'f8', ('lat', 'lon','nv'))
874        newvar[:] = lat_bounds[:]
875   
876# Variable
877ovar = ncobj.variables[opts.varname]
878varinf = variable_inf(ovar)
879varvals = ovar[:]
880
881if varinf.Ndims == 3 and (opts.typen is None or opts.timen is None):
882    print errormsg
883    print '  ' + main + ': curent variable is 3D !!'
884    print '  var dim names:',varinf.dimns
885    print '  var dims lengths:',varinf.dims
886    print '  current variables in file:', list(ncobj.variables)
887    print "    either provide the values for the type (with '-T [typevals]' flag)"
888    print "    either provide the values for the time (with '-I [timevals]' flag)"
889    quit(-1)
890elif varinf.Ndims == 4 and (opts.typen is None or opts.timen is None):
891    print errormsg
892    print '  ' + main + ': curent variable is 4D !!'
893    print '  var dim names:',varinf.dimns
894    print '  var dims lengths:',varinf.dims
895    print '  current variables in file:', list(ncobj.variables)
896    print "    provide the values for the type (with '-T [typevals]' flag)"
897    print "    and provide the values for the time (with '-I [timevals]' flag)"
898    quit(-1)
899
900if opts.range is not None:
901    varmin = retype(opts.range.split(',')[0], varinf.dtype)
902    varmax = retype(opts.range.split(',')[1], varinf.dtype)
903    range2 = retype(0.5*(varmax + varmin), varinf.dtype)
904else:
905    print warnmsg
906    print '  ' + main + ': no range of valid values is provided!!'
907    print '  assuming validity of all values'
908    range2 = retype(0.5*(np.min(varvals) + np.max(varvals)), varinf.dtype)
909
910dimx=varinf.dims[varinf.Ndims-1]
911dimy=varinf.dims[varinf.Ndims-2]
912
913missf = False
914if searchInlist(varinf.attributes,'missing_value'):
915    missf = True
916    missv = retype(ovar.getncattr('missing_value'), varinf.dtype)
917    print warnmsg
918    print '  ' + main + ": variable '" + opts.varname + "' has missing value!!"
919    if type(varvals) != type(ma.array([1, 2, 3])):
920        print warnmsg
921        print '    ' + main + ": variable '" + opts.varname  + "' is not masked!!"
922        print "      masking it using the 'missing_value':", missv
923        varvals = ma.masked_equal(varvals, missv)
924
925# Types
926typevals = None
927if opts.typen is not None:
928    if typen != 'None':
929        if not searchInlist(ncobj.variables, typen):
930            print errormsg
931            print '  ' + main + ": file '" + opts.filen + "' does not have type " +  \
932              "variable '" + typen + "' !!"
933            quit(-1)
934        otype = ncobj.variables[typen]
935        if len(otype.shape) != 2:
936            print warnmsg
937            print '  ' + main + ": file provided types variable '" + typen +         \
938              "' does not have the expected shape (Ntypes, LString)!!"
939            print "    please, remove the vriable (as None,[typedimn]) from the " +  \
940              "'-T [typevals] flag'"
941            quit(-1)
942        else:
943            Ntypes = otype.shape[0]
944            Lstring = otype.shape[1]
945
946        newdim = newncobj.createDimension('Ntypes',Ntypes)
947        newdim = newncobj.createDimension('Lstring',Lstring)
948        typevalues = otype[:]
949        typevals = []
950        for ityp in range(typevalues.shape[0]):
951            string = ''
952            for il in ityp:
953                string = string + ityp[il:il+1]
954            typevals.append(string)
955    else:
956        Ntypes = varinf.dims[varinf.Ndims-3]
957        Lstring = 250
958        newdim = newncobj.createDimension('Ntypes',Ntypes)
959        newdim = newncobj.createDimension('Lstring',Lstring)
960        if opts.typesfile is None:
961            print errormsg
962            print '  ' + main + ": variable has a types dimension '" + typedimn +    \
963              "' but there is not a variable with values !!"
964            print '  provide an ASCII file with the required types values with ' +   \
965              "the flag '-t [ASCII_file]' with the python call"
966            quit(-1)
967        else:
968            typevals = read_ASCIIlist(opts.typesfile)
969            if len(typevals) != Ntypes:
970                print errormsg
971                print '  ' + main + ': number of values:', len(typevals),            \
972                  "from the ASCII file '" + opts.typesfile + "' does not work " +    \
973                  'for the size:', varinf.dims, 'of the current variable!!'
974                quit(-1)
975else:
976    if opts.typesfile is not None:
977        if not os.path.isfile(opts.typesfile):
978            print errormsg
979            print '  ' + main + ":  ASCII types file '" + opts.typesfile +           \
980              "' does not exist !!"
981            quit(-1)
982        typevals = read_ASCIIlist(opts.typesfile)
983        Ntypes = len(typevals)
984        Lstring = 250
985        newdim = newncobj.createDimension('Ntypes',Ntypes)
986        newdim = newncobj.createDimension('Lstring',Lstring)
987
988# Time
989timevals = None
990if opts.timen is not None:
991    if timen != 'None':
992        if not searchInlist(ncobj.variables, timen):
993            print errormsg
994            print '  ' + main + ": file '" + opts.ncfile + "' does not have time " + \
995              " variable '" + timen + "' !!"
996            quit(-1)
997        otime = ncobj.variables[timen]
998        timeinf = variable_inf(otime)
999
1000        dimt = otime.shape[0]
1001        newdim = newncobj.createDimension('time', dimt)
1002        timevals = otime[:]
1003        if timeinf.lname is None or timeinf.units is None or                         \
1004              not searchInlist(timeinf.attribtues, 'calendar'):
1005            print warnmsg
1006            print '  ' + main + ": current time varibale in file '" + timen +        \
1007              "' does not have units!!"
1008            if opts.timeattrs is None:
1009                print "  add flag '-Y [long_name],[timeunits],[calendar]' (@ for " + \
1010                  'spaces) to the python call to use it as the long_name and the ' + \
1011                  'units for the time variable'
1012                quit(-1)
1013            else:
1014                tlname = opts.timeattrs.split(',')[0].replace('@',' ')
1015                tunits = opts.timeattrs.split(',')[1].replace('@',' ')
1016                tcalendar = opts.timeattrs.split(',')[2].replace('@',' ')
1017        else:
1018            tlname = timeinf.lname
1019            tunits = timeinf.units
1020            tcalendar = otime.getncattr('calendar')
1021    else:
1022        dimt = varinf.dims[0]     
1023        newdim = newncobj.createDimension('time', dimt)
1024        timevals = np.arange((dimt), dtype=np.float64)
1025        if opts.timeattrs is None:
1026            print "  add flag '-Y [long_name],[timeunits]' (@ for spaces) to " +     \
1027              'the python call to use it as the long_name and the units for '+       \
1028              ' the time variable'
1029            quit(-1)
1030        else:
1031            tlname = opts.timeattrs.split(',')[0].replace('@',' ')
1032            tunits = opts.timeattrs.split(',')[1].replace('@',' ')
1033            tcalendar = opts.timeattrs.split(',')[2].replace('@',' ')
1034            tu = tunits.split(' ')[0]
1035            if tu[0:5] == 'month':
1036                print warnmsg
1037                print '  ' + main + ": changing time units from '" + tu +            \
1038                  "' to 'days' !!"
1039                tunits = tunits.replace('months', 'days').replace('month','days')
1040                daysmonth = np.array([0., 31., 28., 31., 30., 31., 30., 31., 31.,     \
1041                  30., 31., 30.], dtype=np.float32)
1042                timevals = np.zeros((dimt), dtype=np.float32)
1043                for it in range(1,dimt):
1044                    timevals[it] = np.sum(daysmonth[0:it+1])
1045
1046lons, lats = lonlat2D(ncobj, lonn, latn)
1047
1048# Variable
1049print " Creation of variable '" + opts.varname + "' ..."
1050
1051fillvalt = fillvalue_kind(varinf.dtype, 'std')
1052
1053if varinf.FillValue is None:
1054    print warnmsg
1055    print "  variable does not have 'fill_value' attribute!!"
1056    if opts.oceanval is None:
1057        print warnmsg
1058        print "    you can add flag '-z [oceanval]' to the python call to use it " + \
1059          ' as the current value to mask the ocean values'
1060        oceanval = None
1061    else:
1062        oceanval = retype(opts.oceanval, varinf.dtype)   
1063        print '  using:',oceanval,'as a complementary way to generate the ocean mask'
1064else:
1065    oceanval = retype(opts.oceanval, varinf.dtype)
1066    if oceanval != varinf.FillValue:
1067        print warnmsg
1068        print '  current fill_value:', varinf.FillValue,' and ocean values:',        \
1069          oceanval,' do not coincide!!'
1070        print '  using ocean value to create the mask'
1071
1072slicemask = []
1073for idim in range(varinf.Ndims-2):
1074    slicemask.append(0)
1075slicemask.append(slice(0,varinf.dims[varinf.Ndims-2]))
1076slicemask.append(slice(0,varinf.dims[varinf.Ndims-1]))
1077
1078# Mask
1079print " Creation of 'mask' ..."
1080newvar = newncobj.createVariable('mask', 'i4', ('lat','lon'))
1081maskn = opts.maskn.split('@')[0]
1082if maskn == 'maskfile':
1083    maskf = opts.maskn.split('@')[1]
1084    maskv = opts.maskn.split('@')[2]
1085    maskl = opts.maskn.split('@')[3]
1086    maskL = opts.maskn.split('@')[4]
1087
1088    print warnmsg
1089    print '  ' + main + ": taking mask values from file '" + maskf + "' variable '" +\
1090      maskv + "' !!"
1091    print "    make sure that it is on the same projection and resolution as '" +    \
1092      opts.ncfile + "' !!"
1093
1094    compncobj = NetCDFFile(maskf,'r')
1095    if not searchInlist(compncobj.variables,maskv):
1096        print errormsg
1097        print '  ' + fname + ": file '" + maskf + '" does not have mask variable "' +\
1098          maskn + "' !!"
1099        quit(-1)
1100    if not searchInlist(compncobj.variables,maskl):
1101        print errormsg
1102        print '  ' + fname + ": file '" + maskf + '" does not have lon variable "' + \
1103          maskl + "' !!"
1104        quit(-1)
1105    if not searchInlist(compncobj.variables,maskL):
1106        print errormsg
1107        print '  ' + fname + ": file '" + maskf + '" does not have lat variable "' + \
1108          maskL + "' !!"
1109        quit(-1)
1110
1111    omask = compncobj.variables[maskv]
1112    maskvals = omask[:]
1113    complons, complats = lonlat2D(compncobj, maskl, maskL)
1114
1115# Checking a bit
1116    if maskvals.shape[1] != dimx:
1117        print errormsg
1118        print '    ' + main + ": complementary mask file '" + maskf + "' x-shape ",  \
1119          maskvals.shape[1], "is different than '" + opts.ncfile + "': ", dimx,'!!'
1120        quit(-1)
1121    if maskvals.shape[0] != dimy:
1122        print errormsg
1123        print '    ' + main + ": complementary mask file '" + maskf + "' y-shape ",  \
1124          maskvals.shape[0], "is different than '" + opts.ncfile + "': ", dimy,'!!'
1125        quit(-1)
1126    if not np.all(complons - lons == 0.):
1127        print errormsg
1128        print '    ' + main + ": complementary mask file '" + maskf + "' and '" +    \
1129          opts.ncfile + "' have different longitude projection !!"
1130        quit(-1)
1131    if not np.all(complats - lats == 0.):
1132        print errormsg
1133        print '    ' + main + ": complementary mask file '" + maskf + "' and '" +    \
1134          opts.ncfile + "' have different latitude projection !!"
1135        quit(-1)
1136    compncobj.close()
1137
1138elif maskn == 'maskvarname':
1139    maskv = opts.maskn.split('@')[1]
1140    if not searchInlist(ncobj.variables,maskv):
1141        print errormsg
1142        print '  ' + fname + ": file '" + opts.ncfile + '" does not have variable "'+\
1143          maskn + "' !!"
1144        quit(-1)
1145    else:
1146        omask = ncobj.variables[maskv]
1147        maskvals = omask[:]
1148elif maskn == 'maskbelow':
1149    maskv = retype(opts.maskn.split('@')[1], varinf.dtype)
1150    if missf:
1151        varvalsmask = varvals.filled(range2).reshape(np.prod(varvals.shape[0:        \
1152          varinf.Ndims-2]),dimy,dimx)
1153    else:
1154        varvalsmask = varvals.reshape(np.prod(varvals.shape[0:varinf.Ndims-2]),dimy, \
1155          dimx)
1156    maskvals0 = np.where(np.all(varvalsmask <= maskv, axis=0), 1, 0)
1157    if oceanval is not None:
1158        maskvals00 = np.where(np.all(varvalsmask == oceanval, axis=0), 0, 1)
1159        maskvals0 = maskvals0 * maskvals00
1160#    maskvals = maskvals0[tuple(slicemask)]
1161    maskvals = maskvals0
1162elif maskn == 'maskabove':
1163    maskv = retype(opts.maskn.split('@')[1], varinf.dtype)
1164    if missf:
1165        varvalsmask = varvals.filled(range2).reshape(np.prod(varvals.shape[0:        \
1166          varinf.Ndims-2]),dimy,dimx)
1167    else:
1168        varvalsmask = varvals.reshape(np.prod(varvals.shape[0:varinf.Ndims-2]),dimy, \
1169          dimx)
1170    maskvals0 = np.where(np.all(varvalsmask >= maskv, axis=0), 1, 0)
1171    if oceanval is not None:
1172        maskvals00 = np.where(np.all(varvalsmask == oceanval, axis=0), 0, 1)
1173        maskvals0 = maskvals0 * maskvals00
1174#    maskvals = maskvals0[tuple(slicemask)]
1175    maskvals = maskvals0
1176elif maskn == 'maskrange':
1177    maskv1 = retype(opts.maskn.split('@')[1].split(',')[0], varinf.dtype)
1178    maskv2 = retype(opts.maskn.split('@')[1].split(',')[1], varinf.dtype)
1179    if missf:
1180        varvalsmask = varvals.filled(range2).reshape(np.prod(                        \
1181          varvals.shape[0:varinf.Ndims-2]),dimy,dimx)
1182    else:
1183        varvalsmask = varvals.reshape(np.prod(varvals.shape[0:varinf.Ndims-2]),dimy, \
1184          dimx)
1185    maskvals01 = np.where(np.all(varvalsmask >= maskv1, axis=0), 1, 0)
1186    maskvals02 = np.where(np.all(varvalsmask <= maskv2, axis=0), 1, 0)
1187    maskvals0 = maskvals01 * maskvals02
1188    if oceanval is not None:
1189        maskvals00 = np.where(np.all(varvalsmask == oceanval, axis=0), 0, 1)
1190        maskvals0 = maskvals0 * maskvals00
1191#    maskvals = maskvals0[tuple(slicemask)]
1192    maskvals = maskvals0
1193elif maskn == 'masksumrange':
1194    maskv1 = retype(opts.maskn.split('@')[1].split(',')[0], varinf.dtype)
1195    maskv2 = retype(opts.maskn.split('@')[1].split(',')[1], varinf.dtype)
1196    maskv3 = retype(opts.maskn.split('@')[1].split(',')[2], varinf.dtype)
1197    if missf:
1198        varvalsmask = varvals.filled(range2).reshape(np.prod(varvals.shape[0:        \
1199          varinf.Ndims-2]),dimy,dimx)
1200    else:
1201        varvalsmask = varvals.reshape(np.prod(varvals.shape[0:varinf.Ndims-2]),dimy, \
1202          dimx)
1203    print ' masksumrange: _______'
1204    print '    ', maskv2,'<= varsum <=',maskv1
1205    print '    ', maskv1,'< varsum <=',maskv3
1206    varsum = np.sum(varvals, axis=varinf.Ndims-3)
1207    maskvals1a = np.where(varsum >= maskv2, 1, 0)
1208    maskvals1b = np.where(varsum <= maskv1, 1, 0)
1209    maskvals2a = np.where(varsum > maskv1, 1, 0)
1210    maskvals2b = np.where(varsum <= maskv3, 1, 0)
1211    if varinf.Ndims == 4:
1212        maskvals0 = maskvals1a[0,:,:]*maskvals1b[0,:,:]+maskvals2a[0,:,:]*           \
1213          maskvals2b[0,:,:]
1214    elif varinf.Ndims == 3:
1215        maskvals0 = maskvals1a*maskvals1b+maskvals2a*maskvals2b
1216    else:
1217        print errormsg
1218        print '  ' + main + ": masking by '" + maskn + "' with shape:", varinf.dims, \
1219          'not ready!!'
1220        quit(-1)
1221    if oceanval is not None:
1222        maskvals00 = np.where(np.all(varvalsmask == oceanval, axis=0), 0, 1)
1223        maskvals0 = maskvals0 * maskvals00
1224    if type(maskvals0) == type(ma.masked_array(np.arange(1))):
1225        msk = maskvals0.mask
1226        maskvals = np.where(msk, 0, 1)
1227    else:
1228        maskvals = maskvals0
1229elif maskn == 'maskmissing':
1230    if missf:
1231        maskvals = np.where(varvals.mask[tuple(slicemask)], 0, 1)
1232    else:
1233        print errormsg
1234        print '  ' + main + ": it should be masked following the missing value '" +  \
1235          mask + "' but variable does not have attribute 'missing_value'!!"
1236        quit(-1)
1237else:
1238    print errormsg
1239    print '  ' + main + ": mask type '" + maskn + "' not ready!!"
1240    quit(-1)
1241
1242# This should not be needed, isn't?
1243#if varinf.FillValue is not None:
1244#    msk = msk * varvals.msk
1245#    maskvals = np.where(msk, 0, 1)
1246newvar[:] = maskvals[:]
1247
1248# Dimensions new varibale
1249vardims0 = []
1250for dimn in varinf.dimns:
1251    if dimn == londimn: vardims0.append('lon')
1252    if dimn == latdimn: vardims0.append('lat')
1253    if opts.typen is not None and dimn == typedimn: vardims0.append('Ntypes')
1254    if opts.timen is not None and dimn == timedimn: vardims0.append('time')
1255if not searchInlist(vardims0, 'lon') or not searchInlist(vardims0, 'lat'):
1256    print errormsg
1257    print '  ' + main + ": new variable will not have mandatory 'lon', 'lat' " +     \
1258      "dimensions!!"
1259    print '    check given names of the dimensions:', londimn, ',', latdimn
1260    quit(-1)
1261if len(vardims0) != varinf.Ndims:
1262    print errormsg
1263    print '  ' + main + ': new variable will not have same number of dimensions:',   \
1264      len(vardims0), 'as original one', varinf.Ndims,'!!'
1265    print '    check given names of dimensions lon, lat:', londimn, ',', latdimn
1266    if opts.typen is not None: print '    types:', typedimn
1267    if opts.timen is not None: print '    time:', timedimn
1268    quit(-1)
1269
1270vardims = tuple(vardims0)
1271
1272# Masking variable values according to the mask and the range
1273varmask0 = np.where(maskvals == 0, True, False)
1274varmask = np.zeros(varinf.dims, dtype=bool)
1275varmask[...,:,:] = varmask0
1276if varmin is not None:
1277    if np.min(ovar[:]) < varmin:
1278        Nvals = np.sum(varvals < varmin)
1279        print warnmsg
1280        print '  ' + main + ': masking ', Nvals,'values of variable lower than ',    \
1281         'varmin:', varmin
1282    if np.max(ovar[:]) > varmax:
1283        Nvals = np.sum(varvals > varmax)
1284        print warnmsg
1285        print '  ' + main + ': masking ', Nvals,'values of variable greater than ',  \
1286          'varmax:', varmax
1287
1288    if missf: 
1289        varvalsfill = varvals.filled(varmax + range2)
1290        minrangevarval = ma.masked_less(varvalsfill, varmin)
1291        minrangevarval = np.where(minrangevarval == varmax+range2, missv,            \
1292          minrangevarval)
1293    else:
1294        minrangevarval = ma.masked_less(varvals, varmin)
1295
1296    if missf: 
1297        varvalsfill = varvals.filled(varmin - range2)
1298        maxrangevarval = ma.masked_greater(varvalsfill, varmax)
1299        maxrangevarval = np.where(maxrangevarval == varmax-range2, missv,            \
1300          maxrangevarval)
1301    else:
1302        maxrangevarval = ma.masked_greater(varvals, varmax)
1303
1304    rangevarmask = minrangevarval.mask*maxrangevarval.mask
1305
1306    if missf:
1307        fillvarval = ma.array(varvals.filled(missv), mask=rangevarmask)
1308    else:
1309        fillvarval = ma.array(varvals, mask=rangevarmask)
1310else:
1311    fillvarval = varvals
1312
1313fillvarval = ma.array(fillvarval, mask=varmask)
1314
1315newvar = newncobj.createVariable(opts.varname, nctype(varinf.dtype), vardims,        \
1316  fill_value=fillvalt)
1317newvar[:] = fillvarval[:]
1318
1319if varinf.sname is None:
1320    print warnmsg
1321    print "  variable does not have 'standard_name' attribute!!"
1322    if opts.svname is None:
1323        print errormsg
1324        print "    add flag '-S [standard_name]' to the python call to assign the" + \
1325          ' attribute'
1326        quit(-1)
1327    else:
1328        newattr = set_attributek(newvar, 'standard_name', opts.svname, 'S')
1329else:
1330    newattr = set_attributek(newvar, 'standard_name', varinf.sname, 'S')
1331
1332if varinf.lname is None:
1333    print warnmsg
1334    print "  variable does not have 'long_name' attribute!!"
1335    if opts.lvname is None:
1336        print errormsg
1337        print "    add flag '-L [long_name]' (@ for spaces) to the python call to" + \
1338          'assign the attribute'
1339        quit(-1)
1340    else:
1341        newattr = set_attributek(newvar, 'long_name', opts.lvname.replace('@',' '),  \
1342          'S')
1343else:
1344    newattr = set_attributek(newvar, 'long_name', varinf.lname, 'S')
1345
1346if varinf.units is None:
1347    print warnmsg
1348    print "  variable does not have 'units' attribute!!"
1349    if opts.vunits is None:
1350        print errormsg
1351        print "    add flag '-U [units]' to the python call to assign the attribute"
1352        quit(-1)
1353    else:
1354        newattr = set_attributek(newvar, 'units', opts.vunits, 'S')
1355else:
1356    newattr = set_attributek(newvar, 'units', varinf.units, 'S')
1357
1358newattr = set_attributek(newvar, 'coordinates', 'lon lat', 'S')
1359if opts.range is not None:
1360    newattr = set_attributek(newvar, 'valid_min', varmin, varinf.dtype)
1361    newattr = set_attributek(newvar, 'valid_max', varmax, varinf.dtype)
1362
1363if opts.extrattrs is not None:
1364    if opts.extrattrs.find(',') != -1:
1365        attrs = opts.extrattrs.split(',')
1366    else:
1367        attrs = [opts.extrattrs]
1368    for attr in attrs:
1369        if len(attr.split(':')) != 3:
1370            print errormsg
1371            print '  ' + main + ": wrong number of values:", len(attr.split(':')),   \
1372              'for extra attribute !!'
1373            print "    given values:'" + attr + "' 3 ':' separated are expected"
1374            print '    which gives:', attr.split(':')
1375            quit(-1) 
1376        attrn = attr.split(':')[0]
1377        attrv = attr.split(':')[1]
1378        attrk = attr.split(':')[2]
1379        newattr = set_attributek(newvar, attrn, attrv, attrk)
1380
1381newncobj.sync()
1382# Types variable
1383if typevals is not None:
1384    print " Creation of 'typenames' ..."
1385    newvar = newncobj.createVariable('typenames', 'c', ('Ntypes', 'Lstring'))
1386    writing_str_nc(newvar, typevals, Lstring)
1387
1388# Time variable
1389if timevals is not None:
1390    print " Creation of 'time' ..."
1391    newvar = newncobj.createVariable('time', 'f8', ('time'))
1392    newattr = set_attributek(newvar, 'standard_name', 'time', 'S')
1393    newattr = set_attributek(newvar, 'long_name', tlname, 'S')
1394    newattr = set_attributek(newvar, 'units', tunits, 'S')
1395    newattr = set_attributek(newvar, 'calendar', tcalendar, 'S')
1396    newvar[:] = timevals[:]
1397    newattr = set_attributek(newvar, 'axis', 'T', 'S')
1398    newattr = set_attributek(newvar, '_CoordinateAxisType', 'time', 'S')
1399
1400
1401newncobj.sync()
1402
1403newattr = set_attributek(newncobj, 'creation',"file created with " + main, 'S')
1404newattr = set_attributek(newncobj, 'py_script_version',version, 'S')
1405newattr = set_attributek(newncobj, 'py_script_date',date, 'S')
1406newattr = set_attributek(newncobj, 'py_script_author',"L. Fita", 'S')
1407newattr = set_attributek(newncobj, 'py_script_author_institution',                   \
1408  unicode('Laboratoire de M' + unichr(233)+'t'+unichr(233)+'orologie Dynamique'), 'U')
1409
1410# Getting global attributes
1411ncgattrs = ncobj.ncattrs()
1412for attr in ncgattrs:
1413    attrv = ncobj.getncattr(attr)
1414    attrk = type(attrv)
1415    newattr = set_attributek(newncobj, attr, attrv, attrk)
1416
1417ncobj.close()
1418newncobj.close()
1419print main + ": successfull written of file '" + ofile + "' !!"